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.
8function [U, V, u, v] = init(noise, A, B)
9 U = analys(A + noise);
10 V = analys((B / A) * ones(numel(noise), 1));
11 u = synth(U);
12 v = synth(V);
13end
15function [Un, Vn, u, v] = step(U, V, lam, filt, gx, gy, gz, Vtx, Vty, Vtz, Vpx, Vpy, Vpz, A, B, D1, D2, dt, niter)
16 u = synth(U);
17 v = synth(V);
18 uuv = u .* u .* v;
20 Bu = U + dt * analys(A - (B + 1) * u + uuv);
21 Bv = V + dt * analys(B * u - uuv);
23 Un = Bu ./ (1 + (dt * D1) * lam);
24 Vn = Bv ./ (1 + (dt * D2) * lam);
26 for k = 1:niter
27 % dlap = lap_g - lap_s, evaluated at the current iterate (see
28 % models/schnakenberg.m and docs/richardson-iteration.md for the
29 % derivation).
30 Fu = Un .* filt;
31 Ftu = dtheta(Fu);
32 Fpu = dphi(Fu);
33 dux = Ftu .* Vtx + Fpu .* Vpx;
34 duy = Ftu .* Vty + Fpu .* Vpy;
35 duz = Ftu .* Vtz + Fpu .* Vpz;
36 cux = analys(dux) .* filt;
37 cuy = analys(duy) .* filt;
38 cuz = analys(duz) .* filt;
39 Ftcux = dtheta(cux);
40 Fpcux = dphi(cux);
41 Ftcuy = dtheta(cuy);
42 Fpcuy = dphi(cuy);
43 Ftcuz = dtheta(cuz);
44 Fpcuz = dphi(cuz);
45 lapu = Ftcux .* Vtx + Fpcux .* Vpx;
46 lapu = lapu + Ftcuy .* Vty;
47 lapu = lapu + Fpcuy .* Vpy;
48 lapu = lapu + Ftcuz .* Vtz;
49 lapu = lapu + Fpcuz .* Vpz;
50 dLu = analys(lapu) + lam .* Un;
52 Fv = Vn .* filt;
53 Ftv = dtheta(Fv);
54 Fpv = dphi(Fv);
55 dvx = Ftv .* Vtx + Fpv .* Vpx;
56 dvy = Ftv .* Vty + Fpv .* Vpy;
57 dvz = Ftv .* Vtz + Fpv .* Vpz;
58 cvx = analys(dvx) .* filt;
59 cvy = analys(dvy) .* filt;
60 cvz = analys(dvz) .* filt;
61 Ftcvx = dtheta(cvx);
62 Fpcvx = dphi(cvx);
63 Ftcvy = dtheta(cvy);
64 Fpcvy = dphi(cvy);
65 Ftcvz = dtheta(cvz);
66 Fpcvz = dphi(cvz);
67 lapv = Ftcvx .* Vtx + Fpcvx .* Vpx;
68 lapv = lapv + Ftcvy .* Vty;
69 lapv = lapv + Fpcvy .* Vpy;
70 lapv = lapv + Ftcvz .* Vtz;
71 lapv = lapv + Fpcvz .* Vpz;
72 dLv = analys(lapv) + lam .* Vn;
74 Un = (Bu + (dt * D1) * dLu) ./ (1 + (dt * D1) * lam);
75 Vn = (Bv + (dt * D2) * dLv) ./ (1 + (dt * D2) * lam);
76 end
77end