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)
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)
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)