5
$\begingroup$

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.

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

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}]

Mathematica graphics

sol = u[x] /. First@DSolve[{u''[x] == 6*x, u'[0] == 0, u[1] == 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]

Mathematica graphics

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

enter image description here

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}
 ]
$\endgroup$
5
  • $\begingroup$ which value of h can provide good results, I reduce it but still not good results? $\endgroup$ – user62716 Apr 27 '20 at 21:32
  • $\begingroup$ 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. $\endgroup$ – Nasser Apr 27 '20 at 21:38
  • $\begingroup$ I will check it, many thanks. $\endgroup$ – user62716 Apr 27 '20 at 21:41
  • $\begingroup$ @user62716 fyi, added manipulate to show effect of adding more points. $\endgroup$ – Nasser Apr 27 '20 at 22:53
  • $\begingroup$ 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? $\endgroup$ – user62716 Apr 28 '20 at 11:50
2
$\begingroup$

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[0] - 18 u[1/3] + 9 u[2/3] == 0, 
       -4 + 9 u[1/3] - 18 u[2/3] + 9 u[1] == 0, 
       -((11 u[0])/2) + 9 u[1/3] - 9/2 u[2/3] + u[1] == 0, 
       -1 + u[1] == 0}
(*
{{u[0] -> 0, u[1/3] -> 1/27, u[2/3] -> 8/27, u[1] -> 1}}
 *)

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

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

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

(* Definition of pdetoae isn't included in this post,
   please find it in the link above. *)
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}

enter image description here

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.

$\endgroup$
1
  • $\begingroup$ Dear xzczd many thanks for your help. $\endgroup$ – user62716 Apr 29 '20 at 15:02

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.