Finite Differenctial Element PDE
8.1) Consider the 1D wave equation.
$$\frac{\partial^2 u}{\partial t^2} = v^2 \frac{\partial^2 u}{\partial x^2} $$8.1.a) Write Down the straight foward finite difference approximation.
from (8.5)
Place those in 1D wave equation.
8.1.b) What order approximation is this in time and in space?
Second order in both time(n-1, n, and n+1) and space(j-1, j, and j+1).
8.1.c) Use the von Neumann stability criterion to find the mode amplitudes.
8.1.d) Use this to find a condition on the velocity, time step, and space step for stability (hint: consider the product of the two amplitude solutions).
8.1.e) Do different modes decay at different rates for the stable case?
What is this asking?
8.1.f) Numerically solve the wave equation for the evolution from an initial condition with u = 0 except for one nonzero node, and verify the stability criterion.
import math
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
def update_line(num, data, line):
line.set_data(data[..., :num])
return line,
### Iniital Setups
j = 0
n = 0
dx = 0.2
dt = 0.2
v = 0.9
maxX = 50
t_frame = 500
gamma = 0.1
num_X = int(maxX / dx)
stableF = (v*dt/dx)**2 ## stable factor
print(stableF)
##u[j][n] = (A^n)*math.cos(j*k*dx)
## initialise List
x = []
t0 = []
t1 = []
t2 = []
tnew = []
history_t2_undamped = []
x_Count = int(maxX/dx)
print(x_Count)
for i in range(0, x_Count):
x.append(i)
t0.append(0)
t1.append(0)
t2.append(0)
tnew.append(0)
#t2[1] = 10
#t2[2] = 10
t2[5] = 1
#print(x)
print(t0)
print(t1)
print(t2)
#print(tnew)
# update list t
for n in range (0, t_frame):
for i in range(2,(len(t0)-1)):
term1 = t2[i+1]-2*t2[i]+t2[i-1]
#print (term1)
tnew[i] = stableF*(term1)+2*t2[i]-t1[i]
temp=[]
for j in range(0, len(t0)):
t0[j] = t1[j]
t1[j] = t2[j]
t2[j] = tnew[j]
temp.append(tnew[j])
history_t2_undamped.append(temp)
#print ("history_t2")
#print (history_t2)
## initialise List
x = []
t0 = []
t1 = []
t2 = []
tnew = []
history_t2_damped = []
x_Count = int(maxX/dx)
print(x_Count)
for i in range(0, x_Count):
x.append(i)
t0.append(0)
t1.append(0)
t2.append(0)
tnew.append(0)
t2[5] = 1
for n in range (0, t_frame):
for i in range(2,(len(t0)-1)):
term1 = t2[i+1]-2*t2[i]+t2[i-1]
term2 = t2[i+1]-2*t2[i]+t2[i-1]
tnew[i] = stableF*(term1)+ gamma*(term2)+2*t2[i]-t1[i]
temp=[]
for j in range(0, len(t0)):
t0[j] = t1[j]
t1[j] = t2[j]
t2[j] = tnew[j]
temp.append(tnew[j])
history_t2_damped.append(temp)
plt.plot(x, history_t2[3],'b',c='red',label="theta")
plt.xlabel("x")
plt.ylabel("y")
plt.show()
## simulation_undammped
from matplotlib import animation, rc
from IPython.display import HTML
fig, ax = plt.subplots()
plt.close()
ax.set_xlim(( 0, maxX))
ax.set_ylim((-1, 1))
line, = ax.plot([],[],lw=2)
def init():
line.set_data([], [])
return (line,)
# animation function. This is called sequentially
def animate(i):
X_Axis = np.linspace(0,maxX,x_Count)
Y_Axis = history_t2_undamped[i]
line.set_data(X_Axis, Y_Axis)
#print(i)
return (line,)
anim = animation.FuncAnimation(fig, animate, init_func=init,
frames=t_frame, interval=10, blit=True)
rc('animation', html='jshtml')
anim
## simulation_dammped
from matplotlib import animation, rc
from IPython.display import HTML
fig, ax = plt.subplots()
plt.close()
ax.set_xlim(( 0, maxX))
ax.set_ylim((-1, 1))
line, = ax.plot([],[],lw=2)
def init():
line.set_data([], [])
return (line,)
# animation function. This is called sequentially
def animate(i):
X_Axis = np.linspace(0,maxX,x_Count)
Y_Axis = history_t2_damped[i]
line.set_data(X_Axis, Y_Axis)
#print(i)
return (line,)
anim = animation.FuncAnimation(fig, animate, init_func=init,
frames=t_frame, interval=10, blit=True)
rc('animation', html='jshtml')
anim
8.1.f) If equation is replaced by
$$\frac{\partial^2 u}{\partial t^2} = v^2 \frac{\partial^2 u}{\partial x^2} + \gamma \frac{\partial}{\partial t} \frac{\partial^2 u}{\partial x^2}
$$
assume that
$$u(x,t) = A e^{i(kx-wt)}$$
and find a relationship between $k$ and $w$, and simply it for small $\gamma$. Comment on the relationship to the preceeding question.