Numerically solving system of four PDEs

[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!

 

[Strum]

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

 

[jssdenton]

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.

 

[Strum]

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?

 

[jssdenton]

Sure. Here is a link to a pdf of the article.

http://hopf.chem.brandeis.edu/members_content/yanglingfa/paper/t2.pdf

 

[Strum]

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).

 

[jssdenton]

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!

 

[Strum]

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

 

[Strum]

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).