diff --git a/src/integration.jl b/src/integration.jl index 8ae65cb3..7365894a 100644 --- a/src/integration.jl +++ b/src/integration.jl @@ -24,6 +24,31 @@ function discrete_dynamics(::Type{RK4}, model::AbstractModel, x::SVector, u::SVe x + (k1 + 4k2 + 4k3 + k4)/6 end +function discrete_dynamics(::Type{RK8}, model::AbstractModel, y::SVector{N,T}, u::SVector{M,T},t, dt::T) where {N,M,T} + α = @SVector [ 2/27, 1/9, 1/6, 5/12, .5, 5/6, 1/6, 2/3, 1/3, 1, 0, 1 ]; + β = @SMatrix [ 2/27 0 0 0 0 0 0 0 0 0 0 0 0; + 1/36 1/12 0 0 0 0 0 0 0 0 0 0 0; + 1/24 0 1/8 0 0 0 0 0 0 0 0 0 0; + 5/12 0 -25/16 25/16 0 0 0 0 0 0 0 0 0; + .05 0 0 .25 .2 0 0 0 0 0 0 0 0; + -25/108 0 0 125/108 -65/27 125/54 0 0 0 0 0 0 0; + 31/300 0 0 0 61/225 -2/9 13/900 0 0 0 0 0 0; + 2 0 0 -53/6 704/45 -107/9 67/90 3 0 0 0 0 0; + -91/108 0 0 23/108 -976/135 311/54 -19/60 17/6 -1/12 0 0 0 0; + 2383/4100 0 0 -341/164 4496/1025 -301/82 2133/4100 45/82 45/164 18/41 0 0 0; + 3/205 0 0 0 0 -6/41 -3/205 -3/41 3/41 6/41 0 0 0; + -1777/4100 0 0 -341/164 4496/1025 -289/82 2193/4100 51/82 33/164 12/41 0 1 0] + χ = @SVector[0, 0, 0, 0, 0, 34/105, 9/35, 9/35, 9/280, 9/280, 0, 41/840, 41/840 ]; # Difference between two bottom layers of butcher tableau + + f = y*zeros(eltype(y),1,13) + + f[:,1] = dynamics(model,y,u,t) + for j = 1:12 + f[:,j+1] = dynamics(model, y + dt*f*β[j,:],u,t + α[j]*dt) + end + return y + dt*f*χ +end + ############################################################################################ # EXPLICIT METHODS # diff --git a/src/model.jl b/src/model.jl index a6ddcc0b..95094270 100644 --- a/src/model.jl +++ b/src/model.jl @@ -62,13 +62,16 @@ abstract type Implicit <: QuadratureRule end "Integration rules of the form x′ = f(x,u,x′,u′), where x′,u′ are the states and controls at the next time step." abstract type Explicit <: QuadratureRule end -"Fourth-order Runge-Kutta method with zero-order-old on the controls" +"Eigth-order Runge-Kutta method" +abstract type RK8 <: Implicit end + +"Fourth-order Runge-Kutta method with zero-order-hold on the controls" abstract type RK4 <: Implicit end -"Second-order Runge-Kutta method with zero-order-old on the controls" +"Second-order Runge-Kutta method with zero-order-hold on the controls" abstract type RK3 <: Implicit end -"Second-order Runge-Kutta method with zero-order-old on the controls" +"Second-order Runge-Kutta method with zero-order-hold on the controls" abstract type RK2 <: Implicit end "Third-order Runge-Kutta method with first-order-hold on the controls"