= #CommonsMath_ODE_Lotka-Volterra~Test.jxn

! plots the polpulations r(t) of rabbits and f(t) of foxes

#import org.apache.commons.math3.ode
#import org.apache.commons.math3.ode.nonstiff

slider = @JxnSliderPanel($this)
tend = slider.add( "tend   ", 30., 0., 100. )
a = slider.add( "a", 1., 0.01, 10. )    ! natural growth rate of rabbits
b = slider.add( "b", 0.2, 0.01, 10. )   ! death rate of rabbits due to foxes
c = slider.add( "c", 1., 0.01, 10. )    ! natural death rate of foxes
d = slider.add( "d", 0.04, 0.01, 10. )  ! growth rate of foxes due to rabbits
r0 = slider.add( "r0", 50, 1, 1000 )    ! initial population of rabbits
f0 = slider.add( "f0", 10, 1, 1000 )    ! initial population of foxes

dp853 = @DormandPrince853Integrator( 1.E-8, 100., 0.0001, 0.0001 )
com = @ContinuousOutputModel()
dp853.addStepHandler( com )

! define the ODEs
rdot = "$y[0] * ( a - b * $y[1] );"  ! dr/dt for $y[0] rabbits
fdot = "$y[1] * ( d * $y[0] - c );"  ! df/dt and $y[1] foxes
ode = @JxnODE( $this, { rdot, fdot } )

t0 = 0.
y0 = { r0, f0 }
y  = { 0., 0. }
dp853.integrate( ode, t0, y0, tend, y )

! retrieve and plot output
tt = tend { 0 : 1000 } / 1000
rt = JxnCOMAnalyzer.analyze( com, tt, 0 )  ! r(t) rabbits
ft = JxnCOMAnalyzer.analyze( com, tt, 1 )  ! f(t) foxes
pf1 = plot( tt, rt, ft ).setAutoScale(2)
pf1.setHeadline( "Lotka-Volterra, red: rabbits, blue: foxes").setXYLabels( "t", "r(t), f(t)" )
pf2 = plot(     rt, ft ).setAutoScale(2)
pf2.setHeadline( "Lotka-Volterra: r(t) rabbits, f(t) foxes" ).setXYLabels( "r(t)", "f(t)" )