Finite differencing of the diffusion (heat) equation

Question

I am attempting to code a problem for a meteorology class. Our initial equation was as follows: $$\tfrac{\partial u}{\partial t} = \nu \tfrac{\partial^2 u}{\partial x^2} (*)$$ We were then assigned this finite differencing to program: $$\tfrac{\zeta_{j}^{n+1}-\zeta_{j}^{n-1}} {2\Delta t} = \nu \tfrac{\zeta_{j+1}^{n}-\zeta_{j}^{n+1}-\zeta_{j}^{n-1}+\zeta_{j-1}^{n}}{\Delta x^{2}} (\&)$$
n and j denote time and space respectively.

The problem reads as follows: "Graph (for t = 10000s) one wavelength of the discrete solutions for k=2 [i.e., interval: -$\pi /4$, 3$\pi /4$] and k=4 [i.e., interval: -$\pi /8$, 3$\pi /8$] assuming $\zeta_{0}^{n}=0,\zeta_{j}^{n}=0,\zeta_{j}^{1}=\zeta_{j}^{0}$ for j = 1,......,nx-1, and $\Delta x = L/(nx-1)$, where $L=2\pi /k$. Use $\Delta t = 10s, 1000s, 2000s$, and $5000s.$ How does your analytic solution from solving (*) above compare to the discrete solution a t = 10000s add the analytic solution for t = 10000s to your graph)? Note that there should be two graphs with 5 curves here! What happens to your discrete solution as $\Delta t$ increases?"

I have never seen such a differencing as in the RHS of (&). The LHS is clearly central differencing, yet the right hand side is something I unfamiliar with. It is clearly not forward space, central space, or back space. I am supposed to somehow rearrange (&) to solve for $\zeta_{j}^{n+1}$ and then program the above problem. I do not see how I am supposed to do this, either mathematically or analytically as I do not know of any such differencing methods, and I have spent over 6 hours trying to answer this question.

My apologies if this belongs in the coding forum as I was uncertain whether to post this in math or coding. I am completely unaware of any ways to simultaneously go through two different loops in code, and I have never seen this kind of differencing that combines time and space.

Answer 1

I will give some pointers on how to solve the system numerically.

For simplicity set $\alpha = \left(\frac{\Delta x^2}{2\nu\Delta t}\right)$ and then solve the equation for $\zeta^{n+1}_j$ to find

$$\zeta_{j}^{n+1} = \frac{\zeta_{j+1}^{n}+\zeta_{j-1}^{n}}{1+\alpha} + \frac{\alpha-1}{\alpha + 1}\zeta_{j}^{n-1}$$

This equation tells us that if we know, for any fixed $n$, $\zeta_j^n$ and $\zeta_j^{n-1}$ for all $j$ then the equation above can be used to calculate $\zeta_{j}^{n+1}$ for all $j$.

A simple way to codes this up (there are better memory-saving ways to do this, but this should be good enough here) is to first make an array $\zeta[0:n_t,0:n_x]$ where $n_t$ is the number of time-steps and $n_x$ is the number of grid-cells.

Next we need to input the initial condition $\zeta^1_j = \zeta^0_j = \ldots$ in the entries $\zeta[0,0:n_x]$ and $\zeta[1,0:n_x]$.

Finally we need to propagate the equation to get the solution at all times. This can be best explained by the following pseudo-code

// Loop over time-steps

Loop over $n = 2,3,4,\ldots,n_t$

// Set the boundary condition $\zeta_0^n = \zeta_{n_x}^n = 0$

$\zeta[n,0] = 0$

$\zeta[n,n_x] = 0$

// Calculate the solution for all $j$

Loop over $j=1,2,\ldots, n_x-1$

$\zeta[n,j] = \frac{\zeta[n,j+1]+\zeta[n,j-1]}{1+\alpha} + \frac{\alpha-1}{\alpha + 1}\zeta[n-1,j]$

After this is done you have the solution $\zeta(t,x)$ at the points $(t,x) = \left(n\Delta t, x_{\rm min} + \frac{(x_{\rm max}-x_{\rm min})j}{n_x}\right)$ stored in $\zeta[n,j]$.

Another thing to keep in mind when solving is that a stability analysis give that the numerical scheme is only stable if $\alpha \geq 1$. So if you get nonsensical solutions then I would check if this condition is satisfied for your choice of $\Delta x$, $\Delta t$.

Answer 2

what is on the right hand side is the second derivative approximation of $u$ with respect to $x.$ to see this in one variable case compute $f(1+x) - 2f(1) + f(1-x)$ using taylor series.

$$f(1+h) = f(1) + hf'(1) + \frac{1}{2}h^2f''(1) + \cdots, f(1-h) = f(1) - hf'(1) + \frac{1}{2}h^2f''(1) + \cdots $$

adding the two you get $$f(1+h) + f(1-h) = 2f(1) + h^2f''(1) + \cdots $$

approximating this gets you $f''(1) \simeq \dfrac{f(1+h) - 2f(1) + f(1-h)}{h^2}$

Answer 3

I find it useful to visualize a computational stencil for the algorithm. There are several of these listed on the Wikipedia page for the finite difference method, though, not yours. Yours can be visualized like so:

The dots all represent points in a grid where you want to compute an estimate for your solution. The large dots represent values known at the outset from either an initial or boundary condition. The small dots represent points where we plan to compute the value. At any step, the green disks will represent known values, while the red disk represents the value you will compute, in terms of the known green values using the formula presented by Withner: $$\zeta_{j}^{n+1} = \frac{\zeta_{j+1}^{n}+\zeta_{j-1}^{n}}{1+\alpha} + \frac{\alpha-1}{\alpha + 1}\zeta_{j}^{n-1}.$$ Your loop will essentially move the stencil so that the red disk traces out the row of small dots near the bottom. You'll then increment $j$ to move the stencil up to the next row and repeat.