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

Julia programming

⇠ Back to Blog:Sandbox

Latest revision as of 21:51, 15 March 2021

Here is the same but animated:

Playing with the parameter space:

And the bottom case, animated:

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.">

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:

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

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.

Table of Contents

Crash Course on Scientific Computing

— • 21:51, 15 March 2021

@@ Line 1: / Line 1: @@
-We now turn to the case of coupled equations, in which case the method reads
+<center><wz tip="Still frames for an excitable medium with a=11 and g=2">[[File:Screenshot_20210315_164359.png|400px]]</wz></center>
-$$\vec y'=\vec f(t,\vec y)$$
+Here is the same but animated:
-where
+<center><wz tip="250 iterations of the a=11 g=2 excitable model.">[[File:excitable-medium-a2-b11.gif|400px]]</wz></center>
-$$
+Playing with the parameter space:
-\vec y\equiv
-\begin{pmatrix}
-  y_1\\y_2
-\end{pmatrix}
-\quad\text{and}\quad
-\vec f\equiv
-\begin{pmatrix}
-  f_1\\f_2
-\end{pmatrix}\,.
-$$
-The methods are the same as before just applying directly to the vectors instead. 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="Two (out of the many) configurations:">[[File:Screenshot_20210315_181551.png|400px]]</wz></center>
-\begin{align}
+And the bottom case, 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=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>
-<syntaxhighlight lang="python">
+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:
-function f1(t,y1, y2)
-    α*y1-β*y1*y2
-end
-function f2(t,y1, y2)
-    δ*y1*y2-γ*y2
-end
-</syntaxhighlight>
-Some common initialization for all methods, in particular, we now need two arrays of values to store the evolution of $x$ (in <tt>y1</tt>) and $y$ (in <tt>y2</tt>):
 <syntaxhighlight lang="python">
-tmax=100
+excited = @animate for i ∈ 1:30
-h=.001
+           excite()
-npts=round(Int,tmax/h)
+           plotspace(145, 220, 225, 300)
+       end
-y1=0.0*collect(1:npts);
-y2=0.0*collect(1:npts);
-y1[1]=1;
-y2[1]=.1;
-α=2/3;
-β=4/3;
-γ=1;
-δ=1;
 </syntaxhighlight>
-Euler's method in this case 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>
-$$\vec y_{n+1}=\vec y_n+h\vec f(t_n,\vec y_n)$
+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:
-or, breaking down this equation componentwise:
+<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>
-$$\begin{align}
+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:
-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}$$
-In code:
 <syntaxhighlight lang="python">
-@time for i=1:npts-1
+for i=1:75
-    y1[i+1]=y1[i]+h*f1((i-1)*h,y1[i],y2[i])
+   println(global iteration+=1); excite();
-    y2[i+1]=y2[i]+h*f2((i-1)*h,y1[i],y2[i])
+   Plots.display(plot(space[[100],:]', ylims=(0,a+g), lw=2, legend=false))
 end
 </syntaxhighlight>
-Heun's version reads:
+<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>
-$$\vec y_{n+1}=\vec y_n+{h\over 2}\left(\vec f(t_n,\vec y_n)+\vec f(t_{n+1},\vec y_n+h\vec f(t_i,\vec y_n)\right)$$
-or, breaking down this equation componentwise:
-$$\begin{align}
-  y_{1,n+1} &= y_{1,n} + {h\over2} [f(t_n, y_{1,n}, y_{2,n}) + f(t_n + h, y_{1,n} +  h f(t_n, y_{1,n}, y_{2,n}), y_{2,n})]\\
-  y_{2,n+1} &= y_{2,n} + {h\over2} [f(t_n, y_{1,n}, y_{2,n}) + f(t_n + h, y_{1,n}, y_{2,n} + h f(t_n, y_{1,n}, y_{2,n}))]
-\end{align}$$
-In code; where this time we'll find both more convenient and efficient to define auxiliary quantities, that are furthermore invoked more than once:
-<syntaxhighlight lang="python">
-@time for n=1:npts-1
-    tn=(n-1)*h;
-    f1left=yH1[n]+h*f1(tn,yH1[n],yH2[n])
-    f2left=yH2[n]+h*f2(tn,yH1[n],yH2[n])
-    f1right=f1(tn+h,f1left,f2left)
-    f2right=f2(tn+h,f1left,f2left)
-    fH=(h/2)*(f1left+f1right)
-    yH1[n+1]=yH1[n]+fH
-    yH2[n+1]=yH2[n]+fH
-end
-</syntaxhighlight>
-RK4's version reads:
-Results can be shown in two ways:
-<syntaxhighlight lang="python">
-plot([[y1[i] for i=1:npts], [y2[i] for i=1:npts]])
-</syntaxhighlight>
-<syntaxhighlight lang="python">
-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}}