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/sand 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 matplotlibmatplotlib.use('pdf')import os import matplotlib.pyplot as pltimport numpy as npfrom tqdm import tqdmfrom mpl_toolkits.mplot3d import Axes3Dfrom matplotlib import cmimport matplotlib.pyplot as pltfrom matplotlib.colors import Normalizec = 1.#3.*(10.**8.)#z spaceNgridz = 400 #number of intervals in z-axiszmax = L = 5.dz = zmax/Ngridz z = np.linspace(0.0,zmax,Ngridz+1) #timeNgridt = 400 #number of intervals in t-axistmax = L/cdt = tmax/Ngridtt = 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[1]-t[0]#    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 kTTdef dt2(X,k,i,z,t,kX,TYPE):     Approximation for second time derivative - CENTRAL DIFFERENCE    ht = t[1]-t[0]    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 kTTEp = 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 bartotal = 100.progressz = (total)/(len(z))k = 0while 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).realnewpath = 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/GwPCtOutput of backward difference approximation: http://imgur.com/a/66jemAs 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.
Using backward vs central finite difference approximation
finite difference;python;wave propagation
Higher order methods often have a smaller radius of convergence, i.e., they require smaller time steps. In your context, this means that they require a smaller CFL number, often significantly smaller than 1. Did you take this into account?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.To put things differently: If you increase your wave speed by a factor of 10, then everything happens ten times as fast, and it only makes sense that you then have to choose a time step ten times smaller than before.