# Numerically solving system of four PDEs

• jssdenton

#### jssdenton

Hi Forum,

I'm trying to use Mathematica to graphically explore a system of four PDEs, as defined in Yang et al. (2002).
Spatial Resonances and Superposition Patterns in a Reaction-Diffusion Model
with Interacting Turing Modes. Physical Review Letters 88(20). The equations are:

$\frac{\partial x_i}{\partial t} = D_{x_i} \nabla^2 x_i + \alpha(x_j - x_i) + f(x_i, y_i)$
$\frac{\partial y_i}{\partial t} = D_{y_i} \nabla^2 y_i + \alpha(y_j - y_i) + g(x_i, y_i)$

where

$f(x, y) = a - (1+b)x + x^2y$
$g(x,y) = bx - x^2y$

The reactive chemicals, x and y, are defined by subscripts i and j to determine which of the two layers they are in. Each reactive chemical has a diffusion coefficient, D.

I'm having trouble using NDSolve to generate the plots, however. I set up the constants:

a = 3;
myDu1 = 1.31;
myDv1 = 9.87;
myDu2 = 34;
myDv2 = 344.9;
(* bcrit=(1.+a*Sqrt[myDu/myDv])^2 *)
b = 6
x1min = -16*Pi;
x1max = 16*Pi;
y1min = x1min;
y1max = x1max; x2min = x1min;
x2max = x1max;
y2min = x1min;
y2max = x1max;
tmax = 1000;
f[x_, y_, d_, c_] := a + Sin[x*y] + Sin[d*c];
g[x_, y_, d_, c_] := b/a + Sin[x*y] + Sin[d*c];\.01

and then define the function for NDSolve:

(* Setup in NDSolve . Set \[Alpha] = \[Beta] = 1.0 *)

my2DBrusselator3 = NDSolve[
{

D[u1[x, y, d, c, t], t] ==
myDu1*(Laplacian[u1[x, y, d, c, t], {x, y, d, c}]) + (u2[x, y, d,
c, t] - u1[x, y, d, c, t]) + a - (b + 1)*u1[x, y, d, c, t] +
u1[x, y, d, c, t]^2*v1[x, y, d, c, t],

D[v1[x, y, d, c, t], t] ==
myDv1*(Laplacian[v1[x, y, d, c, t], {x, y, d, c}]) + (v2[x,
y, d, c, t] - v1[x, y, d, c, t]) + b*u1[x, y, d, c, t] -
u1[x, y, d, c, t]^2*v1[x, y, d, c, t],\.01

D[u2[x, y, d, c, t], t] ==
myDu2*(Laplacian[u2[x, y, d, c, t], {x, y, d, c}]) + (u1[x, y, d,
c, t] - u2[x, y, d, c, t]) + a - (b + 1)*u2[x, y, d, c, t] +
u2[x, y, d, c, t]^2*v2[x, y, d, c, t],

D[v2[x, y, d, c, t], t] ==
myDv2*(Laplacian[v2[x, y, d, c, t], {x, y, d, c}]) + (v1[x,
y, d, c, t] - v2[x, y, d, c, t]) + b*u1[x, y, d, c, t] -
u1[x, y, d, c, t]^2*v1[x, y, d, c, t],

u1[x1min, y, d, c, t] == f[x1min, y, d, c],
u1[x1max, y, d, c, t] == f[x1max, y, d, c],
u1[x, y1min, d, c, t] == f[x, y1min, d, c],
u1[x, y1max, d, c, t] == f[x, y1max, d, c],
u1[x, y, x2min, c, t] == f[x, y, x2min, c],
u1[x, y, x2max, c, t] == f[x, y, x2max, c],
u1[x, y, d, y2min, t] == f[x, y, d, y2min],
u1[x, y, d, y2max, t] == f[x, y, d, y2max],
v1[x, y, d, c, 0] == g[x, y, d, c],
v1[x1min, y, d, c, t] == g[x1min, y, d, c],
v1[x1max, y, d, c, t] == g[x1max, y, d, c],
v1[x, y1min, d, c, t] == g[x, y1min, d, c],
v1[x, y1max, d, c, t] == g[x, y1max, d, c],
v1[x, y, x2min, c, t] == g[x, y, x2min, c],
v1[x, y, x2max, c, t] == g[x, y, x2max, c],
v1[x, y, d, y2min, t] == g[x, y, d, y2min],
v1[x, y, d, y2max, t] == g[x, y, d, y2max],
u2[x1min, y, d, c, t] == f[x1min, y, d, c],
u2[x1max, y, d, c, t] == f[x1max, y, d, c],
u2[x, y1min, d, c, t] == f[x, y1min, d, c],
u2[x, y1max, d, c, t] == f[x, y1max, d, c],
u2[x, y, x2min, c, t] = f[x, y, x2min, c],
u2[x, y, x2max, c, t] = f[x, y, x2max, c],
u2[x, y, d, y2min, t] = f[x, y, d, y2min],
u2[x, y, d, y2max, t] = f[x, y, d, y2max],
v2[x, y, d, c, 0] == g[x, y, d, c],
v2[x1min, y, d, c, t] == g[x1min, y, d, c],
v2[x1max, y, d, c, t] == g[x1max, y, d, c],
v2[x, y1min, d, c, t] == g[x, y1min, d, c],
v2[x, y1max, d, c, t] == g[x, y1max, d, c],
v2[x, y, x2min, c, t] == g[x, y, x2min, c],
v2[x, y, x2max, c, t] == g[x, y, x2max, c],
v2[x, y, d, y2min, t] == g[x, y, d, y2min],
v2[x, y, d, y2max, t] == g[x, y, d, y2max]},

{u1, v1, u2, v2},

{x, x1min, x1max},
{y, y1min, y1max},
{d, x2min, x2max},
{c, y2min, y2max},
{t, 0, tmax},
PrecisionGoal -> 2]

