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)?
Compile
some of your functions. $\endgroup$Compile
is not going to be a gain. $\endgroup$f
, so overlooked that issue... Hopefully OP's true problem will not / can be made not to suffer from that. $\endgroup$