# Using backward vs central finite difference approximation

I am solving the simple 2nd-order wave equation:

$$\frac {\partial ^2 E}{\partial t^2} = c^2 \frac {\partial ^2 E}{\partial z^2}$$

Over a domain of (in SI units):

$z = [0,L=5]$m, $t = [0,t_{max} = 5]$s, $c = 1$ m/s

and boundary/initial conditions:

$$E(z=0) = E(z=L) = 0$$

$$E(t=0) = -E(t=t_{max}) = sin(\frac{\pi z}{L})$$

I know the analytic solution, but am trying to solve it numerically. (The numerical and analytic solutions are compared to test accuracy.) I just had a couple questions:

1) When I solve the wave equation by applying a central difference approximation of order 2 on the second derivative of time, the code works perfectly fine. Although when I apply backward finite difference approximations of higher orders for the second derivative in time, my solution diverges. Is there any particular reason why using a backward difference approximation would be worse than a central difference approximation of the same order? I was under the impression using a higher order backward difference method would give me higher accuracy, but it appears to not work at all. I can supply the code if needed, although it's a fairly basic implementation.

2) I am able to solve the equation above when using a central difference approximation and setting $c = 1$ m/s, but if I use $c = 3 \times 10^8$ m/s, the solution diverges. This makes sense as the Courant number ($C_o = c \frac {\Delta t}{\Delta x}$) is much greater than one. But I am curious: are there any possible finite difference schemes that can be applied to numerically solve the equation if $C_o >> 1$? I have read a bit about BDF methods to help solve stiff equations, but in this case, the only problem is that a more precise grid is needed, correct? I had also heard about spectral methods, but am still trying to learn how to implement them and if they would be applicable/accurate enough if other (nonlinear) terms were included in the wave equation.

Edit:

After further testing, it appears the output is actually slightly off even without the backward differences. The code:

import matplotlib
matplotlib.use('pdf')
import os
import matplotlib.pyplot as plt
import numpy as np
from tqdm import tqdm
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
import matplotlib.pyplot as plt
from matplotlib.colors import Normalize

c = 1.#3.*(10.**8.)

#z space
Ngridz = 400 #number of intervals in z-axis
zmax = L = 5.
dz = zmax/Ngridz
z = np.linspace(0.0,zmax,Ngridz+1)

#time
Ngridt = 400 #number of intervals in t-axis
tmax = L/c
dt = tmax/Ngridt
t = np.linspace(0.0,tmax,Ngridt+1)

#def dt2(X,k,i,z,t,kX,TYPE):
#    """Approximation for second time derivative""" - BACKWARD DIFFERENCE
#    ht = t-t
#    if TYPE == 'CONJ':
#        if i == 0:
#            kTT = np.conj(X[k,i]-2.*X[k,i-1]+X[k,i-2])/(ht**2.)
#        elif i == 1:
#            kTT = np.conj(X[k,i-1]-2.*X[k,i]-X[k,i+1])/(ht**2.)
#        elif i == 2:
#            kTT = np.conj(-(X[k,i])+2.*(X[k,i-1])-(X[k,i-2]))/(ht**2.)
#        elif i == 3:
#            kTT = np.conj(-2.*(X[k,i])+5.*(X[k,i-1])-
#            4.*(X[k,i-2])+(X[k,i-3]))/(ht**2.)
#        elif i == 4:
#            kTT = np.conj((-35/12.)*(X[k,i])+(26./3.)*(X[k,i-1])-
#            (19./2.)*(X[k,i-2])+(14./3.)*(X[k,i-3])-
#            (11./12.)*(X[k,i-4]))/(ht**2.)
#        elif i == 5:
#            kTT = np.conj((-15./4.)*(X[k,i])+(77./6.)*(X[k,i-1])-
#            (107./6.)*(X[k,i-2])+13.*(X[k,i-3])-
#            (61./12.)*(X[k,i-4])+(5./6.)*(X[k,i-5]))/(ht**2.)
#        else:
#            kTT = np.conj((-203./45.)*(X[k,i])+(87./5.)*(X[k,i-1])-
#            (117./4.)*(X[k,i-2])+(254./9.)*(X[k,i-3])-
#            (33./2.)*(X[k,i-4])+(27./5.)*(X[k,i-5])-
#            (137./180.)*(X[k,i-6]))/(ht**2.)
#    else:
#        if i == 0:
#            kTT = (X[k,i]-2.*X[k,i-1]+X[k,i-2])/(ht**2.)
#        elif i == 1:
#            kTT = (X[k,i-1]-2.*X[k,i]-X[k,i+1])/(ht**2.)
#        elif i == 2:
#            kTT = (-(X[k,i])+2.*(X[k,i-1])-(X[k,i-2]))/(ht**2.)
#        elif i == 3:
#            kTT = (-2.*(X[k,i])+5.*(X[k,i-1])-
#            4.*(X[k,i-2])+(X[k,i-3]))/(ht**2.)
#        elif i == 4:
#            kTT = ((-35/12.)*(X[k,i])+(26./3.)*(X[k,i-1])-
#            (19./2.)*(X[k,i-2])+(14./3.)*(X[k,i-3])-
#            (11./12.)*(X[k,i-4]))/(ht**2.)
#        elif i == 5:
#            kTT = ((-15./4.)*(X[k,i])+(77./6.)*(X[k,i-1])-
#            (107./6.)*(X[k,i-2])+13.*(X[k,i-3])-
#            (61./12.)*(X[k,i-4])+(5./6.)*(X[k,i-5]))/(ht**2.)
#        else:
#            kTT = ((-203./45.)*(X[k,i])+(87./5.)*(X[k,i-1])-
#            (117./4.)*(X[k,i-2])+(254./9.)*(X[k,i-3])-
#            (33./2.)*(X[k,i-4])+(27./5.)*(X[k,i-5])-
#            (137./180.)*(X[k,i-6]))/(ht**2.)
#    return kTT

