1% Brusselator reaction-diffusion on a closed surface.
2%
3% du/dt = D1*lap_g(u) + A - (B+1)*u + u^2*v
4% dv/dt = D2*lap_g(v) + B*u - u^2*v
5%
6% Same scheme as models/schnakenberg.m, including the grouped transforms:
7% [a, b] = synth(x, y) runs the group as batched Legendre dispatches.
9function [U, V, u, v] = init(noise, A, B)
10 [U, V] = analys(A + noise, (B / A) * ones(numel(noise), 1));
11 [u, v] = synth(U, V);
12end
14function [Un, Vn, u, v] = step(U, V, lam, filt, gx, gy, gz, p1, p2, q2, r, A, B, D1, D2, dt, niter)
15 [u, v] = synth(U, V);
16 uuv = u .* u .* v;
18 ru = A - (B + 1) * u + uuv;
19 rv = B * u - uuv;
20 [Ru, Rv] = analys(ru, rv);
21 Bu = U + dt * Ru;
22 Bv = V + dt * Rv;
24 Un = Bu ./ (1 + (dt * D1) * lam);
25 Vn = Bv ./ (1 + (dt * D2) * lam);
27 for k = 1:niter
28 % dlap = lap_g - lap_s, evaluated at the current iterate in flux form
29 % (see models/schnakenberg.m, docs/richardson-iteration.md and
30 % docs/reduced-transforms.md for the derivation and the ordering).
31 Fu = Un .* filt;
32 Fv = Vn .* filt;
33 vtu = dthetac(Fu);
34 vpu = dphic(Fu);
35 vtv = dthetac(Fv);
36 vpv = dphic(Fv);
37 [Ftu, Fpu, Ftv, Fpv] = synth(vtu, vpu, vtv, vpv);
38 Pu = p1 .* Ftu + p2 .* Fpu;
39 Qu = p2 .* Ftu + q2 .* Fpu;
40 Pv = p1 .* Ftv + p2 .* Fpv;
41 Qv = p2 .* Ftv + q2 .* Fpv;
42 [PAu, QAu, PAv, QAv] = analys(Pu, Qu, Pv, Qv);
43 Pcu = PAu .* filt;
44 Qcu = QAu .* filt;
45 Pcv = PAv .* filt;
46 Qcv = QAv .* filt;
47 scu = dthetac(Pcu) + dphic(Qcu);
48 scv = dthetac(Pcv) + dphic(Qcv);
49 [Lu, Lv] = synth(scu, scv);
50 lapu = r .* Lu;
51 lapv = r .* Lv;
52 [LAu, LAv] = analys(lapu, lapv);
53 dLu = LAu + lam .* Un;
54 dLv = LAv + lam .* Vn;
56 Un = (Bu + (dt * D1) * dLu) ./ (1 + (dt * D1) * lam);
57 Vn = (Bv + (dt * D2) * dLv) ./ (1 + (dt * D2) * lam);
58 end
59end