class Spitzy::Ode

Numerically solves an initial value problem of the form

Attributes

fx[R]

value of f at every point in x

method[R]

the numerical scheme applied

mx[R]

number of points (i.e. length of x)

tol[R]

the error tolerance of the method (only applicable for methods with automatic step size adjustment)

u[R]

the numerical solution as an array

x[R]

array of the points in the domain at which the numerical solution was evaluated

Public Class Methods

new(xrange:, dx:, yini:, tol: 1e-2, maxiter: 1e6, method: :dopri, &f) click to toggle source

Constructor for all solver routines for the initial value problem:

  • dy/dx = f(x,y), xmin < x < xmax,

  • y(xmin) = yini

It initializes the parameters and solves the equation using one of the methods: Dormand-Prince, Forward Euler, Adams-Bashforth order 2.

Arguments

  • xrange - An array of the form [xmin, xmax]. The interval on which the solution will be evaluated.

  • dx - Determines the grid of x values. The length of the x step for methods without automatic step size adjustment, or the maximum step size for methods with automatic step size adjustment.

  • yini - Initial value for y, i.e. y(xmin).

  • tol - The desired error tolerance of the method (only for methods with automatic step size correction).

  • maxiter - The maximal number of performed iterations for methods with automatic step size adjustment.

  • method - The numerical scheme used to solve the ODE. Possible values are: :dopri (Dormand-Prince), :euler (Forward Euler), :ab2 (Adams-Bashforth of order 2). :dopri performs an automatic step size adjustment and error estimation in order to keep the error of the numerical solution below tol. :euler and :ab2 operate on a fixed step size dx, and cannot estimate the error of the numerical solution. Currently, the default is :dopri (Dormand-Prince is currently also the default method in MATLAB and GNU Octave's ode45 solver and is the default choice for the Simulink's model explorer solver).

  • +&f+ - The right hand side of the differential equation which must be supplied as a proc object. It is a function of x and y, where x should be the first argument, and y the second.

Usage

# Dormand-Prince example 1
f = proc { |t,y| 1.0 + y/t + (y/t)**2 }
dopri_sol = Spitzy::Ode.new(xrange: [1.0,4.0], dx: 0.1, yini: 0.0, &f) 
puts "The numerical solution at x=4 is u(4)=#{dopri_sol.u[-1]}"

# Dormand-Prince example 2
f = proc { |t,y| -2.0 * y + Math::exp(-2.0 * (t - 6.0)**2) }
dopri_sol = Spitzy::Ode.new(xrange: [0.0,10.0], dx: 1.5, yini: 1.0, 
                            tol: 1e-6, maxiter: 1e6, &f) 
# Plot dopri_sol.x agains dopri_sol.u to observe step size adaptivity

# Euler example 1
f = proc { |t,y| 1.0 + y/t + (y/t)**2 }
euler_sol = Spitzy::Ode.new(xrange: [1.0,4.0], dx: 0.01,
                            yini: 0.0, method: :euler, &f) 

# Adams-Bashforth example 1
f = proc { |t,y| -2.0*t*y }
ab2_sol = Spitzy::Ode.new(xrange: [0.0,4.0], dx: 0.1, 
                          yini: 1.0, method: :ab2, &f) 
exact_sol = ab2_sol.x.map { |tt| Math::exp(-(tt**2)) }
maxerror1 = exact_sol1.each_with_index.map {|n,i| n - ab2_sol1.u[i] }.max.abs
puts "Number of x steps: #{ab2_sol1.mx}"
puts "Error: #{maxerror1}maxerror2}"
puts "Error ratio: #{maxerror1/maxerror2}"
# File lib/spitzy/ode.rb, line 84
def initialize(xrange:, dx:, yini:, tol: 1e-2, maxiter: 1e6, method: :dopri, &f)

  raise(ArgumentError, "Expected xrange to be an array of length 2") unless xrange.length == 2
  @dx = dx
  @xmin = xrange[0]
  @xmax = xrange[1]
  @yini = yini
  @maxiter = maxiter
  @tol = tol
  @f = f 
  @x = [] # Stores the x grid
  @fx = [] # Stores the values of f at every point in x
  @u = [] # Stores numerical solution

  @method = method
  case @method
    when :euler then euler
    when :ab2 then ab2
    when :am3 then am3
    when :dopri then dopri
    else raise(ArgumentError, "#{@method} is not a valid method.")
  end
end

Public Instance Methods

f(x0, y0) click to toggle source

Evaluate the function f(x,y) at x=x0 and y=y0

Arguments

  • x0 - A floating point number

  • y0 - A floating point number representing y(x0)

# File lib/spitzy/ode.rb, line 116
def f(x0, y0)
  @f.call(x0, y0)
end

Private Instance Methods

ab2() click to toggle source

Solve the initial value problem

* dy/dx = f(x,y), xmin < x < xmax, 
* y(xmin) = yini

with the Adams-Bashforth method of order 2, given by the explicit formula u = u + dx/2 * (3*f - f). Since the method requires the last two functional values for the approximation of the next functional value, a Runge-Kutta method of order 2 (Heun’s method) is applied to approximate the functional value at the second time step.

# File lib/spitzy/ode.rb, line 152
def ab2 
  @x = (@xmin..@xmax).step(@dx).to_a # x steps
  @x << @xmax if @x.last < @xmax
  @mx = @x.length # Number of x steps

  @u[0] = @yini
  @fx[0] = self.f(@x[0], @u[0])
  # Runge-Kutta of order 2 for the second time step
  @u[1] = @u[0] + @dx/2.0 * (@fx[0] + self.f(@x[1], @u[0] + @dx*@fx[0]))
  @fx[1] = self.f(@x[1], @u[1])

  1.upto(@mx-2) do |n|
    @u[n+1] = @u[n] + @dx/2.0 * (3.0*@fx[n] - @fx[n-1])
    @fx[n+1] = self.f(@x[n+1], @u[n+1])
  end

  #tol and maxiter not relevant for the this method
  @tol = nil
  @maxiter = nil
