From c3d210a756b3f99f7adadc8f2d35f19952ca73c9 Mon Sep 17 00:00:00 2001 From: inivekin Date: Sat, 20 Mar 2021 10:23:11 +0800 Subject: [PATCH] Adds runge-kutta implementation for n dimensions --- extra/math/runge-kutta/authors.txt | 1 + .../math/runge-kutta/examples/examples.factor | 52 +++++++++++++++ .../math/runge-kutta/runge-kutta-docs.factor | 13 ++++ .../math/runge-kutta/runge-kutta-tests.factor | 23 +++++++ extra/math/runge-kutta/runge-kutta.factor | 65 +++++++++++++++++++ extra/math/runge-kutta/summary.txt | 1 + extra/math/runge-kutta/tags.txt | 1 + 7 files changed, 156 insertions(+) create mode 100644 extra/math/runge-kutta/authors.txt create mode 100644 extra/math/runge-kutta/examples/examples.factor create mode 100644 extra/math/runge-kutta/runge-kutta-docs.factor create mode 100644 extra/math/runge-kutta/runge-kutta-tests.factor create mode 100755 extra/math/runge-kutta/runge-kutta.factor create mode 100644 extra/math/runge-kutta/summary.txt create mode 100644 extra/math/runge-kutta/tags.txt diff --git a/extra/math/runge-kutta/authors.txt b/extra/math/runge-kutta/authors.txt new file mode 100644 index 0000000000..17bb1288dc --- /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 index 0000000000..b79c559d44 --- /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 * - ; + +: ( -- 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 + { 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 * + * ; + +:: ( 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 { 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 index 0000000000..66c82b6c2c --- /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: +{ $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 index 0000000000..e171004e23 --- /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 index 0000000000..77b0cda6dc --- /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 >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 + +: ( 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 index 0000000000..82fff83cf1 --- /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 index 0000000000..ede10ab61b --- /dev/null +++ b/extra/math/runge-kutta/tags.txt @@ -0,0 +1 @@ +math -- 2.34.1