Illustrates integrating systems of ordinary differential equations (ODE's) in IronPython.
import numerics
from Extreme.Mathematics import *
from Extreme.Mathematics.Calculus.OrdinaryDifferentialEquations import *
#/ Illustrates integrating systems of ordinary differential equations
#/ (ODE's) using classes in the
#/ Extreme.Mathematics.Calculus.OrdinaryDifferentialEquations
#/ namespace of the Extreme Optimization Mathematics Library for .NET.
#
# The ClassicRungeKuttaIntegrator class implements the
# well-known 4th order fixed step Runge-Kutta method.
rk4 = ClassicRungeKuttaIntegrator()
# The differential equation is expressed in terms of a
# DifferentialFunction delegate. This is a function that
# takes a double (time value) and two Vectors (y value and
# return value) as arguments.
#
# The Lorentz function below defines the differential function
# for the Lorentz attractor. Note that dy may be None on input:
def Lorentz(t, y, dy):
if dy == None:
dy = Vector.Create(3)
sigma = 10.0
beta = 8.0 / 3.0
rho = 28.0
dy[0] = sigma * (y[1] - y[0])
dy[1] = y[0] * (rho - y[2]) - y[1]
dy[2] = y[0] * y[1] - beta * y[2]
return dy
rk4.DifferentialFunction = Lorentz
# To perform the computations, we need to set the initial time...
rk4.InitialTime = 0.0
# and the initial value.
rk4.InitialValue = Vector.Create(1.0, 0.0, 0.0)
# The Runge-Kutta integrator also requires a step size:
rk4.InitialStepsize = 0.1
print "Classic 4th order Runge-Kutta"
for i in range(0,6):
t = 0.2 * i
# The Integrate method performs the integration.
# It takes as many steps as necessary up to
# the specified time and returns the result:
y = rk4.Integrate(t)
# The IterationsNeeded always shows the number of steps
# needed to arrive at the final time.
print "{0:.2f}: {1:18.14f} ({2} steps)".format(t, y, rk4.IterationsNeeded)
# The RungeKuttaFehlbergIntegrator class implements a variable-step
# Runge-Kutta method. This method chooses the stepsize, and so
# is generally more reliable.
rkf45 = RungeKuttaFehlbergIntegrator()
rkf45.InitialTime = 0.0
rkf45.InitialValue = Vector.Create(1.0, 0.0, 0.0)
rkf45.DifferentialFunction = Lorentz
rkf45.AbsoluteTolerance = 1e-8
print "Classic 4/5th order Runge-Kutta-Fehlberg"
for i in range(0,6):
t = 0.2 * i
y = rkf45.Integrate(t)
print "{0:.2f}: {1:18.14f} ({2} steps)".format(t, y, rkf45.IterationsNeeded)
# The CVODE integrator, part of the SUNDIALS suite of ODE solvers, # is the most advanced of the ODE integrators.
cvode = CvodeIntegrator()
cvode.InitialTime = 0.0
cvode.InitialValue = Vector.Create(1.0, 0.0, 0.0)
cvode.DifferentialFunction = Lorentz
cvode.AbsoluteTolerance = 1e-8
print "CVODE (multistep Adams-Moulton)"
for i in range(0,6):
t = 0.2 * i
y = cvode.Integrate(t)
print "{0:.2f}: {1:18.14f} ({2} steps)".format(t, y, cvode.IterationsNeeded)
#
# Other properties and methods
#
# The IntegrateSingleStep method takes just a single step
# in the direction of the specified final time:
cvode.IntegrateSingleStep(2.0)
# The CurrentTime property returns the corresponding time value.
print "t after single step:", cvode.CurrentTime
# CurrentValue returns the corresponding value:
print "Value at this t: {0:.14f}".format(cvode.CurrentValue)
# The IntegrateMultipleSteps method performs the integration
# until at least the final time, and returns the last
# value that was computed:
cvode.IntegrateMultipleSteps(2.0)
print "t after multiple steps:", cvode.CurrentTime
#
# Stiff systems
#
print "\nStiff systems"
# The CVODE integrator is the only ODE integrator that can
# handle stiff problems well. The following example uses
# an equation for the size of a flame. The smaller
# the initial size, the more stiff the equation is.
# We compare the performance of the variable step Runge-Kutta method
# and the CVODE integrator:
# The differential function for the flame expansion problem:
def Flame(t, y, dy):
if dy == None:
dy = Vector.Create(1)
dy[0] = y[0] * y[0] * (1 - y[0])
return dy
delta = 0.0001
tFinal = 2 / delta
rkf45 = RungeKuttaFehlbergIntegrator()
rkf45.InitialTime = 0.0
rkf45.InitialValue = Vector.Create(delta)
rkf45.DifferentialFunction = Flame
print "Classic 4/5th order Runge-Kutta-Fehlberg"
for i in range(0,11):
t = i * 0.1 * tFinal
y = rkf45.Integrate(t)
print "{0:8}: {1:18.14f} ({2} steps)".format(t, y, rkf45.IterationsNeeded)
# The CVODE integrator will use a special (implicit) method
# for stiff problems. To select this method, pass
# EdeKind.Stiff as an argument to the constructor.
cvode = CvodeIntegrator(OdeKind.Stiff)
cvode.InitialTime = 0.0
cvode.InitialValue = Vector.Create(delta)
cvode.DifferentialFunction = Flame
# When solving stiff problems, a Jacobian for the system must be
# supplied. The Jacobian is the matrix of partial derivatives of each
# equation in the system with respect to each variable in the system.
# <paramref name="J"/> may be <see langword="null"/>
# on input.</remarks>
def FlameJacobian(t, y, dy, J):
if J == None:
J = Matrix.Create(1, 1)
J[0, 0] = (2 - 3 * y[0]) * y[0]
return J
cvode.SetDenseJacobian(FlameJacobian)
# Notice how the CVODE integrator takes a lot fewer steps, # and is also more accurate than the variable-step
# Runge-Kutta method.
print "CVODE (Stiff - Backward Differentiation Formula)"
for i in range(0,11):
t = i * 0.1 * tFinal
y = cvode.Integrate(t)
print "{0:8}: {1:18.14f} ({2} steps)".format(t, y, cvode.IterationsNeeded)