# 2nd Order finite difference for 1D wave equation matlab issue

I'm trying verify that a 2nd order finite difference in space and time approximation of the 1D wave equation is really 2nd order. My Matlab implementation tells me otherwise - I'm not sure of what I've done incorrectly.

For simplicity I've decided to use the smooth solution $u(x,t)=\cos{(ct)}\cos{x}$. The PDE has periodic boundary conditions:

$$u_{tt} = c^2u_{xx},~~~u(x,0)=\cos{x},~~~u_t(x,0)=0,~~~u(-2\pi,t)=u(2\pi,t).$$ I've taken $c=1$, and the domain to be $[-2\pi,2\pi]$.

The finite difference scheme I'm using is centered second order in space and time: $$\frac{u(x_j,t_{i+1})-2u(x_j,t_{i})+u(x_j,t_{i-1})}{(\Delta t)^2} = c^2\frac{u(x_{j+1},t_{i})-2u(x_{j},t_{i})+u(x_{j-1},t_{i})}{(\Delta x)^2}.$$

To approximate the initial condition $u_t(x,0)=0$, I tried both the first order forward approximation $$\frac{u(x_j,t_{j+1})-u(x_j,t_{j})}{\Delta t}=0,$$ and the second order centered approximation $$\frac{u(x_j,t_{j+1})-u(x_j,t_{j-1})}{2\Delta t}=0.$$ When I fix $\Delta t$ and refine $\Delta x$, both of these approximations to $u_t(x,0)=0$ seem to give an order of approximately 1.5 according to my Matlab code. The largest error looks like it occurs near the boundaries $x=\pm 2\pi$. Can anyone see what I'm doing incorrectly?

a = -2*pi;
b =  2*pi;

NN = 2.^(3:9);
for jj = 1 : length(NN)
n = NN(jj);               %Number of grid points
dx = (b-a)/(n-1);         %Spatial mesh spacing
x = a : dx : b;           %Mesh in x direction
c = 1;
u = @(x,t) cos(c*t).*cos(x);
u0 = u(x,0.0)';

e = ones(n,1);
A = spdiags([e -2*e e], -1:1, n, n);   %Laplacian matrix
A(1,end) = 1;  A(end,1) = 1;           %Periodic boundary conditions

T = 1;                                 %Final simulation time
dt = 1e-3;                             %Temporal mesh spacing
nsteps = (floor(T/dt));                %Number if time steps

%u_old = u0;  %First order approxiation of u_t(x,0)=0.
%u_now = u0;

%Second order approxiation of u_t(x,0)=0.
u_old = u0;
u_now = u0 + 0.5*(dt*c/dx)^2*(A* u0);  %(approximates u(x,dt))

for ii = 2 : nsteps  %First time step was directly above
u_new = 2*u_now - u_old + (dt*c/dx)^2*(A* u_now);
u_old = u_now;
u_now = u_new;
end

LI(jj) =  norm(u_new - u(x,T)',inf);
L2(jj) = dx*norm(u_new - u(x,T)',2);
if ( jj > 1 )  %Look at order of convergence
log ( L2(jj-1)/L2(jj) ) / log ( 2 )
end
end

figure(1)
semilogy( L2,'-')
figure(2)
plot(x,u_new,'ro-',x,u(x,T),'--')

• Is this a homework question? – Kirill Oct 24 '14 at 20:26
• Not a homework question – user107904 Oct 24 '14 at 20:54
• Just a guess, but when you calculate the errors since you're comparing vectors of different lengths you have to divide by the norm of the exact solution. I know in my experience not doing so usually leads to a reduction of 0.5 in the rate of convergence. – Lukas Bystricky Oct 24 '14 at 21:33
• Hmm, when I change the L2 error to the relative L2 error I get that the order is approximately 1. – user107904 Oct 24 '14 at 23:18

The boundary condition $y(-2\pi)=y(2\pi)$ means the points $\pm2\pi$ are identified. So the points to the left and right of $2\pi$ are $2\pi-h$ and $-2\pi+h$. However, your mesh has both (non-identical) points $\pm2\pi$ included in it, as separate points, which looks completely suspicious to me. Then the point to the right of $2\pi$ is $-2\pi$, which can't be right as it's the same point.