Numerically solving system of four PDEs

  • #1
5
0
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:

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

Answers and Replies

  • #2
105
4
What is c and d? ( yeah not much help :/ )
 
  • #3
5
0
As I understand it, there should be 4 equations, with [itex] x_1, y_1, x_2, y_2 [/itex]. I was using x, y, d, and c to represent those functions because I was worried about bungling subscripts.
 
  • #4
105
4
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?
 
  • #6
105
4
As I understand it, there should be 4 equations, with [itex] x_1, y_1, x_2, y_2 [/itex]. 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).
 
  • #7
5
0
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!
 
  • #8
105
4
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
 
  • #9
105
4
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).
 

Related Threads on Numerically solving system of four PDEs

  • Last Post
Replies
1
Views
2K
  • Last Post
Replies
2
Views
5K
  • Last Post
Replies
14
Views
4K
  • Last Post
Replies
1
Views
709
Replies
1
Views
818
  • Last Post
Replies
3
Views
6K
  • Last Post
Replies
1
Views
2K
  • Last Post
Replies
1
Views
3K
  • Last Post
Replies
0
Views
2K
  • Last Post
Replies
3
Views
3K
Top