ODE integration example

This tutorial shows how to perform the numerical integration of an ordinary differential equation (ODE) using DA objects as the state components' type. It also demonstrates how to extract the state transition matrix (STM) from the polynomial expansion of the dynamics flow.

In this case, the state vector represents an object in orbit around a central body and subject only to the gravitational pull of the latter. Its motion is described by the differential equations for the Kepler problem expressed in Cartesian coordinates.

Install dependencies

Make sure the required packages are installed

using Pkg
Pkg.add("DACE")

Using DACE

Load the required modules

using OrdinaryDiffEq
using DACE

Define the parameters and intial conditions for a circular orbit with a normalized radius equal to one

μ = 1.0
x0 = [1.0, 0.0, 0.0, 0.0, 1.0, 0.0]

6-element Vector{Float64}:
 1.0
 0.0
 0.0
 0.0
 1.0
 0.0

Set the integration time span equal to one revolution

t0 = 0.0
tf = 2π

6.283185307179586

Define the equations of motion for the resticted two-body problem

function kepler_ode!(du, u, μ, _)
    du[1:3] .= u[4:6]
    du[4:6] .= -u[1:3] * μ ./ sum(u[1:3].^2)^(3/2)
end

kepler_ode! (generic function with 1 method)

Compute the nominal solution

prob = ODEProblem(kepler_ode!, x0, [t0, tf], μ)
sol = solve(prob, Vern9(), abstol=1e-12, reltol=1e-12)

xf = sol.u[end]

6-element Vector{Float64}:
  0.9999999999999963
 -4.963822535880929e-14
  0.0
  4.3859021831616145e-14
  0.9999999999999949
  0.0

Initialize DACE to compute second-order expansions

DACE.init(2,6)

Define the initial conditions as a DA vector

dx_dace = [DA(i,1.0) for i in 1:6]
x0_dace = x0 .+ dx_dace

println(x0_dace)

DACE.DAAllocated[     I  COEFFICIENT              ORDER EXPONENTS
     1    1.0000000000000000e+00   0   0  0  0  0  0  0
     2    1.0000000000000000e+00   1   1  0  0  0  0  0
------------------------------------------------
,      I  COEFFICIENT              ORDER EXPONENTS
     1    1.0000000000000000e+00   1   0  1  0  0  0  0
------------------------------------------------
,      I  COEFFICIENT              ORDER EXPONENTS
     1    1.0000000000000000e+00   1   0  0  1  0  0  0
------------------------------------------------
,      I  COEFFICIENT              ORDER EXPONENTS
     1    1.0000000000000000e+00   1   0  0  0  1  0  0
------------------------------------------------
,      I  COEFFICIENT              ORDER EXPONENTS
     1    1.0000000000000000e+00   0   0  0  0  0  0  0
     2    1.0000000000000000e+00   1   0  0  0  0  1  0
------------------------------------------------
,      I  COEFFICIENT              ORDER EXPONENTS
     1    1.0000000000000000e+00   1   0  0  0  0  0  1
------------------------------------------------
]

Compute the polynomial expansion of the dynamics flow

sol_dace = solve(remake(prob, u0=x0_dace), Vern9(), abstol=1e-12, reltol=1e-12)

xf_dace = sol_dace.u[end]
xf_cons = DACE.cons.(xf_dace)

6-element Vector{Float64}:
  0.9999999999999967
 -6.053491041768666e-14
  0.0
  5.355438315035599e-14
  0.9999999999999952
  0.0

Verify the (nominal) final state

println("Max abs error on final state: " * string(maximum(abs.(xf_cons - xf))))

Max abs error on final state: 1.0896685058877372e-14

Extract the state transition matrix (STM)

stm_dace = DACE.jacobian(xf_dace)
stm_cons = DACE.cons.(stm_dace)

6×6 Matrix{Float64}:
   1.0          -1.83384e-13  0.0          -2.40419e-13   -7.38964e-13   0.0
 -18.8496        1.0          0.0           7.24656e-13  -18.8496        0.0
   0.0           0.0          1.0           0.0            0.0          -6.05349e-14
  18.8496       -5.46722e-13  0.0           1.0           18.8496        0.0
  -1.10489e-12   1.54696e-13  0.0           9.8001e-14     1.0           0.0
   0.0           0.0          5.35544e-14   0.0            0.0           1.0

This page was generated using Literate.jl.