2
$\begingroup$

I have solved 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 found u(0), u(1/3) and u(2/3)

I construct the linear system $A\,u = 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}.

I want to solve the equation with different types of boundary conditions. Please what happen if the boundary conditions change to

  1. u(0) = 0, u'(1) = 3,
  2. u'(0) = 0, u'(1) = 3

I need to construct A and b.

$\endgroup$
0

1 Answer 1

4
$\begingroup$

case 1. u(0) = 0, u'(1) = 3.

enter image description here

The first step is to solve for $u_{n+1}$. From equation (1) $$ u_{n+1}=2h\alpha+u_{n-1} $$ Substituting this into (2) gives the equation for the last point \begin{align} \frac{u_{n-1}-2u_{n}+u_{n+1}}{h^{2}} & =f_{n}\nonumber\\ \frac{u_{n-1}-2u_{n}+\left( 2h\alpha+u_{n-1}\right) }{h^{2}} & =f_{n}\nonumber\\ \frac{2u_{n-1}-2u_{n}+2h\alpha}{h^{2}} & =f_{n}\nonumber\\ 2u_{n-1}-2u_{n} & =h^{2}f_{n}-2h\alpha\tag{3} \end{align} Therefore, the equations are: For the first node $u_{1}=\beta$, for the second node at $i=2$ and using $\frac{u_{i-1}-2u_{i}+u_{i+1}}{h^{2}}=f_{i}$ gives \begin{align*} \frac{u_{1}-2u_{2}+u_{3}}{h^{2}} & =f_{2}\\ u_{1}-2u_{2}+u_{3} & =h^{2}f_{2} \end{align*} And for the third node \begin{align*} \frac{u_{2}-2u_{3}+u_{4}}{h^{2}} & =f_{3}\\ u_{2}-2u_{3}+u_{4} & =h^{2}f_{3} \end{align*} And so on until node $i=n$ which is (3) $$ 2u_{n-1}-2u_{n}=h^{2}f_{n}-2h\alpha $$ Putting these matrix form gives \begin{align*} \begin{pmatrix} 1 & 0 & 0 & \cdots & 0 & 0 & 0\\ 1 & -2 & 1 & 0 & \cdots & 0 & 0\\ 0 & 1 & -2 & 1 & 0 & \cdots & 0\\ 0 & 0 & 1 & -2 & 1 & 0 & 0\\ 0 & 0 & 0 & 1 & \ddots & 0 & 0\\ 0 & 0 & 0 & \cdots & 1 & -2 & 1\\ 0 & 0 & 0 & 0 & 0 & 2 & -2 \end{pmatrix} \begin{pmatrix} u_{1}\\ u_{2}\\ u_{3}\\ \vdots\\ u_{n-2}\\ u_{n-1}\\ u_{n} \end{pmatrix} & = \begin{pmatrix} \beta\\ h^{2}f\left( x_{2}\right) \\ h^{2}f\left( x_{3}\right) \\ \vdots\\ h^{2}f\left( x_{n-2}\right) \\ h^{2}f\left( x_{n-1}\right) \\ h^{2}f\left( x_{n}\right) -2h\alpha \end{pmatrix} \\ Au & =b \end{align*}

enter image description here

code

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, 1]] = 1;
      A[[1, 2]] = 0;
      A[[-1, -2]] = 2;
      A
      ];
makeB[n_, h_, force_, leftBC_, rightBC_] := Module[{b, i},
      b = Table[0, {i, n}];
      Do[

    b[[i]] = 
     If[i == 1, leftBC, 
      If[i < n, f[(i - 1)*h]*h^2, (h^2*f[(i - 1)*h] - 2 h rightBC) ]
            ]
        , {i, 1, n}
        ];
      b
      ];
