5
$\begingroup$

Update: The question may seem too broad, but the answer provided by belisarius responds perfectly to the question, as the use of $ListConvolve$ and $NestList$ meets the desired criteria in in addition to the speed performance.

I'm currently simulating an infection (traditional Volterra infection agent simulation) in Mathematica. A very complicated but complete model can be found in the Wolfram blog

The set up is the following:

  • We have a $n\times n$ matrix that represent the world. In each box we have three numbers: infected, susceptible and dead.
  • The equations that govern the model are the following:

$S^{n+1}_{ij}= S^n_{ij} - \beta \left( I^n_{ij} S^n_{ij}+S^n_{i+1j}I^n_{i+1j}+S^n_{i-1j} I^n_{i-1j}+S^n_{ij-1} I^n_{ij-1}+S^n_{ij+1} I^n_{ij+1} \right) $ $I^{n+1}_{ij}= I^n_{ij} + \beta \left( I^n_{ij} S^n_{ij}+S^n_{i+1j}I^n_{i+1j}+S^n_{i-1j} I^n_{i-1j}+S^n_{ij-1} I^n_{ij-1}+S^n_{ij+1} I^n_{ij+1} \right) $

$D^{n+1}=S^0-S^{n+1}-I^{n+1}$

(With the initial conditions $I^0$ and $S^0$ to choose according pleases).

As you can see, for each time interval (n) you have to go through the matrix (i,j) in order to calculate (n+1), so it's a recursive thing. As you can see in the Wolfram blog, they use a combination of $Table$ (or $ParallelTable$) and recursive functions declared as $f[n,i,j]:=f[n-1,i,j]+...$.

The question is the following:

I want a functional-programming style approach to do this in order to be elegant and speed up the calculations. I was wondering if ListConvolve (as you can find in my other question ( Functional programming approach to avoid traditional loops ) can be used with NestList somehow.

The ideal features of the approach would be:

  1. No loops.
  2. Functional-programming style.
  3. Speed (Real time evaluation ?¿ )
  4. Elegant code.

As you can see, the core of the question is how to avoid recursive and procedural programing functions with loops that work very bad in Mathematica in order to speed up things. Of course the following is allowed:

  • Compiled functions.
  • Parallel evaluation.
  • Magic.

Thank you all for your time and i hope you will find this question interesting enough.

$\endgroup$
  • $\begingroup$ What's your trial? $\endgroup$ – xyz Mar 4 '15 at 1:15
  • $\begingroup$ Trial? You mean my current attempt? What i have done uses procedural programming stuff and is not useful. $\endgroup$ – Dargor Mar 4 '15 at 1:24
7
$\begingroup$
beta = .1;
m = 50;
SeedRandom[42];
{s, i} = RandomReal[{0, 1}, {2, m, m}];
newSI[{s_, i_}] := Clip[{s - #, i + #}] &[beta ListConvolve[CrossMatrix[1], ArrayPad[s i, 1]]]
ListAnimate[gr = GraphicsRow /@ Map[MatrixPlot, NestList[newSI, {s, i}, 50], {2}]]

enter image description here

$\endgroup$
  • $\begingroup$ This is sooooo good. I'm currently learning so, would be too much to ask if you can break-down the code? Sorry for the bothering and thank you so much! $\endgroup$ – Dargor Mar 4 '15 at 2:25
  • $\begingroup$ @PabloGalindoSalgado La mejor forma de aprender Mathematica es ejecutando una a una las funciones en tu máquina. No tengo problema en explicar algun detalle, pero funciones como ListConvolve[] me llevarían páginas ... $\endgroup$ – Dr. belisarius Mar 4 '15 at 2:28
  • $\begingroup$ Perfecto, me pondré con ello inmediatamente! ¿Mi deficiente inglés es lo que ha delatado mi idioma nativo? $\endgroup$ – Dargor Mar 4 '15 at 2:30
  • $\begingroup$ @PabloGalindoSalgado Nope. Just your name ;) $\endgroup$ – Dr. belisarius Mar 4 '15 at 2:47
  • $\begingroup$ Ok, code understood. Just some remark: As the population must remain between 0 and 1 i suggest to use Clip[{s - beta #, i + beta #}] to avoid numerical differences below 0 and up 1 to overflow. $\endgroup$ – Dargor Mar 4 '15 at 2:55

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.