7

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.

Dargor
  • 1,369
  • 9
  • 17

1 Answers1

9
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

Dr. belisarius
  • 115,881
  • 13
  • 203
  • 453
  • 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!

    – Dargor Mar 04 '15 at 02:25
  • @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 ... – Dr. belisarius Mar 04 '15 at 02:28
  • Perfecto, me pondré con ello inmediatamente! ¿Mi deficiente inglés es lo que ha delatado mi idioma nativo? – Dargor Mar 04 '15 at 02:30
  • @PabloGalindoSalgado Nope. Just your name ;) – Dr. belisarius Mar 04 '15 at 02:47
  • 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. – Dargor Mar 04 '15 at 02:55
  • @PabloGalindoSalgado Done. But I don't think that is a good way to "normalize" results – Dr. belisarius Mar 05 '15 at 17:44
  • It is not normalization, is a problem of time discretization. In a continuous case, the solution goles assimptotically to 0, but in this case the finite difference allow the solución to go under 0 and then to overflow. – Dargor Mar 05 '15 at 18:30
  • @PabloGalindoSalgado I would just ensure that (without deaths) the total pop in each cell remains constant and proportionally distributed to the I[n] and S[n] coefficients. Although perhaps the non-neg limit ought to be enforced. – Dr. belisarius Mar 05 '15 at 18:46