PDF (letter size)

## Notes on some numerical schemes

May 26, 2022   Compiled on January 29, 2024 at 3:04am

### 1 Introduction

These are notes to describe some numerical schemes.

#### 1.1 Centered diﬀerence for 1D wave PDE

Here we want to solve $$u_{tt}=c^{2}u_{xx}$$ in 1D ﬁnite domain $$0<x<L$$ and $$t>0$$, with boundary conditions $$u\left ( 0,t\right ) =f\left ( t\right ) ,u\left ( L,t\right ) =g\left ( t\right )$$ and initial position $$u\left ( x,0\right ) =\alpha \left ( x\right )$$ and initial velocity $$\frac {\partial u\left ( x,0\right ) }{\partial t}=\beta \left ( x\right )$$.

Centered diﬀerence is used.

At each internal node $$j$$ we have the following ﬁnite diﬀerence represenation of $$u_{tt}=c^{2}u_{xx}$$

To handle initial conditions, initial velocity is used to solve for $$u_{j}^{-1}$$

This gives all the information needed to ﬁnd the matrices to use. Let $$k=\Delta t$$. From Eq(1) \begin {align*} \frac {u_{j}^{1}-u_{j}^{-1}}{2k} & =\alpha \\ u_{j}^{-1} & =u_{j}^{1}-2k\alpha \end {align*}

Substituting this in Eq(2) gives\begin {align*} \frac {\left ( u_{j}^{1}-2k\alpha \right ) -2u_{j}^{0}+u_{j}^{1}}{k^{2}} & =c^{2}\frac {u_{j-1}^{0}-2u_{j}^{0}+u_{j+1}^{0}}{h^{2}}\\ 2u_{j}^{1} & =\frac {k^{2}c^{2}}{h^{2}}\left ( u_{j-1}^{0}-2u_{j}^{0}+u_{j+1}^{0}\right ) +2u_{j}^{0}+2k\alpha \\ u_{j}^{1} & =\frac {1}{2}\frac {k^{2}c^{2}}{h^{2}}\left ( u_{j-1}^{0}-2u_{j}^{0}+u_{j+1}^{0}\right ) +u_{j}^{0}+k\alpha \end {align*}

Therefore for $$n=1$$ only and for $$j=1\cdots N$$ where $$N$$ is number of nodes$\begin {pmatrix} u_{1}^{1}\\ u_{2}^{1}\\ u_{3}^{1}\\ u_{4}^{1}\\ u_{5}^{1}\end {pmatrix} =\frac {1}{2}\frac {k^{2}c^{2}}{h^{2}}\begin {pmatrix} 1 & 0 & 0 & 0 & 0\\ 1 & -2 & 1 & 0 & 0\\ 0 & 1 & -2 & 1 & 0\\ 0 & 0 & 1 & -2 & 1\\ 0 & 0 & 0 & 0 & 1 \end {pmatrix}\begin {pmatrix} u_{1}^{0}\\ u_{2}^{0}\\ u_{3}^{0}\\ u_{4}^{0}\\ u_{5}^{0}\end {pmatrix} +\begin {pmatrix} u_{1}^{0}\\ u_{2}^{0}\\ u_{3}^{0}\\ u_{4}^{0}\\ u_{5}^{0}\end {pmatrix} +k\alpha$ Where $$\begin {pmatrix} u_{1}^{0}\\ u_{2}^{0}\\ u_{3}^{0}\\ u_{4}^{0}\\ u_{5}^{0}\end {pmatrix}$$ is known and comes from boundary and initial conditions. $$u_{1}^{0}$$ is left B.C. and $$u_{N}^{0}$$ comes from right B.C. and $$u_{2}^{0}\cdots u_{N-1}^{0}$$ comes from initial conditions $$u\left ( x,0\right )$$. Now, for $$n=2$$ or higher times\begin {align*} \frac {u_{j}^{n-1}-2u_{j}^{n}+u_{j}^{n+1}}{k^{2}} & =c^{2}\frac {u_{j-1}^{n}-2u_{j}^{n}+u_{j+1}^{n}}{h^{2}}\\ u_{j}^{n-1}-2u_{j}^{n}+u_{j}^{n+1} & =\frac {k^{2}c^{2}}{h^{2}}\left ( u_{j-1}^{n}-2u_{j}^{n}+u_{j+1}^{n}\right ) \\ u_{j}^{n+1} & =\frac {k^{2}c^{2}}{h^{2}}\left ( u_{j-1}^{n}-2u_{j}^{n}+u_{j+1}^{n}\right ) +2u_{j}^{n}-u_{j}^{n-1} \end {align*}

