PDE Toolbox reaction-diffusion model collapses to zero, but FDM works

I’m trying to simulate a two–species reaction–diffusion system (modified Thomas model, from Liu & Maini) in MATLAB.

The equations are:

with parameters

delta = 6;
D = 0.45;
alpha = 0.899;
beta = -0.91;
r2 = 2;
r3 = 3.5;
gamma = -alpha;

Finite difference code (works)

When I discretise with a simple finite–difference scheme (explicit Euler, Neumann BCs), I see the expected pattern formation:

%Uconc = deltat*(D*delta*LaplaceU + alpha*Ucurrent + Vcurrent - r2*Ucurrent.*Vcurrent...
%    - alpha*r3*Ucurrent.*Vcurrent.^2) + Ucurrent;
%Vconc = deltat*(delta*LaplaceV + gamma*Ucurrent + beta*Vcurrent + r2*Ucurrent.*Vcurrent...
%    +alpha*r3*Ucurrent.*Vcurrent.^2) + Vcurrent;

Domain: [0,50]×[0,50], mesh 50×50, deltat=0.02.

This produces spots/stripes as expected (matches paper results).

PDE Toolbox code (fails)

I tried to rewrite the same model using the PDE Toolbox:

model = createpde(2);
L = 200;  % square domain
R1 = [3 4 0 L L 0 0 0 L L]';
geometryFromEdges(model, decsg(R1));

generateMesh(model, 'Hmax', 10);

applyBoundaryCondition(model,'neumann', ...
   'Edge', 1:model.Geometry.NumEdges,'g',[0;0],'q',[0;0]);

specifyCoefficients(model, ...
   'm',0, 'd',1, ...
   'c',[D*delta; delta], ...
   'a',[0;0], ...
   'f',@(region,state)[
       alpha*state.u(1,:) + state.u(2,:) ...
       - r2*state.u(1,:).*state.u(2,:) ...
       - alpha*r3*state.u(1,:).*(state.u(2,:)).^2;
       gamma*state.u(1,:) + beta*state.u(2,:) ...
       + r2*state.u(1,:).*state.u(2,:) ...
       + alpha*r3*state.u(1,:).*(state.u(2,:)).^2 ] );

setInitialConditions(model,@(r)[2*rand(1,numel(r.x))-1; ...
                                 2*rand(1,numel(r.x))-1]);

tlist = linspace(0,1000,500);
results = solvepde(model,tlist);
finalU = results.NodalSolution(:,1,end);
finalV = results.NodalSolution(:,2,end);
figure(1); clf;
pdeplot(model, 'XYData', finalU, 'Mesh', 'off');
axis equal; title(sprintf('u at t = 1000'));

But in this FEM version, the solution always decays to the trivial homogeneous state (U=V=0), never forming patterns.

Can anyone spot some reasons that my setup in the PDE toolbox is not correct?

[Torsten]

All numerical methods rely on differentiable data. Since you start with purely random initial conditions, I wouldn't trust in any of the two solutions - neither with your FDM nor with the PDE Toolbox.

What kind of "patterns" in the solutions for u and v do you expect to see - especially since you start purely random ?

[Heath]

@Torsten thanks for respoinding! The patterns expected are the same as the one I have attatched (this is the concetration in U, and V is approximately the inverse of this pattern). This was generated using the FDM and uses the same style of random initial condition as I am using for PDE toolbox FEM. The patterns in the FDM reliably form this style with approximately uniform spots forming. This also agrees with the results found in this paper 212.pdf which the work is based off. For these parameter values Turing patterns form and the dominant wavelength from linear analysis can be found and shows that this spotted pattern with an almost uniform spacing is expected.

My reasoning for trying to get the FEM working is that I want to try and integrate the on a more complex domain than a square. I have started with the square as this is a simple test case to see if the FEM method is working as expecrted (which it currently isn't).

[Torsten]

Your mesh is very coarse - you should at least use a finer mesh to reproduce such small local patterns. Although I'm still not convinced that numerical solvers can give senseful results for random initial data ... I'm surprised the FEM solver converges at all.

Answer 1

[Ankita]

Hi Heath,

It appears that the finite-difference code generates the expected patterns, while the PDE Toolbox version does not. One possible factor could be how the reaction terms are specified in the PDE Toolbox.

In some setups, all reaction terms (both linear and nonlinear) may be placed in the 'f' parameter. However, within the PDE Toolbox, linear reaction terms (those involving only U or V, such as alpha*U, beta*V, V, or gamma*U) are typically included in the 'a' matrix, which can be a 2x2 matrix for a two-equation system. Nonlinear terms are generally placed in 'f'. Adjusting the assignment of linear and nonlinear terms in this way might help the PDE Toolbox simulation to behave more like the FD scheme and potentially allow pattern formation to emerge. For further details on specifying coefficients in the PDE Toolbox, can refer to the official documentation: https://www.mathworks.com/help/pde/ug/pde.pdemodel.specifycoefficients.html

Another aspect to consider is the initial conditions. If these are too close to zero or too homogeneous, and the system is not sufficiently perturbed, patterns may not develop. Introducing some random perturbation to the initial conditions could be beneficial.

It may also help to match the domain size and mesh density to the FD code, and using a finer mesh might improve accuracy.

Hopefully, some of these suggestions will help achieving the desired pattern.