0af3386Reaction-diffusion on spherical-harmonic surfacesJeremy Magland 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% lap_g is the Laplace-Beltrami operator of the surface, which the embedding
7% (gx, gy, gz) determines. The reaction is explicit on the grid; diffusion is
8% implicit, which on a general surface is no longer a divide — see `step`.
9%
10% Provided by the caller: synth/analys (the transforms), lam = l(l+1) per
11% coefficient, gx/gy/gz (the surface, on the grid) and Gx/Gy/Gz (the same
12% surface as spherical-harmonic coefficients), noise (the seeded perturbation),
13% niter (iterations of the implicit solve), and the parameters. Grid fields are
14% npts x 1; spectral fields are real 2 x nlm (row 1 real part, row 2
15% imaginary), so no complex arithmetic is needed. Each function returns the new
16% spectral state followed by the grid fields to display.
18function [U, V, u, v] = init(noise, a, b)
19 us = a + b;
20 vs = b / (us * us);
21 U = analys(us + noise);
22 V = analys(vs * ones(numel(noise), 1));
23 u = synth(U);
24 v = synth(V);
25end
27function [Un, Vn, u, v] = step(U, V, lam, gx, gy, gz, a, b, D1, D2, dt, niter)
28 u = synth(U);
29 v = synth(V);
30 uuv = u .* u .* v;
32 % Explicit reaction. Bu, Bv are the right-hand side of the implicit
33 % diffusion solve (I - dt*D*lap_g) Unew = B.
34 Bu = U + dt * analys(a - u + uuv);
35 Bv = V + dt * analys(b - uuv);
37 % Split the surface Laplacian as lap_g = lap_s + dlap, where lap_s is the
38 % round-sphere one. lap_s is diagonal in spherical-harmonic space with
39 % eigenvalues -lam, so (I - dt*D*lap_s) inverts in a single divide and the
40 % whole geometry sits inside dlap. That turns the implicit solve into
41 %
42 % Unew = (Bu + dt*D1*dlap(Unew)) ./ (1 + dt*D1*lam)
43 %
44 % which the loop below iterates from the round-sphere answer: preconditioned
45 % Richardson, with the exactly invertible round-sphere operator as the
46 % preconditioner. It converges while dt*D*dlap stays small against
47 % (I - dt*D*lap_s), which is what would keep the cost to a few transforms.
48 Un = Bu ./ (1 + (dt * D1) * lam);
49 Vn = Bv ./ (1 + (dt * D2) * lam);
51 for k = 1:niter
52 % ---- placeholder: dlap = lap_g - lap_s --------------------------------
53 % The geometry enters HERE and nowhere else. What belongs here is the
54 % Laplace-Beltrami operator of the embedding minus the round-sphere one,
55 % which needs the induced metric (from derivatives of gx, gy, gz) and
56 % surface derivatives of the field — transforms this project does not have
57 % yet. See "The geometry is not in the operator yet" in the README.
58 %
59 % Until then dlap is identically zero. Not a stand-in that happens to be
60 % small: it is exactly the round sphere, so this loop provably changes
61 % nothing, every iterate equals the line above, and the scheme is bit for
62 % bit the one turing-sphere runs. The surface is drawn, not solved on.
63 dLu = 0 * Un;
64 dLv = 0 * Vn;
65 % ----------------------------------------------------------------------
66 Un = (Bu + (dt * D1) * dLu) ./ (1 + (dt * D1) * lam);
67 Vn = (Bv + (dt * D2) * dLv) ./ (1 + (dt * D2) * lam);
68 end
69end