# Finite difference method for 1D Poisson equation

I want to solve the following 1D Poisson equation using finite difference method:

$$u'' = 6 x,\ u' (0) = 0,\ u (1) = 1$$

where $$h=1/3$$ i.e I need to find $$u(0)$$, $$u(1/3)$$ and $$u(2/3)$$.

I construct the linear system $$Au=b$$, where

      A = {{-2, 2, 0}, {1, -2, 1}, {0, 1, -2}};
b = {0, 2/9, -5/9};
LinearSolve[A, b]


I got {1/9, 1/9, 1/3} which is wrong it should be {0, 1/27, 8/27, since the exact $$u=x^3$$, I don't know what is the mistake?

Please what happen if the boundary conditions change to

b)

$$u(0)=0,\ u'(1)=3$$

c)

$$u'(0)=0,\ u'(1)=3$$?

Thanks.

• No error. {1/9, 1/9, 1/3} is correct. I just solved it using FDM and got same answer. Apr 27 '20 at 21:21
• Dear Nasser, thanks for your response but the error very big? why? Apr 27 '20 at 21:25

I solved it using FDM and the answer is correct.

ClearAll[u, x];
h = 1/3;
eq1 = -2 u0 + 2 u1 == 0;
eq2 = u0 - 2 u1 + u2 == (6*h)*h^2;
eq3 =  u1 - 2 u2 == (6*2 h - 1/h^2)*h^2;
pts = Solve[{eq1, eq2, eq3}, {u0, u1, u2}] sol = u[x] /. First@DSolve[{u''[x] == 6*x, u' == 0, u == 1}, u[x], x]
p1 = Plot[sol, {x, 0, 1}];
p2 = ListPlot[{{0, 1/9}, {h, 1/9}, {2 h, 1/3}, {1, 1}},  PlotStyle -> Red];
Show[p1, p2] The error is large, since $$h$$ is large. With more points, it will improve.

Here is a quick hack to show the effect of adding more points to FDM makeA[n_] := Module[{A, i, j},
A = Table[0, {i, n}, {j, n}];
Do[
Do[
A[[i, j]] = If[i == j, -2, If[i == j + 1 || i == j - 1, 1, 0]],
{j, 1, n}
],
{i, 1, n}
];
A[[1, 2]] = 2;
A
];
makeB[n_, h_, f_] := Module[{b, i},
b = Table[0, {i, n}];
Do[
b[[i]] = If[i == 1, 0,
If[i < n, f[(i - 1)*h]*h^2, (f[(i - 1)*h] - 1/h^2)*h^2]
]
, {i, 1, n}
];
b
];
f[x_] := 6*x;(*RHS of ode*)
Manipulate[
Module[{h, A, b, sol, solN, p1, p2, x},
h = 1/(nPoints - 1);
A = makeA[nPoints - 1];
b = makeB[nPoints - 1, h, f];
sol = LinearSolve[A, b];
solN = Table[{n*h, sol[[n + 1]]}, {n, 0, nPoints - 2}];
AppendTo[solN, {1, 1}];
p1 = Plot[x^3, {x, 0, 1}];
p2 = ListLinePlot[solN, PlotStyle -> Red, Mesh -> All];
Grid[{
{Row[{" h = ", NumberForm[N@h, {5, 4}]}]},
{Show[p1, p2,
PlotLabel -> "Red is numerical, Blue is exact solution",
GridLines -> Automatic,
GridLinesStyle -> LightGray, ImageSize -> 400]}}
]
],
{{nPoints, 3, "How many points?"}, 3, 20, 1, Appearance -> "Labeled"},
TrackedSymbols :> {nPoints}
]

• which value of h can provide good results, I reduce it but still not good results? Apr 27 '20 at 21:32
• You need to use some norm measure. Then you specify the error you are willing to accept and use the h which gives norm measure you want. There are 3 basic types of norm measures. max-norm, 1-norm or 2-norm. I stole this from this page it has more information on this. Basically, reduce $h$ (i.e. add more points) until your are happy with the amount of error you have compared to exact solution. Apr 27 '20 at 21:38
• I will check it, many thanks. Apr 27 '20 at 21:41
• @user62716 fyi, added manipulate to show effect of adding more points. Apr 27 '20 at 22:53
• Dear Nasser, many thanks for your great work and help, please what happen if the boundary conditions change to b) u(0)=0,u'(1)=3,c)u'(0)=0,u'(1)=3? Apr 28 '20 at 11:50

I got {1/9, 1/9, 1/3} which is wrong it should be 0, 1/27, 8/27, since the exact $$u=x^3$$.

As illustrated by Nasser, {1/9, 1/9, 1/3} is not wrong. If you want to obtain {0, 1/27, 8/27}, a possible solution (not sure if it's the only solution) is to use a higher order difference formula to approximate the b.c. at $$x=0$$. The formula you've used at $$x=0$$ is the first order forward difference formula:

$$f'(0)\approx \frac{f(h)-f(0)}{h}$$

If we use the third order one-sided difference formula instead:

$$f'(0)\approx-\frac{11 f(0)}{6 h}+\frac{3 f(h)}{h}-\frac{3 f (2 h)}{2 h}+\frac{f (3 h)}{3 h}$$

Solve@{-2 + 9 u - 18 u[1/3] + 9 u[2/3] == 0,
-4 + 9 u[1/3] - 18 u[2/3] + 9 u == 0,
-((11 u)/2) + 9 u[1/3] - 9/2 u[2/3] + u == 0,
-1 + u == 0}
(*
{{u -> 0, u[1/3] -> 1/27, u[2/3] -> 8/27, u -> 1}}
*)


BTW these difference equations can be generated easily with my pdetoae:

eq = u''[x] == 6 x;
bc = {u'  == 0, u  == 1};

grid = Range[0, 1, 1/3];
difforder = 3;

(* Definition of pdetoae isn't included in this post,
ptoafunc = pdetoae[u[x], grid, difforder];

ae = ptoafunc[eq][[2 ;; -2]]
aebc = ptoafunc@bc

Solve@Flatten@{ae, aebc}


If you want the $$A$$ and $$b$$, just use CoefficientArrays:

{barray, marray} = CoefficientArrays[Flatten@{ae, aebc}, u /@ grid];
MatrixForm /@ {barray, marray} LinearSolve[marray, -barray]
(* {0, 1/27, 8/27, 1} *)


With the help of pdetoae, dealing with b.c.s in (b) and (c) is trivial, so I'd like to omit the code here.

• Dear xzczd many thanks for your help. Apr 29 '20 at 15:02