def dt2(X,k,i,z,t,kX,TYPE):
"""Approximation for second time derivative""" - CENTRAL DIFFERENCE
ht = t-t
if TYPE == 'CONJ':
if i == 0:
kTT = np.conj(X[k,i]-2.*X[k,i+1]+X[k,i+2])/(ht**2.)
else:
kTT = np.conj(X[k,i-1]-2.*X[k,i]-X[k,i+1])/(ht**2.)
else:
if i == 0:
kTT = (X[k,i]-2.*X[k,i+1]+X[k,i+2])/(ht**2.)
else:
kTT = (X[k,i-1]-2.*X[k,i]+X[k,i+1])/(ht**2.)
return kTT

Ep = np.zeros((len(z),len(t)),dtype=np.complex_)
TEST = np.zeros((len(z),len(t)),dtype=np.complex_)

progress = tqdm(total=100.) #this provides a progress bar
total = 100.
progressz = (total)/(len(z))

k = 0
while k < (len(z) - 1):
progress.update(progressz)

EpM = Ep

hz = z[k+1]-z[k] #hz is positive

i = 0
while i < (len(t) - 1):
ht = t[i+1] - t[i]

EpM[0,i] = 0. #at z = 0
#        EpM[Ngridz-1,:] = 0.
EpM[k,0] = np.sin(np.pi*z[k])
#        EpM[k+1,0] = np.sin(np.pi*z[k+1])
EpM[k,Ngridt-1] = -np.sin(np.pi*z[k])

EpM[k+1,i] = (-EpM[k,i-1] + 2.*EpM[k,i] + ((hz/(c))**2.)*dt2(EpM,k,i,z,t,0.,'x'))
#((hz/(c*ht))**2.)*(EpM[k,i+1] - 2.*EpM[k,i] + EpM[k,i-1]))

#        print 'k='+str(k)+',i='+str(i)+',dE='+str(EpM[k+1,i]-EpM[k,i])

EpM[0,i] = 0. #at z = 0
#        EpM[Ngridz-1,:] = 0.
EpM[k,0] = np.sin(np.pi*z[k])
#        EpM[k+1,0] = np.sin(np.pi*z[k+1])
EpM[k,Ngridt-1] = -np.sin(np.pi*z[k])

TEST[k,i] = np.sin(np.pi*z[k])*np.cos(np.pi*t[i])

i = i + 1

Ep[k,:] = EpM[k,:]

k = k + 1

if k == (len(z)-1):
progress.update(progressz)

Ereal = (Ep).real

newpath = r'test_wave_equation'
if not os.path.exists(newpath):
os.makedirs(newpath)

