# How to increase speed of iterations with nonlinear functions

Let me ask about part of a problem which I have in solving of a nonlinear DE (numerically). I am working with vectors of length $xN$. I have an initial vector $v$, a given $xN\times xN$ matrix $A$ and a given nonlinear function $f:\mathbb{R}^2\to\mathbb{R}$. Next, I need to construct the following two-dimensional list of vectors (I do not use here Mathematica's brackets just to simplify notations) using the rule $$w_{n,m}[k]=\begin{cases} v[k], & \text{if } n=0 \text{ or } m=1,\\ w_{n,m-1}[k]+f\Bigl[(A.w_{n-1,m})[k],w_{n-1,m}[k] \Bigr], &\text{otherwise}. \end{cases}$$ Here $k=1,\ldots,xN$ are vector indexes, $n=0,\ldots nN$ are iteration numbers, and $m=1,\ldots,tN$ are 'time'-steps. After all $nN$ iterations are done, I have to save the last iteration $$u_m[k]=w_{nN,m}[k], \qquad m=1,\ldots,tN, \quad k=1,\ldots,xN,$$ and I have to redefine the initial condition by the value of the last iteration at the last 'moment of time' $$v[k]=u_{tN}[k], \qquad k=1,\ldots,xN,$$ and repeat the same procedure to get now values for vectors $u_{tN+m}$, $m=1,\ldots,tN$, and so on until I find all vectors $u_{j}$ for $j=1,\ldots, tN\cdot tS$, for some chosen natural $tS$.

My realization is the following (here the matrix $A$ and the initial $v$ are random, whereas I need a particular; on the other hand, for some $f$ one can have an overflow, thus I use very simple $f$ just to show the slowness even in this case):

ClearAll["Global*"];
xN = 20;
tN = 50;
tS = 30;
nN = 20;
A = RandomReal[NormalDistribution[0, 1], {xN, xN}];
v = RandomReal[NormalDistribution[0, 1], {xN}];
f[p_?NumericQ, r_?NumericQ] := p + r;
rhs[vec_] := rhs[vec] = MapThread[f, {A.vec, vec}];
Do[
w[n_, m_] := w[n, m] = If[n == 0 || m == 1, v, w[n, m - 1] + rhs[w[n - 1, m]]];
Do[u[m + (i - 1) tN] = w[nN, m], {m, tN, 1, -1}];
v = w[nN, tN];
Clear[w],
{i, 1, tS}]; // AbsoluteTiming


It works, but for $xN=20$ it spends 18 sec, whereas for $xN=30$ almost 40 sec; unfortunately, I need $xN=200$ and less trivial $f$. (Note that the code above is very sensitive to $f$, if $f(p,r)=0.01(p+r)$ then everything is ten times faster, hence something is bad in the code).

Could you suggest me an improvement for the code (I am very newbie in Mathematica)?

• You can try Compile some of your functions. – Silvia Nov 15 '15 at 10:20
• Could you advice more precisely? – Dmitri Nov 15 '15 at 14:30
• @Silvia The numbers early on exceed machine double ranges so Compile is not going to be a gain. – Daniel Lichtblau Nov 15 '15 at 17:13
• @DanielLichtblau You're right! I actually didn't test on OP's f, so overlooked that issue... Hopefully OP's true problem will not / can be made not to suffer from that. – Silvia Nov 15 '15 at 17:29

Note your iteration rule can be reduced to an explicit formula:

$$\boldsymbol{\mathrm{w}}_{n,m}=\boldsymbol{v}+\sum_{k=2}^m f(\boldsymbol{\mathrm{A}}\cdot \boldsymbol{\mathrm{w}}_{n-1,k},\boldsymbol{\mathrm{w}}_{n-1,k})\mathrm{,\;\;(for\;}n>1\mathrm{)}$$

So we can calculate the $n$-th row of $\boldsymbol{\mathrm{w}}$ in one time by performing Accumulateon a list derived from the $(n-1)$-th row of $\boldsymbol{\mathrm{w}}$. By Nest-ing the calculation, we get the $\boldsymbol{\mathrm{u}}_m$.

Code is straightforward. Note I increased xN, tN, nN, and used a nonlinear $f(\alpha,\beta):=\exp(-\alpha^2)+\sin(\beta)$ here:

uFunc = Compile[{
{v, _Real, 1},  {A, _Real, 2},
{tN, _Integer}, {nN, _Integer}
},
Nest[
Accumulate[
Join[{v},
Function[{α, β},
Exp[-α^2] + Sin[β]
][#.Transpose[A], #] & @ Rest[#]
]
] &,
Table[v, tN], (* 1st row of w *)
nN + 1
]
]

xN = 200;
tN = 500;
nN = 200;
A = RandomReal[NormalDistribution[0, 1], {xN, xN}];
v = RandomReal[NormalDistribution[0, 1], {xN}];

AbsoluteTiming[uFunc[v, A, tN, nN];][[1]]


2.44941

• In Table[v,tN] the brackets are missed: Table[v, {tN}]. – Dmitri Nov 15 '15 at 18:47
• @Dmitri They are the same. – Silvia Nov 16 '15 at 1:49
• Are you in v10? in v9 and v8 Table[v, mi] breaks the code. – xzczd Nov 28 '15 at 3:25
• @xzczd I didn't notice that! Compare the first usage information of Table` in v10's doc and in v9's. I think they changed the syntax without notifying it at the top-right corner and it's a doc bug. – Silvia Nov 28 '15 at 15:27
• @Silvia Could you please also look here? This is a continuation of the question: I need to calculate your code not for tN=500, but say for tN=5, then the last row of the result needs to be saved and considered as a new value for vector v, and so on .. 50 times. I've realised this with a NestedList but I am still looking for an optimisation. – Dmitri Dec 1 '15 at 10:11