Adds runge-kutta implementation for n dimensions

diff --git a/extra/math/runge-kutta/authors.txt b/extra/math/runge-kutta/authors.txt
new file mode 100644 (file)
index 0000000..17bb128
--- /dev/null
+++ b/extra/math/runge-kutta/authors.txt
@@ -0,0 +1 @@
+Kevin Cope

diff --git a/extra/math/runge-kutta/examples/examples.factor b/extra/math/runge-kutta/examples/examples.factor
new file mode 100644 (file)
index 0000000..b79c559
--- /dev/null
+++ b/extra/math/runge-kutta/examples/examples.factor
@@ -0,0 +1,52 @@
+USING: kernel io accessors arrays math.runge-kutta sequences math math.matrices 
+       ui.gadgets ui.gadgets.charts ui.gadgets.panes ui.gadgets.charts.lines ui.theme ;
+IN: math.runge-kutta.examples
+
+: lorenz-dx/dt ( tx..n -- dx )
+    rest first2
+    swap - 10 * ;
+
+: lorenz-dy/dt ( tx..n -- dy )
+    rest first3
+    28 swap - [ swap ] dip * swap - ;
+
+: lorenz-dz/dt ( tx..n -- dz )
+    rest first3
+    [ * ] dip 8/3 * - ;
+
+: <lorenz> ( -- dx..n/dt delta tx..n t-limit )
+    { [ lorenz-dx/dt ] [ lorenz-dy/dt ] [ lorenz-dz/dt ] } 0.01 { 0 0 1 21/20 } 150 ;
+
+
+: lorenz. ( -- )
+    chart new { { -20 20 } { -20 20 } } >>axes
+    line new link-color >>color
+             <lorenz> <runge-kutta-4> { 0 3 } cols-except >>data
+    add-gadget
+    gadget. ;
+
+
+:: rf-dx/dt ( tx..n gamma -- dx )
+    tx..n rest first3 :> ( x y z )
+    y z 1 - x sq + * gamma x * + ;
+
+:: rf-dy/dt ( tx..n gamma -- dy )
+    tx..n rest first3 :> ( x y z )
+    x 3 z * 1 + x sq - * gamma y * + ;
+
+:: rf-dz/dt ( tx..n alpha -- dz )
+    tx..n rest first3 :> ( x y z )
+    2 neg z * alpha x y * + * ;
+
+:: <rabinovich-fabrikant> ( gamma alpha -- dx..n/dt delta tx..n t-limit )
+    gamma '[ _ rf-dx/dt ] gamma '[ _ rf-dy/dt ] alpha '[ _ rf-dz/dt ]
+    3array
+    0.01 { 0 -1 0 0.5 } 150 ;
+
+
+: rabinovich-fabrikant. ( -- )
+    chart new { { -2 2 } { -2 2 } } >>axes
+    line new link-color >>color
+             0.1 0.14 <rabinovich-fabrikant> <runge-kutta-4> { 0 3 } cols-except >>data
+             add-gadget
+    gadget. ;

diff --git a/extra/math/runge-kutta/runge-kutta-docs.factor b/extra/math/runge-kutta/runge-kutta-docs.factor
new file mode 100644 (file)
index 0000000..66c82b6
--- /dev/null
+++ b/extra/math/runge-kutta/runge-kutta-docs.factor
@@ -0,0 +1,13 @@
+USING: help.markup help.syntax ;
+IN: math.runge-kutta 
+
+HELP: <runge-kutta-4>
+{ $description "Simple runge-kutta implementation for generating 4th-order approximated solutions for a set of first order differential equations" }
+{ $examples
+    "A lorenz attractor is a popular system to model with this: "
+    { $code "USING: runge-kutta runge-kutta.examples ;" "lorenz." }
+    "note that the produced chart is a 2 dimensional representation of a 3 dimensional solution. "
+    "Similarly, the rabinovich-fabrikant system (stable alpha-gamma limit cycle): "
+    { $code "USING: runge-kutta runge-kutta.examples ;" "rabinovich-fabrikant." }
+} ;
+

