Consider the 1D wave equation
$\frac{\partial^2 u}{\partial t^2} = v^2 \frac{\partial^2 u}{\partial x^2}$
(a) Write down the straightforward finite-difference approximation.
Referring to Equation (8.20) in the textbook, we can easily derive the straightforward approximation:
$\frac{\partial^2 u}{\partial t^2} = \frac{u_j^{n+1}-2u_j^n+u_j^{n-1}}{(\Delta t)^2}, \; \frac{\partial^2 u}{\partial x^2} = \frac{u_{j+1}^n-2u_j^n+u_{j-1}^n}{(\Delta x)^2}$
$\frac{u_j^{n+1}-2u_j^n+u_j^{n-1}}{(\Delta t)^2} = v^2 \left(\frac{u_{j+1}^n-2u_j^n+u_{j-1}^n}{(\Delta x)^2}\right)$
$u_j^{n+1} = 2u_j^n - u_j^{n-1} + v^2 \frac{(\Delta t)^2}{(\Delta x)^2}[u_{j+1}^n - 2u_j^n + u_{j-1}^n]$
(b) What order approximation is this in time and space?
We can clearly see that this is a second order approximation in both space and time.
(c) Use the von Neumann stability criterion to find the mode amplitudes.
$A = 2 - \frac{1}{A} + v^2 \frac{(\Delta t)^2}{(\Delta x)^2}[e^{ik \Delta x} - 2 + e^{-ik \Delta x}]$
$= 2 - \frac{1}{A} + v^2 \frac{(\Delta t)^2}{(\Delta x)^2} 2 [\text{cos}(k \Delta x) - 1]$
$= 2 - \frac{1}{A} - 4v^2 \frac{(\Delta t)^2}{(\Delta x)^2}\text{sin}^2\frac{k \Delta x}{2}$
$A^2 + \left[4v^2 \frac{(\Delta t)^2}{(\Delta x)^2}\text{sin}^2\frac{k \Delta x}{2} - 2\right]A + 1 = 0$
We will substitute the coefficient for the first order $A$ to be defined as $2b$, which gives us:
$A^2 + 2bA + 1 = 0$
$A = -b \pm \sqrt{b^2 - 1}$
(d) Use this to find a condition on the velocity, time step, and space step for stability (hint: consider the product of the two amplitude solutions).
Multiplying the two $A$ solutions give us:
$b^2 - (b^2 - 1) = 1$
In order for this to be true, either both solutions for $A$ have to be $1$, or one of the solutions will be larger than $1$, and the other will be smaller than $1$, which would make it unstable. For both of the solutions to have a magnitude of $1$, that means they must be complex conjugates:
$A = -b \pm i \sqrt{1 - b^2}$
$|A| = b^2 + (1 - b^2) = 1$
Which requires the term inside of the square root to be positive, so we have:
$1 - b^2 \ge 0$
Substituting back in our coefficient for $2b$ from part (c), we get:
$1 - \left(2v^2 \frac{(\Delta t)^2}{(\Delta x)^2} \text{sin}^2 \frac{k \Delta x}{2} - 1\right) \ge 0$
Since the sine term will become a $1$, we only need to worry about the term with velocity, time, and position:
$v \frac{\Delta t}{\Delta x} \le 1$
This is the Courant condition.
(e) Do different modes decay at different rates for the stable case?
No, because as we can see from part (d), $A$ is independent of $k$, so there is no damping.
(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.
See Neil's str.html example.
(g) If the equation is replaced by
$\frac{\partial^2 u}{\partial t^2} = v^2 \frac{\partial^2 u}{\partial x^2} + \gamma \frac{\partial}{\partial t} \frac{\partial^2 u}{\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 $\gamma$. Comment on the relationship to the preceeding question.
$-\omega^2 = -v^2k^2 + i \gamma k^2 \omega$
$\omega^2 + i \gamma k^2 \omega - v^2k^2 = 0$
$\omega = \frac{1}{2}\left(-i \gamma k^2 \pm \sqrt{-\gamma^2k^4 + 4v^2k^2}\right)$
$= -\frac{1}{2}i\gamma k^2 \pm vk \sqrt{1 - \frac{\gamma^2 k^2}{4v^2}}$
$\approx -\frac{1}{2}i\gamma k^2 \pm vk\left(1 - \frac{\gamma^2 k^2}{8v^2}\right)$
$\approx -\frac{1}{2}i\gamma k^2 \pm vk$
These solutions represent travelling waves in either direction that are exponentially damped in time.
(h) Repeat the numerical solution of the wave equation with the same initial conditions, but include the damping term.
Again, see Neil's str.html example.
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?
We will plug into equation (8.21), which gives us:
$\frac{2D\Delta t}{(\Delta x)^2} = \frac{2 \times 1 \times \Delta t}{1^2} \le 1 \implies \Delta t \le \frac{1}{2}$
See Neil's diffexp.html example.
(b) Now repeat this using implicit finite differences and compare the stability.
See Neil's diffimp.html example.
Use ADI to solve a 2D diffusion problem on a lattice, starting with randomly seeded values.
See Neil's diffadi.html example.
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.
See Neil's sor.html example.