TheCanadian
- 361
- 13
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:
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?
$$ \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:
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""" #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?