In Matrix form$\begin {pmatrix} u_{1}^{n+1}\\ u_{2}^{n+1}\\ u_{3}^{n+1}\\ u_{4}^{n+1}\\ u_{5}^{n+1}\end {pmatrix} =\frac {k^{2}c^{2}}{h^{2}}\begin {pmatrix} 1 & 0 & 0 & 0 & 0\\ 1 & -2 & 1 & 0 & 0\\ 0 & 1 & -2 & 1 & 0\\ 0 & 0 & 1 & -2 & 1\\ 0 & 0 & 0 & 0 & 1 \end {pmatrix}\begin {pmatrix} u_{1}^{n}\\ u_{2}^{n}\\ u_{3}^{n}\\ u_{4}^{n}\\ u_{5}^{n}\end {pmatrix} +2\begin {pmatrix} u_{1}^{n}\\ u_{2}^{n}\\ u_{3}^{n}\\ u_{4}^{n}\\ u_{5}^{n}\end {pmatrix} -\begin {pmatrix} u_{1}^{n-1}\\ u_{2}^{n-1}\\ u_{3}^{n-1}\\ u_{4}^{n-1}\\ u_{5}^{n-1}\end {pmatrix}$ So to ﬁnd $$u_{j}^{n+1}$$ we need to know the last time step solution and also the solution for the step before that.

Small Mathematica was written to implement the above scheme. Here is screen show of the GUI for one example.

The Mathematica notebook of the above is here

The following is an animation of $$u_{tt}=4u_{xx}$$ with ﬁxed ends, and string length $$L=5$$ with initial position given by $$8x\frac {\left ( L-x\right ) ^{2}}{L^{3}}$$ and zero initial velocity. This uses $$\Delta t=0.02$$ seconds and $$30$$ grid points over the length $$5$$.

 Pause Play Restart Step forward Step back

The following is listing of the source code

Finite difference for 1D wave PDE
By Nasser M. Abbasi. Version May 8, 2020

L = 5;
leftBC[x_, t_] := 0;(*t^2;*)
rightBC[x_, t_] := 0;(*t^2+25;*)
initialPosition[x_] := 8 x*(5 - x)^2/5^3; (*x^2;*)
initialVelocity = 0;

AccountingForm[v, f, NumberSigns -> {"-", " "},
(*these 2 functions thanks to xzczd*)
numberForm[a_List, n_] := numberForm[#, n] & /@ a;

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, -1]] = 1;
A[[-1, -2]] = 0;
A];

makeInitialU[nPoints_, L_, h_, leftBC_, rightBC_, initialPosition_] :=
Module[{u, j, t = 0},
u = Table[0, {j, nPoints}];
Do[
u[[j]] =
If[j == 1, leftBC[0, 0],
If[j == nPoints, rightBC[L, 0], initialPosition[(j - 1)*h]]],
{j, 1, nPoints}
];
u
];

makePlot[currentTime_, showMMA_, grid_, currentU_, u_, opt_, opt1_,
yRangeMin_, yRangeMax_, solN_, showMatrix_, k_, c_, h_, A_,
initialVelocity_] := Module[{},

Grid[{
{Dynamic@If[showMMA,
Show[
ListLinePlot[Transpose[{grid, u}], Evaluate[opt]],
Plot[solN[x, currentTime], {x, 0, 5}, Evaluate[opt1]],
PlotRange -> {{0, 5}, {-yRangeMin, yRangeMax}}
],
ListLinePlot[Transpose[{grid, u}],
Evaluate@
Join[opt, {PlotRange -> {{0, 5}, {-yRangeMin, yRangeMax}}}]
]
]
},
{Dynamic@If[showMatrix,
Row[{"U = ", NumberForm[k^2*c^2/2*h^2], " ", MatrixForm[A],
" . ", MatrixForm[numberForm[u, {5, 4}]], " + ",
MatrixForm[numberForm[u, {5, 4}]],
If[initialVelocity != 0, Row[{" + ", k*initialVelocity}]],
" = ", MatrixForm[numberForm[currentU, {5, 4}]]}]
,
"No matrix display"
]}
}, Spacings -> {1, 1}, Frame -> True]
];