But I get errors about boundary values being specifiable for only 1 independent variable. Could somebody help me figure out how to plot this system? It should look like Figure 3c from Yang et al. (2002): a mixture of small and large dots, of two size classes. I am a biologist, and have never had a course in PDEs!

Thanks!

What is c and d? ( yeah not much help :/ )

As I understand it, there should be 4 equations, with $x_1, y_1, x_2, y_2$. I was using x, y, d, and c to represent those functions because I was worried about bungling subscripts.

When I read your equations I thought ## x_{i} = x_{i}(t,x,y) ##, such that the ## x_{i} ## are functions of the spatial coordinates ## (x,y) ##. If so then you have 4 equation descretized on a 2d lattice, in which case you only need two coordinates. Perhaps this is not the case? Is it impossible to find a free copy of the article?

As I understand it, there should be 4 equations, with $x_1, y_1, x_2, y_2$. I was using x, y, d, and c to represent those functions because I was worried about bungling subscripts.

If you use x y d and c to represent the functions ## x_1, x_2, y_1, y_2 ##, then what is u1 u2 v1 and v2? As it looks now you laplacian terms seems to be the 4 dimensional laplacian, while it should only be 2( the spatial dimension).

I think I see what you are saying, but aren't there two sets of spatial coordinates in 2D--the x,y layer for the first system, and the x,y layer for the second system? I am not sure how to formulate the D[] functions, Laplacian, or initial/boundary conditions with this in mind.

I had tried to extend some code I found for a two-equation system, which used u[x,y,t] and v[x,y,t] with f[x_, y_] and g[x_, y_], but certainly bungled it!

No you are supposed to solve in 1 2d system. I don't usually use mathematica for these problems so I can't help you much with sýntax. However, the coupling terms in the v2 equation looks wrong. Some 1 index -> 2 index.
Do you know the steady states for the system? And do you know what they are actually plotting in the figure 3 in the article?

edit: Found the answer to the last question. It is just x1 they are plotting

Also do you know if they have made an expository article? It seems like the x1 and x2 takes care of the small and large scale structures, while the y fields are some inverse fields( regions of large value in x1 is a region of low value in y1 and so on).