Spectral Method for Wave Equation using Finite Differences

Clear["Global`*"]

Variable coefficient wave equation (Trefethen’s example p6.m)

Solve numerically the differential equation

((∂u(x,t))/(∂t))=c(x)((∂u(x,t))/(∂x)),              c(x) =(1/5)+(sin^2)(x-1)

for x ∈ [0,2 π], t > 0 with the following initial condition

u(x,0)=(e^(-100(((x-1))^2)))

and periodic  boundary conditions

u(0,t)=u(2π,t)

Note that the initial function is not periodic but it is so close to zero at the ends of the interval that it can be regarded as periodic in practice.

(* Grid, variable coefficient, and initial data *)
n=128;
h=2.0*Pi/n;
x=Table[h*i,{i,1,n}];
t=0;
dt=h/40;
c=0.2+Sin[x-1]^2;
v=Exp[-100*(x-1)^2];
(* For the leap-frog method below, we need another initial function for time -dt *)
vold=Exp[-100*(x-0.2*dt-1)^2]; (* c(x) at x = 1 is close to -1/5 *)
ListPlot[{Transpose[{x,c}],Transpose[{x,v}],Transpose[{x,vold}]}, PlotStyle->{PointSize[0.015]},Joined->True,Mesh->All,PlotLegends->{"c(x)","u(x,0)","u(x,-dt)"}]

(* Setting times for plotting and time step accordingly *)
tmax=8;
tplot=0.1;
plotgap=Round[tplot/dt];
dt=tplot/plotgap;
nplots=Round[tmax/tplot];
data=Join[{v},Table[ConstantArray[0,n],{nplots}]];
tdata={t};

(* Determine the differentiation matrix for finite differences *)
getSparseFDMatrix[n_,p_,h_]:=Module[{v,pf2,matrix},
v=ConstantArray[0,n];
pf2=Factorial[p]^2;
Do[
v[[s+1]]=(-1)^(s+1)*pf2/(s*Factorial[p-s]*Factorial[p+s]);
v[[n-s+1]]=-v[[s+1]],
{s,1,p}
];
matrix=SparseArray[Flatten[{Table[Band[{1,i}]->v[[i]],{i,1,n}],
Table[Band[{i,1}]->v[[n+2-i]],{i,2,n}]}],{n,n}];
Return[matrix/h]];

(* Here, we use 5-point central difference method, you can try to change p to increase accuracy *)
DM=getSparseFDMatrix[n,2,h];

(* Time evolution of the wave *)
Timing[For[i=1,i≤nplots,i++,
For[j=1,j≤plotgap,j++,t=t+dt;
(* Time-stepping by the leap-frog formula *)
vnew=DM.v;
vnew=vold-2.0*dt*c*vnew;
vold=v;
v=vnew;
];
data[[i+1]]=v;
AppendTo[tdata,t];
];
]

ListPlot[Table[Transpose[{x,data[[it,All]]}],{it,1,81,10}],PlotRange->All,Joined->True,AxesLabel->{"x","u(x,t)"},PlotLabel->"Wave function at various times starting on the left"]

ListPlot3D[Flatten[Table[{x[[j]],tdata[[i]],data[[i,j]]},{i,1,Length[tdata]},{j,1,Length[x]}],1],PlotRange->{{0,2*Pi},{0,tmax},{-1,5}},AxesLabel->{"x","t","u"},Mesh->{nplots}]