f[x_] := 6*x;(*RHS of ode*)
Manipulate[
  Module[{h, A, b, sol, solN, p1, p2, x, leftBC, rightBC},
    h = 1/(nPoints - 1);
    leftBC = 0;
    rightBC = 3;
    A = makeA[nPoints];
    b = makeB[nPoints, h, f, leftBC, rightBC];
    sol = LinearSolve[A, b];
    solN = Table[{n*h, sol[[n + 1]]}, {n, 0, nPoints - 1}];

    p1 = Plot[x^3, {x, 0, 1}, AxesOrigin -> {0, 0}]; (*exact solution*)

    p2 = ListLinePlot[solN, PlotStyle -> Red, Mesh -> All];

    Grid[{
            {Row[{" h = ", NumberForm[N@h, {5, 4}]}]},
            {Row[{MatrixForm[A], MatrixForm[Array[u, nPoints]] , 
                        " = ", MatrixForm[N@b] }]},
            {Show[p1, p2,                      
                 PlotLabel -> "Red is numerical, Blue is exact solution",
                 GridLines -> Automatic,
                 GridLinesStyle -> LightGray, ImageSize -> 400
                 ]
             }
          }, Frame -> All, Spacings -> {1, 2}
       ]
    ],
  {{nPoints, 3, "How many points?"}, 3, 8, 1, Appearance -> "Labeled"},
  TrackedSymbols :> {nPoints}
  ]

case 2 u'(0) = 0, u'(1) = 3

enter image description here

The first step is to solve for $u_{n+1}$. From equation (1) $$ u_{n+1}=2h\alpha+u_{n-1} $$ Substituting this into (2) gives the equation for the last point \begin{align} \frac{u_{n-1}-2u_{n}+u_{n+1}}{h^{2}} & =f_{n}\nonumber\\ \frac{u_{n-1}-2u_{n}+\left( 2h\alpha+u_{n-1}\right) }{h^{2}} & =f_{n}\nonumber\\ \frac{2u_{n-1}-2u_{n}+2h\alpha}{h^{2}} & =f_{n}\nonumber\\ 2u_{n-1}-2u_{n} & =h^{2}f_{n}-2h\alpha\tag{3} \end{align}

Similary we solve for $u_{0}$. From equation (3) $$ u_{0}=u_{2}-2h\beta $$ Substituting this into (4) gives the equation for the first point \begin{align} \frac{u_{0}-2u_{1}+u_{2}}{h^{2}} & =f_{0}\nonumber\\ \frac{\left( u_{2}-2h\beta\right) -2u_{1}+u_{2}}{h^{2}} & =f_{0}\nonumber\\ \frac{2u_{2}-2h\beta-2u_{1}}{h^{2}} & =f_{0}\nonumber\\ 2u_{2}-2u_{1} & =f_{0}h^{2}+2h\beta\tag{3} \end{align}

Therefore, the equations are: For the first node $$ 2u_{2}-2u_{1}=f_{0}h^{2}+2h\beta $$ For the second node at $i=2$ and using $\frac{u_{i-1}-2u_{i}+u_{i+1}}{h^{2} }=f_{i}$ gives \begin{align*} \frac{u_{1}-2u_{2}+u_{3}}{h^{2}} & =f_{2}\\ u_{1}-2u_{2}+u_{3} & =h^{2}f_{2} \end{align*} And for the third node \begin{align*} \frac{u_{2}-2u_{3}+u_{4}}{h^{2}} & =f_{3}\\ u_{2}-2u_{3}+u_{4} & =h^{2}f_{3} \end{align*} And so on until node $i=n$ which is (3) $$ 2u_{n-1}-2u_{n}=h^{2}f_{n}-2h\alpha $$ Putting these matrix form gives \begin{align*} \begin{pmatrix} -2 & 2 & 0 & \cdots & 0 & 0 & 0\\ 1 & -2 & 1 & 0 & \cdots & 0 & 0\\ 0 & 1 & -2 & 1 & 0 & \cdots & 0\\ 0 & 0 & 1 & -2 & 1 & 0 & 0\\ 0 & 0 & 0 & 1 & \ddots & 0 & 0\\ 0 & 0 & 0 & \cdots & 1 & -2 & 1\\ 0 & 0 & 0 & 0 & 0 & 2 & -2 \end{pmatrix} \begin{pmatrix} u_{1}\\ u_{2}\\ u_{3}\\ \vdots\\ u_{n-2}\\ u_{n-1}\\ u_{n} \end{pmatrix} & = \begin{pmatrix} f_{0}h^{2}+2h\beta\\ h^{2}f\left( x_{2}\right) \\ h^{2}f\left( x_{3}\right) \\ \vdots\\ h^{2}f\left( x_{n-2}\right) \\ h^{2}f\left( x_{n-1}\right) \\ h^{2}f\left( x_{n}\right) -2h\alpha \end{pmatrix} \\ Au & =b \end{align*}

