Trying to generate Turing patterns for Brusselator equations

[roastedwater]

so my project is on reaction diffusion equations in 2d. i have been asked to reproduce known and published work to start off with and the system is the brusselator equations with diffusion terms. to put it simply my program doesn't work. now the parameter values for the equations are taken from existing work so that's no issue. i am putting up the code i have written. its a forward difference method. i haven't been able to code implicit or adi methods for nonlinear problems. can someone please look at the code and tell me where i am goin wrong? I am a real amateur when it comes to computing.

 

[jackmell]

That's a mess dude. The syntax looks all wrong for Mathematica and even so, don't code a package. Just code a straight mathematica notebook. This is what I suggest. First just show us what Brusselator you're trying to solve. Let me try:

[tex]\frac{\partial u}{\partial t}=u^2 v-2u+1/4\nabla^2 u[/tex]

[tex]\frac{\partial v}{\partial t}=u-u^2 v+1/4 \nabla^2 v [/tex]

If that's not it, then edit it and show us and give the values of all the parameters as well as the initial and boundary conditions. Then finally, scrap that hard-code approach and just use the build-in numerical DE solver in Mathematica called NDSolve. Now, often it takes some effort to code PDEs in NDSolve but it is much, much easier to work with than hard-coding a numerical algorithm. If you're not familiar with NDSolve, then just start by using it with a simple ODE: just type NDSolve in a Mathematica notebook, put the cursor over it and hit F1 to get help on it. Study the exampls. Then do a coupled system of ODEs, solve it, plot them, do a few more, then do say a heat equation, get that one working and plot the solution using Plot3D, then move up to a coupled system of PDEs that you know the solution of, then finally, try and code the brusselator.

I've never coded a brusselator before but here's my first attempt in Mathematica with just a bunch of sine and cosine for initial and boundary values. It's not much but it's a first attempt and there are some warning messages reported by Mathematica:

myBrusselator = NDSolve[
   {D[u[x, y, t], t] == 
     u[x, y, t]^2*v[x, y, t] - 
      2*u[x, y, t] + (1/4)*
       (D[u[x, y, t], x, x] + 
        D[u[x, y, t], y, y]), 
    D[v[x, y, t], t] == u[x, y, t] - 
      u[x, y, t]^2*v[x, y, t] + 
      (1/4)*(D[v[x, y, t], x, x] + 
        D[v[x, y, t], y, y]), 
    u[x, y, 0] == Sin[-(Pi*x + 2*Pi*y)], 
    u[0, y, t] == Sin[-(2*Pi*y + 3*t)], 
    u[2, y, t] == Sin[-(2*Pi + 2*Pi*y + 
        t/2)], u[x, 0, t] == 
     Sin[-(Pi*x + t)], u[x, 2, t] == 
     Sin[-(Pi*x + 4*Pi - t)], 
    v[x, y, 0] == Cos[3*Pi*x + Pi*y], 
    v[0, y, t] == Cos[Pi*y + t/2], 
    v[2, y, t] == Cos[6*Pi + Pi*y + t], 
    v[x, 0, t] == Cos[3*Pi*x + t/4], 
    v[x, 2, t] == Cos[3*Pi*x + 2*Pi - 
       t]}, {u, v}, {x, 0, 2}, 
   {y, 0, 2}, {t, 0, 1}]
   
   
