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;

    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 ) 

semilogy( L2,'-')
  • $\begingroup$ Is this a homework question? $\endgroup$
    – Kirill
    Oct 24 '14 at 20:26
  • $\begingroup$ Not a homework question $\endgroup$
    – user107904
    Oct 24 '14 at 20:54
  • $\begingroup$ 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. $\endgroup$ Oct 24 '14 at 21:33
  • $\begingroup$ Hmm, when I change the L2 error to the relative L2 error I get that the order is approximately 1. $\endgroup$
    – user107904
    Oct 24 '14 at 23:18

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.

  • $\begingroup$ It looks like this worked. Thanks. Changing my spatial grid to x = x(1:end-1) seemed to do the trick. $\endgroup$
    – user107904
    Oct 24 '14 at 23:24

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.