end
dopri() click to toggle source

Solve the initial value problem

  • dy/dx = f(x,y), xmin < x < xmax,

  • y(xmin) = yini

with the Dormand-Prince method.

This method automatically adapts the step size in order to keep the error of the numerical solution below the tolerance level tol. However, the step size is not allowed to exceed the specified maximal step size dx (dx is also the initial step size). The algorithm throws an exception if it fails to compute the numerical solution on the given x domain in less than or equal to maxiter iterations.

Reference

  • J R Dormand, P J Prince, “A family of embedded Runge-Kutta formulae”

  • E Hairer, S P Norsett and G Wanner, “Solving Differential Equations I, Nonstiff Problems”, p. 176.

# File lib/spitzy/ode.rb, line 190
def dopri
  a21 = 1.0/5.0
  a31 = 3.0/40.0
  a32 = 9.0/40.0
  a41 = 44.0/45.0
  a42 = -56.0/15.0
  a43 = 32.0/9.0
  a51 = 19372.0/6561.0
  a52 = -25360.0/2187.0
  a53 = 64448.0/6561.0
  a54 = -212.0/729.0
  a61 = 9017.0/3168.0
  a62 = -355.0/33.0
  a63 = 46732.0/5247.0
  a64 = 49.0/176.0
  a65 = -5103.0/18656.0
  a71 = 35.0/384.0
  a72 = 0.0
  a73 = 500.0/1113.0
  a74 = 125.0/192.0
  a75 = -2187.0/6784.0
  a76 = 11.0/84.0

  c2 = 1.0 / 5.0
  c3 = 3.0 / 10.0
  c4 = 4.0 / 5.0
  c5 = 8.0 / 9.0
  c6 = 1.0
  c7 = 1.0

  b1order5 = 35.0/384.0
  b2order5 = 0.0
  b3order5 = 500.0/1113.0
  b4order5 = 125.0/192.0
  b5order5 = -2187.0/6784.0
  b6order5 = 11.0/84.0
  b7order5 = 0.0

  b1order4 = 5179.0/57600.0
  b2order4 = 0.0
  b3order4 = 7571.0/16695.0
  b4order4 = 393.0/640.0
  b5order4 = -92097.0/339200.0
  b6order4 = 187.0/2100.0
  b7order4 = 1.0/40.0

  @x[0] = @xmin
  @u[0] = @yini
  @fx[0] = self.f(@x[0], @u[0])
  h = @dx 
  i = 0

  0.upto(@maxiter) do |iter|
     # Compute the function values
     k1 = @fx[i] 
     k2 = self.f(@x[i] + c2*h, @u[i] + h*(a21*k1))
     k3 = self.f(@x[i] + c3*h, @u[i] + h*(a31*k1+a32*k2))
     k4 = self.f(@x[i] + c4*h, @u[i] + h*(a41*k1+a42*k2+a43*k3))
     k5 = self.f(@x[i] + c5*h, @u[i] + h*(a51*k1+a52*k2+a53*k3+a54*k4))
     k6 = self.f(@x[i] + h, @u[i] + h*(a61*k1+a62*k2+a63*k3+a64*k4+a65*k5))
     k7 = self.f(@x[i] + h, @u[i] + h*(a71*k1+a72*k2+a73*k3+a74*k4+a75*k5+a76*k6))

     error = (b1order5 - b1order4)*k1 + (b3order5 - b3order4)*k3 + (b4order5 - b4order4)*k4 + 
       (b5order5 - b5order4)*k5 + (b6order5 - b6order4)*k6 + (b7order5 - b7order4)*k7
     error = error.abs

     # error control
     if error < @tol then
        @x[i+1] = @x[i] + h
        @u[i+1] = @u[i] + h * (b1order5*k1 + b3order5*k3 + b4order5*k4 + b5order5*k5 + b6order5*k6)
        @fx[i+1] = self.f(@x[i+1], @u[i+1])
        i = i+1
     end

     delta = 0.84 * (@tol / error)**0.2
     if delta <= 0.1 then
        h = h * 0.1
     elsif delta >= 4.0 then
        h = h * 4.0
     else 
        h = delta * h
     end

     # set h to the user specified maximal allowed value
     h = @dx if h > @dx 

     if @x[i] >= @xmax then
       break
     elsif @x[i] + h > @xmax then
       h = @xmax - @x[i]
     end
  end

  @mx = @x.length # Number of x steps

  raise(RuntimeError, "Maximal number of iterations reached 
        before evaluation of the solution on the entire x interval 
        was completed (try to increase maxiter or use a different method") if @x.last < @xmax
end
euler() click to toggle source

Solve the initial value problem

* dy/dx = f(x,y), xmin < x < xmax, 
* y(xmin) = yini

with the forward Euler scheme u = u + dx * f.

# File lib/spitzy/ode.rb, line 128
def euler 
  @x = (@xmin..@xmax).step(@dx).to_a # x steps
  @x << @xmax if @x.last < @xmax
  @mx = @x.length # Number of x steps
  @u[0] = @yini
  @fx[0] = self.f(@x[0], @u[0])
  0.upto(@mx-2) do |n|
    @u[n+1] = @u[n] + @dx * @fx[n]
    @fx[n+1] = self.f(@x[n+1], @u[n+1])
  end
  #tol and maxiter not relevant for the Euler method
  @tol = nil
  @maxiter = nil
end