[genius] Thu Aug 14 00:54:57 2014 Jiri (George) Lebl <jirka 5z com>

[George Lebl <jirka src gnome org>]

commit 76ab135d44267717594cc059b3cd0b2b9bca07f9
Author: Jiri (George) Lebl <jiri lebl gmail com>
Date:   Thu Aug 14 00:55:00 2014 -0500

    Thu Aug 14 00:54:57 2014  Jiri (George) Lebl <jirka 5z com>
    
        * lib/equation_solving/diffyqs.gel: slight optimization for
          RungeKutta and RungeKuttaFull

 ChangeLog                        |    5 +++++
 lib/equation_solving/diffeqs.gel |   24 ++++++++++++++----------
 2 files changed, 19 insertions(+), 10 deletions(-)
---
diff --git a/ChangeLog b/ChangeLog
index 733e6bd..88384a2 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,8 @@
+Thu Aug 14 00:54:57 2014  Jiri (George) Lebl <jirka 5z com>
+
+       * lib/equation_solving/diffyqs.gel: slight optimization for
+         RungeKutta and RungeKuttaFull
+
 Wed Aug 13 18:41:27 2014  Jiri (George) Lebl <jirka 5z com>
 
        * src/graphing.c: fix using a single point in LinePlotDrawPoints if
diff --git a/lib/equation_solving/diffeqs.gel b/lib/equation_solving/diffeqs.gel
index 6a25599..184049d 100644
--- a/lib/equation_solving/diffeqs.gel
+++ b/lib/equation_solving/diffeqs.gel
@@ -66,13 +66,15 @@ function RungeKutta(f,x0,y0,x1,n) = (
        x := float(x0);
        y := float(y0);
        for k = 1 to n do (
-               k1 := f(x,y);
-               k2 := f(x + hover2,y + k1*hover2);
-               k3 := f(x + hover2,y + k2*hover2);
+               #k1 := f(x,y);
+               #k2 := f(x + hover2,y + k1*hover2);
+               #k3 := f(x + hover2,y + k2*hover2);
+               k3 := f(x + hover2,y + (k2 := f(x + hover2,y + (k1:=f(x,y))*hover2))*hover2);
                increment x by h;
-               k4 := f(x,y + k3*h);
+               #k4 := f(x,y + k3*h);
+               #y := y + (k1 + 2*(k2 + k3) + k4)*h/6.0
                # y could be a vector, so do not use increment
-               y := y + (1/6)*(k1 + 2*k2 + 2*k3 + k4)*h
+               y := y + (k1 + 2*(k2 + k3) + f(x,y + k3*h))*h/6.0
        );
        y
 )
@@ -94,13 +96,15 @@ function RungeKuttaFull(f,x0,y0,x1,n) = (
        out@(1,1) := x := float(x0);
        out@(1,2) := y := float(y0);
        for k = 2 to n+1 do (
-               k1 := f(x,y);
-               k2 := f(x + hover2,y + k1*hover2);
-               k3 := f(x + hover2,y + k2*hover2);
+               #k1 := f(x,y);
+               #k2 := f(x + hover2,y + k1*hover2);
+               #k3 := f(x + hover2,y + k2*hover2);
+               k3 := f(x + hover2,y + (k2 := f(x + hover2,y + (k1:=f(x,y))*hover2))*hover2);
                increment x by h;
-               k4 := f(x,y + k3*h);
+               #k4 := f(x,y + k3*h);
+               #y := y + (k1 + 2*(k2 + k3) + k4)*h/6.0;
                # y could be a vector, so do not use increment
-               out@(k,2) := y := y + (1/6)*(k1 + 2*k2 + 2*k3 + k4)*h;
+               out@(k,2) := y := y + (k1 + 2*(k2 + k3) + f(x,y + k3*h))*h/6.0;
                out@(k,1) := x;
        );
        out