The analytical solution for $u^{\prime\prime}\left( x\right) =6x$ with $u^{\prime}\left( 0\right) =0,u^{\prime}\left( 1\right) =3$ is not unique. It is $x^{3}+C$. The constant $C$ is arbitrary and an infinite number of solutions exist. A solution exist which is up to an arbitrary additive constant. To select a constant for the purpose of the numerical analysis, the constant is found to give the analytical solution a zero mean which is done by solving

\begin{align*} \int_{0}^{1}\left( x^{3}+C\right) dx & =0\\ \left[ \frac{x^{4}}{4}+Cx\right] _{0}^{1} & =0\\ \frac{1}{4}+C & =0\\ C & =-\frac{1}{4} \end{align*}

Hence the solution $u\left( x\right) =x^{3}-\frac{1}{4}$ is used.

enter image description here

code

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[[-1, -2]] = 2;
      A
      ];
makeB[n_, h_, force_, leftBC_, rightBC_] := Module[{b, i},
      b = Table[0, {i, n}];
      Do[              
         b[[i]] = If[i == 1, f[0]*h^2 + 2*h*leftBC, 
           If[i < n, f[(i - 1)*h]*h^2, (h^2*f[(i - 1)*h] - 2 h rightBC) ]
              ]
        , {i, 1, n}
        ];
      b
      ];
f[x_] := 6*x;(*RHS of ode*)
Manipulate[
  Module[{h, A, b, sol, solN, p1, p2, x, leftBC, rightBC, 
   normalizationConstant},
    h       = 1/(nPoints - 1);
    leftBC  = 0;
    rightBC = 3;
    A       = makeA[nPoints];
    b       = makeB[nPoints, h, f, leftBC, rightBC];
    sol     = LinearSolve[A, b];
    solN    = Table[{n*h, sol[[n + 1]]}, {n, 0, nPoints - 1}];
    normalizationConstant = -1/4;
    solN[[All, 2]] =        solN[[All, 2]] - Mean[solN[[All, 2]]]; (*To match normalization *)

    p1 = Plot[x^3 + normalizationConstant, {x, 0, 1}, 
    AxesOrigin -> {0, 0}, 
    PlotRange -> {Automatic, {-.4, 1}}]; (*exact solution*)
    p2 = ListLinePlot[solN, PlotStyle -> Red, Mesh -> All];
    Grid[{
            {Row[{" h = ", NumberForm[N@h, {5, 4}]}]},
            {Row[{MatrixForm[A], MatrixForm[Array[u, nPoints]] , 
               " = ", MatrixForm[N@b] }]},
            {Show[p1, p2, 
                PlotLabel -> "Red is numerical, Blue is exact solution",
                GridLines -> Automatic,
                GridLinesStyle -> LightGray, ImageSize -> 400
                ]
             }
          }, Frame -> All, Spacings -> {1, 2}
       ]
    ],
  {{nPoints, 3, "How many points?"}, 3, 20, 1, 
  Appearance -> "Labeled"},
  TrackedSymbols :> {nPoints}
  ]

Update, answering comment

Hard to answer this in comment as I have to show large code.

solve u''=2, u'(0)=0,u'(1)=2, the exact now x^2+c, c=-1/3,.. the error is big, why

The error is not big at all. as you add more nodes it goes down as expected. you must made mistake somewhere changing the above code for the new ode you are now asking about. You have to change f and change boundary conditions and change the normalization constant. Here is the Manipulate for the above ODE

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[[-1, -2]] = 2;
      A
      ];