DynamicModule[{solN, lastU, currentU, currentTime = 0, A, h,
showMatrix = False,
showMMA = False, k = 0.02, nPoints = 30, maxTime = 4,
yRangeMax = 1.2, yRangeMin = 1.2,
opt, opt1, pde, ic, bc, grid, g = 0, u, x, t, nextU, c = 4,
state = "STOP", tick = False},

opt = {PlotStyle -> Red, AxesOrigin -> {0, 0}, Mesh -> All,
MeshStyle -> {Blue, PointSize[0.01]},
ImageSize -> 400, ImagePadding -> 10, ImageMargins -> 10};
opt1 = {PlotStyle -> Blue, AxesOrigin -> {0, 0}, ImageSize -> 400,
ImagePadding -> 10, ImageMargins -> 10};

Dynamic[
tick;
If[currentTime == 0,
A = makeA[nPoints];
h = L/(nPoints - 1);
lastU =
N@makeInitialU[nPoints, L, h, leftBC, rightBC, initialPosition];
currentU =
0.5 (c^2*k^2)/h^2*(A . lastU) + lastU + (k*initialVelocity);
currentU[[1]] = leftBC[0, k];
currentU[[-1]] = rightBC[L, k];
pde = D[u[x, t], {t, 2}] == c ^2 D[u[x, t], {x, 2}];
ic = {u[x, 0] == initialPosition[x],
Derivative[0, 1][u][x, 0] == initialVelocity};
bc = {u[0, t] == leftBC[0, t], u[L, t] == rightBC[L, 0]};
solN =
Quiet@NDSolveValue[{pde, ic, bc}, u, {x, 0, 5}, {t, 0, maxTime}];
grid = Range[0, L, h];
g = makePlot[currentTime, showMMA, grid, currentU, lastU, opt,
opt1, yRangeMin, yRangeMax, solN, showMatrix, k, c, h, A,
initialVelocity];
If[state == "RUN" || state == "STEP",
If[(currentTime + k) <= maxTime,
currentTime = currentTime + k
,
state == "STOP"
]
]
,
If[state != "STOP",
nextU = (c^2*k^2)/h^2*A . currentU + 2 currentU - lastU;
nextU[[1]] = leftBC[0, currentTime];
nextU[[-1]] = rightBC[L, currentTime];

g = makePlot[currentTime, showMMA, grid, currentU, nextU, opt,
opt1, yRangeMin, yRangeMax, solN, showMatrix, k, c, h, A,
initialVelocity];

If[state == "RUN" || state == "STEP",
If[(currentTime + k) <= maxTime,
currentTime = currentTime + k
]
];
If[state == "STEP", state = "STOP"];
lastU = currentU;
currentU = nextU
]
];

Row[{Grid[{
{Row[{Button[
Text@Style["run", 12], {currentTime = 0; state = "RUN"},
ImageSize -> {60, 40}],

Button[Text@Style["stop", 12], {state = "STOP"},
ImageSize -> {60, 40}],

Button[Text@Style["step", 12], {state = "STEP"},
ImageSize -> {60, 40}],

Button[Text@Style["reset", 12], {currentTime = 0;
state = "STOP"}, ImageSize -> {60, 40}]}]
},
{Row[{"Show matrix", Spacer[3],
Checkbox[
Dynamic[showMatrix, {showMatrix = #;
tick = Not[tick]} &]]}]},
{Row[{"Show Mathematica solution", Spacer[3],
Checkbox[
Dynamic[showMMA, {showMMA = #; tick = Not[tick]} &]]}]},
{Row[{"Number of grid points? ",
Manipulator[
Dynamic[nPoints, {nPoints = #; currentTime = 0;
state = "STOP"} &], {3, 50, 1}, ImageSize -> Tiny],
Dynamic[nPoints]}]},
{Row[{"Wave speed (c) ? ",
Manipulator[
Dynamic[c, {c = #; currentTime = 0;
state = "STOP"} &], {0.01, 5, 0.01}, ImageSize -> Tiny],
Dynamic[c]}]},
{Row[{"Time step? (delT) ? ",
Manipulator[
Dynamic[k, {k = #; currentTime = 0;
state = "STOP"} &], {0.001, 0.05, 0.01},
ImageSize -> Tiny], Dynamic[k]}]},
{Row[{"max run time ?",
Manipulator[
Dynamic[maxTime, {maxTime = #; currentTime = 0;
state = "STOP"} &], {0, 5, 0.01}, ImageSize -> Tiny],
Dynamic[maxTime]}]},
{Row[{"yRangeMax ?",
Manipulator[
Dynamic[yRangeMax, {yRangeMax = #; tick = Not[tick]} &], {1,
30, 0.01}, ImageSize -> Small], Dynamic[yRangeMax]}]},
{Row[{"yRangeMin ?",
Manipulator[
Dynamic[yRangeMin, {yRangeMin = #; tick = Not[tick]} &], {1,
30, 0.01}, ImageSize -> Small], Dynamic[yRangeMin]}]}
}, Alignment -> Left, Spacings -> {1, 1}, Frame -> All
], g}
]
,
ContinuousAction -> False,
TrackedSymbols :> {currentTime, state, tick}
]

]