# Numerical Solver for 1D Wave Equation¶

In [ ]:
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation

def wave_eq_solver(dx, dt, c, X, T):
C2 = (c*dt/dx)**2
num_t = int(T / dt)
num_x = int(X / dx)
u = np.zeros((num_t,num_x))

# set initial values - only single node
u[0,2] = 1.0
u[1,2] = 1.0

#set initial pulse shape to something more interesting
#     def gaussian(x):
#         val = np.exp(-(x**2)/0.25)
#         if val<.001:
#             return 0.0
#         else:
#             return val

#     for j in range(0,num_t):
#         u[0,j]=gaussian(j*dx)
#         u[1,j]=u[0,j]

# run the solver
for n in range(1,num_t-1):
for j in range(1,num_x-1):
u[n+1,j] = 2*u[n,j]-u[n-1,j]+C2*(u[n,j+1]+u[n,j-1]-2*u[n,j])

return u

if __name__ == "__main__":
dx = 0.1
dt = 0.01
T = 2
X = 10
c = 1.0

u = wave_eq_solver(dx, dt, c, X, T)
print u

fig = plt.figure()
plts = []             # get ready to populate this list the Line artists to be plotted
plt.hold("off")
num_t = int(T / dt)
for i in range(num_t):
p, = plt.plot(u[i,:], 'k')   # this is how you'd plot a single line...
plts.append( [p] )           # ... but save the line artist for the animation
ani = animation.ArtistAnimation(fig, plts, interval=50, repeat_delay=3000)   # run the animation
ani.save('wave.mp4')    # optionally save it to a file
plt.show()

In [ ]:
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation

def wave_eq_solver(dx, dt, c, gamma, X, T):
#     C2 = (c*dt/dx)**2
num_t = int(T / dt) + 1
num_x = int(X / dx) + 1
u = np.zeros((num_t,num_x))
print num_t
print num_x

alpha = ((c**2)*(dt**2) / (dx**2))
beta = (gamma * dt / (dx**2))

# set initial values - only single node
u[0,2] = 1.0
u[1,2] = 1.0

#set initial pulse shape
#     def gaussian(x):
#         val = np.exp(-(x**2)/0.25)
#         if val<.001:
#             return 0.0
#         else:
#             return val

#     for j in range(0,num_t):
#         u[0,j]=gaussian(j*dx)
#         u[1,j]=u[0,j]

# run the solver with the damping term included

for n in range(1,num_t - 1):
for j in range(1,num_x - 1):
u[n+1][j] = alpha*(u[n][j+1] - 2*u[n][j] + u[n][j-1]) + beta*(u[n][j+1] - 2*u[n][j] + u[n][j-1] - u[n-1][j+1] + 2*u[n-1][j] - u[n-1][j-1]) + 2*u[n][j] - u[n-1][j]

return u

if __name__ == "__main__":
dx = 0.1
dt = 0.01
T = 2
X = 10
c = 1.0
gamma = 0.1

u = wave_eq_solver(dx, dt, c, gamma, X, T)

print u

fig = plt.figure()
plts = []             # get ready to populate this list the Line artists to be plotted
plt.hold("off")
num_t = int(T / dt)
print num_t
for i in range(num_t):
p, = plt.plot(u[i,:], 'k')   # this is how you'd plot a single line...
plts.append( [p] )           # ... but save the line artist for the animation
ani = animation.ArtistAnimation(fig, plts, interval=50, repeat_delay=3000)   # run the animation
ani.save('wave_damping.mp4')    # optionally save it to a file


# 1D Diffusion Solver for 500 Point Lattice - Explicit Method¶

In [ ]:
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation

def diffusion_eq_solver(dx, dt, D, X, T):
num_t = int(T / dt) + 1
num_x = int(X / dx) + 1
u = np.zeros((num_t,num_x))

# set initial values - only single node at central site
u[0,num_x/2] = 1.0
u[1,num_x/2] = 1.0

print ((D*dt)/(dx**2))

# run the solver with the damping term included
for n in range(1,num_t - 1):
for j in range(1,num_x - 1):
u[n+1][j] = u[n][j] + (((D*dt)/(dx**2))*(u[n][j+1] - 2.0*(u[n][j]) + u[n][j-1]))

