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

Question

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),'--')

Answer 1

This is only a guess, I haven't run your code.

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.