In [5]:
%matplotlib inline
In [4]:
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import animation, rc
from IPython.display import HTML
(a) Consider the 1D wave equation $\quad \frac{\partial^2 u}{\partial t^2} = v^2 \frac{\partial^2 u}{\partial x^2}$¶
$ \text{First we construct the finite difference approximation for a second derivative}\\ f^{''} = \lim_{h\to\infty}\frac{f^{'}(x)-f^{'}(x-h)}{h}\quad \text{where:}\quad f^{'}(x) = \frac{f(x+h)-f(x)}{h}\quad f^{'}(x-h) = \frac{f(x)-f(x-h)}{h}\\ \therefore \: f^{''} = \lim_{h\to\infty}\frac{1}{h}\frac{(f(x+h)-f(x))-(f(x)-f(x-h))}{h} \implies f^{''} = \lim_{h\to\infty}\frac{(f(x+h)-2f(x)+f(x-h)}{h^2}\\ \text{thus the finite difference approximation of:} \quad \frac{\partial^2 u}{\partial t^2} = v^2 \frac{\partial^2 u}{\partial x^2} \quad \text{can be written as}\\ \frac{u(t+\Delta t)-2u(t)+u(t-\Delta t)}{\Delta t^2} = v^2 \frac{u(x+\Delta x)-2u(x)+u(x-\Delta x)}{\Delta x^2} \quad $
(b) What order approximation is this in time and in space?¶
This approximation is correct to first order in both dimensions
(c) Use the von Neumann stability criterion to find the mode amplitudes.¶
$ \frac{1}{\Delta t^2}(U_{j}^{n+1}-2U_{j}^{n}+U_{j}^{n-1}) = \frac{v^2}{\Delta x^2}(U_{j+1}^{n}-2U_{j}^{n}+U_{j-1}^{n}) \quad \text{using the ansatz:} \quad A^n e^{\:ijk\Delta x} \quad\\ \frac{e^{\: ijk\Delta x}}{\Delta t^2}(A-2+\frac{1}{A}) = \frac{v^2}{\Delta x^2}(e^{\: i(j+1)k\Delta x}-2e^{\: ijk\Delta x}+ e^{\: i(j-1)k\Delta x}) \implies A-2+\frac{1}{A} = \frac{v^2 \Delta t^2}{\Delta x^2}(2 cos(k \Delta x) + 1)(e^{\: ik\Delta x}+1+e^{\: -ik\Delta x})\\ A-2+\frac{1}{A} = \frac{v^2 \Delta t^2}{\Delta x^2}(2 cos(k \Delta x) + 1) \quad \text{define:}\quad u = \frac{v^2 \Delta t^2}{\Delta x^2}(2 cos(k \Delta x) + 1)\\ A^2 - 2A +1 = uA \implies = A^2-(u+2)A+1= 0 \\ A = \frac{1}{2} (-(u+2)\pm \sqrt{u^2+4u} ) \implies \frac{1}{2} (-(u+2)\pm \sqrt{u^2+4u} ) <= 1\\ \therefore -u \pm sqrt{u^2+4u} >= 4 \implies u <= -4\\ \frac{v^2 \Delta t^2}{\Delta x^2}(2 cos(k \Delta x) + 1) <= -4 \quad \text{min} \quad cos(x) = -1\\ $
In [3]:
def wave(dt,dx,v,steps,length):
u = np.zeros(shape=(steps+1,length))
u[0,25] = 1
u[1,25] = 1
s = (((v**2)*(dt**2))/(dx**2))
for i in range(1,steps-1):
u[i+1,1:-1] = 2*(u[i,1:-1])-(u[i-1,1:-1])+s*(u[i,:-2]-2*u[i,1:-1]+u[i,2:])
return(u)
In [4]:
def wave_damp(dt,dx,v,gamma,steps,length):
u = np.zeros(shape=(steps+1,length))
u[0,25] = 1
u[1,25] = 1
s = ((((v**2)*(dt**2))+(dt*gamma))/(dx**2))
for i in range(1,steps-1):
#u[i+1,1:-1] = 2*(u[i,1:-1])-(u[i-1,1:-1])+(((v**2)*(dt**2))/(dx**2))*(u[i,:-2]-2*u[i,1:-1]+u[i,2:])+((gamma*dt)/(dx**2))*(u[i-1,:-2]-2*u[i-1,1:-1]+u[i-1,2:])
u[i+1,1:-1] = 2*(u[i,1:-1])-(u[i-1,1:-1])+s*(u[i,:-2]-2*u[i,1:-1]+u[i,2:])
return(u)
In [5]:
def run(data):
x = np.linspace(-1, 1, 51)
y1, y2 = data
line[0].set_data(x, y1)
line[1].set_data(x, y2)
return line
In [6]:
def data():
tstep = 0
while tstep < 1000:
tstep+=1
y1 = u1[tstep,:]
y2 = u2[tstep,:]
yield y1, y2
In [7]:
u2 = wave(0.004,0.00401,1,1000,51)
u1 = wave(0.004,0.004,1,1000,51)
# wave(dt,dx,v,timestep,numpoints)
In [8]:
fig1, ((ax1, ax2)) = plt.subplots(nrows=2, ncols=1,)
line1, = ax1.plot([], [], lw=4, color = 'black')
line2, = ax2.plot([], [], lw=4, color = 'orange')
line = [line1, line2]
ax1.set_ylim(-1, 1)
ax1.set_xlim(-1, 1)
ax2.set_ylim(-1, 1)
ax2.set_xlim(-1, 1)
ax1.set_title("ratio of $v^2 \Delta t^2$ to $\Delta x^2 = 1.0$")
ax2.set_title("ratio of $v^2 \Delta t^2$ to $\Delta x^2 = 0.997$")
plt.tight_layout()
plt.close()
anim = animation.FuncAnimation(fig1, run, data, blit=True, interval=50, repeat=True)
HTML(anim.to_html5_video())
Out[8]:
In [9]:
u3 = wave_damp(1,1,1,0.0001,1000,51)
u4 = wave_damp(1,1,1,0.0001,1000,51)
# wave(dt,dx,v,timestep,numpoints)
In [10]:
def data2():
tstep = 0
while tstep < 10000:
tstep+=1
y1 = u3[tstep,:]
y2 = u4[tstep,:]
yield y1, y2
In [11]:
fig, ((ax1, ax2)) = plt.subplots(nrows=2, ncols=1,)
line1, = ax1.plot([], [], lw=4, color = 'black')
line2, = ax2.plot([], [], lw=4, color = 'orange')
line = [line1, line2]
ax1.set_ylim(-1, 1)
ax1.set_xlim(-1, 1)
ax2.set_ylim(-1, 1)
ax2.set_xlim(-1, 1)
ax1.set_title("ratio of $v^2 \Delta t^2$ to $\Delta x^2 = 1.0$")
ax2.set_title("ratio of $v^2 \Delta t^2$ to $\Delta x^2 = 0.997$")
plt.tight_layout()
plt.close()
anim = animation.FuncAnimation(fig, run, data2, blit=True, interval=50,repeat=True)
HTML(anim.to_html5_video())
Out[11]:
In [314]:
def diffuse(D,dx,dt,steps,length,start):
u = np.zeros(shape=(steps+1,length))
u[0,start] = 1
u[1,start] = 1
for i in range(1,steps-1):
u[i+1,1:-1] = (u[i,1:-1])+((D*dt)/(dx**2))*(u[i,:-2]-2*u[i,1:-1]+u[i,2:])
return(u)
In [317]:
def rundiff(data):
x = np.linspace(-1, 1, 500)
y1, y2, y3 = data
line[0].set_data(x, y1)
line[1].set_data(x, y2)
line[2].set_data(x, y3)
return line
In [318]:
def data3():
tstep = 0
while tstep < 1000:
tstep+=1
y1 = d1[tstep,:]
y2 = d2[tstep,:]
y3 = d3[tstep,:]
yield y1, y2, y3
In [315]:
d1 = diffuse(0.1,1,1,1000,500,250)
d2 = diffuse(0.1,1,0.5,1000,500,250)
d3 = diffuse(0.1,1,0.1,1000,500,250)
In [322]:
fig, (ax1, ax2, ax3) = plt.subplots(1,3)
line1, = ax1.plot([], [], lw=4, color = 'black')
line2, = ax2.plot([], [], lw=4, color = 'orange')
line3, = ax3.plot([], [], lw=4, color = 'red')
line = [line1, line2, line3]
ax1.set_ylim(-1, 1)
ax1.set_xlim(-0.1, 0.1)
ax2.set_ylim(-1, 1)
ax2.set_xlim(-0.1, 0.1)
ax3.set_ylim(-1, 1)
ax3.set_xlim(-0.1, 0.1)
st = fig.suptitle("Explicit method",fontsize=14)
ax1.set_title("$\Delta t = 1$")
ax2.set_title("$\Delta t = 0.5$")
ax3.set_title("$\Delta t = 0.5$")
plt.tight_layout()
st.set_y(0.95)
fig.subplots_adjust(top=0.85)
plt.close()
anim = animation.FuncAnimation(fig, rundiff, data3, blit=True, interval=100,repeat=True)
HTML(anim.to_html5_video())
Out[322]:
In [ ]:
def triagonal(n,a,b,c):
A = a*np.roll(np.eye(n),1,axis = 0)
A[0,:] = 0
C = c*np.roll(np.eye(n),-1,axis = 0)
C[n-1,:] = 0
B = b*np.eye(n)
return A+B+C
In [279]:
def thomas(n,alpha,y):
# n,a,b,c are integers and d is a vector of the initial values
x = np.zeros((n,1))
b = (1+2*alpha)*np.ones((n,1))
a = -alpha*np.ones((n,1))
c = -alpha*np.ones((n,1))
c[0,0] = 0
a[n-1,0] = 0
x = np.zeros((n,1))
# forward pass
c[0] = c[0]/b[0]
y[0] = y[0]/b[0]
for k in range(1,n):
c[k] = c[k]/(b[k]-a[k]*c[k-1])
y[k] = (y[k]-a[k]*y[k-1])/(b[k]-a[k]*c[k-1])
# backward pass
x[n-1] = y[n-1]
for k in range(n-2,-1,-1):
x[k] = y[k] - c[k]*x[k+1]
return x
In [345]:
def runthomas(n,alpha):
x = np.zeros((1000,n))
y = np.zeros((n,1))
y[int(np.floor((n-1)/2))] = 1
i = 0
while i < 1000:
temp = thomas(n,alpha,y)
y = temp
x[i,:] = np.transpose(temp)
i += 1
return x
In [350]:
n = 500
alpha1 = 0.1
alpha2 = 0.1*0.5
alpha3 = 0.1*0.1
d1 = runthomas(n,alpha1)
d2 = runthomas(n,alpha2)
d3 = runthomas(n,alpha3)
In [351]:
fig, (ax1, ax2, ax3) = plt.subplots(1,3)
line1, = ax1.plot([], [], lw=4, color = 'black')
line2, = ax2.plot([], [], lw=4, color = 'orange')
line3, = ax3.plot([], [], lw=4, color = 'red')
line = [line1, line2, line3]
ax1.set_ylim(-1, 1)
ax1.set_xlim(-0.1, 0.1)
ax2.set_ylim(-1, 1)
ax2.set_xlim(-0.1, 0.1)
ax3.set_ylim(-1, 1)
ax3.set_xlim(-0.1, 0.1)
st = fig.suptitle("Implicit method", fontsize=14)
ax1.set_title("$\Delta t = 1$")
ax2.set_title("$\Delta t = 0.5$")
ax3.set_title("$\Delta t = 0.5$")
plt.tight_layout()
st.set_y(0.95)
fig.subplots_adjust(top=0.85)
plt.close()
anim = animation.FuncAnimation(fig, rundiff, data3, blit=True, interval=100,repeat=True)
HTML(anim.to_html5_video())
Out[351]:
In [ ]:
fig = plt.figure
plt.contourf(X,Y,Z,100,alpha=1,cmap=plt.get_cmap('RdBu'))
plt.xlim(-.5,0.5)
plt.ylim(-.5,0.5)
plt.colorbar()
In [ ]:
d1 = diffuse(0.1,1,1,1000,500,200)
In [ ]:
x = y = np.arange(-1, 1, 2/500)
X, Y = np.meshgrid(x, y)
a, b = np.meshgrid(d1[i,:],d1[i,:])
Z = a*b
In [19]:
d = diffuse(0.1,1,1,1000,500,250)
fig = plt.figure()
plt.xlim(220,280)
plt.ylim(220,280)
def f(j):
x, y = np.meshgrid(d[j,:],d[j,:])
return x * y
j = 0
image = []
for i in range(100):
j += 5
im = plt.imshow(f(j), animated=True,alpha=1,cmap=plt.get_cmap('RdBu'))
image.append([im])
anim = animation.ArtistAnimation(fig, image, interval=50, blit=True, repeat = True)
plt.close()
HTML(anim.to_html5_video())
Out[19]:
In [34]:
def diff(dt,dx,D,steps,length,start):
u = np.zeros(shape=(steps+1,length))
u[0,start] = 1
u[1,start] = 1
for i in range(1,steps-1):
u[i+1,1:-1] = (u[i,1:-1])+((D*dt)/(dx**2))*(u[i,:-2]-2*u[i,1:-1]+u[i,2:])
return(u)
d1 = diff(0.1,1,1,1000,500,266)
d2 = diff(0.1,1,1,1000,500,252)
d3 = diff(0.1,1,1,1000,500,215)
d4 = diff(0.1,1,1,1000,500,233)
fig = plt.figure()
plt.xlim(220,280)
plt.ylim(220,280)
def f(i):
x1, y1 = np.meshgrid(d1[j,:],d1[j,:])
x2, y2 = np.meshgrid(d2[j,:],d2[j,:])
x3, y3 = np.meshgrid(d3[j,:],d3[j,:])
x4, y4 = np.meshgrid(d4[j,:],d4[j,:])
z1 = x1 * y1
z2 = x2 * y2
z3 = x3 * y3
z4 = x4 * y4
return z1+z2+z3+z4
j = 0
image = []
for i in range(100):
j += 5
im = plt.imshow(f(j), animated=True,alpha=1,cmap=plt.get_cmap('RdBu'))
image.append([im])
anim = animation.ArtistAnimation(fig, image, interval=50, blit=True, repeat = True)
plt.close()
HTML(anim.to_html5_video())
Out[34]:
In [60]:
In [313]:
fig = plt.figure
plt.plot(x)
plt.show()
In [ ]:
def thomas(n,alpha,y):
# n,a,b,c are integers and d is a vector of the initial values
x = np.zeros((n,1))
b = (1+2*alpha)*np.ones((n,1))
a = -alpha*np.ones((n,1))
c = -alpha*np.ones((n,1))
c[0,0] = 0
a[n-1,0] = 0
x = np.zeros((n,1))
# forward pass
c[0] = c[0]/b[0]
y[0] = y[0]/b[0]
for k in range(1,n):
c[k] = c[k]/(b[k]-a[k]*c[k-1])
y[k] = (y[k]-a[k]*y[k-1])/(b[k]-a[k]*c[k-1])
# backward pass
x[n-1] = y[n-1]
for k in range(n-2,-1,-1):
x[k] = y[k] - c[k]*x[k+1]
return x
In [ ]:
def runADI(n,alpha):
x = np.zeros((1000,n))
y = np.zeros((n,1))
y[int(np.floor((n-1)/2))] = 1
i = 0
while i < 1000:
temp = thomas(n,alpha,y)
y = temp
x[i,:] = np.transpose(temp)
i += 1
return x
In [ ]:
n = 500
alpha1 = 0.1
alpha2 = 0.1*0.5
alpha3 = 0.1*0.1
d1 = runthomas(n,alpha1)
d2 = runthomas(n,alpha2)
d3 = runthomas(n,alpha3)