Differential Equations QuickStart Sample (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)