Goal: predict whether each represented mode decays, changes sign or grows. Use A = 0.1 to connect to the same field as the earlier chapters.
An exact solution to compare with
Take the one-dimensional heat equation on a periodic interval,
\[ u_t=D u_{xx},\qquad D>0,\qquad u(x,0)=1+A\cos(kx). \]
On \([0,2\pi)\), \(k\) is an integer. The mean remains 1 and the oscillation decays:
\[ u(x,t)=1+A e^{-Dk^2t}\cos(kx). \]
Substitution into the PDE checks the result independently of a numerical solver. Shorter wavelengths decay faster because the rate is proportional to \(k^2\).
Space becomes a matrix, then time advances it
The centred second-difference Laplacian is
\[ (L_hu)_j=\frac{u_{j-1}-2u_j+u_{j+1}}{h^2}. \]
On a Fourier mode its eigenvalue is
\[ \lambda_{\mathrm{FD2}}(\theta) =-\frac{4}{h^2}\sin^2(\theta/2),\qquad \theta=kh. \]
A Fourier spectral second derivative instead uses \(\lambda_{\mathrm{F}}=-k^2\). Unlike the real first derivative discussed in the previous lesson, the Nyquist cosine has a nonzero second derivative. For an even grid, its \(k^2=(\pi/h)^2\) participates in the diffusion stability bound.
Each mode now obeys an ordinary differential equation, \(a'=D\lambda a\). Forward Euler replaces it by \(a_{n+1}=g a_n\), where \(g=1+D\Delta t\lambda\). Writing \(r=D\Delta t/h^2\) gives
\[ g_{\mathrm{FD2}}=1-4r\sin^2(\theta/2),\qquad g_{\mathrm{F}}=1-r\theta^2. \]
The exact PDE gain over one time step is \(g_{\mathrm{exact}}=\exp(-r\theta^2)\).
Stability is a condition on every represented mode
For nonpositive real eigenvalues, forward Euler requires \(|g|\leq1\). The worst mode on an even periodic 1D grid gives
\[ r\leq\frac12\quad\text{for FD2},\qquad r\leq\frac{2}{\pi^2}\approx0.203\quad\text{for Fourier differentiation}. \]
A more accurate spatial derivative can therefore require a smaller explicit time step. Stability also does not mean faithful decay: \(g=-1\) changes sign forever instead of damping. For strict decay of every nonconstant mode, stay inside the bound.
Try to break the method
The first plot shows every grid wavenumber’s one-step gain. Curves outside the horizontal band \([-1,1]\) reveal a forward-Euler instability. The second follows your selected mode through repeated steps. Its vertical coordinate is \(\log_{10}|a_n/a_0|\): a value of \(-3\) means an amplitude one thousandth of its initial value. Signs are reported separately through the one-step gains.
Interactive experiment / periodic 1D heat equation
Selected mode: gFD = 0.700, gFourier = 0.630. FD2 grid: stable; Fourier grid: stable for forward Euler in 1D. A negative g flips sign each step; |g| > 1 grows. A grid can be unstable even when the selected mode decays. Interactive controls appear when the experiment is ready.
Analytic modal symbols evaluated in the browser, not OpenPFC solver output. Exact zero amplitude is displayed at −300 on the logarithmic plot. With N fixed, changing r changes Δt and the final physical time.
For the highest grid frequency, \(\theta=\pi\):
| \(r\) | FD2 + Euler gain | Fourier + Euler gain | What to expect |
|---|---|---|---|
| 0.15 | 0.40 | −0.48 | Both bounded; Fourier’s sign alternates |
| 0.35 | −0.40 | −2.45 | FD2 bounded; Fourier grows |
| 0.70 | −1.80 | −5.91 | Both grow |
The experiment can show a decaying selected mode while another mode on the grid is unstable. Stability is about all modes, including errors introduced by roundoff or other operations, rather than only the initial signal you intended to put there.
Change the time integrator
For this constant-coefficient linear equation, Fourier transformation diagonalizes the spatial operator. Instead of forward Euler, apply
\[ \widehat u_k^{\,n+1} =\exp(-Dk^2\Delta t)\widehat u_k^{\,n}. \]
This is exact evolution of each represented Fourier mode over the step. It has no forward-Euler diffusion restriction for this equation. The exponential curve in the experiment is that exact single-mode reference. It does not mean arbitrary nonlinear equations can be advanced exactly or without time-step accuracy constraints. Exponential time differencing uses this idea for a chosen linear part and approximates the remainder.
The mean is the \(k=0\) mode. Both Euler symbols and the exponential have \(g(0)=1\), so they preserve the mean in exact arithmetic. This is useful to check, but it does not establish that nonzero modes decay at the right rate.
Checkpoint: a decaying mode on an unstable grid
At r = 0.35, can a low Fourier mode decay while the full Fourier + Euler grid is unstable?
Yes. For a low mode, \(\theta=kh\) is small and \(|1-r\theta^2|<1\). At the highest mode, \(\theta=\pi\) and \(g=1-0.35\pi^2\approx-2.45\), outside the stability interval. Errors in that mode can grow even if the intended initial mode decays. The stability bound must cover every represented mode. Reducing \(r\) below \(2/\pi^2\) restores the 1D Fourier + Euler bound; then an error study is still needed.
Take the experiment back to OpenPFC
The diffusion starter uses a three-dimensional periodic FD2 grid and forward Euler. Its input validation requires \(D\Delta t/h^2\leq1/6\) for its unit grid spacing. Our browser lesson is one-dimensional, with a bound of \(1/2\); do not copy that bound into the 3D program. In \(d\) dimensions with equal spacing, the most negative FD2 Laplacian eigenvalue is \(-4d/h^2\).
Try these controlled comparisons:
- In the browser, select the middle preset. Explain why Fourier’s derivative accuracy does not prevent an unstable Euler update.
- In the starter, halve \(\Delta t\) and double the step count. Compare at the same final time, using the analytic decay rate as your reference.
- Keep a fixed physical domain while refining the grid. What happens to the explicit stability limit? What changes if you accidentally keep unit spacing and enlarge the domain instead?
Continue with the C++ examples, the numerics limits and Trefethen’s finite-difference and spectral text for further study of accuracy, stability and convergence.