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, jhat, 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 % Mean-J preconditioned solve (see models/schnakenberg.m), then iterate
39 % the geometric correction.
40 lamJ = lam ./ jhat;
41 Un = Bu ./ (1 + (dt * D1) * lamJ);
42 Vn = Bv ./ (1 + (dt * D2) * lamJ);
44 for k = 1:niter
45 % dlap = lap_g - lap_s, evaluated at the current iterate (Algorithm 3 of
46 % evolving_surface/notes/algos.tex): surface gradient of the field,
47 % contracted through the inverse metric quantities Vt*/Vp*; each
48 % Cartesian component re-analysed and differentiated again; recombined
49 % into the surface divergence. lam.*Un adds back -lap_s(Un), since lam
50 % holds +l(l+1). filt zeroes the top two degrees, where the theta/phi
51 % derivative recurrences cannot exactly represent a derivative. See
52 % docs/richardson-iteration.md.
53 Fu = Un .* filt;
54 Ftu = dtheta(Fu);
55 Fpu = dphi(Fu);
56 dux = Ftu .* Vtx + Fpu .* Vpx;
57 duy = Ftu .* Vty + Fpu .* Vpy;
58 duz = Ftu .* Vtz + Fpu .* Vpz;
59 cux = analys(dux) .* filt;
60 cuy = analys(duy) .* filt;
61 cuz = analys(duz) .* filt;
62 Ftcux = dtheta(cux);
63 Fpcux = dphi(cux);
64 Ftcuy = dtheta(cuy);
65 Fpcuy = dphi(cuy);
66 Ftcuz = dtheta(cuz);
67 Fpcuz = dphi(cuz);
68 lapu = Ftcux .* Vtx + Fpcux .* Vpx;
69 lapu = lapu + Ftcuy .* Vty;
70 lapu = lapu + Fpcuy .* Vpy;
71 lapu = lapu + Ftcuz .* Vtz;
72 lapu = lapu + Fpcuz .* Vpz;
73 dLu = (analys(lapu) + lamJ .* Un) .* filt;
75 Fv = Vn .* filt;
76 Ftv = dtheta(Fv);
77 Fpv = dphi(Fv);
78 dvx = Ftv .* Vtx + Fpv .* Vpx;
79 dvy = Ftv .* Vty + Fpv .* Vpy;
80 dvz = Ftv .* Vtz + Fpv .* Vpz;
81 cvx = analys(dvx) .* filt;
82 cvy = analys(dvy) .* filt;
83 cvz = analys(dvz) .* filt;
84 Ftcvx = dtheta(cvx);
85 Fpcvx = dphi(cvx);
86 Ftcvy = dtheta(cvy);
87 Fpcvy = dphi(cvy);
88 Ftcvz = dtheta(cvz);
89 Fpcvz = dphi(cvz);
90 lapv = Ftcvx .* Vtx + Fpcvx .* Vpx;
91 lapv = lapv + Ftcvy .* Vty;
92 lapv = lapv + Fpcvy .* Vpy;
93 lapv = lapv + Ftcvz .* Vtz;
94 lapv = lapv + Fpcvz .* Vpz;
95 dLv = (analys(lapv) + lamJ .* Vn) .* filt;
97 Un = (Bu + (dt * D1) * dLu) ./ (1 + (dt * D1) * lamJ);
98 Vn = (Bv + (dt * D2) * dLv) ./ (1 + (dt * D2) * lamJ);
99 end
100end