Animate[DensityPlot[
   Evaluate[v[x, y, t]] /. 
    myBrusselator, {x, 0, 2}, {y, 0, 2}, 
   ColorFunction -> "CMYKColors", 
   PlotPoints -> 35, MeshFunctions -> 
    {#3 & , #3 & }, Mesh -> 25, 
   MeshStyle -> {Black, Dashed}], 
  {t, 0, 1}, AnimationRate -> 0.01]

Maybe someone can improve it to illustrate better Turing patterns.

 

[roastedwater]

thank you so much jackmell for your response. but there's a msicommunication here. the code is not a mathematica code. its a MATLAB code. and are you sure i can generate turing patterns using the initial conditions like the one you have set? correct me if i am wrong but i thought the initial conditions had to be small perturbations about steady state while you're taking the initial conditions to be trigonometric functions. again i am not sure about this and will definitely try the conditions you've taken. the brusselator equations i am using are as follows.

∂u/∂t= Du *∇^2 (u)+A+(B-1)*u+ u^2 *v
∂v/∂t= Dv *∇^2 (v)+bu- u^2 *v

the steady state in absence of diffusion is given by

u = A
v = B/A

the turing patterns are generated for a particular range of Du, Dv, A and B. the set of values i have used is taken from an example where these same values are shown to give the patterns. in this case a patterns of stripes is what should come.
In the code i have used d1 for Du and d2 for Dv. the code is messy as i have tried a lot things in combination and i didnt want to type in everytime. i was hoping one of the combinations would click. i am very unsure about how to set the boundary conditions for this problem. this code works perfectly well for the simple heat equation but not for the reaction diffusion equations. i used the variable delta for the perturbation about steady state. also p and q are intermediate values of the concentrations u and v. since the equations are coupled i store the updated values in different variables and then overwrite u and v. i have also tried an adi code that i borrowed off the web but which promptly blew up in 3 iterations. i wish i could just use a package for this but my guide is very particular about writing my own code.

 

[jackmell]

Ok. Sorry about thinking it was Mathematica code and I didn't realize you were required to write the code directly. Also, I'm not sure that code I wrote can generate Turing patterns although it does generate turing-looking patterns when I run it. I'm not experienced with the brusselator but pretty good with DEs in general and believe in "roughing it in" with something even if it's wrong, then cleaning it up to make it work so that's why I coded that above. Let me first make sure I have the equations correct. Is it:

[tex]\frac{\partial u}{\partial t}=A+(B-1)u+u^2 v+D_u \left(\frac{\partial^2 u}{\partial x^2}+\frac{\partial ^2 u}{\partial y^2}\right)[/tex]

[tex]\frac{\partial v}{\partial t}=bu-u^2 v+D_v \left(\frac{\partial^2 v}{\partial x^2}+\frac{\partial ^2 v}{\partial y^2}\right)[/tex]

So assume that's right. What values of A, B b, [itex]D_u[/itex] and [itex]D_v[/itex] are you using?

Also, you said you're not sure about the initial and boundary conditions. Just to get started, use any reasonable initial and boundary values and keep "playing" with them until you get something that works. For example, just run it in the unit square 0<x<1 and 0<y<1 for some initial profile for some time period say 20 seconds.

 

[roastedwater]

there was a typo in the equations i sent you. the second equation has the term Bu not bu. there is no b in the equations.
the equations are

du/dt = Du*(d2u/dx2 + d2u/dy2) + A - (B + 1)*u + v*u^2
dv/dt = Dv*(d2v/dx2 + d2v/dy2) + B*u - v*u^2

A = 3,
B = 9,
Du = 5;
Dv = 12;

The pattern I'm supposed to get is stripes. i think in the solved example is done over a 50 by 50 square grid. the stripe thickness is about 5 or so..
Also, the boundary conditions are periodic boundary conditions or no flux boundary conditions. initial conditions if i am not wrong are steady state conditions with slight perturbation. the perturbation can be random.
the steady state is given by
u = A, v = B/A
any small perturbation will do. the example i looked at had random perturbations. i am trying to attach the document I am following but there's some problem with the proxy on my university campus. if you open wiki page for reaction diffusion equation you will see the results for another set of equations called fitzhugh-nagumo equations. the initial condition is noisy. small random perturbations about steady state.

 

[Filip Larsen]

I haven't followed the discussion here, but in case you need parameters for stable Turing stripes, you may want to have a look at [1], particularly the paper Stability of Turing patterns in the Brusselator model by Pena and Perez-Garcia which looks interesting.


[1] http://www.google.com/search?q=Turing+patterns+in+the+Brusselator+model

Later: I just noticed the mentioned paper actually references a paper my professor wrote based on the results from my master thesis. Nice to know someone out there had a use for that

 

[roastedwater]

hello filip. thank you so much for your response. in fact i have read this particular paper once before. though i got quite lost halfway through the paper! its quite mathematically involved! i haven't even reached a stage where i can analyse equations in such a comprehensive manner. i am still trying to write a few codes to simulate behaviour by just taking some given values. i suppose that such an analysis is what will be required of me after this stage is over. the work i doing is also towards a master's thesis and if your master's thesis involved work of such complexity you have my deepest respect! :) i would love to read your thesis if there was any way you could send it across to me. an email maybe? please do respond.

 

[Filip Larsen]

i would love to read your thesis if there was any way you could send it across to me.

Unfortunately the only thing I have left of it now it now is a paper copy (and a stack of pattern printouts), both the LaTeX and C++ source files has been lost since then. Also, if I recall correctly, I analysed only the Gray-Scott and Lengyel-Epstein model, not the Brusselator. I did get my solver library to produce both Turing stripes, Turing hexagons, Hopf-oscillation, as well as mixed-mode states using finite differencing and ADI methods in an object oriented framework (though details are very blurry for me, even after a quick skim of my old report) so if you have specific questions to that I may be able to dig out some old and dusty pointers, but probably not much more than that. Sorry.

 

[roastedwater]

ahh.. that's too bad.. i was hoping to read some work at a masters level on this topic.. i don't know anyone in my university who works on this infact.. hm.. i am not sure where the problem is occurring but i have attached a code to this thread. the first post. i wrote it for the explicit method. the code for adi blew up. i found one off the web and tried to modify it for this but it didnt work. i am unsure as to what the initial conditions are and what the boundary conditions are. the document i am following says periodic or no flux boundary conditions. i don't know how to code that though. i have a feeling this is the problem. though again it could well be that i made some fundamental errors in the code. could you see the code and tell me where i am making any mistakes? that would be of great help! its a MATLAB code.

 

[Filip Larsen]

First some general pointers you probably already know:

Part of the process of making a thesis is to get an opportunity to tackle some of the complex issues that comes with such a thing, in this case including analyzing why the code you write does not work as intended. There may be many convoluted issues that muddles the picture, so you may want to take smaller steps and ensure that every step along the way is something you understand. If you are completely baffled it may help to back off a bit and try something more simple first. And underlying all this is of course the need to research (and learn) the various domains you are utilizing.

Regarding boundary conditions:

The specific choice of condition, in combination with the choice of differencing approximation for the Laplacian, gives the differencing expression for points on boundary, that is, it determines the algebraic equations you end up having to solve. For instance, if c₀, c₁, c₂, ... is a sample line with c₀ being on the boundary, a second order central approximation for the Laplacian at c₁ is (c₀-2c₁+c₂)/h², with h being the sample distance. With a boundary condition you have an approximation for c₀ as a function of other samples in the domain, say c₁ and c₂, which means that c₀ is no longer a free variable in your algebraic equations but can be "moved to the right hand side". With other approximations get other equations.

Periodic conditions are when the, usually rectangular, reactor state space is "placed" on a torus so that the state space "wraps around" to the other end of that particular dimension. For instance c_-1 = c_m and c_m+1 = c₀, where m is the right-most sample. This "geometric" relationship is then used when you write out the differencing operators. In connection with ADI, which gives a tridiagonal systems of equations of order N, periodic boundary conditions means the system now becomes a cyclic tridiagonal system of order N which can be reduced to two tridiagonal systems of order N-1.

No-flux boundary condition is a Neuman condition with zero flux, that is the condition c'₀ = (dc/dx)₀ = 0. In our example, that would give a tree-point approximation of c₀ = (4/3)c₁-(1/3)c₂, if I read my old notes correctly.

Regarding initial conditions:

There are probably many ways to make a sensible initialization. If I recall correctly, I used some kind of random concentrations (within a range, of course) initially. Later when examining the stability of patterns, starting with the pattern overlaid with a bit of noise seems like a good idea.

Note, that starting from a purely homogeneous concentration (c_i = c_init for all i) probably is going to produce "degenerate" results without spatial patterns.

 

[jackmell]

Guys, I messed up with this ok by jumping the gun and going directly for the Turing patterns without first studying the dynamics with a simpler model. I'd like to try and correct my mistake please as jumping the gun like that is really not the way to do math. First I should point out the main feature of a Turing instability as I understand it now:

Certain reaction-diffusion systems can bifurcate: perturbed from their equilibrium state the system will either return to this state or settle into a new state. This is determined by a bifurcation parameter which in the case of the brusselator is the value of b. For b<b_c, perturbations around the equilibrium point return to the equilibrium point but for b>b_c, the state evolves into a new state. For the brusselator, the bifurcation point is given by the expression:

[tex]b_c=(1+a\sqrt{D_u/D_v})^2[/tex]

The value of [tex]b_c[/tex] can be determined by linearizing the system similar to the method used to linearize coupled ODE systems. I was not aware this could be done with PDEs until now. See for example:

http://www.teemuleppanen.net/content/science/turing/brusselator/instructions.pdf

With that in mind, I'd like to scale-back my attempt at this with a 1-D model:[tex]\frac{\partial u}{\partial t}=8 \frac{\partial^2 u}{\partial x^2}+4.5-(b+1)u+u^2 v[/tex]

[tex]\frac{\partial v}{\partial t}=16 \frac{\partial^2 v}{\partial x^2}+bu-u^2v[/tex]

with [itex]b_c\approx 17.5[/itex]

Also, I'm beginning to believe the boundary conditions are not so important so in my first attempt to better understand this system below, I keep the boundaries at the equilibrium points [itex]u_0=a[/itex] and [itex]v_0=b/a[/itex]. For the initial conditions, I just use a sine wave oscillating around the equilibrium points of both u and v:

myDu = 8; 
myDv = 16; 
a = 4.5; 
bcrit = (1 + a*Sqrt[myDu/myDv])^2
b = 17; 
xmin = -6*Pi; 
xmax = 6*Pi; 
tmax = 100; 
u0 = a; 
v0 = b/a; 

myBrusselator = NDSolve[{D[u[x, t], t] == myDu*D[u[x, t], x, x] + 
      a - (b + 1)*u[x, t] + u[x, t]^2*v[x, t], 
    D[v[x, t], t] == myDv*D[v[x, t], x, x] + b*u[x, t] - 
      u[x, t]^2*v[x, t], u[x, 0] == a + Sin[2*x], u[xmin, t] == a, 
    u[xmax, t] == a, v[x, 0] == b/a + Sin[2*x], v[xmin, t] == b/a, 
    v[xmax, t] == b/a}, {u, v}, {x, xmin, xmax}, {t, 0, tmax}, 
   MaxStepSize -> 0.1]
   
Plot3D[u[x, t] /. myBrusselator, {x, xmin, xmax}, {t, 0, tmax}, 
  PlotRange -> All, PlotPoints -> 50]

In the first u-plot below, b=17, just below the bifurcation point. Note how the sine-wave perturbation is quickly quenched as the value of u(x,t) returns to it's equilibrium point a=4.5. In the second plot I set b=18, just above the critical point and note how u(x,t) settles into a new state although at this point I don't know if u(x,t) would eventually return to it's equilibrium point at u(x,t)=4.5.

 

[jackmell]

In case anyone is interested, here's the link to download the (matlab) code for that exercise I cited above:

http://www.teemuleppanen.net/content/science/turing/turing_learning.shtml

Here's my second attempt at a 2D model using Mathematica and I stress that to really grasp this thing, I really should be focusing on the underlying dynamics but for now, just kinda' playing around with it:

[tex]
\frac{\partial u}{\partial t}=5\left(\frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 u}{\partial y^2}\right)+3-(b+1)u+u^2 v,\quad -6\pi\leq \{x,y\}\leq 6\pi,\quad 0\leq t\leq 10
[/tex]

