A finite system of nonlinear ODEs
Goal: derive the time update and verify its order at fixed spatial cutoff. After the cubic Galerkin projection, each retained coefficient obeys
\[\dot c_k=L_kc_k+\mathcal N_k(c),\qquad L_k=1-Dk^2,\quad \mathcal N_k(c)=-\widehat{u_K^3}_k.\]
The reaction’s linear term \(+u\) is included in \(L\); the cubic belongs to \(\mathcal N\). A different split changes the numerical method and its accuracy properties. Variation of constants gives
\[c(t+h)=e^{hL}c(t)+\int_0^h e^{(h-s)L}\mathcal N(c(t+s))\,ds.\]
The diagonal exponential is cheap in Fourier space. The remaining integral still depends on a changing nonlinear field.
First freeze the nonlinear rate
Taking the nonlinear rate at the start of the step gives ETD1:
\[a=e^{hL}c^n+h\varphi_1(hL)\mathcal N(c^n),\qquad \varphi_1(z)=\frac{e^z-1}{z}.\]
This is exact for the chosen linear equation when \(\mathcal N=0\). For nonlinear dynamics it is generally first order in time. Accurately computing \(e^{hL}\) does not make the frozen cubic exact.
Add a second nonlinear evaluation
Use the predictor \(a\) above, evaluate \(\mathcal N(a)\), and linearly interpolate the nonlinear rate across the step. Integrating that interpolant gives exponential Heun, a second-order ETD method:
\[c^{n+1}=a+h\varphi_2(hL)\bigl[\mathcal N(a)-\mathcal N(c^n)\bigr],\] \[\varphi_2(z)=\frac{e^z-1-z}{z^2}.\]
There are two nonlinear evaluations per step. Each needs the same alias-free product strategy. This method is ETD2; it is not ETDRK4. Kassam and Trefethen’s paper provides a deeper comparison of exponential and other integrators for stiff PDEs, including higher-order methods.
The apparent singularities are removable:
\[\varphi_1(0)=1,\qquad\varphi_2(0)=\tfrac12,\] \[\varphi_1(z)=1+\tfrac z2+\tfrac{z^2}{6}+\cdots,\qquad \varphi_2(z)=\tfrac12+\tfrac z6+\tfrac{z^2}{24}+\cdots.\]
Use expm1 and short Taylor expansions near zero. Directly subtracting nearly equal exponentials can destroy the coefficient accuracy. With \(D=1\), modes \(k=\pm1\) have \(L_k=0\) even though the PDE remains nonlinear.
Refine time without changing the spatial problem
The experiment uses the course’s initial field, \(K=7\), domain length \(2\pi\) and final time \(T=1\). It computes the cubic through direct coefficient convolution. Its reference uses classical RK4 with 2,048 small steps at the same K. Changing the number of ETD steps holds both the final time and the finite-dimensional spatial model fixed.
Interactive experiment / fixed Galerkin system, variable time step
Fixed K=7, domain length 2π and final time T=1. Current errors: ETD1 3.694e-3, ETD2 7.387e-5. Last-refinement orders 1.00 and 2.01. RK4 reference uses the same retained Galerkin system. Interactive controls appear when the experiment is ready.
Link to this setupDownload these values (CSV)
Direct-convolution browser model, not OpenPFC output. The reported field error is measured at 256 evenly spaced points. CSV contains both profiles and refinement values. The RK4 reference is a numerical solution of the finite Galerkin ODE, not an exact PDE solution.
Halving \(h\) should eventually reduce ETD1’s error by about 2 and ETD2’s by about 4. Those orders apply in the asymptotic regime for this smooth problem. They do not prove sufficient spatial resolution. The browser reference is also checked by step doubling in the mathematical tests; the native verifier independently repeats that check for its own RK4 oracle.
Stability, energy and accuracy ask different questions
The continuum equation is a gradient flow of
\[E[u]=\int_0^{2\pi}\left[\frac D2u_x^2+\frac{(u^2-1)^2}{4}\right]dx,\] \[\frac{dE}{dt}=-\int_0^{2\pi}u_t^2\,dx\leq0.\]
The exactly integrated Galerkin ODE also inherits this identity because its projection is orthogonal. ETD2 is a time approximation; the formula above does not guarantee that every finite ETD2 step decreases energy. Check observed energy and time convergence separately. Its nonlinear part remains explicit, so removing the linear diffusion restriction does not give unrestricted accuracy or unconditional nonlinear stability.
The mean follows \(d\langle u\rangle/dt=\langle u-u^3\rangle\). A mean change is expected here. Reusing the diffusion course’s mean-conservation assertion would test the wrong equation.
Checkpoint: what happens when L = 0?
For a mode with Lₖ = 0, does the ETD2 correction disappear?
No. Take the continuous limits of the phi functions. The predictor is \(a_k=c_k^n+h\mathcal N_k(c^n)\) and the corrected mode is \(c_k^{n+1}=c_k^n+\tfrac h2[\mathcal N_k(c^n)+\mathcal N_k(a)]\). That is the explicit trapezoidal/Heun update for that component. The nonlinear rates still couple it to all other retained modes.