← All NMath Code Examples
Imports System
Imports System.Collections.Generic
Imports System.Diagnostics
Imports CenterSpace.NMath.Core
Namespace CenterSpace.NMath.Examples.VisualBasic
<summary>
A .NET Example in Visual Basic showing how to use the RungeKutta45OdeSolver
to solve a nonstiff set of equations describing the motion of a rigid body
without external forces:
y1= y2*y3, y1(0) = 0
y2= -y1*y3, y2(0) = 1
y3= -0.51*y1*y2, y3(0) = 1
</summary>
Module RungeKutta45OdeSolverExample
Sub Main()
Console.WriteLine()
Simple example solving the system of differential equations of
motion for a rigid body with no external forces. The equations
are defined above by the function Rigid().
Console.WriteLine("Simple Solve -----------------------------------")
Console.WriteLine()
SimpleSolve()
Console.WriteLine(Environment.NewLine)
A mass matrix and an output function are options that may
be specified for the solver. The output function is a callback
that is invoked at initialization, each integration step and
after the conclusion of the last integration step. Thus the
solution process can be monitored, and even terminated Imports
the output function.
Console.WriteLine("Output Function and Mass Matrix Solve ---------")
Console.WriteLine()
WithOutputFunctionAndMassMatrix()
Console.WriteLine(Environment.NewLine)
The above examples use the NMath class RungeKutta45OdeSolver. This
class uses an adaptive algorithm to figure out the optimal time step
size at each integration step. NMath also supplies a class,
RungeKutta5OdeSolver that does NOT use an adaptive step size
algorithm, allowing the user to specify the time values for
each integration step.
Console.WriteLine("Non-adaptive Step Size -------------------------")
Console.WriteLine()
NonAdaptiveStepSizeExample()
Console.WriteLine()
Console.WriteLine("Press Enter Key")
Console.Read()
End Sub
Sub SimpleSolve()
Construct the solver.
Dim Solver As New RungeKutta45OdeSolver()
Construct the time span vector. If this vector contains exactly
two points, the solver interprets these to be the initial and
final time values. Step size and function output points are
provided automatically by the solver. Here the initial time
value t0 is 0.0 and the final time value is 12.0.
Dim TimeSpan As New DoubleVector(0.0, 12.0)
Initial y vector.
Dim y0 As New DoubleVector(0.0, 1.0, 1.0)
Construct solver options. Here we set the absolute and relative tolerances to use.
At the ith integration step the error, e(i) for the estimated solution
y(i) satisfies
e(i) <= max(RelativeTolerance * Math.Abs(y(i)), AbsoluteTolerance(i))
The solver can increase the number of output points by a specified factor
called Refine (useful for creating smooth plots). The default value is
4. Here we set the Refine value to 1, meaning we do not wish any
additional output points to be added by the solver.
Dim SolverOptions As New RungeKutta45OdeSolver.Options()
SolverOptions.AbsoluteTolerance = New DoubleVector(0.0001, 0.0001, 0.00001)
SolverOptions.RelativeTolerance = 0.0001
SolverOptions.Refine = 1
Construct the delegate representing our system of differential equations...
Dim odeFunction As New Func(Of Double, DoubleVector, DoubleVector)(AddressOf Rigid)
...and solve. The solution is returned as a key/value pair. The first Keyelement of the pair is
the time span vector, the second Valueelement of the pair is the corresponding solution values.
That is, if the computed solution function is y then
y(soln.Key(i)) = soln.Value(i)
Dim Soln As RungeKutta45OdeSolver.Solution(Of DoubleMatrix) = Solver.Solve(odeFunction, TimeSpan, y0, SolverOptions)
Console.WriteLine("T = " & Soln.T.ToString("G5"))
Console.WriteLine("Y = ")
Console.WriteLine(Soln.Y.ToTabDelimited("G5"))
End Sub
<summary>
Function describing the system of differential equations.
</summary>
<param name="t">Time parameter.</param>
<param name="y">State vector.</param>
<returns>The vector of values of the derivatives of the functions
at a specified time t and function values y:
y1= y2*y3,
y2= -y1*y3,
y3= -0.51*y1*y2
</returns>
Function Rigid(ByVal T As Double, ByVal Y As DoubleVector) As DoubleVector
Dim dy As New DoubleVector(3)
dy(0) = Y(1) * Y(2)
dy(1) = -Y(0) * Y(2)
dy(2) = -0.51 * Y(0) * Y(1)
Return dy
End Function
<summary>
The output function to be called by the ODE solver at each integration step.
If this function returns true the solver continues with the next step, if
it returns false the solver is halted.
This output function just outputs to the console upon each invocation.
</summary>
<param name="t">The time value at the current integration step. If flag is equal to Initialize,
this will be the initial time value t0.</param>
<param name="y">The calculated function value at the current integration step. If flag is equal to Initialize,
this will be the initial function value y0.</param>
<param name="flag">Flag indicating what stage the solver is at:
Initialize - before solving begins. y and t values are the problems initial values.
IntegrationStep - just finished an integration step. y and t values are the values
calculated at that step.
Done - just finished last step. y and t values are the last values in the returned
solution.</param>
<returns>true if the solver is to proceed with the calculation, false forces the solver
to stop.</returns>
Function OutputFunction(ByVal T As DoubleVector, ByVal Y As DoubleMatrix, ByVal Flag As RungeKutta45OdeSolver.OutputFunctionFlag) As Boolean
If Flag = RungeKutta45OdeSolver.OutputFunctionFlag.Initialize Then
Console.WriteLine("Output function : Initialize")
Return True
End If
If Flag = RungeKutta45OdeSolver.OutputFunctionFlag.IntegrationStep Then
Console.WriteLine("Output function Integration step:")
Console.WriteLine("t = " + T.ToString("G5"))
Console.WriteLine("y = ")
Console.WriteLine(Y.ToTabDelimited("G5"))
Return True
End If
If Flag = RungeKutta45OdeSolver.OutputFunctionFlag.Done Then
Console.WriteLine("Output function: Done")
Return True
End If
Console.WriteLine("Output function: Unknown flag")
Return False
End Function
Sub WithOutputFunctionAndMassMatrix()
Construct the solver.
Dim Solver As New RungeKutta45OdeSolver()
Construct the time span vector. If this vector contains exactly
two points, the solver interprets these to be the initial and
final time values. Step size and function output points are
provided automatically by the solver. Here the initial time
value t0 is 0.0 and the final time value is 12.0.
Dim TimeSpan As New DoubleVector(0.0, 12.0)
Initial y vector.
Dim y0 As New DoubleVector(0.0, 1.0, 1.0)
Construct solver options. Here we set the absolute and relative tolerances to use.
At the ith integration step the error, e(i) for the estimated solution
y(i) satisfies
e(i) <= max(RelativeTolerance*Math.Abs(y(i)), AbsoluteTolerance(i))
Dim SolverOptions As New RungeKutta45OdeSolver.Options()
SolverOptions.AbsoluteTolerance = New DoubleVector(0.0001, 0.0001, 0.00001)
SolverOptions.OutputFunction = New Func(Of DoubleVector, DoubleMatrix, RungeKutta45OdeSolver.OutputFunctionFlag, Boolean)(AddressOf OutputFunction)
SolverOptions.RelativeTolerance = 0.0001
Construct the delegate representing our system of differential equations...
Dim odeFunction As New Func(Of Double, DoubleVector, DoubleVector)(AddressOf Rigid)
...and solve. The solution is returned as a key/value pair. The first Keyelement of the pair is
the time span vector, the second Valueelement of the pair is the corresponding solution values.
That is, if the computed solution function is y then
y(soln.Key(i)) = soln.Value(i)
Dim Soln As RungeKutta45OdeSolver.Solution(Of DoubleMatrix) = Solver.Solve(odeFunction, TimeSpan, y0, SolverOptions)
Print out a few values
Console.WriteLine("Number of time values = " & Soln.T.Length)
For I As Integer = 0 To Soln.T.Length - 1
If (I Mod 7 = 0) Then
Console.WriteLine("y({0}) = {1}", Soln.T(I).ToString("G5"), Soln.Y.Row(I).ToString("G5"))
End If
Next
We can also specify a mass matrix in the solver options...
Dim M As New DoubleMatrix("3x3[1 2 1 2 1 3 9 4 1]")
SolverOptions.MassMatrix = M
Soln = Solver.Solve(odeFunction, TimeSpan, y0, SolverOptions)
Print out a few values for the mass matrix solution.
Console.WriteLine(Environment.NewLine & "With constant mass matrix:")
Console.WriteLine("Number of time values = " & Soln.T.Length)
For I = 0 To Soln.T.Length - 1
If (I Mod 7 = 0) Then
Console.WriteLine("y({0}) = {1}", Soln.T(I).ToString("G5"), Soln.Y.Row(I).ToString("G5"))
End If
Next
You can also solve with default values for all the options
Soln = Solver.Solve(odeFunction, TimeSpan, y0)
End Sub
Sub NonAdaptiveStepSizeExample()
Construct the solver.
Dim Solver As New RungeKutta5OdeSolver()
Construct the time span vector. Imports the non-adaptive solve
function, integration steps will be performed only at the points
in the time span vector.
Dim TimeSpan As New DoubleVector(85, 0, 0.15)
Initial y vector.
Dim y0 As New DoubleVector(0.0, 1.0, 1.0)
Construct the delegate representing our system of differential equations...
Dim odeFunction As New Func(Of Double, DoubleVector, DoubleVector)(AddressOf Rigid)
...and solve. Each row of the returned matrix contains the calculated function values at the
corresponding point in the input time span vector.
Dim Soln As DoubleMatrix = Solver.Solve(odeFunction, TimeSpan, y0)
Print out a few values
Console.WriteLine("Number of time values = " & TimeSpan.Length)
For I As Integer = 0 To TimeSpan.Length - 1
If (I Mod 7 = 0) Then
Console.WriteLine("y({0}) = {1}", TimeSpan(I).ToString("G5"), Soln.Row(I).ToString("G5"))
End If
Next
Solve with a constant mass matrix
Dim M As New DoubleMatrix("3x3[1 2 1 2 1 3 9 4 1]")
Dim SolverOptions As New RungeKutta5OdeSolver.Options()
SolverOptions.MassMatrix = M
Soln = Solver.Solve(odeFunction, TimeSpan, y0, SolverOptions)
Print out a few values for the mass matrix solution.
Console.WriteLine(Environment.NewLine & "With constant mass matrix:")
Console.WriteLine("Number of time values = " & TimeSpan.Length)
For I As Integer = 0 To TimeSpan.Length - 1
If (I Mod 7 = 0) Then
Console.WriteLine("y({0}) = {1}", TimeSpan(I).ToString("G5"), Soln.Row(I).ToString("G5"))
End If
Next
Solve with an output function that is a class member function.
Dim outputFunctionObject As New SimpleNonAdaptiveOutput()
SolverOptions.OutputFunction = New Func(Of Double, DoubleVector, RungeKutta45OdeSolver.OutputFunctionFlag, Boolean)(AddressOf outputFunctionObject.OutputFunction)
Soln = Solver.Solve(odeFunction, TimeSpan, y0, SolverOptions)
End Sub
End Module
Class containing an example output function for the Dormand Prince ODE non-adaptive solve
function. The class does nothing but accumulate the solutions at each solver
integration stop.
Class SimpleNonAdaptiveOutput
Gets and sets a list of the time values.
Private _T As List(Of Double) = Nothing
Public Property T As List(Of Double)
Get
Return _T
End Get
Set(ByVal value As List(Of Double))
_T = value
End Set
End Property
Gets and sets a list of the calculated function values.
Private _Y As List(Of DoubleVector) = Nothing
Public Property Y As List(Of DoubleVector)
Get
Return _Y
End Get
Set(ByVal value As List(Of DoubleVector))
_Y = value
End Set
End Property
Gets and sets a flag indicating whether or not the output
function/class has been initialized.
Private _Init As Boolean = False
Public Property IsInitialized As Boolean
Get
Return _Init
End Get
Set(ByVal value As Boolean)
_Init = value
End Set
End Property
Flag indicating the output function is done being invoked.
Private _IsDone As Boolean = False
Public Property IsDone As Boolean
Get
Return _IsDone
End Get
Set(ByVal value As Boolean)
_IsDone = value
End Set
End Property
Constructs and initializes a SimpleNonAdaptiveOutput object.
Public Sub SimpleNonAdaptiveOutput()
IsInitialized = False
IsDone = False
T = Nothing
Y = Nothing
End Sub
<summary>
The output function to be called by the ODE solver at each integration step. If this function returns true
the solver continues with the next step, if it returns false the solver is halted.
</summary>
<param name="Time">The time value at the current integration step. If flag is equal to Initialize,
this will be the initial time value t0.</param>
<param name="YValue">The calculated function value at the current integration step. If flag is equal to Initialize,
this will be the initial function value y0.</param>
<param name="Flag">Flag indicating what stage the solver is at:
Initialize - before solving begins. y and t values are the problems initial values.
IntegrationStep - just finished an integration step. y and t values are the values
calculated at that step.
Done - just finished last step. y and t values are the last values in the returned
solution.</param>
<returns>true if the solver is to proceed with the calculation, false forces the solver
to stop.</returns>
Public Function OutputFunction(ByVal Time As Double, ByVal YValue As DoubleVector, ByVal Flag As RungeKutta45OdeSolver.OutputFunctionFlag) As Boolean
If (Flag = RungeKutta45OdeSolver.OutputFunctionFlag.Initialize) Then
T = New List(Of Double)()
T.Add(Time)
Y = New List(Of DoubleVector)()
Y.Add(YValue)
IsInitialized = True
IsDone = False
Return True
End If
If (Flag = RungeKutta45OdeSolver.OutputFunctionFlag.Done) Then
IsDone = True
T.Add(Time)
Y.Add(YValue)
Return True
End If
If (Flag = RungeKutta45OdeSolver.OutputFunctionFlag.IntegrationStep) Then
T.Add(Time)
Y.Add(YValue)
Return True
End If
Throw New InvalidArgumentException("Unknown output function flag: " + Flag)
End Function
End Class
End Namespace
← All NMath Code Examples