makeB[n_, h_, force_, leftBC_, rightBC_] := Module[{b, i},
      b = Table[0, {i, n}];
      Do[                     b[[i]] = 
        If[i == 1, f[0]*h^2 + 2*h*leftBC, 
        If[i < n, f[(i - 1)*h]*h^2, (h^2*f[(i - 1)*h] - 2 h rightBC) ]
            ]
        , {i, 1, n}
        ];
      b
      ];
f[x_] := 2;(*RHS of ode*)
Manipulate[
  Module[{h, A, b, sol, solN, p1, p2, x, leftBC, rightBC, 
   normalizationConstant},
    h = 1/(nPoints - 1);
    leftBC = 0;
    rightBC = 2;
    A = makeA[nPoints];
    b = makeB[nPoints, h, f, leftBC, rightBC];
    sol = LinearSolve[A, b];
    solN = Table[{n*h, sol[[n + 1]]}, {n, 0, nPoints - 1}];
    normalizationConstant = -1/3;
    solN[[All, 2]] = solN[[All, 2]] - Mean[solN[[All, 2]]]; (*To match normalization *)

    p1 = Plot[x^2 + normalizationConstant, {x, 0, 1}, 
       AxesOrigin -> {0, 0}, 
       PlotRange -> {Automatic, {-.4, 1}}]; (*exact solution*)
       p2 = ListLinePlot[solN, PlotStyle -> Red, Mesh -> All];
    Grid[{
            {Row[{" h = ", NumberForm[N@h, {5, 4}]}]},
            {Row[{MatrixForm[A], MatrixForm[Array[u, nPoints]] , 
                  " = ", MatrixForm[N@b] }]},
            {Show[p1, p2, 

      PlotLabel -> "Red is numerical, Blue is exact solution",
                 GridLines -> Automatic,
                 GridLinesStyle -> LightGray, ImageSize -> 400
                 ]
             }
          }, Frame -> All, Spacings -> {1, 2}
       ]
    ],
  {{nPoints, 3, "How many points?"}, 3, 20, 1, 
  Appearance -> "Labeled"},
  TrackedSymbols :> {nPoints}
  ]

If you run this, you'll see the error is small and gets smaller with larger nodes.

$\endgroup$
15
  • $\begingroup$ Thank you very much Dr.Nasser, I highly appreciate your help. Just for other case which is Neumann B.C. I think it is not uniquely solvable due to the derivative, i.e. if u=x^3 solution then x^3+c also solution right? $\endgroup$
    – user62716
    Commented Apr 29, 2020 at 12:24
  • $\begingroup$ @user62716 Yes, The analytical solution for case 2 is not unique. The constant $c$ is arbitrary and an infinite number of solutions exist. A solution exist which is up to an arbitrary additive constant. To select a constant for the purpose of the numerical analysis, the constant is selected to give the analytical solution a zero mean. i.e. solve $\int_0^1 (x^3+c) \, dx=0$ for $c$. This gives $c=-\frac{1}{4}$ and use this $c$ to compare the numerical solution with the analytical solution now by fixing $c$. If it is still not clear I can add case 2 also later on. $\endgroup$
    – Nasser
    Commented Apr 29, 2020 at 17:31
  • $\begingroup$ Dear Dr. Nasser, many thanks for your help and top clarification, I will be happy to see Neumann case 2, but I am sorry I took too much time from you. $\endgroup$
    – user62716
    Commented Apr 29, 2020 at 21:17
  • $\begingroup$ Many thanks Dr.Nasser for your great help today I highly appreciate that. Just minor correction in your great well-written part in both case in third term in vector b, I thing it is h^2 f(x3) not h^2 f(x2). Best wishes and regards, I will move to 2D Laplace or maybe 3D with different types of boundary conditions. $\endgroup$
    – user62716
    Commented Apr 29, 2020 at 22:41
  • $\begingroup$ @user62716 thanks. There were latex cut/paste error. Corrected. $\endgroup$
    – Nasser
    Commented Apr 29, 2020 at 22:43

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