This wikilog article is a draft, it was not published yet.

Julia programming

⇠ Back to Blog:Sandbox
m
m
 
(23 intermediate revisions by one user not shown)
Line 1: Line 1:
We now turn to an interesting ''coupled'' system of ODE, namely, the [https://en.wikipedia.org/wiki/Lotka%E2%80%93Volterra_equations Volterra-Lotka model of prey-predator] interactions. This is modelled with the pair of coupled ODE:
+
<center><wz tip="Still frames for an excitable medium with a=11 and g=2">[[File:Screenshot_20210315_164359.png|400px]]</wz></center>
  
\begin{align}
+
Here is the same but animated:
\frac{dx}{dt} &= \alpha x - \beta x y, \\
+
\frac{dy}{dt} &= \delta x y - \gamma y,
+
\end{align}
+
  
The functions are implemented as:
+
<center><wz tip="250 iterations of the a=11 g=2 excitable model.">[[File:excitable-medium-a2-b11.gif|400px]]</wz></center>
  
<syntaxhighlight lang="python">
+
Playing with the parameter space:
function f1(t,y1, y2)
+
    α*y1-β*y1*y2
+
end
+
function f2(t,y1, y2)
+
    δ*y1*y2-γ*y2
+
end
+
</syntaxhighlight>
+
  
Euler's method in this case reads:
+
<center><wz tip="Two (out of the many) configurations:">[[File:Screenshot_20210315_181551.png|400px]]</wz></center>
  
$$\begin{align}
+
And the bottom case, animated:
y_{1,n+1}&=y_{1,n}+hf_1(t_n,y_{1,n},y_{2,n})\\
+
y_{2,n+1}&=y_{2,n}+hf_2(t_n,y_{1,n},y_{2,n})
+
\end{align}$$
+
  
Heun's version reads:
+
<center><wz tip="250 iterations of the a=10 g=20 excitable model with a lot of initially infected cells in various stages, preventing extinction.">[[File:excitable-medium-a10-b20.gif|400px]]</wz></center>
 +
 
 +
It is interesting to focus on the surviving spot. We can zoom on this area and see what is going on there. The cells have locked into this stable pattern, involving, apparently, two cells only:
  
 
<syntaxhighlight lang="python">
 
<syntaxhighlight lang="python">
@time for n=1:npts-1
+
excited = @animate for i ∈ 1:30
    y1[n+1]=y1[n]+(h/2)*(f1((n-1)*h,y1[n],y2[n])+f1(n*h,y1[n]+h*f1((n-1)*h,y1[n],y2[n]),y2[n]))
+
          excite()
    y2[n+1]=y2[n]+(h/2)*(f2((n-1)*h,y1[n],y2[n])+f2(n*h,y1[n],y2[n]+h*f2((n-1)*h,y1[n],y2[n])))
+
          plotspace(145, 220, 225, 300)
end
+
      end
 
</syntaxhighlight>
 
</syntaxhighlight>
  
RK4's version reads:
+
<center><wz tip="A stable 2-cells pattern that survives in the g>a configuration by locking two neighbours with a period a+b=30.">[[File:stable-2cells-pattern.gif|400px]]</wz></center>
  
 +
It is interesting as well to change parameters during the simulation. This locks the patterns into some dislocated sequences which give rise to apparent cycles of stationary-evolution in space, which, however, are due to structures that formed under previous conditions of evolution. For instance, the following was obtained by iterating 50 iterations with parameters <tt>a=50 and g=25</tt> then changing <tt>g=5</tt> for 50 iterations and coming back to the initial <tt>g=25</tt> and cycling for ever:
  
 +
<center><wz tip="Strange oscillations in the wake of changing parameters in mid-flight. This is an interesting optical effect due to dephasing of the wavefronts caused by dislocations in the structure following various types of evolutions.">[[File:strange-oscillations.gif|400px]]</wz></center>
  
<syntaxhighlight lang="python">
+
There clearly appears to have two types of oscillations, one seemingly stationary in space whereas the over travels as wavefronts, that get periodically frozen. This is, however, merely an optical artifact, that can be well understood by looking at the state of the cell's evolution in time. Here is such a cut obtained in a similar configuration as the one displayed as a full density plot:
npts=10000
+
h=.001
+
 
+
y1=0.0*collect(1:npts);
+
y2=0.0*collect(1:npts);
+
 
+
y1[1]=1;
+
y2[1]=.1;
+
 
+
α=2/3;
+
β=4/3;
+
γ=1;
+
δ=1;
+
 
+
@time for i=1:npts-1
+
    y1[i+1]=y1[i]+h*f1((i-1)*h,y1[i],y2[i])
+
    y2[i+1]=y2[i]+h*f2((i-1)*h,y1[i],y2[i])
+
    end
+
</syntaxhighlight>
+
  
 
<syntaxhighlight lang="python">
 
<syntaxhighlight lang="python">
plot([[y1[i] for i=1:npts], [y2[i] for i=1:npts]])
+
for i=1:75
 +
  println(global iteration+=1); excite();
 +
  Plots.display(plot(space[[100],:]', ylims=(0,a+g), lw=2, legend=false))
 +
end
 
</syntaxhighlight>
 
</syntaxhighlight>
  
<syntaxhighlight lang="python">
+
<center><wz tip="Slice through the cellular automaton showing two types of oscillations, as a result of dislocations: drifting and collapsing.">[[File:wavefronts-excitable-medium-cut.gif|400px]]</wz></center>
plot([(y1[i], y2[i]) for i=1:npts])
+
</syntaxhighlight>
+
 
+
 
+
  
http://tutorial.math.lamar.edu/Classes/DE/Bernoulli.aspx
+
One can see how the dislocations cause this alternances of drifting or collapsing evolution of the cell states.
  
http://calculuslab.deltacollege.edu/ODE/7-C-3/7-C-3-h.html
+
Clearly there are many variations and configurations to explore, this is barely scratching the surface. Instead of delving further into this model, we turn to the most famous and also more impressive Hodgepodge machine.
  
Backward Euler method, or Implicit Euler method,
+
{{WLP6}}

Latest revision as of 21:51, 15 March 2021

Screenshot 20210315 164359.png

Here is the same but animated:

Excitable-medium-a2-b11.gif

Playing with the parameter space:

Screenshot 20210315 181551.png

And the bottom case, animated:

Excitable-medium-a10-b20.gif

It is interesting to focus on the surviving spot. We can zoom on this area and see what is going on there. The cells have locked into this stable pattern, involving, apparently, two cells only:

excited = @animate for i ∈ 1:30
           excite()
           plotspace(145, 220, 225, 300)
       end
a configuration by locking two neighbours with a period a+b=30.">Stable-2cells-pattern.gif

It is interesting as well to change parameters during the simulation. This locks the patterns into some dislocated sequences which give rise to apparent cycles of stationary-evolution in space, which, however, are due to structures that formed under previous conditions of evolution. For instance, the following was obtained by iterating 50 iterations with parameters a=50 and g=25 then changing g=5 for 50 iterations and coming back to the initial g=25 and cycling for ever:

Strange-oscillations.gif

There clearly appears to have two types of oscillations, one seemingly stationary in space whereas the over travels as wavefronts, that get periodically frozen. This is, however, merely an optical artifact, that can be well understood by looking at the state of the cell's evolution in time. Here is such a cut obtained in a similar configuration as the one displayed as a full density plot:

for i=1:75
   println(global iteration+=1); excite();
   Plots.display(plot(space[[100],:]', ylims=(0,a+g), lw=2, legend=false))
end
Wavefronts-excitable-medium-cut.gif

One can see how the dislocations cause this alternances of drifting or collapsing evolution of the cell states.

Clearly there are many variations and configurations to explore, this is barely scratching the surface. Instead of delving further into this model, we turn to the most famous and also more impressive Hodgepodge machine.