# Finite difference method for 1D wave equation

I want to solve the following 1D wave equation:

utt = uxx with t > 0, 0 <= x <= 5 and ic = u(x, 0) = x^2 and du/dt(x, 0) = 0 and u(0, t) = t^2, u(5, t) = t^2 + 25.

ClearAll["Global*"]

heqn = D[u[x, t], {t, 2}] == D[u[x, t], {x, 2}];
ic = {u[x, 0] == (x^2), Derivative[0, 1][u][x, 0] == 0};
bc = {u[0, t] == t^2, u[5, t] == t^2 + 25};

sol = DSolve[{heqn, ic, bc }, u[x, t], {x, t}]}]


I got the exact solution u = t^2 + x^2.

Now I implement the finite difference method:

utt = ui^(n + 1) - 2 ui^n + ui^(n - 1)/delta t;
uxx = ui + 1^n - 2 ui^n + ui - 1^(n - 1)/delta x,


then I got:

ui^(n+1) = -ui^(n - 1) + 2 ui^n + r^2[ui + 1^n - 2 ui^n + ui - 1^n]      (* 1 *)


where r = delta t/delta x;

Now using (1), I want to construct A, b and get u at x = 1, 2, 3, 4; t = 0.25, 0.5, 0.75, 1.

I do not know how to get the required values numerically.

## Introduction

This uses implicit finite difference method. Using standard centered difference scheme for both time and space. To make it more general, this solves $$u_{tt} = c^2 u_{xx}$$ for any initial and boundary conditions and any wave speed $$c$$. It also shows the Mathematica solution (in blue) to compare against the FDM solution in red (with the dots on it).

The more grid point is used, the more accurate the solution becomes. Will first show couple of demos. The first for the problem in the question

$$u_{xx} = u_{tt}$$ with initial position of the string $$u(x,0)=x^2$$ and zero initial velocity. The boundary conditions $$u(0,t)=t^2,u(5,t)=t^2+25$$. Normally string is fixed on both ends. Running this for 1 second, with 6 grid points, using $$\Delta t=0.01$$ shows Second example is a fixed string on both ends with higher wave speed. $$u_{tt} = 4 u_{xx}$$ with fixed left and right ends, and initial position $$u(x,0)= 8 x+(L-x)^2/L^3$$ where $$L=5$$ is the length. The length is always fixed at 5 in this version. Zero initial velocity also. Mathematica's solution above is more accurate because the time step used in FDM is large $$0.03$$ and only $$13$$ points are used. Making the time step smaller makes it more accurate but will take longer to run.

## Short description of the scheme used

Centered difference is used. As follows To handle initial conditions, initial velocity is used to solve for $$u^{-1}_j$$ This gives all the information needed to find 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 find $$u_{j}^{n+1}$$ we need to know the last time step solution and also the solution for the step before that.

It is the above $$A$$ and solution vectors which is displayed below the plot.

## Code

Edit these lines to change initial and boundary conditions. These are for example 1 above

L = 5;
leftBC[x_, t_] := t^2;
rightBC[x_, t_] := t^2 + 25;
initialPosition[x_] := x^2;
initialVelocity := 0;


These are for example 2 above (the fixed string)

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;


These are helper functions

padIt1[v_, f_List] :=
AccountingForm[v, f, NumberSigns -> {"-", "+"},
NumberPadding -> {"0", "0"}, SignPadding -> True];
(*these 2 functions thanks to xzczd*)
numberForm[a_List, n_] := numberForm[#, n] & /@ a;
numberForm[a_, n_] := padIt1[a, n];

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[{
{Row[{"time ", NumberForm[Dynamic@currentTime, {4, 2}]}]},
{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]
];


This is the DynamicModule

DynamicModule[{solN, lastU, currentU, currentTime = 0, A, h,
showMatrix = True,
showMMA = True, k = 0.01, nPoints = 6, maxTime = 1, yRangeMax = 30,
yRangeMin = 2,
opt, opt1, pde, ic, bc, grid, g = 0, u, x, t, nextU, c = 1,
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[] = 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[] = 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,
Checkbox[
Dynamic[showMatrix, {showMatrix = #;
tick = Not[tick]} &]]}]},
{Row[{"Show Mathematica solution", Spacer,
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}
]

]


## Notebook download

It might be hard to copy all the code above. Here is link to the notebook itself which might make it easier.

• Wow! An impressively comprehensive answer. – MarcoB May 8 '20 at 13:50
• Dear Dr.Nasser, many thanks for your perfect answer, please is it possible to send me the math file, I got the code but I want the above steps, if possible...Best regards. – user62716 May 8 '20 at 14:56
• @user62716 Please find the pdf file here – Nasser May 8 '20 at 19:04
• Dear Dr.Nasser, many thanks for your great and perfect work....Best regards – user62716 May 8 '20 at 19:14

This question should be brilliantly answered by Mathematica documentation. Have a close look at:

The Numerical Method of Lines

This is an introduction to Mathematica NDSolve'FiniteDifferenceDerivative` and has several examples starting with the heat equation and the asked 1D wave equation of this very question. It shows how to do the tables and lists and matrices and how to use Mathematic for that.

For example:  and all that is needed to solve the task numerically. The documentation goes even further to give formulas for error calculation.

This is an example how straight forward that is: f is general. Order is 2 and too general, h can be the finite difference both in time and space.

• Dear user2432923, thank you for your response, I need to construct the matrix A and vector b to get the solution not use the builtin functions.....Best regards – user62716 May 6 '20 at 19:53