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)