1% Schnakenberg reaction-diffusion on a closed surface.
2%
3% du/dt = D1*lap_g(u) + a - u + u^2*v
4% dv/dt = D2*lap_g(v) + b - u^2*v
5%
6% Explicit reaction, implicit diffusion (IMEX Euler): the step forms the
7% right-hand side B of the linear system (I - dt*D*lap_g) Unew = B, and hands
8% the solve to solvers/richardson.m — which applies the operator's geometric
9% part through lib/dlap.m. Trying a different solver for the same operator is
10% a change to these two call lines. Grid fields are npts x 1; spectral fields
11% are real 2 x nlm. See docs/richardson-iteration.md.
13function [U, V, u, v] = init(noise, a, b)
14 us = a + b;
15 vs = b / (us * us);
16 U = analys(us + noise);
17 V = analys(vs * ones(numel(noise), 1));
18 u = synth(U);
19 v = synth(V);
20end
22function [Un, Vn, u, v] = step(U, V, lam, filt, gx, gy, gz, Vtx, Vty, Vtz, Vpx, Vpy, Vpz, a, b, D1, D2, dt, niter)
23 u = synth(U);
24 v = synth(V);
25 uuv = u .* u .* v;
27 % Right-hand side of the implicit solve (I - dt*D*lap_g) Unew = B.
28 Bu = U + dt * analys(a - u + uuv);
29 Bv = V + dt * analys(b - uuv);
31 % The species diffuse independently, so each gets its own solve.
32 Un = richardson(Bu, dt * D1, lam, filt, Vtx, Vty, Vtz, Vpx, Vpy, Vpz, niter);
33 Vn = richardson(Bv, dt * D2, lam, filt, Vtx, Vty, Vtz, Vpx, Vpy, Vpz, niter);
34end