- #1

- 5

- 0

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:

[itex] \frac{\partial x_i}{\partial t} = D_{x_i} \nabla^2 x_i + \alpha(x_j - x_i) + f(x_i, y_i) [/itex]

[itex] \frac{\partial y_i}{\partial t} = D_{y_i} \nabla^2 y_i + \alpha(y_j - y_i) + g(x_i, y_i) [/itex]

where

[itex] f(x, y) = a - (1+b)x + x^2y [/itex]

[itex] g(x,y) = bx - x^2y [/itex]

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!