[tex]
\frac{\partial v}{\partial t}=12\left(\frac{\partial^2 v}{\partial x^2}+\frac{\partial^2 v}{\partial y^2}\right)+bu-u^2v
[/tex]

[tex]\begin{align*}
u(x,y,0)&=a+\sin(x y)&\quad v(x,y,0)&=b/a+\sin(x y) \\
u(-6\pi,y,t)&=a+\sin(-6\pi y) &\quad v(-6\pi,y,t)&=b/a+\sin(-6\pi y)
\\
u(6\pi,y,t)&=a+\sin(6\pi y) &\quad v(6\pi,y,t)&=b/a+\sin(6\pi y)\\
u(x,-6\pi,t)&=a+\sin(-6\pi x) &\quad v(x,-6\pi,t)&=b/a+\sin(-6\pi x)\\
u(x,6\pi,t)&=a+\sin(6\pi x) &\quad v(x,6\pi,t)&=b/a+\sin(6\pi x)
\end{align*}
[/tex]

Then [itex]b_c\approx 8.6[/itex]

In the first density plot of u(x,y,10) below, b=9. In the second density plot of u(x,y,10), b=12 using the following code although Mathematica had difficulty generating all this data.

a = 3; 
myDu = 5; 
myDv = 12; 
bcrit = (1. + a*Sqrt[myDu/myDv])^2
b = 9
xmin = -6*Pi; 
xmax = 6*Pi; 
ymin = xmin; 
ymax = xmax; 
tmax = 10; 
f[x_, y_] := a + Sin[x*y]; 
g[x_, y_] := b/a + Sin[x*y]; 

