1
$\begingroup$

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[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 kTT

def 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 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: https://i.stack.imgur.com/NianM.jpg

Output of backward difference approximation: https://i.stack.imgur.com/Lvmxv.jpg

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.

Improve this question
$\endgroup$
2
  • $\begingroup$ Which approximation do you use for the spatial derivative ? $\endgroup$
    – Coriolis
    Oct 6, 2016 at 9:17
  • $\begingroup$ @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/… $\endgroup$
    – Mathews24
    Oct 6, 2016 at 19:12

1 Answer 1

Reset to default
6
$\begingroup$

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.

Improve this answer
$\endgroup$
3
  • $\begingroup$ Thank you for the explanation. That clears things up a lot. I am new to using numerical methods and was not aware that higher order methods have a smaller radius of converge (actually, I was under the impression the opposite was true). Does radius of convergence mean anything beside step-size? Also, would you happen to have any references for reading about the radius of convergence of different methods (e.g. Runge-Kutta), and in particular, on the backward difference approximations? $\endgroup$
    – Mathews24
    Oct 7, 2016 at 1:09
  • $\begingroup$ Wow, I'm still so surprised at how little of that I knew. I earlier read Gear's method would be helpful as it was supposed to be stable, but didn't realize this conveyed nothing about its accuracy. $\endgroup$
    – Mathews24
    Oct 7, 2016 at 1:17
  • $\begingroup$ Radius of convergence is certainly related to the maximal time step. It's been a long time since I learned about that, but I'm pretty sure the standard book about ODE solvers, Hairer and Wanner, will have a good description. $\endgroup$
    – Wolfgang Bangerth
    Oct 8, 2016 at 11:34

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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