import numpy as np
import matplotlib.pyplot as plt
from matplotlib import animation
import sympy as sym
sym.init_printing(use_latex='mathjax')
%matplotlib inline
from IPython.display import HTML
The straightforward finite difference approximation for second order differentials is $$ \frac{u(x+\Delta x, t) - 2u(x, t) + u(x-\Delta x, t)}{\Delta x^2} = \frac{\partial^2 u}{\partial x^2}\bigg|_{x, t} + \mathcal{O}[\Delta x^2] $$ So using the simplification $u(j\Delta x, n\Delta t) = u^n_j$ the 2D wave equation becomes $$ \frac{u_j^{n+1} - 2u^n_j + u_j^{n-1}}{\Delta t^2} = v^2\frac{u_{j+1}^n - 2u^n_j + u_{j-1}^n}{\Delta x^2}$$ which can be solved from initial conditions $u^n_j$ knowing the first step and past time and let $\alpha = \frac{v\Delta t}{\Delta x}$, $$ u^{n+1}_j = -u_j^{n-1} + 2(1 - \alpha^2)u_j^n + \alpha^2(u_{j+1}^n + u_{j-1}^n) $$
(b)¶
what order approximation is this in time and in space?
This approximation is correct to first order in both space and time.
(c)¶
Use the von Neumann stability criterion to find the mode amplitudes.
Von Neumann stability analysis assumes an oscillatory dependence on space and a exponential dependence on time with the ansatz $$ u^n_j = A(k)^ne^{ijk\Delta x} $$ Then plugging this in to our equation and assuming a starting point and letting n be any integer gives $$ A^2e^{ijk\Delta x} = -e^{ijk\Delta x} + 2(1-\alpha^2)Ae^{ijk\Delta x} + \alpha^2\big(Ae^{i(j+1)k\Delta x} + Ae^{i(j-1)k\Delta x}\big) $$ becomes a quadratic equation in A $$ A^2 - 2(1-2\alpha^2 + 2\alpha^2\cos k\Delta x)A + 1 = 0$$ Let $\beta = 1-2\alpha^2 + 2\alpha^2\cos k\Delta x = 1-2\alpha^2\sin^2\big(\frac{k\Delta x}{2}\big)$
A, a, k, x = sym.symbols("A alpha k Delta")
B = sym.symbols("beta")
ans = sym.solve(-A**2 + 2*B*A -1, A )
ans
therefore, $$ A = \beta \pm \sqrt{\beta^2 -1}$$ where $$ \beta = 1-2\bigg(\frac{v\Delta x}{\Delta x}\bigg)^2\sin^2\bigg(\frac{k\Delta x}{2}\bigg) $$
(d)¶
Use this to find a condition on the velocity, time step, and space step for stability.
For stability, $|A| < 1$ and we have two different modes that must be met, as expected for a second order equation. Examining $A$, $\beta$ needs to be <1 to keep the positive case less than 1. Therefore, $A$ will be complex with magnitude $$ |A|^2 = \beta^2 + 1 - \beta^2 = 1 $$ $$ -1 \leq 1-2\bigg(\frac{v\Delta x}{\Delta x}\bigg)^2\sin^2\big(\frac{k\Delta x}{2}\big) \leq 1 $$ Which simplifies to the Courant-Fredrichs-Levy criterion $$ \frac{|v|\Delta t}{\Delta x} \leq 1$$
(e)¶
Do different modes decay at different rates for the stable case?
(f)¶
Numerically solve the wave equation for the evolution from an initial condition with $u=0$ except for one nonzero node, and verify the stability criterion.
dx = 1
dt = 0.1
v = 1
max_t = 500 # number of frames to be generated
max_x = 50 # number of points to be drawn
init_pos = 1
# u is an 2 x max_t array - [[position],[time]]
u_acc = np.zeros((max_t, max_x))
u_acc[0][1] = 1
i = 0
for n in range(1, max_t-1):
for j in range(1, max_x-1):
t_prev = u_acc[n-1][j]
x_next = u_acc[n][j+1]
x_prev = u_acc[n][j-1]
curr_state = u_acc[n][j]
next_u = -t_prev + 2*(1-((v*dt)/dx)**2)*curr_state + (v*dt/dx)**2*(x_next + x_prev)
u_acc[n+1][j] = next_u
plt.plot(np.arange(max_x), u_acc[40])
fig = plt.figure()
ax = plt.axes(xlim=(0, max_x), ylim=(-6, 6))
line, = ax.plot([], [], lw=2)
def init():
line.set_data([], [])
return line,
def animate(i):
x = np.arange(max_x)
y = u_acc[i]
line.set_data(x, y)
return line,
anim = animation.FuncAnimation(fig, animate, init_func=init,
frames=max_t, interval=20)
HTML(anim.to_html5_video())
# dt = 0.1
# v = 1
# dx = 1
# max_t = 200
# max_x = 50
HTML(anim.to_html5_video())
# dt = 0.1
# v = 1
# dx = 1
# max_t = 500
# max_x = 50
(g) If the equation is replaced by $$ \frac{\partial^2 u}{\partial t^2} = v^2\frac{\partial^2 u}{\partial x^2} + \lambda\frac{\partial}{\partial t}\frac{\partial^2u}{\partial x^2}$$ assume that $$ u(x,t) = Ae^{i(kx-\omega t)}$$ and find a relationship between $k$ and $\omega$, and simplify it for small $\lambda$. Comment on the relationship to the preceeding question.
(h) Repeat the numerical solution of the wave equation with the same initial conditions, but include the damping term.
8.2¶
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 = \Delta x = 1$, and use fixed boundary conditions set equal to zero.
(a) Use the explicit finite difference scheme, and look at the behavior for $\Delta t = 1, 0.5$, and 0.1. What step size is required by the Courant condition?
(b) Now repeat this using implicit finite differences and compare the stability.
8.3¶
Use ADI to solve a 2D diffusion problem on lattice, starting with randomly seeded values.
8.4¶
Use SOR to solve Laplace's equation in 2D, with boundary conditions $u_{j, 1} = u_{1, k} = 0, u_{N,k} = -1, u_{j,N}=1$, and explore how the convergence rate depends on $\alpha$, and how the best choice for $\alpha$ depends on the lattice size.