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

Question

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]

Answer 1

"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}
*)

Answer 2

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