diff --git a/week5/README.MD b/week5/README.MD new file mode 100644 index 0000000000000000000000000000000000000000..24bc2ce7329164b22219cd3b7025a123a00847ec --- /dev/null +++ b/week5/README.MD @@ -0,0 +1,24 @@ +# Week 5 - Partial Differential Equations + +### 9.1 + +The wave equation with and without damping for a dirac delta and sinusoidal initial conditions. + + + + + +### 9.2 + +Increasing dt's above (and below) the courant condition, resulting in numerical instability + + + + + +### 9.4 + +The 2D laplace equation using successive over-relaxation with fixed boundary conditions (and a fixed dirac-delta) + + + \ No newline at end of file diff --git a/week5/img/dt0_1.gif b/week5/img/dt0_1.gif new file mode 100644 index 0000000000000000000000000000000000000000..c045bf1185467cc2c03c184d6ce77c857aea962a Binary files /dev/null and b/week5/img/dt0_1.gif differ diff --git a/week5/img/dt0_5.gif b/week5/img/dt0_5.gif new file mode 100644 index 0000000000000000000000000000000000000000..22d41ed8a698a262bb0376b3dd221e848c94a218 Binary files /dev/null and b/week5/img/dt0_5.gif differ diff --git a/week5/img/dt1_0.gif b/week5/img/dt1_0.gif new file mode 100644 index 0000000000000000000000000000000000000000..33808d3eee36d57263c370e69cebaedca855a6e5 Binary files /dev/null and b/week5/img/dt1_0.gif differ diff --git a/week5/img/laplace_2D_dirac.gif b/week5/img/laplace_2D_dirac.gif new file mode 100644 index 0000000000000000000000000000000000000000..12bfb29922621a0d02032d2efb0f5e5e285285f5 Binary files /dev/null and b/week5/img/laplace_2D_dirac.gif differ diff --git a/week5/img/laplace_2D_edges.gif b/week5/img/laplace_2D_edges.gif new file mode 100644 index 0000000000000000000000000000000000000000..9cabe0b91b8d20fa25c060277d01c5d06d6734f7 Binary files /dev/null and b/week5/img/laplace_2D_edges.gif differ diff --git a/week5/img/wave_dirac.gif b/week5/img/wave_dirac.gif new file mode 100644 index 0000000000000000000000000000000000000000..2d4f1d04f7eaadeedde8345d81f9b3fd2ae2cec5 Binary files /dev/null and b/week5/img/wave_dirac.gif differ diff --git a/week5/img/wave_no_damp.gif b/week5/img/wave_no_damp.gif new file mode 100644 index 0000000000000000000000000000000000000000..d32588d088d640e6862db1962f0db51de781e7ba Binary files /dev/null and b/week5/img/wave_no_damp.gif differ diff --git a/week5/img/wave_yes_damp.gif b/week5/img/wave_yes_damp.gif new file mode 100644 index 0000000000000000000000000000000000000000..3b3974dc2f410950f0cfe2bb220829a2f4e7946c Binary files /dev/null and b/week5/img/wave_yes_damp.gif differ diff --git a/week5/week5_1.py b/week5/week5_1.py new file mode 100644 index 0000000000000000000000000000000000000000..8382fbba48bf60244ad9632a18c0391dec078f7c --- /dev/null +++ b/week5/week5_1.py @@ -0,0 +1,83 @@ +import numpy as np +import matplotlib.pyplot as plt +import matplotlib.pyplot as plt +import matplotlib.animation as animation +from matplotlib.animation import FuncAnimation, writers + +''' +6.1) +The wave equation suggests that where you have a large curvature (rate of change of the slope) you will have a large change with respect to time, proportional to c. +https://www.youtube.com/watch?v=lHGMjfCPZik +Recusion relation taken from here: https://wiki.seg.org/wiki/Solving_the_wave_equation_in_1D +y(t+dt) = 2*y(t) - y(t-dt) + r*(y_right - 2*y_self + y_left) +where r = (c*dt/dx)^2 +''' + +class tensioned_string(): + + def __init__(self, x, y0, c): + self.x, self.y, self.y0, self.c = np.copy(x), self.pad_array(y0), self.pad_array(y0), c + self.y_prev = np.copy(self.y0) + + def pad_array(self, arr): # pads the array with an additional first and last element, which will allow us to convolve the edges (these will share the boundary condition) + return np.concatenate((arr[:1], arr, arr[-1:])) + + def increment(self, dt): + ''' increment by dt ''' + r = (self.c * dt / dx)**2 # you can use np.gradient for non-constant dx + temp = np.copy(self.y) + # now apply the recursion relation + self.y[1:-1] = 2*self.y[1:-1] - self.y_prev[1:-1] + r*(self.y[2:] - 2*self.y[1:-1] + self.y[:-2]) + # with damping (assume dt is built into gamma because it was going unstable) + # self.y[1:-1] -= 0.005*(self.y[1:-1] - temp[1:-1]) + self.y_prev = temp + + # boundary condition of clamped end points. We accept that they've changed above, but now we will just reset them to their original values + self.y[[0,1,-2,-2]] = self.y0[[0,1,-2,-2]] + +L = 100 +dx = 0.5 +x = np.linspace(0, 100, int(L/dx)) # string is 100 units long, with 128 points +y = np.empty_like(x) + +d0 = 0.1 # initial displacement of the string + +# Initial Conditions - sinusoidal ripple +for i in range(len(x)): + y[i] = d0*(np.sin((2*np.pi)*x[i]/10) / (10 + (x[i] - 10)**2)) + +# Initial Conditions - plucked string +# d0_location = 0.7 # % where along the string the initial displacement occurs +# y[64] = d0 +# set everything between our "pluck" location to be linear to that point +# y[np.where(x <= d0_location*x[-1])] = d0/(d0_location*x[-1]) * x[np.where(x <= d0_location*x[-1])] +# y[np.where(x > d0_location*x[-1])] = -d0/((1.0 - d0_location)*x[-1]) * (x[np.where(x > d0_location*x[-1])] - x[-1]) + +c = 5 # speed of waves +dt = 0.01 # time step +total_time = 300 +fps = 20 + +string_1 = tensioned_string(x, y, c) + +fig,ax = plt.subplots() +ax.plot(x, y) +line, = ax.plot(x,y) +ax.set_ylim([-1.5*max(y),1.5*max(y)]) # set our axes limits to be slighty past the intial displacement + +def init(): + return line, + +def update(i): + string_1.increment(dt) + line.set_ydata(string_1.y[1:-1]) # update the data + return line, + +anim = animation.FuncAnimation(fig, update, init_func=init, frames=int(600), blit=True, interval = 0.05, repeat=True) +# plt.show() + +my_dpi = 96 +f = r"c://Users/david/Desktop/NMM/week5/img/wave_no_damp.gif" +fps2 = 20 +writervideo = animation.FFMpegWriter(fps=fps2, extra_args=['-vcodec', 'libx264']) +anim.save(f, writer='imagemagick', fps=fps2, dpi=my_dpi) \ No newline at end of file diff --git a/week5/week5_2.py b/week5/week5_2.py new file mode 100644 index 0000000000000000000000000000000000000000..a506547b7abfaa88616a044a7d5b43e05a4d0d83 --- /dev/null +++ b/week5/week5_2.py @@ -0,0 +1,70 @@ +import numpy as np +import matplotlib.pyplot as plt +import matplotlib.pyplot as plt +import matplotlib.animation as animation + +''' +6.2) (Main Question) +Write a program to solve a 1D diffusion problem on a lattice of 500 sites, with an initial condition of zero +at all the sites, except the central site which starts at the value 1.0. Take D = ∆x = 1, and use fixed boundary conditions set equal to zero. +D = ∆x = 1 +''' + +lattice_num = 500 +total_time = 50 +dt = 0.1 + +state = np.zeros(lattice_num) +state[lattice_num // 2] = 1 # set the middle site to 1, // rounds down +all_states = np.zeros((int(total_time/dt),lattice_num)) +all_states[0] = np.copy(state) +diffusion = 1 +dx = 1 + +def time_step(last_state, diffusion): + new_state = np.copy(last_state) # deep copy the last state + alpha = diffusion * dt / (dx ** 2) # alpha is the potentially time/space dependent diffusion constant + filt = np.array([alpha, 1 - 2*alpha, alpha]) # convolution filter for diffusion + new_state[1:-1] = np.convolve(last_state, filt, mode='valid') # convolution + return new_state + +for t in range(total_time): + state = time_step(state, diffusion) + # input() + # print(state) + all_states[t] = np.copy(state) + +my_dpi = 96 +fig = plt.figure(figsize=(400/my_dpi, 400/my_dpi), dpi=my_dpi) +ax = fig.add_subplot(1, 1, 1) +ax.set_xlim(0, lattice_num) +ax.set_ylim(-0.5, 1) +ax.set_title(f'dt: {dt: .2f} s', fontdict=None) +ax.grid() + +x = np.arange(0,lattice_num) +plot = ax.plot(x, all_states[0]) + +def animate(i): + plot[0].set_data(x, all_states[i]) + return plot + +anim = animation.FuncAnimation(fig, animate, range(total_time), blit=False, interval=50, repeat=True) +plt.show() + +# f = r"c://Users/david/Desktop/NMM/week5/img/dt1_0.gif" +# writervideo = animation.FFMpegWriter(fps=30, extra_args=['-vcodec', 'libx264']) +# anim.save(f, writer='imagemagick', fps=30, dpi=my_dpi) + +''' +6.2) A +Use the explicit finite difference scheme, and look at the behavior for ∆t = 1, 0.5, and 0.1. +What step size is required by the Courant condition? + +See the 3 gifs at each dt, which show the "boom" at dt=1, the start of stability at dt=0.5, and the acceptable solution at dt=0.1. +Courant condition: dt <= dx^2 / (2*diffusion) +See answer of 0.5 below: +''' + +dt_max = dx**2 / (2*diffusion) +print(dt_max) \ No newline at end of file diff --git a/week5/week5_3.py b/week5/week5_3.py new file mode 100644 index 0000000000000000000000000000000000000000..93e172c4441f959fbafec48891067a4aa2e70463 --- /dev/null +++ b/week5/week5_3.py @@ -0,0 +1,22 @@ +import numpy as np +import matplotlib.pyplot as plt +import matplotlib.pyplot as plt +import matplotlib.animation as animation + +''' +6.3) +Use ADI to solve a 2D diffusion problem on a lattice, starting with randomly seeded values. +''' + +# ADI = Alternating Direction Implicit -- https://en.wikipedia.org/wiki/Alternating-direction_implicit_method +lattice_x_num = 50 +lattice_y_num = 50 +total_time = 100 +dt = 0.1 + +state = np.zeros((lattice_y_num, lattice_x_num)) +state[lattice_y_num // 2][lattice_x_num // 2] = 1 # set the middle site to 1, // rounds down +all_states = np.zeros((int(total_time/dt),lattice_num)) +all_states[0] = np.copy(state) +diffusion = 1 +dx = 1 \ No newline at end of file diff --git a/week5/week5_4.py b/week5/week5_4.py new file mode 100644 index 0000000000000000000000000000000000000000..d6b5505218aff88adc13225d68bf3ba5f5f85d0f --- /dev/null +++ b/week5/week5_4.py @@ -0,0 +1,64 @@ +import numpy as np +import matplotlib.pyplot as plt +from scipy import signal +import matplotlib.animation as animation +from matplotlib import cm + +''' +6.4) +Use SOR to solve Laplace's equation in 2D, with boundary conditions +uj,1 = u1,k = 0, uN,k = -1,uj,N = 1 +and explore how the convergence rate depends on alpha, and how the best choice for alpha depends on the lattice size. + +Equation 9.69 with rho = 0 (forcing function?), therefore our last term is zero +''' + +# SOR = Successive Over-Relaxation -- https://en.wikipedia.org/wiki/Successive_over-relaxation + +def SOR_step(u): + # we are literally just going to implement equation 9.69 here, starting with the first term + new_u = (1 - alpha) * u + # enforce boundary or initial conditions + new_u[0, :], new_u[:, 0], new_u[-1, :], new_u[:, -1] = 0, 0, 1, -1 + # new_u[lattice_size//2, lattice_size//2] = 1 # set the middle point to 1 + + + # Do it with a scipy convolution for speed + filt = np.array([[ 0, alpha/4, 0], + [alpha/4, 0, alpha/4], + [ 0, alpha/4, 0]]) + new_u = signal.convolve2d(u, filt, boundary='fill', mode='same') # convolve it up + new_u[0, :], new_u[:, 0], new_u[-1, :], new_u[:, -1] = 0, 0, 1, -1 + # new_u[lattice_size // 2, lattice_size //2] = 1 + return new_u + +lattice_size = 50 +total_time = 100 +dt = 0.05 +alpha = 1 +fps = 30 + +u = np.zeros((lattice_size, lattice_size)) + +fig = plt.figure() +ax = fig.add_subplot(111, projection='3d') + +xs = np.linspace(-1, 1, len(u)) +ys = np.linspace(-1, 1, len(u)) +X, Y = np.meshgrid(xs, ys) +u = SOR_step(u) + +surf = ax.plot_surface(X, Y, u, rstride=2, cstride=2 ,cmap=cm.jet) # plot a 3d surface plot + +def update(i): + global u + u = SOR_step(u) + ax.cla() # fairly certain that ax.cla() is very slow, but can't find the update function + ax.plot_surface(X, Y, u, cmap = 'viridis', edgecolor = 'none') + +anim = animation.FuncAnimation(fig, update, total_time, interval=1000/fps, blit = False) +# plt.show() + +f = r"c://Users/david/Desktop/NMM/week5/img/laplace_2D_edges.gif" +writervideo = animation.FFMpegWriter(fps=fps/2, extra_args=['-vcodec', 'libx264']) +anim.save(f, writer='imagemagick', fps=fps/3, dpi=96)