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 implicit solve
7% splits lap_g = lap_s + dlap: the round-sphere part lap_s is diagonal in
8% spherical-harmonic space (eigenvalues -lam), and the loop iterates the
9% geometric correction dlap from that exact solve. Grid fields are npts x 1;
10% spectral fields are real 2 x nlm. See docs/richardson-iteration.md.
11%
12% This is the 12-transform reference: the correction evaluates lap_g in
13% Cartesian-gradient form (Algorithm 4 of evolving_surface/notes/algos.tex),
14% carrying grad_g(u) as three ambient components through the inverse metric
15% quantities Vt*/Vp*. models/schnakenberg.m computes the same operator in
16% flux form with 6 transforms per species per iteration
17% (docs/reduced-transforms.md); this variant is kept live
18% for A/B comparison, in the app and in the tests.
20function [U, V, u, v] = init(noise, a, b)
21 us = a + b;
22 vs = b / (us * us);
23 U = analys(us + noise);
24 V = analys(vs * ones(numel(noise), 1));
25 u = synth(U);
26 v = synth(V);
27end
29function [Un, Vn, u, v] = step(U, V, lam, filt, gx, gy, gz, Vtx, Vty, Vtz, Vpx, Vpy, Vpz, a, b, D1, D2, dt, niter)
30 u = synth(U);
31 v = synth(V);
32 uuv = u .* u .* v;
34 % Right-hand side of the implicit solve (I - dt*D*lap_g) Unew = B.
35 Bu = U + dt * analys(a - u + uuv);
36 Bv = V + dt * analys(b - uuv);
38 % Round-sphere solve, then iterate the geometric correction.
39 Un = Bu ./ (1 + (dt * D1) * lam);
40 Vn = Bv ./ (1 + (dt * D2) * lam);
42 for k = 1:niter
43 % dlap = lap_g - lap_s, evaluated at the current iterate (Algorithm 3 of
44 % evolving_surface/notes/algos.tex): surface gradient of the field,
45 % contracted through the inverse metric quantities Vt*/Vp*; each
46 % Cartesian component re-analysed and differentiated again; recombined
47 % into the surface divergence. lam.*Un adds back -lap_s(Un), since lam
48 % holds +l(l+1). filt zeroes the top two degrees, where the theta/phi
49 % derivative recurrences cannot exactly represent a derivative. See
50 % docs/richardson-iteration.md.
51 Fu = Un .* filt;
52 Ftu = dtheta(Fu);
53 Fpu = dphi(Fu);
54 dux = Ftu .* Vtx + Fpu .* Vpx;
55 duy = Ftu .* Vty + Fpu .* Vpy;
56 duz = Ftu .* Vtz + Fpu .* Vpz;
57 cux = analys(dux) .* filt;
58 cuy = analys(duy) .* filt;
59 cuz = analys(duz) .* filt;
60 Ftcux = dtheta(cux);
61 Fpcux = dphi(cux);
62 Ftcuy = dtheta(cuy);
63 Fpcuy = dphi(cuy);
64 Ftcuz = dtheta(cuz);
65 Fpcuz = dphi(cuz);
66 lapu = Ftcux .* Vtx + Fpcux .* Vpx;
67 lapu = lapu + Ftcuy .* Vty;
68 lapu = lapu + Fpcuy .* Vpy;
69 lapu = lapu + Ftcuz .* Vtz;
70 lapu = lapu + Fpcuz .* Vpz;
71 dLu = analys(lapu) + lam .* Un;
73 Fv = Vn .* filt;
74 Ftv = dtheta(Fv);
75 Fpv = dphi(Fv);
76 dvx = Ftv .* Vtx + Fpv .* Vpx;
77 dvy = Ftv .* Vty + Fpv .* Vpy;
78 dvz = Ftv .* Vtz + Fpv .* Vpz;
79 cvx = analys(dvx) .* filt;
80 cvy = analys(dvy) .* filt;
81 cvz = analys(dvz) .* filt;
82 Ftcvx = dtheta(cvx);
83 Fpcvx = dphi(cvx);
84 Ftcvy = dtheta(cvy);
85 Fpcvy = dphi(cvy);
86 Ftcvz = dtheta(cvz);
87 Fpcvz = dphi(cvz);
88 lapv = Ftcvx .* Vtx + Fpcvx .* Vpx;
89 lapv = lapv + Ftcvy .* Vty;
90 lapv = lapv + Fpcvy .* Vpy;
91 lapv = lapv + Ftcvz .* Vtz;
92 lapv = lapv + Fpcvz .* Vpz;
93 dLv = analys(lapv) + lam .* Vn;
95 Un = (Bu + (dt * D1) * dLu) ./ (1 + (dt * D1) * lam);
96 Vn = (Bv + (dt * D2) * dLv) ./ (1 + (dt * D2) * lam);
97 end
98end