5
$\begingroup$

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)?

$\endgroup$
  • 1
    $\begingroup$ You can try Compile some of your functions. $\endgroup$ – Silvia Nov 15 '15 at 10:20
  • $\begingroup$ Could you advice more precisely? $\endgroup$ – Dmitri Nov 15 '15 at 14:30
  • 2
    $\begingroup$ @Silvia The numbers early on exceed machine double ranges so Compile is not going to be a gain. $\endgroup$ – Daniel Lichtblau Nov 15 '15 at 17:13
  • $\begingroup$ @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. $\endgroup$ – Silvia Nov 15 '15 at 17:29
3
$\begingroup$

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

$\endgroup$
  • $\begingroup$ In Table[v,tN] the brackets are missed: Table[v, {tN}]. $\endgroup$ – Dmitri Nov 15 '15 at 18:47
  • $\begingroup$ @Dmitri They are the same. $\endgroup$ – Silvia Nov 16 '15 at 1:49
  • 1
    $\begingroup$ Are you in v10? in v9 and v8 Table[v, mi] breaks the code. $\endgroup$ – xzczd Nov 28 '15 at 3:25
  • 2
    $\begingroup$ @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. $\endgroup$ – Silvia Nov 28 '15 at 15:27
  • $\begingroup$ @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. $\endgroup$ – Dmitri Dec 1 '15 at 10:11

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.