Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Numerically solving system of four PDEs

  1. Jan 21, 2015 #1
    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!
     
  2. jcsd
  3. Jan 24, 2015 #2
    What is c and d? ( yeah not much help :/ )
     
  4. Jan 24, 2015 #3
    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.
     
  5. Jan 25, 2015 #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. Jan 25, 2015 #5
  7. Jan 25, 2015 #6
    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).
     
  8. Jan 25, 2015 #7
    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!
     
  9. Jan 26, 2015 #8
    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
     
  10. Jan 26, 2015 #9
    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).
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook




Similar Discussions: Numerically solving system of four PDEs
  1. System of PDEs (Replies: 2)

  2. Solving PDEs (Replies: 3)

  3. Solving PDEs? (Replies: 1)

Loading...