Differential Equations in IronPython QuickStart Sample

Illustrates integrating systems of ordinary differential equations (ODE’s) in IronPython.

This sample is also available in: C#, Visual Basic, F#.

Overview

This QuickStart sample demonstrates how to solve ordinary differential equations (ODEs) using Numerics.NET’s ODE solver capabilities.

The sample illustrates three different ODE solvers:

The classic 4th order Runge-Kutta method with fixed step size
The Runge-Kutta-Fehlberg method (RKF45) with adaptive step size control
The CVODE solver from the SUNDIALS suite, which can efficiently handle both non-stiff and stiff problems

The sample uses these solvers on two test problems:

The Lorenz attractor - a nonlinear system of three coupled equations that exhibits chaotic behavior. This demonstrates basic usage of the solvers on a non-stiff system.
A flame expansion problem - a single nonlinear equation that becomes increasingly stiff as the initial condition approaches zero. This shows how CVODE’s implicit methods can efficiently handle stiff problems where explicit methods struggle.

The code shows how to:

Set up and configure different ODE solvers
Define systems of differential equations using delegate functions
Provide Jacobian matrices for stiff problems
Control integration accuracy through step size and tolerance settings
Access solver statistics and intermediate results

The code

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 Extreme Numerics.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)