plt.figure(1)
fig, ax = plt.subplots(figsize=(20, 20))
plt.subplot(221)
plt.plot(t,Ereal[0,:],'k:',linewidth = 1.5)
plt.ylabel('Normalized E.real')
plt.subplot(222)
plt.plot(t,Ereal[int(Ngridz*0.33),:],'k:',linewidth = 1.5)
plt.subplot(223)
plt.plot(t,Ereal[int(Ngridz*0.66),:],'k:',linewidth = 1.5)
plt.subplot(224)
plt.plot(t,Ereal[int(Ngridz),:],'k:',linewidth = 1.5)
plt.xlabel(r" t (s)")
plt.savefig(str(newpath)+'/E.real_vs_t.pdf')
#plt.show()

plt.figure(2)
fig, ax = plt.subplots(figsize=(20, 20))
plt.subplot(221)
plt.plot(z,Ereal[:,0],'k:',linewidth = 1.5)
plt.ylabel('Normalized E.real')
plt.subplot(222)
plt.plot(z,Ereal[:,int(Ngridt*0.33)],'k:',linewidth = 1.5)
plt.subplot(223)
plt.plot(z,Ereal[:,int(Ngridt*0.66)],'k:',linewidth = 1.5)
plt.subplot(224)
plt.plot(z,Ereal[:,Ngridt],'k:',linewidth = 1.5)
plt.xlabel(r" z (m)")
plt.savefig(str(newpath)+'/E.real_vs_z.pdf')

plt.figure(3)
fig, ax = plt.subplots(figsize=(20, 20))
plt.subplot(221)
plt.plot(t,TEST[0,:],'k:',linewidth = 1.5)
plt.ylabel('True E')
plt.subplot(222)
plt.plot(t,TEST[int(Ngridz*0.33),:],'k:',linewidth = 1.5)
plt.subplot(223)
plt.plot(t,TEST[int(Ngridz*0.66),:],'k:',linewidth = 1.5)
plt.subplot(224)
plt.plot(t,TEST[int(Ngridz),:],'k:',linewidth = 1.5)
plt.xlabel(r" t (s)")
plt.savefig(str(newpath)+'/E.real_vs_t.pdf')
#plt.show()

plt.figure(4)
fig, ax = plt.subplots(figsize=(20, 20))
plt.subplot(221)
plt.plot(z,TEST[:,0],'k:',linewidth = 1.5)
plt.ylabel('Normalized E.real')
plt.subplot(222)
plt.plot(z,TEST[:,int(Ngridt*0.33)],'k:',linewidth = 1.5)
plt.subplot(223)
plt.plot(z,TEST[:,int(Ngridt*0.66)],'k:',linewidth = 1.5)
plt.subplot(224)
plt.plot(z,TEST[:,Ngridt],'k:',linewidth = 1.5)
plt.xlabel(r" z (m)")
plt.savefig(str(newpath)+'/E.real_vs_z.pdf')


Output: http://imgur.com/a/GwPCt

Output of backward difference approximation: http://imgur.com/a/66jem

As you may see, my output does not seem to match the analytic solution, and this discrepancy only appears to be present initially when using the central difference approximation. As the time- and z-components advance, it appears my solution does converge to the right answer. I have tried increasing the number of grid points (i.e. decreasing interval size) and debugging but I don't seem to see any obvious errors. Any ideas on what may be causing the problem? I have commented out the function for the backward difference approximation of the second time derivative, but when I uncomment it and use it instead of the central difference function for dt2, the solution diverges and the original problem still remains.

• Which approximation do you use for the spatial derivative ? – Coriolis Oct 6 '16 at 9:17
• @Coriolis I have edited the question to include the code and output. The first link shows E vs t and E vs z; one is the analytic solution plotted, while the other is my numerical solution plotted when using central differences. The second link is the output plotted when I use backward difference approximations for the second time derivative. The coefficients for the second derivative with accuracy up to 5 using backward differences (which I used) are here: en.wikipedia.org/wiki/… – Mathews24 Oct 6 '16 at 19:12

As for the question of whether there are methods that can deal with CFL numbers much larger than one: This really doesn't make any sense. The CFL number is defined as the ratio of the time step divided by the time it takes information to cross one cell. If you wanted to have CFL $\gg$ 1, that would mean that within one time step, information travels far more than one cell. Because within this time interval, the signal may interact with boundaries or with itself, you can not expect that you can accurately resolve the solution if your CFL number is large. In other words, there are of course methods that remain stable, but they will not be accurate if the CFL number is much greater than one.