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()
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
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
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
Some trouble with this, to discuss in class
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()