my2DBrusselator = NDSolve[{D[u[x, y, t], t] == 
     myDu*(D[u[x, y, t], x, x] + D[u[x, y, t], y, y]) + a - 
      (b + 1)*u[x, y, t] + u[x, y, t]^2*v[x, y, t], 
    D[v[x, y, t], t] == myDv*(D[v[x, y, t], x, x] + 
        D[v[x, y, t], y, y]) + b*u[x, y, t] - u[x, y, t]^2*v[x, y, t], 
    u[x, y, 0] == f[x, y], u[xmin, y, t] == f[xmin, y], 
    u[xmax, y, t] == f[xmax, y], u[x, ymin, t] == f[x, ymin], 
    u[x, ymax, t] == f[x, ymax], v[x, y, 0] == g[x, y], 
    v[xmin, y, t] == g[xmin, y], v[xmax, y, t] == g[xmax, y], 
    v[x, ymin, t] == g[x, ymin], v[x, ymax, t] == g[x, ymax]}, {u, v}, 
   {x, xmin, xmax}, {y, ymin, ymax}, {t, 0, tmax}]

DensityPlot[
 u[x, y, tmax] /. my2DBrusselator, {x, xmin, xmax}, {y, ymin, ymax}, 
 PlotPoints -> 50]

 

[jackmell]

