Testing code to solve 2nd order wave equation

[TheCanadian]

As an exercise, I am trying to solve the 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]=0,t) = 0 ##

## E(z[1],t) = \sin(\pi z[1]) \cos(\pi t) ##

## E(z,t[0]=0) = sin(\frac{\pi z}{L}) ##

I know the analytic solution, which is:

## E(z,t) = \sin(\pi z) \cos(\pi t) ##

but am trying to solve it numerically. (The numerical and analytic solutions are compared to test accuracy.) Here is my very simple code, which applies a central difference approach:

    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""" #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_)
    EpM = 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)
     
        hz = z[k+1]-z[k] #hz is positive
     
        i = 1
        while i < (len(t) - 1):
            ht = t[i+1] - t[i]
             
            EpM[0,i] = 0.
            EpM[0,i+1] = 0.
    #        EpM[k,Ngridt-1] = (np.cos(np.pi*t[Ngridt-1]))*(np.sin(np.pi*z[k]))
            EpM[1,i] = (np.cos(np.pi*t[i]))*(np.sin(np.pi*z[1]))
    #        EpM[Ngridz-1,i] = (np.cos(np.pi*t[i]))*(np.sin(np.pi*z[Ngridz-1]))
         
            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]))
             
            EpM[0,i] = 0.
            EpM[0,i+1] = 0.
    #        EpM[k,Ngridt-1] = (np.cos(np.pi*t[Ngridt-1]))*(np.sin(np.pi*z[k]))
            EpM[1,i] = (np.cos(np.pi*t[i]))*(np.sin(np.pi*z[1]))
    #        EpM[Ngridz-1,i] = (np.cos(np.pi*t[i]))*(np.sin(np.pi*z[Ngridz-1]))      
         
            TEST[k,i] = np.sin(np.pi*z[k])*np.cos(np.pi*t[i])
         
            i = i + 1
         
        Ep[k,:] = EpM[k,:]
        Ep[k+1,:] = EpM[k+1,:]
       
        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,label='z = 0')
    plt.ylabel('Numerical E')
    plt.legend()
    plt.subplot(222)
    plt.plot(t,Ereal[int(Ngridz*0.33),:],'k:',linewidth = 1.5,label='z = 0.33*zmax')
    plt.legend()
    plt.subplot(223)
    plt.plot(t,Ereal[int(Ngridz*0.66),:],'k:',linewidth = 1.5,label='z = 0.66*zmax')
    plt.legend()
    plt.subplot(224)
    plt.plot(t,Ereal[int(Ngridz),:],'k:',linewidth = 1.5,label='z = zmax')
    plt.xlabel(r" t (s)")
    plt.legend()
    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,label='t = 0')
    plt.ylabel('Numerical E')
    plt.legend()
    plt.subplot(222)
    plt.plot(z,Ereal[:,int(Ngridt*0.33)],'k:',linewidth = 1.5,label='t = 0.33*tmax')
    plt.legend()
    plt.subplot(223)
    plt.plot(z,Ereal[:,int(Ngridt*0.66)],'k:',linewidth = 1.5,label='t = 0.66*tmax')
    plt.legend()
    plt.subplot(224)
    plt.plot(z,Ereal[:,Ngridt],'k:',linewidth = 1.5,label='t = tmax')
    plt.xlabel(r" z (m)")
    plt.legend()
    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,label='z = 0')
    plt.ylabel('True E')
    plt.legend()
    plt.subplot(222)
    plt.plot(t,TEST[int(Ngridz*0.33),:],'k:',linewidth = 1.5,label='z = 0.33*zmax')
    plt.legend()
    plt.subplot(223)
    plt.plot(t,TEST[int(Ngridz*0.66),:],'k:',linewidth = 1.5,label='z = 0.66*zmax')
    plt.legend()
    plt.subplot(224)
    plt.plot(t,TEST[int(Ngridz),:],'k:',linewidth = 1.5,label='z = zmax')
    plt.xlabel(r" t (s)")
    plt.legend()
    plt.savefig(str(newpath)+'/E.true_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,label='t = 0')
    plt.ylabel('True E')
    plt.legend()
    plt.subplot(222)
    plt.plot(z,TEST[:,int(Ngridt*0.33)],'k:',linewidth = 1.5,label='t = 0.33*tmax')
    plt.legend()
    plt.subplot(223)
    plt.plot(z,TEST[:,int(Ngridt*0.66)],'k:',linewidth = 1.5,label='t = 0.66*tmax')
    plt.legend()
    plt.subplot(224)
    plt.plot(z,TEST[:,Ngridt],'k:',linewidth = 1.5,label='t = tmax')
    plt.xlabel(r" z (m)")
    plt.legend()
    plt.savefig(str(newpath)+'/E.true_vs_z.pdf')

Here is the output of my code: http://imgur.com/a/z1geu

It seems as though my numerical result is converging to the analytic solution initially, but it appears to go to 0 afterwards. When I uncomment the boundary conditions for the end of my grid (which are commented out in the above code), the output shows the opposite behaviour where my numerical output is initially 0 but begins to converge to the analytic solution afterwards.

I've tried applying the boundary/initial conditions from the analytic solution to interior points in my numerical code, decreasing the step size, and trying higher order difference approximations, but none of these seem to work. I've continued looking through the code but don't seem to be able to find any mistakes. I've been going through it for a couple days now and the error seems very simple, but I'm missing it. Any ideas?

 

[lofgran]

As a fellow scientist, I understand the frustration of trying to solve a problem and not getting the desired results. After looking through your code, I have a few suggestions that may help you troubleshoot the issue:

1. Check your boundary conditions: It's possible that there is a mistake in how you are implementing the boundary conditions. Double check the equations and make sure they are applied correctly at the boundaries.

2. Try a different numerical method: While central difference is a commonly used approach, it may not be the most suitable for this specific problem. You could try using a different method, such as finite difference or finite element, and see if that produces better results.

3. Check your time step: Make sure your time step is small enough to accurately capture the dynamics of the system. If it is too large, it may lead to instability or inaccuracies in the solution.

4. Use a smaller domain: Instead of using a large domain, try using a smaller one and see if that produces better results. This could help pinpoint where the issue is occurring.

5. Debug your code: Use print statements or a debugger to check the values of your variables at different points in your code. This can help you identify where the error is occurring and fix it.

I hope these suggestions help you in solving your problem. Good luck!