VB.NET Runge Kutta Solver Example

Imports System
Imports System.Collections.Generic

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

Namespace RungeKuttaExample

  Module Program

    ' Example illustrating the use of the RungeKuttaSolver class for solving
    ' first order initial value differential equations.
  
    Sub Main()

      ' Consider the first order initial value problem
      ' y' = (3 - 4*y)/2*x, y(1) = -4
      ' which has solution
      ' y = 0.75 - 19.0/(4*x^2)
      ' To solve this differential equation using the 
      ' Runge Kutta method, we first construct a FirstOrderInitialValueProblem
      ' object which contains delegate for computing y' from x and y, and
      ' the initial values.
      Dim F As New NMathFunctions.DoubleBinaryFunction(AddressOf MyFunction)
    
      Dim X0 As Double = 1.0
      Dim y0 As Double = -4.0
      Dim Prob As New FirstOrderInitialValueProblem(F, X0, y0)

      ' The Runge Kutta method returns the solution function, f, as a set of tabulated values
      ' xi and f(xi). The initial x-value, x0 is the x0 of the initial condition. Number of
      ' specified points and delta determine the remaining xi as
      ' xi = x0 + i*delta
      Dim NumPoints As Integer = 3000
      Dim Delta As Double = 0.0005

      ' Construct a first order Runge Kutta solver (also known as Eulers method) and solve
      ' our first order initial values differential equation.
      Dim Solver As New RungeKuttaSolver(NumPoints, Delta, RungeKuttaSolver.SolverOrder.First)
      Dim TabulatedValues As KeyValuePair(Of Double(), Double()) = Solver.Solve(Prob)

      ' Construct a delegate for the actual solution function for comparison.
      Dim ActualSolutionFunction As New NMathFunctions.DoubleUnaryFunction(AddressOf myFunction2)

      ' The solution to the differential equation is returned as a key-value pair. The
      ' key is an array of doubles containing the x-values and the value is an array
      ' containing the y = f(x) values. Below we compare the computed function values
      ' to the know function values by computing the maximum absolute value of their
      ' differences (l-infinity norm of the difference).
      Dim XValues As New DoubleVector(TabulatedValues.Key)
      Dim calculatedSolutionYvalues As New DoubleVector(TabulatedValues.Value)
      Dim actualSolutionYvalues = XValues.Apply(ActualSolutionFunction)
      Dim MaxDiff As Double = (calculatedSolutionYvalues - actualSolutionYvalues).InfinityNorm()
      Console.WriteLine("Maximum error using first order Runge Kutta = " & MaxDiff)

      ' Solve the same problem with a second order Runge Kutta method and compute
      ' the error. Note the error decreases.
      solver.Order = RungeKuttaSolver.SolverOrder.Second
      tabulatedValues = solver.Solve(prob)
      calculatedSolutionYvalues = New DoubleVector(tabulatedValues.Value)
      maxDiff = (calculatedSolutionYvalues - actualSolutionYvalues).InfinityNorm()
      Console.WriteLine("Maximum error using second order Runge Kutta = " & MaxDiff)

      ' You can also specify a TabulatedFunction object into which to place the
      ' solution. We specify a NaturalCubicSpline function and use a fourth order
      ' Runge Kutta to solve the same problem. Note the solution improves.
      Dim SolutionFunction As TabulatedFunction = New NaturalCubicSpline()
      solver.Order = RungeKuttaSolver.SolverOrder.Fourth
      Solver.Solve(Prob, SolutionFunction)
      MaxDiff = 0

      Dim I As Integer
      For I = 0 To SolutionFunction.NumberOfTabulatedValues - 1
        Dim absDiff As Double = Math.Abs(SolutionFunction.GetY(I) - ActualSolutionFunction(SolutionFunction.GetX(I)))
        MaxDiff = Math.Max(MaxDiff, absDiff)
      Next

      Console.WriteLine("Maximum error using fourth order Runge Kutta = " & MaxDiff)

      Console.WriteLine()
      Console.WriteLine("Press Enter Key")
      Console.Read()

    End Sub

    Private Function MyFunction(ByVal x As Double, ByVal y As Double) As Double
      Return (3.0 - 4.0 * y) / (2.0 * x)
    End Function

    Private Function MyFunction2(ByVal x As Double) As Double
      Return 0.75 - (19.0 / (4.0 * x * x))
    End Function

  End Module
End Namespace