Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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]
share|improve this question
    
In addition to @rasher's answer: I tried to compile some information on memoization here - you might find it useful. –  Leonid Shifrin Apr 14 at 23:07

2 Answers 2

"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}
*)
share|improve this answer
    
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 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
share|improve this answer
1  
+1. This reminds me of this one. Might even be considered a dupe. –  Leonid Shifrin Apr 15 at 1:05
    
@Leonid I missed that one. +1 on your answer there. –  Mr.Wizard Apr 15 at 1:53

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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