I decided to resort to (direct) finite differences to help run this brusselator since Mathematica's NDSolve runs into problems primarily because it saves all the interim iterations which quickly exhaust memory for any decent run. Maybe there is a way to avoid that but I don't know how at present. With that in mind, I contacted Teemu Leppanen, the author of the exercise I cited above, and received the C source code. After I made the required modifications and edited the output file to be compatable with Mathematica's list format, I used ListDensityPlot and I was able to finally achieve what I would consider some resolution to this problem in the form of brusselator stripes and dots. The parameters I used was:

a=4.5
b=6.75 (stripes), 7.5 (dots)
D_u=2
D_v=16

and I ran the simulator 50,000 iterations to obtain the plots below.

RoastedWater, if you're still interested, I can help you edit Teemu's code to obtain the results above but maybe you've moved on. Anyway, it was a fun problem for me as I've been interested in biological evolution for a very long time. :)

 

[roastedwater]

hello jackmell. i wasnt able to log on till now as i had my semester end exams for the last couple of weeks. i would definitely like to know how the code in c works. i have never used c. the languages i know are fortran, basic and matlab. however, i followed up your links and obtained the code from the site you mentioned. it was very helpful. thank you very much. but i still am not able to get the same results that you have obtained (referring to your latest post, the stripes and the spots). i wrote a much neater code in MATLAB though after going through the c code. i will upload that code and the functions i have used in it. one issue that did get sorted out was the boundary conditions. i also tried ndsolve as you did but my mathematica says out of memory for any tmax over 50. also, the code obtained from the link you posted had a few gaps. i suppose since it was a code given out as an exercise the author wanted the students to fill it out and complete it. i will nonetheless attach the fresh code i wrote and we could maybe compare the two.

 

[roastedwater]

this is in continuation to the previous post. i attached the plots i obtained as well as the laplacian for the code i used. its taken from leppanen's code of course. i seem to have to gave obtained the reverse of what you did. i got stripes where you got spots ans spots where you got stripes.

 

[jackmell]

