1
$\begingroup$

I am trying to do a recursive numerical method to solve an ODE in Mathematica. Whatever I try, I cannot get the computation time down. I need to do over 1000 steps, but I cannot get even 20 steps within a reasonable amount of time. I have been reading and attempting "Memoization" and using nest, but they aren't doing well because I have three values. Can someone help me get better time on this? One of the recursive functions I have is as follows:

    Timing[
         h = .01;
         f[x_] := Dot[{{0, 1}, {-1, 0}}, x];
         wp[0] := {0, 1};
         wp[1] = wp[0] + h*f[wp[0]];
         wp[2] = wp[1] + h*f[wp[1]];
         wp[n_] := wp[n - 1] + (h/12) (23*f[wp[n - 1]] - 16*f[wp[n - 2]] + 5 f{wp[n - 3]]);
          ]
$\endgroup$

1 Answer 1

1
$\begingroup$

Memoization is the way to go here. Your problem is that when Mathematica sees wp[1000], it will rewrite it using your rule for wp[n_], and then it will try to evaluate wp[999], wp[998], and wp[997]. But each of those evaluations depend on, e.g., wp[996], and so Mathematica will end up evaluating wp[996] three times. This ends up exploding combinatorially: think about how many times Mathematica would evaluate wp[500], for example.

You need to tell Mathematica that once it evaluates wp[996], or any value of wp, then it's done: any time it sees wp[996] again, it should spit out the same value without thinking. This is what =, rather than := is for: when Mathematica sees =, it evaluates the right hand side once, and saves it as the value of the left hand side, never going through a long evaluation again.

So in summary: to use this "memoization" procedure on a rule like f[x_] := rhs, just write f[x_] := f[x] = rhs. The first time Mathematica sees f[x], it will evaluate rhs, then save that value as f[x] and use that value any time it sees f[x] in the future.

Timing[h = .01;
 f[x_] := Dot[{{0, 1}, {-1, 0}}, x];
 wp[0] := {0, 1};
 wp[1] = wp[0] + h*f[wp[0]];
 wp[2] = wp[1] + h*f[wp[1]];
 wp[n_] := 
  wp[n] = wp[
     n - 1] + (h/12) (23*f[wp[n - 1]] - 16*f[wp[n - 2]] + 
       5 f[wp[n - 3]]);
 wp[1000]
 ]

gives {0.03125, {-0.544073, -0.839153}}.

$\endgroup$
1
  • $\begingroup$ You have saved my life $\endgroup$ Oct 26, 2020 at 0:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.