: lorenz-dy/dt ( tx..n -- dy )
rest first3
- 28 swap - [ swap ] dip * swap - ;
+ 28 swap - swapd * swap - ;
: lorenz-dz/dt ( tx..n -- dz )
rest first3
:: rf-dz/dt ( tx..n alpha -- dz )
tx..n rest first3 :> ( x y z )
- 2 neg z * alpha x y * + * ;
+ -2 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 ]
runge-kutta-2-transform ;
: runge-kutta-4-transform ( rk3 tx..n delta -- tx..n' delta )
- [ [ swap ] dip [ v*n ] keep prefix v+ ] keep ;
+ [ swapd [ 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
+ swapd 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
+ spin
[ { [ ]
[ runge-kutta-2-transform ]
[ runge-kutta-3-transform ]
call( -- rk1 rk2 rk3 rk4 rk4 ) drop 4array
! make array of zeroes of appropriate length to reduce into
- dup first length>> 0 <repetition> >array
+ dup first length>> 0 <array>
! reduce the results to the estimated change for the timestep
[ v+ ] reduce