2
$\begingroup$

I have this simple Do Loop that, for some reason, takes a really long time for $h<\frac{1}{15}$. Why is this so? Is there a more efficient way in programming such a recursion?

h = 1/20;
v1[0] = 1;
v2[0] = 0;
Do[v1[j + h] = h v2[j] + v1[j]; 
 v2[j + h] = h (v2[j] + 1/v1[j]) + v2[j], {j, 0, 2, h}]
$\endgroup$
3
  • 2
    $\begingroup$ Yes, there is a simple solution. Replace exact numbers with approximate numbers, and it will be significantly faster. This means: replace v1[0] = 1; v2[0] = 0; with v1[0] = 1.; v2[0] = 0.; $\endgroup$
    – Domen
    Aug 4 at 12:37
  • 1
    $\begingroup$ To elaborate on @Domen's comment, your code uses exact numbers, and with say h = 1/15, your v2[2] is a rational number where numerator and denominator have millions of digits. $\endgroup$
    – user293787
    Aug 4 at 12:41
  • 1
    $\begingroup$ Are you trying to solve a differential equation? If yes, then a built-in method will be much more effective: NDSolve[{v1'[j] == v2[j], v2'[j] == v2[j] + 1/v1[j], v1[0] == 1, v2[0] == 0}, {v1, v2}, {j, 0, 2}] or even just NDSolve[{v1''[j] == v1'[j] + 1/v1[j], v1[0] == 1, v1'[0] == 0}, v1, {j, 0, 2}] $\endgroup$
    – Roman
    Aug 4 at 13:17

1 Answer 1

7
$\begingroup$

As @Domen says, use machine-precision numbers instead of exact rationals. Apart from this, I'd comment that there is no need to pre-compute the numbers in a Do loop: it is more Mathematica-style to define a memoizing recursion and then access the functions $v_1(j)$ and $v_2(j)$ randomly (without having to be explicit about the order of their evaluation):

h = 1/20;
Clear[v1, v2];
v1[0] = 1.;
v2[0] = 0.;
v1[j_ /; Divisible[j, h] && j > 0] := v1[j] = h v2[j - h] + v1[j - h];
v2[j_ /; Divisible[j, h] && j > 0] := v2[j] = h (v2[j - h] + 1/v1[j - h]) + v2[j - h];

ListPlot[Transpose@Table[{{j, v1[j]}, {j, v2[j]}}, {j, 0, 2, h}], PlotLegends -> {"v1", "v2"}]

enter image description here

The NDSolve curves in the above plot give the limit $h\to0^+$ and were calculated with

NDSolve[{v1''[j] == v1'[j] + 1/v1[j], v1[0] == 1, v1'[0] == 0}, v1, {j, 0, 2}]

Random access without explicit in-order precomputation:

v1[3]
(*    12.2335    *)

v2[7]
(*    629.188    *)
$\endgroup$

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.