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[[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}
]
]
Notebook download
It might be hard to copy all the code above. Here is link to the notebook itself which might make it easier.
