Data Analysis Mathematics Linear Algebra Statistics
New Version 8.1!

Supports .NET 6.0. Try it for free with our fully functional 60-day trial version. QuickStart Samples

# Differential Equations QuickStart Sample (IronPython)

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 = sigma * (y - y)
dy = y * (rho - y) - y
dy = y * y - beta * y
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

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 = y * y * (1 - y)
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) * y

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)```