Differential Equations in IronPython QuickStart Sample
Illustrates integrating systems of ordinary differential equations (ODE’s) in IronPython.
View this sample in: C# Visual Basic F#
```Python 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 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) ```