# Why are my recursive functions running so slowly? And how do I fix them?

I am working with making numerical solvers for a modeling class I am in, and I want to use the internal recursion handling that is built into mathematcia, but there seems to be a problem. NDSolve gives me very rapid results that match up with these, but they are so slow. Euler's method for 20 steps is taking nearly a minute to complete, and Runge-Kutta simply sits there for forever if I go above 6 steps. I know they should be faster, and I could do them as for loops easily enough, but I want to know how to make these work as quickly as they should be. I am sure that I am just missing something.

(*Euler's Method*)
a = 20;
h1 = 2./(a);
f1[t_, x_] := 2 x - 2 t^2 - 3;
t1[k_] := t1[k - 1] + h1;
y1[k_] := y1[k - 1] + h1*f1[t1[k - 1], y1[k - 1]];
t1[0] = 0;
y1[0] = 2;
y1[a]

(*Runge-Kutta*)
b = 6;
h2 = 2/b;
K1[k_] := h2*f2[t2[k], y2[k]];
K2[k_] := h2*f2[t2[k] + 1/2*h2, y2[k] + 1./2*K1[k]];
K3[k_] := h2*f2[t2[k] + 1/2*h2, y2[k] + 1/2*K2[k]];
K4[k_] := h2*f2[t2[k] + h2, y2[k] + K3[k]];
f2[t_, y_] := 2 y - 2 t^2 - 3;
t2[k_] := t2[k - 1] + h2;
y2[k_] := y2[k - 1] + 1/6 (K1[k - 1] + 2*K2[k - 1] + 2*K3[k - 1] + K4[k - 1]);
t2[0] = 0;
y2[0] = 2;
y2[b]

• In addition to @rasher's answer: I tried to compile some information on memoization here - you might find it useful. – Leonid Shifrin Apr 14 '14 at 23:07

"Memoize" your recursion (if you search the docs for "recursion", you'll see a hit for Functions That Remember Values They Have Found. Read it if the concept is new to you. )

(*Euler's Method*)
Timing[
a = 15;
h1 = 2./(a);
f1[t_, x_] := 2 x - 2 t^2 - 3;
t1[k_] := t1[k - 1] + h1;
y1[k_] := y1[k - 1] + h1*f1[t1[k - 1], y1[k - 1]];
t1[0] = 0;
y1[0] = 2;
y1[a]
]

Timing[
ClearAll[f1, t1, y1];
a = 15;
h1 = 2./(a);
f1[t_, x_] := f1[t, x] = 2 x - 2 t^2 - 3;
t1[k_] := t1[k] = t1[k - 1] + h1;
y1[k_] := y1[k] = y1[k - 1] + h1*f1[t1[k - 1], y1[k - 1]];
t1[0] = 0;
y1[0] = 2;
y1[a]
]

(*
{1.513210, 5.75543}
{0., 5.75543}
*)

• Thank you so much, this is exactly what I needed. The sims are running so much faster, and now I can see that I have an error I have to track down in my runge-kutta method, since it doesn't actually converge! Fun times today. – Adam Apr 16 '14 at 14:32

The reason your code is slow is that as written it will retrace the entire recursive tree at every step of the evaluation. One way around this is "memoization" as rasher showed. However, this is still not optimal as it stores a large number of definitions that are not required for the calculation. Instead, when possible, it is better to realize that your operation can be written iteratively because there are a fixed number of values used at each step of the computation. Here is how I would write it:

a = 20;
h1 = 2/a;

f1[t_, y_] := 2 y - 2 t^2 - 3;

f2[{t_, y_}] := {t + h1, y + h1*f1[t, y]};

Last @ Nest[f2, {0, 2}, a]

6.13312


This uses Nest to repeatedly apply a transformation function (f2) to a list of starting values: {0, 2}. These are the t, y values used by the f2 function and they represent the values of your y1 and t1 functions at each step.

Encapsulated as a stand-alone function, with f1 and f2 functions merged:

func[a_Integer?NonNegative] :=
Module[{h1, fn},
h1 = 2/a;
fn[{t_, y_}] := {t + h1, y + h1 (2 y - 2 t^2 - 3)};
Last @ Nest[fn, {0, 2}, a]
]


Unlike a recursive definition, even with memoization, this code is almost without an upper bound. By this I mean that if you try a = 300 with rasher's code (and default Mathematica settings) you will get $RecursionLimit::reclim: errors. If you set $RecursionLimit = Infinity and try a = 50000 you will crash Mathematica (at least version 7.) However func can be applied without problem:

func[50000]

7.99893


Or even larger:

func[5000000]

7.99999

• +1. This reminds me of this one. Might even be considered a dupe. – Leonid Shifrin Apr 15 '14 at 1:05
• @Leonid I missed that one. +1 on your answer there. – Mr.Wizard Apr 15 '14 at 1:53