I am trying to solve a 1-D 4th order differential equation using finte difference method. I am implementing the algorithm in maple and was able to get approximate solutions for second order differential equations. Now I stepped up to 4th order, since this is the actual type of equation that I am trying to solve, but I am failing to find numeric solutions.
Below you can find the maple code. Here is a rough description of what I am doing in my FDM implementation:
Right now, I am trying to solve the following boundary problem:
$$ 1000 f(x) + f''''(x) = 1 \qquad f(0)=0 \qquad f'(0)=0 \qquad f''(L)=0 \qquad f'''(L)=0 $$
I discretize the differential equation using backwards difference method ($d$ being the step size):
$$ 1000 f_i + \frac{f_i - 4 f_{i - 1} + 6 f_{i - 2} - 4 f_{i - 3} + f_{i - 4}}{d^4}= 1 \quad i \in [1, n] $$
Now I can build the matrix for the lineare system of equations. Next, I discretize the boundary conditions:
$$ f(0)=0 \rightarrow f_1 = 0 $$ $$ f'(0)=0 \rightarrow \frac{f_{i+1}-f_i}{d} = 0 \rightarrow f_{2}-f_1 = 0 \quad \text{(forward difference here)} $$
$$ f''(L)=0 \rightarrow \frac{f_i - 2 f_{i - 1} + f_{i - 2}}{d^2} = 0 \rightarrow f_n - 2 f_{n - 1} + f_{n - 2} = 0 $$
$$ f'''(L)=0 \rightarrow \frac{f_i - 3 f_{i - 1} + 3 f_{i - 2} - f_{i - 3}}{d^3} = 0 \rightarrow f_n - 3 f_{n - 1} + 3 f_{n - 2} - f_{n - 3} = 0 $$
Now I replace first two and the last two equations of my liner system of equations with the previous four respectively.
Then I define the right hand side vector and solve the system of equations. As you can see in the picture, the numeric solution is not even close to the precise solution. It even seems like the last two boundary condtions don't even have an influence on the numeric solution, so I am not sure if it is implemented correctly. I would really appreciate your help!
The maple code:
restart
with(LinearAlgebra):
with(plots):
param := [L=1]
#define the actual differential equation and solve it
dgl_real:= 1000*g(x) + diff(g(x), x, x, x, x)=1;
res_real := dsolve([eval(dgl_real,param), g(0)=0,D(g)(0)=0, D(D(g))(L)=0, D(D(D(g)))(L)=0], g(x)):
#define the left hand side of the descritized differential equation
dgl:=1000*f[i] + ddddf
#approximation of f'(x)
df := (f[i]-f[i-1])/d
#approximation of f''(x)
ddf := (f[i]-2*f[i-1]+f[i-2])/d^2
#approximation of f'''(x)
dddf := (f[i]-3*f[i-1]+3*f[i-2]-f[i-3])/d^3
#approximation of f''''(x)
ddddf := (f[i]-4*f[i-1]+6*f[i-2]-4*f[i-3]+f[i-4])/d^4
#left hand side of boundary condition f'(0) (left hand side)
BC1 := eval(df, [i=1])
#left hand side of boundary condition f''(L) (left hand side)
BC2 := eval(ddf,[i=n]);
#left hand side of boundary condition f'''(L)
BC3 := eval(dddf,[i=n]);
n := 100: #number of points
d := L/(n-1): #step size
#define the values for x
X:=Vector(n):
for j from 1 to n do
x[j] := (j-1)*d;
X[j] := (j-1)*d;#used for plotting
end do:
A := Matrix(n, n):
#for k from 1 to n do
#A[1, k] := coeff(BC1, f[k]);
#end do:
#boundary conditions on left end: f(0), f'(0)
A[1,1]:=1:
A[2,1]:=-1/d:
A[2,2]:=1/d:
# Populate the matrix
for j from 3 to (n-2) do
term := eval(dgl, [i=j]);
for k from 1 to n do
A[j, k] := coeff(term, f[k]);
end do;
end do:
#boundary condtions on the right end: f''(L), f'''(L)
for k from 1 to n do
A[n-1, k] := coeff(BC2, f[k]);
end do:
for k from 1 to n do
A[n, k] := coeff(BC3, f[k]);
end do:
A[n-5..n, n-5..n]
Q := Vector(n):
for j from 1 to n do
Q[j] := 1;
end do:
Mat:=eval(A, param):
Q[1] := 0:
Q[2] := 0:
Q[n-1] := 0:
Q[n] := 0:
u := MatrixInverse(Mat).Q:
plot1:=pointplot(eval(X, param), u):
plot2:=plot(rhs(eval(res_real,param)), x=0..eval(L,param)):
display(plot1, plot2)
Follows a MATHEMATICA script which formally is quite readable for a programmer
Here in blue the finite difference solution, and in red the exact solution.