Ok. Good you got the results you were looking for. Maybe I have the stripes and spots reversed. I barely have a handle on this at present but plan to work on it more. Also, I'm interested in verifying the results. How do we know after 50 thousand Euler iterations, the final results satisfy the PDE? I know we get results qualitatively similar to published results but how do we know those results accurately reflect the evolution of these equations? In Mathematica, it's very easy to back-substitute NDSolve data into the PDE and directly verify the results. May take a little more effort to do so with the Teemu data. I'll try and figure out how to do that and plot the discrepancy.

 

[roastedwater]

ah.. that's an interesting question. i think since the finite difference scheme is stable, there would be an error but it will still largely satisfy the original pdes, the stability condition being D*ht/hx/hx being less that 0.25. this can vary with the boundary conditions but for this set of boundary conditions it stable. of course stability does not guarantee accuracy but it shows that the solution has a good chance of being accurate. but then again, these are first order methods and the errors would be larger than say a second or a third order scheme. infact, the book, numerical recipes by teukolsky warns that programmers resorting to higher order methods may not always get accurate results as this might introduce errors in the form of oscillations resulting from higher harmonics.

 

[jackmell]

I stored a range of (time) iterates for both u and v so that I can back-substitute the solution into the PDE system. I ran the simulation and waited for the last 10 iterates. Since dt=0.01 in this particular run and I ran 5000 iterates, I saved a time slice from 49.90 to 49.99. I then ran InterpolationFunction on both the u data and v data at t=49.95, and computed mydifFunction=(lhs-rhs) of the first PDE and plotted the results. The plot below suggests the results are accurate to about 5 decimal places although I'm not yet quite confident of this.

myuData = << "c:\\Documents and \
Settings\\Owner\\My \
Documents\\ufile.dat"; 

myvData = << "c:\\Documents and \
Settings\\Owner\\My \
Documents\\vfile.dat"; 

myu = Interpolation[myuData]
myv = Interpolation[myvData]

aval = 4.5; 
bval = 6.75; 

mydifFunction[x2_, y2_] := 
  D[myu[x, y, t], t] - 
    2*(D[myu[x, y, t], x, x] + 
      D[myu[x, y, t], y, y]) - aval + 
    (bval + 1)*myu[x, y, t] - 
    myu[x, y, t]^2*myv[x, y, t] /. 
   {x -> x2, y -> y2, t -> 49.95}

mydifTable = Table[{xval, yval, 
     Evaluate[mydifFunction[xval, 
       yval]]}, {xval, 1, 99, 10}, 
    {yval, 1, 99, 10}]; 

ListPlot3D[Flatten[mydifTable, 1], 
  PlotRange -> All]

 

[Pythagorean]

Interesting thread so far. I went back to the motivation of turing patterns, and found it in Turing's paper on morphogenesis in organisms. From there, I found this paper on Turing instability that I'm still working through.

http://www.math.lsa.umich.edu/~rauch/256/turingexample.pdf

 

[jackmell]

Hi Pythagorean. Thanks for that. Here's the link to his paper "Chemical Basis of Morphogenesis". Looks like a paper I would be interested in studying:

http://www.dna.caltech.edu/courses/cs191/paperscs191/turing.pdf

 

[roastedwater]

classic paper. that was the first paper i was asked to read as a part of my project. the amount of work contained in it is enormous! they didn't have computers to help them with numerics at the time. simply phenomenal!

 

[dennisb]

Hi guys. I've tried your code. How can it work on Octave?

 

[jackmell]

Dennis, what you're asking is a little awkward because first, this is an old thread and usually better to just start a new one and if necessary, make a reference to this one. Also, Octave? Is that a programming language? I suppose so. But this forum is mostly for the math although I did run ram-rod with code in here. What I would suggest you do, unless someone helps you here is to go down to the "Math and Science Software" sub forum below, start a new thread with a title such as:

Octave: How to code Brusselator

Then describe what you want to do and maybe make a cross reference to this thread.

Think you'll get a better chance at help that way.

 

[dennisb]

thank you jackmell. I am sorry. it seems that my question is a little odd. let me rephrase that.

my problem is the title of this thread. I am trying to generate turing patterns for brusselator equations, but with a different programming language which is Octave. It is quite similar to Matlab. I have my equations but i don't know how to implement it.

 

[nahumxhdez]

code in mathematica

Podrian proporcionarme el codigo de patrones de turing en mathematica


Soy estudiante de preparatoria, y estoy interesado en la modelacion de sistemas biologicos