diff --git a/extra/math/runge-kutta/runge-kutta-tests.factor b/extra/math/runge-kutta/runge-kutta-tests.factor
new file mode 100644 (file)
index 0000000..e171004
--- /dev/null
+++ b/extra/math/runge-kutta/runge-kutta-tests.factor
@@ -0,0 +1,23 @@
+USING: math math.runge-kutta sequences tools.test ;
+
+! 1-dimensional exponential increase of x with t
+! initial state at { t x } of { 0 1 } with time step of 1
+! x is always increasing by x so dx/dt is just the x value from the initial state
+! gradient at beginning of interval
+{ { 1 } } [ { 0 1 } 1 { [ second ] } (runge-kutta) ] unit-test
+! gradient at midpoint of interval using previous estimation
+{ { 3 } } [ { 1 } { 0 1 } 1 runge-kutta-2-transform { [ second ] } (runge-kutta) ] unit-test
+! again gradient at midpoint of interval using previous estimation
+{ { 5 } } [ { 3 } { 0 1 } 1 runge-kutta-3-transform { [ second ] } (runge-kutta) ] unit-test
+! gradient at end of interval using previous estimation
+{ { 6 } } [ { 5 } { 0 1 } 1 runge-kutta-4-transform { [ second ] } (runge-kutta) ] unit-test
+
+! the summation of vectors for each stage, divided by 6, added to inital
+{ 7/2 } [ { 6 5 3 1 } 0 [ + ] reduce 6 / 1 + ] unit-test
+{ { 1 7/2 } } [ { [ second ] } 1 { 0 1 } (runge-kutta-4) ] unit-test
+
+! alters time dimension by delta
+{ { 1 1 2 3 } } [ 1 { 0 1 2 3 } increment-time ] unit-test
+
+! alters spatial dimensions by approximation
+{ { 1 3/2 7/3 13/4 } } [ { 1/2 1/3 1/4 } { 1 1 2 3 } increment-state-by-approximation ] unit-test

diff --git a/extra/math/runge-kutta/runge-kutta.factor b/extra/math/runge-kutta/runge-kutta.factor
new file mode 100755 (executable)
index 0000000..77b0cda
--- /dev/null
+++ b/extra/math/runge-kutta/runge-kutta.factor
@@ -0,0 +1,65 @@
+USING: kernel accessors fry combinators sequences sequences.generalizations
+       arrays math math.vectors math.functions ;
+IN: math.runge-kutta
+
+: runge-kutta-2-transform ( rk1 tx..n delta -- tx..n' delta )
+    [ swap [ [ v*n ] keep prefix 1/2 v*n ] dip v+ 2 v*n ] keep ;
+
+: runge-kutta-3-transform ( rk2 tx..n delta -- tx..n' delta )
+    runge-kutta-2-transform ;
+
+: runge-kutta-4-transform ( rk3 tx..n delta -- tx..n' delta )
+    [ [ swap ] dip [ v*n ] keep prefix v+ ] keep ;
+
+: (runge-kutta) ( delta tx..n dx..n/dt -- rk )
+    [ swap ] dip dup length>> [ cleave ] dip narray swap v*n
+    ; inline
+
+: runge-kutta-differentials ( dx..n/dt -- seq )
+    '[ _ (runge-kutta) ] ;
+
+: runge-kutta-transforms ( tx..n delta dx..n/dt -- seq )
+    -rot swap
+    [ { [ ]
+      [ runge-kutta-2-transform ]
+      [ runge-kutta-3-transform ]
+      [ runge-kutta-4-transform ] } ] dip
+    '[ _ runge-kutta-differentials compose ] map
+    [ 2curry ] 2with map ;
+
+: increment-time ( delta tx..n -- t+dtx..n )
+    [ swap [ 0 ] 2dip '[ _ + ] change-nth ] keep ; inline
+
+: increment-state-by-approximation ( rk4 t+dtx..n -- t+dtx'..n' )
+    swap 0 prefix v+ ;
+
+: (runge-kutta-4) ( dx..n/dt delta tx..n -- tx..n' )
+    [
+        ! set up the set of 4 equations
+        ! NOTE differential functions and timestep are curried
+        runge-kutta-transforms [ [ dup ] compose ] map
+
+        ! using concat instead causes slow down by an order of magnitude
+        [ ] [ compose ] reduce
+
+        ! this will always produce 4 outputs with a duplication of the last result
+        ! NOTE the dup is kept in the transform so that function
+        ! can be used in higher-order estimation
+        call( -- rk1 rk2 rk3 rk4 rk4 ) drop 4array 
+
+        ! make array of zeroes of appropriate length to reduce into
+        dup first length>> 0 <repetition> >array
+
+        ! reduce the results to the estimated change for the timestep
+        [ v+ ] reduce
+        1/6 v*n
+    ] 2keep
+    increment-time
+    increment-state-by-approximation
+    ;
+
+: time-limit-predicate ( t-limit -- quot: ( tx..n -- ? ) )
+    '[ dup first _ <= ] ; inline
+
+: <runge-kutta-4> ( dxn..n/dt delta initial-state t-limit -- seq )
+    time-limit-predicate [ [ (runge-kutta-4) ] [ 3drop f ] if ] compose 2with follow ; inline

diff --git a/extra/math/runge-kutta/summary.txt b/extra/math/runge-kutta/summary.txt
new file mode 100644 (file)
index 0000000..82fff83
--- /dev/null
+++ b/extra/math/runge-kutta/summary.txt
@@ -0,0 +1 @@
+runge-kutta 4-stage implementation for n dimensions

diff --git a/extra/math/runge-kutta/tags.txt b/extra/math/runge-kutta/tags.txt
new file mode 100644 (file)
index 0000000..ede10ab
--- /dev/null
+++ b/extra/math/runge-kutta/tags.txt
@@ -0,0 +1 @@
+math