Finite difference for 1D wave PDE

By Nasser M. Abbasi. Version May 8, 2020

Example 1

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

Example 2

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

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 ", padIt1[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]
];

Main module

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+2currentU-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}
]

]