return u

if __name__ == "__main__":
dx = 1
dt = 0.1
T = 20
X = 500
D = 1

u = diffusion_eq_solver(dx, dt, D, X, T)

print u

fig = plt.figure()
plts = []             # get ready to populate this list the Line artists to be plotted
plt.hold("off")
num_t = int(T / dt)
for i in range(num_t):
p, = plt.plot(u[i,:])   # this is how you'd plot a single line...
plts.append( [p] )           # ... but save the line artist for the animation
ani = animation.ArtistAnimation(fig, plts, interval=50, repeat_delay=3000)   # run the animation
ani.save('diffusion_1D.mp4')    # optionally save it to a file


# 1D Diffusion Solver for 500 Point Lattice - Implicit Method¶

In [ ]:
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation

def gaussian_inv(M):
# replace with gaussian eliminator
return np.linalg.inv(M)

def tridiag(a, b, c, k1=-1, k2=0, k3=1):
return np.diag(a, k1) + np.diag(b, k2) + np.diag(c, k3)

def diffusion_solver_implicit(dx, dt, D, X, T):
num_x = int((X / dx) + 1)
num_t = int((T / dt) + 1)
alpha = D * dt / (dx**2)
# left and right diag
l_diag = [-1.0*alpha] * (num_x - 1)
# center diag
c_diag = [1 + 2*alpha] * num_x
M = tridiag(l_diag, c_diag, l_diag)
M_inv = gaussian_inv(M)
u = np.zeros((num_t,num_x))

# set initial values - only single node at central site
u[0,num_x/2] = 1.0
u[1,num_x/2] = 1.0

# run the solver
for n in range(1,num_t - 1):
for j in range(1,num_x - 1):
#u[n+1][j] = u[n][j] + (D*dt)/(dx**2) * (u[n+1][j+1] - 2*u[n+1][j] + u[n+1][j-1])
u[n+1][j] = u[n][j] + (D*dt)/(dx**2) * ( np.dot(M_inv[j+1],u[n]) - 2*(np.dot(M_inv[j],u[n])) + np.dot(M_inv[j-1],u[n]) )

return u

if __name__ == "__main__":
dx = 1
dt = 0.1
T = 20
X = 500
D = 1

u = diffusion_solver_implicit(dx, dt, D, X, T)
print u

fig = plt.figure()
plts = []             # get ready to populate this list the Line artists to be plotted
plt.hold("off")
num_t = int(T / dt)
for i in range(num_t):
p, = plt.plot(u[i,:], 'k')   # this is how you'd plot a single line...
plts.append( [p] )           # ... but save the line artist for the animation
ani = animation.ArtistAnimation(fig, plts, interval=50, repeat_delay=3000)   # run the animation
ani.save('diffusion_1D_implicit.mp4')    # optionally save it to a file


# 2D Diffusion Lattice Problem - Solve with ADI¶

Some trouble with this, to discuss in class

# 2D Laplace Equation - Solve with SOR¶

In [ ]:
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation

def laplace_eq_solver(num_j, num_k, num_steps, alpha):
images = []
u = np.zeros((num_j,num_k))

# set initial values
u[num_j-1:] = -1
u[:num_k-1] = 1

# run the solver with the damping term included
for n in range(1, num_steps - 1):
for j in range(1,num_j - 1):
for  k in range(1,num_k - 1):
u[j][k] = ((1-alpha)*u[j][k]) + 0.25*alpha*(u[j+1][k] + u[j-1][k] + u[j][k+1] + u[j][k-1])
images.append(u)

print (images[-1] == images[-5])
return u, images

if __name__ == "__main__":
num_j = 50
num_k = 50
num_steps = 100
alpha = 0.3

u,images = laplace_eq_solver(num_j, num_k, num_steps, alpha)

print u
print len(images)

fig = plt.figure()
plts = []
plt.hold("off")
for s in range(num_steps - 2):
p = plt.imshow(images[s])
plts.append( [p] )
ani = animation.ArtistAnimation(fig, plts, interval=50, repeat_delay=3000)   # run the animation
ani.save('laplace_SOR.mp4')
plt.show()


In [ ]: