Reaction Diffusion in two-dimensional space
Diffusion-induced instability

Schnakenberg kinetics (No diffusion, or spatially uniform solution)

Schnakenberg kinetics

fn1 = k1 - k2 u + k3 u^2 v;
fn2 = k4 - k3 u^2 v;

para = {k1->0.2, k2->1, k3->1, k4->0.5};

{fn1N, fn2N} = {fn1, fn2}/.para

g1=ContourPlot[fn1N==0,{u,0,3},{v,0,3},ContourStyle->RGBColor[0,0,1]]

g2=ContourPlot[fn2N==0,{u,0,3},{v,0,3},ContourStyle->RGBColor[1,0,0]]

Show[g1,g2]

fn1N

fn2N

FindRoot[{fn1N==0,fn2N==0},{u,0.5},{v,0.9}]

Solve[{fn1==0,fn2==0},{u,v}]

deq ={u'[t]==k1 - k2 u[t] + k3 u[t]^2 v[t],
	v'[t]==k4 - k3 u[t]^2 v[t],u[0]==u0,v[0]==v0}

para

deqN = deq/.{u0->0,v0->0}/.para

tEnd = 200;
sol = NDSolve[deqN,{u[t],v[t]},{t,0,tEnd}]

gPhase=ParametricPlot[ Evaluate[{u[t],v[t]}/.sol],{t,0,tEnd},PlotRange->{{0,3},{0,3}},AspectRatio->1,AxesLabel->{"u","v"}]

Plot[ Evaluate[{u[t],v[t]}/.sol],{t,0,tEnd},PlotRange->All,AxesLabel->{"t","u, v"}]

Show[g1,g2,gPhase]

2

Schnakenberg kinetics + Diffusion

Schnakenberg kinetics + Diffusion -- -> TURING PATTERN!!

SetDirectory["/Users/takasu/Library/Developer/Xcode/DerivedData/reaction_diffusion-drcxqrxbavdoubdxqelreislhgkm/Build/Products/Debug/"]

size1 = 101;
size2 = 51;

datau = ReadList["RD-2dim-schnakenberg-u",{Real,Real,Real}];
datav = ReadList["RD-2dim-schnakenberg-v",{Real,Real,Real}];

datau = Partition[datau,size1*size2];
datav = Partition[datav,size1*size2];

Length[datau]
Length[datav]

ListPlot3D[datau[[1]],AxesLabel->{"x","y","u"}]

ListPlot3D[datav[[1]],AxesLabel->{"x","y","v"}]

ListPlot3D[Last[datau],PlotRange->All,AxesLabel->{"x","y","u"}]

ListPlot3D[Last[datav],PlotRange->All,AxesLabel->{"x","y","v"}]

ListContourPlot[Last[datau],AspectRatio->0.5]

(* 3次元の図 ListPlot3D は描画に時間がかかるので、かわりに2次元の図　ListContourPlot で表す　*)

glist={};
Do[
g=ListContourPlot[ datau[[i]] ,PlotRange->{0,2},AspectRatio->0.5,Axes->True];
AppendTo[glist,g];
PrintTemporary[i],{i,1,Length[datau],1}
]

ListAnimate[glist]

Bi-stable dynamics (No diffusion, or spatially uniform solution)

Schnakenberg kinetics

fn1 = u - u^3 - v;
fn2 = epsilon(u - a1 v - a0);

para = {epsilon->0.05, a0->-0.1, a1->2.0};

{fn1N, fn2N} = {fn1, fn2}/.para

g1=ContourPlot[fn1N==0,{u,-1,1},{v,-1,1},ContourStyle->RGBColor[0,0,1]]

g2=ContourPlot[fn2N==0,{u,-1,1},{v,-1,1},ContourStyle->RGBColor[1,0,0]]

Show[g1,g2]

fn1N

fn2N

FindRoot[{fn1N==0,fn2N==0},{u,0.5},{v,0.9}]

FindRoot[{fn1N==0,fn2N==0},{u,-0.5},{v,-0.4}]

deq ={u'[t]==u[t]-u[t]^3-v[t],
	v'[t]==epsilon(u[t]-a1 v[t]-a0),u[0]==u0,v[0]==v0}

para

deqN1 = deq/.{u0->-0.5,v0->0}/.para

tEnd = 200;
sol1 = NDSolve[deqN1,{u[t],v[t]},{t,0,tEnd}]

gPhase1=ParametricPlot[ Evaluate[{u[t],v[t]}/.sol1],{t,0,tEnd},PlotRange->{{-1,1},{-1,1}},AspectRatio->1,AxesLabel->{"u","v"}]

Plot[ Evaluate[{u[t],v[t]}/.sol1],{t,0,tEnd},PlotRange->All,AxesLabel->{"t","u, v"}]

deqN2 = deq/.{u0->0.5,v0->0}/.para

tEnd = 200;
sol2 = NDSolve[deqN2,{u[t],v[t]},{t,0,tEnd}]

gPhase2=ParametricPlot[ Evaluate[{u[t],v[t]}/.sol2],{t,0,tEnd},PlotRange->{{-1,1},{-1,1}},AspectRatio->1,AxesLabel->{"u","v"}]

Plot[ Evaluate[{u[t],v[t]}/.sol2],{t,0,tEnd},PlotRange->All,AxesLabel->{"t","u, v"}]

Show[g1,g2,gPhase1,gPhase2]

2

Bi-stable dynamics + Diffusion

Schnakenberg kinetics + Diffusion -- -> TURING PATTERN!!

SetDirectory["/Users/takasu/Library/Developer/Xcode/DerivedData/reaction_diffusion-drcxqrxbavdoubdxqelreislhgkm/Build/Products/Debug/"]

size1 = 101;
size2 = 51;

datau = ReadList["RD-2dim-laby-u",{Real,Real,Real}];
datav = ReadList["RD-2dim-laby-v",{Real,Real,Real}];

datau = Partition[datau,size1*size2];
datav = Partition[datav,size1*size2];

Length[datau]
Length[datav]

ListPlot3D[datau[[1]],AxesLabel->{"x","y","u"}]

ListPlot3D[datav[[1]],AxesLabel->{"x","y","v"}]

ListPlot3D[Last[datau],PlotRange->All,AxesLabel->{"x","y","u"}]

ListPlot3D[Last[datav],PlotRange->All,AxesLabel->{"x","y","v"}]

ListContourPlot[Last[datau],AspectRatio->0.5]

(* 3次元の図 ListPlot3D は描画に時間がかかるので、かわりに2次元の図　ListContourPlot で表す　*)

glist={};
Do[
g=ListContourPlot[ datau[[i]] ,PlotRange->{0,2},AspectRatio->0.5,Axes->True];
AppendTo[glist,g];
PrintTemporary[i],{i,1,Length[datau],1}
]

ListAnimate[glist]