VB Runge Kutta Solver Example

← All NMath Code Examples

 

Imports System
Imports System.Collections.Generic
Imports System.Diagnostics

Imports CenterSpace.NMath.Core
Imports CenterSpace.NMath.Analysis


Namespace CenterSpace.NMath.Analysis.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 'Key' element of the pair is
      ' the time span vector, the second 'Value' element 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 'Key' element of the pair is
      ' the time span vector, the second 'Value' element 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
Top