Creo que seria bueno comenzar con algo en mathematica

thank you


My e-mail is nahumxhdez@gmail.com

 

[nahumxhdez]

Hi! guys

They could give me the code of Turing patterns in mathematica


I'm a high school student, and I'm interested in the modeling of biological systems


My email is nahumxhdez@gmail.com
I think it would be good to start with something in mathematica

 

[jackmell]

I'd like if I may, post a description of the Mathematica code I used to generate the Turing patterns above with a little bit of explanation:

We study:

$$\frac{\partial u}{\partial t}=A+(B-1)u+u^2 v+D_u \left(\frac{\partial^2 u}{\partial x^2}+\frac{\partial ^2 u}{\partial y^2}\right)=F(u,v,\nabla_u^2)$$

$$\frac{\partial v}{\partial t}=bu-u^2 v+D_v \left(\frac{\partial^2 v}{\partial x^2}+\frac{\partial ^2 v}{\partial y^2}\right)=G(u,v,\nabla_v^2)$$

with periodic boundary conditions, that is, over a torus. We use Euler's method:

$$u_{i+1}=u_i+\Delta t F(u_i,v_i,\nabla_{u_i}^2)$$
$$v_{i+1}=v_i+\Delta t G(u_i,v_i,\nabla_{v_i}^2)$$

which for stability reasons we let [itex]\Delta x=\Delta y=1[/itex]. See VonNeumann stability analysis for more info about this: http://en.wikipedia.org/wiki/Von_Neumann_stability_analysis).

The initial values for u(x,y,0) and v(x,y,0) are random points. For example, u is initialized with the code:

[tex]u=a+0.3*\text{RandomReal}[\{-0.5,0.5\},\{n,n\}][/tex]

For example, use a rectangular grid 64x64 and allow the edges to ``wrap around'', that is the last column of the array has for computational purposes, on it's right side, the first column of the array and likewise, the last row of the array has below it, the first row of the array, then we can implement the periodic boundary conditions. With a little work studying the RotateLeft and RotateRight functions in Mathematica, it's not difficult to show the forward-difference approximation for [itex]\nabla_u^2[/itex] with this ``wrap-around'' geometry is:

[tex]\text{RotateLeft}[u,\{1,0\}]+\text{RotateLeft}[u,\{0,1\}]+\text{RotateRight}[u,\{1,0\}]+\text{RotateRight}[u,\{0,1\}]-4*u[/tex]

In the code below, the variables u and v are 2D arrays so all the arithmetic is being done in array arithmetic. Anyone interested should be able to simply cut-and-paste the code into Mathematica.

n = 64; 
a = 4.5; 
(* b=7.5 for stripes and 6.75 for dots *)
b = 7.5; 
du = 2; 
dv = 16; 
dt = 0.01; 
totaliter = 10000; 
u = a + 0.3*RandomReal[{-0.5, 0.5}, {n, n}]; 
v = b/a + 0.3*RandomReal[{-0.5, 0.5}, {n, n}]; 
cf = Compile[{{uIn, _Real, 2}, {vIn, _Real, 2}, {aIn, _Real}, {bIn, _Real}, {duIn, _Real}, 
     {dvIn, _Real}, {dtIn, _Real}, {iterationsIn, _Integer}}, 
    Block[{u = uIn, v = vIn, lap, dt = dtIn, k = bIn + 1, 
      kern = N[{{0, 1, 0}, {1, -4, 1}, {0, 1, 0}}], du = duIn, dv = dvIn}, 
     Do[lap = RotateLeft[u, {1, 0}] + RotateLeft[u, {0, 1}] + RotateRight[u, {1, 0}] + 
          RotateRight[u, {0, 1}] - 4*u; u = u + dt*(du*lap + a - u*(k - v*u)); 
        lap = RotateLeft[v, {1, 0}] + RotateLeft[v, {0, 1}] + RotateRight[v, {1, 0}] + 
          RotateRight[v, {0, 1}] - 4*v; v = v + dt*(dv*lap + u*(b - v*u)); , {iterationsIn}]; 
      {u, v}]]; 
Timing[c1 = cf[u, v, a, b, du, dv, dt, totaliter]; ]
ListDensityPlot[c1[[1]]]