Skip to content

Instantly share code, notes, and snippets.

@ChrisRackauckas

ChrisRackauckas/differential_equations_taylor_integration_benchmark.md

Last active August 31, 2020 05:17
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save ChrisRackauckas/1301b23aa12ad83de7138d8e41d64dd6 to your computer and use it in GitHub Desktop.
Save ChrisRackauckas/1301b23aa12ad83de7138d8e41d64dd6 to your computer and use it in GitHub Desktop.
Download ZIP
DifferentialEquations.jl Runge-Kutta integrators vs Taylor Integration Benchmarks
Raw
differential_equations_taylor_integration_benchmark.md

DifferentialEquations.jl Runge-Kutta Integrators vs Taylor Integration

This is against a Taylor integration method that utilizes taylor-mode automatic differentiation in a specifically optimized way for ODE integration via its @taylorize macro (for more details, see the documentation). This demonstration is done on the classic Pleiades n-body problem, a non-stiff ODE since the Taylor integration methods are only for non-stiff equations.

Raw
runge_kutta_vs_taylor.jl
using OrdinaryDiffEq, DiffEqDevTools, TaylorIntegration
function f(du,u,p,t)
@inbounds begin
x = view(u,1:7) # x
y = view(u,8:14) # y
v = view(u,15:21) # x′
w = view(u,22:28) # y′
du[1:7] .= v
du[8:14].= w
for i in 15:28
du[i] = zero(u[1])
end
for i=1:7,j=1:7
if i != j
r = ((x[i]-x[j])^2 + (y[i] - y[j])^2)^(3/2)
du[14+i] += j*(x[j] - x[i])/r
du[21+i] += j*(y[j] - y[i])/r
end
end
end
end
prob = ODEProblem(f,[3.0,3.0,-1.0,-3.0,2.0,-2.0,2.0,3.0,-3.0,2.0,0,0,-4.0,4.0,0,0,0,0,0,1.75,-1.5,0,0,0,-1.25,1,0,0],(0.0,3.0))
abstols = 1.0 ./ 10.0 .^ (6:9)
reltols = 1.0 ./ 10.0 .^ (3:6);
sol = solve(prob,Vern8(),abstol=1/10^12,reltol=1/10^10,maxiters=1000000)
test_sol = TestSolution(sol);
using Plots; gr()
setups = [Dict(:alg=>Tsit5())
Dict(:alg=>Vern6())
Dict(:alg=>Vern7())
Dict(:alg=>Vern9())
Dict(:alg=>TaylorMethod(8))
]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;appxsol=test_sol,numruns=100,maxiters=10000,verbose=false)
plot(wp)
# Possibly not the fastest way, so expand the ODE to allow for its `@taylorize` to work and get full speed
using ModelingToolkit
sys = modelingtoolkitize(prob)
clipboard(ODEFunctionExpr(sys))
@inbounds @taylorize function taylor_f(var"##MTIIPVar#313", var"##MTKArg#309", var"##MTKArg#310", var"##MTKArg#311")
(x₁, x₂, x₃, x₄, x₅, x₆, x₇, x₈, x₉, x₁₀, x₁₁, x₁₂, x₁₃, x₁₄, x₁₅, x₁₆, x₁₇, x₁₈, x₁₉, x₂₀, x₂₁, x₂₂, x₂₃, x₂₄, x₂₅, x₂₆, x₂₇, x₂₈, t) = (var"##MTKArg#309"[1], var"##MTKArg#309"[2], var"##MTKArg#309"[3], var"##MTKArg#309"[4], var"##MTKArg#309"[5], var"##MTKArg#309"[6], var"##MTKArg#309"[7], var"##MTKArg#309"[8], var"##MTKArg#309"[9], var"##MTKArg#309"[10], var"##MTKArg#309"[11], var"##MTKArg#309"[12], var"##MTKArg#309"[13], var"##MTKArg#309"[14], var"##MTKArg#309"[15], var"##MTKArg#309"[16], var"##MTKArg#309"[17], var"##MTKArg#309"[18], var"##MTKArg#309"[19], var"##MTKArg#309"[20], var"##MTKArg#309"[21], var"##MTKArg#309"[22], var"##MTKArg#309"[23], var"##MTKArg#309"[24], var"##MTKArg#309"[25], var"##MTKArg#309"[26], var"##MTKArg#309"[27], var"##MTKArg#309"[28], var"##MTKArg#311")
var"##MTIIPVar#313"[1] = x₁₅
var"##MTIIPVar#313"[2] = x₁₆
var"##MTIIPVar#313"[3] = x₁₇
var"##MTIIPVar#313"[4] = x₁₈
var"##MTIIPVar#313"[5] = x₁₉
var"##MTIIPVar#313"[6] = x₂₀
var"##MTIIPVar#313"[7] = x₂₁
var"##MTIIPVar#313"[8] = x₂₂
var"##MTIIPVar#313"[9] = x₂₃
var"##MTIIPVar#313"[10] = x₂₄
var"##MTIIPVar#313"[11] = x₂₅
var"##MTIIPVar#313"[12] = x₂₆
var"##MTIIPVar#313"[13] = x₂₇
var"##MTIIPVar#313"[14] = x₂₈
var"##MTIIPVar#313"[15] = (getproperty(Base, :+))((getproperty(Base, :+))((getproperty(Base, :+))((getproperty(Base, :+))((getproperty(Base, :+))((getproperty(Base, :+))(0, (getproperty(Base, :/))((getproperty(Base, :*))(2, (getproperty(Base, :-))(x₂, x₁)), (getproperty(Base, :^))((getproperty(Base, :+))((getproperty(Base, :^))((getproperty(Base, :-))(x₁, x₂), 2), (getproperty(Base, :^))((getproperty(Base, :-))(x₈, x₉), 2)), 1.5))), (getproperty(Base, :/))((getproperty(Base, :*))(3, (getproperty(Base, :-))(x₃, x₁)), (getproperty(Base, :^))((getproperty(Base, :+))((getproperty(Base, :^))((getproperty(Base, :-))(x₁, x₃), 2), (getproperty(Base, :^))((getproperty(Base, :-))(x₈, x₁₀), 2)), 1.5))), (getproperty(Base, :/))((getproperty(Base, :*))(4, (getproperty(Base, :-))(x₄, x₁)), (getproperty(Base, :^))((getproperty(Base, :+))((getproperty(Base, :^))((getproperty(Base, :-))(x₁, x₄), 2), (getproperty(Base, :^))((getproperty(Base, :-))(x₈, x₁₁), 2)), 1.5))), (getproperty(Base, :/))((getproperty(Base, :*))(5, (getproperty(Base, :-))(x₅, x₁)), (getproperty(Base, :^))((getproperty(Base, :+))((getproperty(Base, :^))((getproperty(Base, :-))(x₁, x₅), 2), (getproperty(Base, :^))((getproperty(Base, :-))(x₈, x₁₂), 2)), 1.5))), (getproperty(Base, :/))((getproperty(Base, :*))(6, (getproperty(Base, :-))(x₆, x₁)), (getproperty(Base, :^))((getproperty(Base, :+))((getproperty(Base, :^))((getproperty(Base, :-))(x₁, x₆), 2), (getproperty(Base, :^))((getproperty(Base, :-))(x₈, x₁₃), 2)), 1.5))), (getproperty(Base, :/))((getproperty(Base, :*))(7, (getproperty(Base, :-))(x₇, x₁)), (getproperty(Base, :^))((getproperty(Base, :+))((getproperty(Base, :^))((getproperty(Base, :-))(x₁, x₇), 2), (getproperty(Base, :^))((getproperty(Base, :-))(x₈, x₁₄), 2)), 1.5)))
var"##MTIIPVar#313"[16] = (getproperty(Base, :+))((getproperty(Base, :+))((getproperty(Base, :+))((getproperty(Base, :+))((getproperty(Base, :+))((getproperty(Base, :+))(0, (getproperty(Base, :/))((getproperty(Base, :*))(1, (getproperty(Base, :-))(x₁, x₂)), (getproperty(Base, :^))((getproperty(Base, :+))((getproperty(Base, :^))((getproperty(Base, :-))(x₂, x₁), 2), (getproperty(Base, :^))((getproperty(Base, :-))(x₉, x₈), 2)), 1.5))), (getproperty(Base, :/))((getproperty(Base, :*))(3, (getproperty(Base, :-))(x₃, x₂)), (getproperty(Base, :^))((getproperty(Base, :+))((getproperty(Base, :^))((getproperty(Base, :-))(x₂, x₃), 2), (getproperty(Base, :^))((getproperty(Base, :-))(x₉, x₁₀), 2)), 1.5))), (getproperty(Base, :/))((getproperty(Base, :*))(4, (getproperty(Base, :-))(x₄, x₂)), (getproperty(Base, :^))((getproperty(Base, :+))((getproperty(Base, :^))((getproperty(Base, :-))(x₂, x₄), 2), (getproperty(Base, :^))((getproperty(Base, :-))(x₉, x₁₁), 2)), 1.5))), (getproperty(Base, :/))((getproperty(Base, :*))(5, (getproperty(Base, :-))(x₅, x₂)), (getproperty(Base, :^))((getproperty(Base, :+))((getproperty(Base, :^))((getproperty(Base, :-))(x₂, x₅), 2), (getproperty(Base, :^))((getproperty(Base, :-))(x₉, x₁₂), 2)), 1.5))), (getproperty(Base, :/))((getproperty(Base, :*))(6, (getproperty(Base, :-))(x₆, x₂)), (getproperty(Base, :^))((getproperty(Base, :+))((getproperty(Base, :^))((getproperty(Base, :-))(x₂, x₆), 2), (getproperty(Base, :^))((getproperty(Base, :-))(x₉, x₁₃), 2)), 1.5))), (getproperty(Base, :/))((getproperty(Base, :*))(7, (getproperty(Base, :-))(x₇, x₂)), (getproperty(Base, :^))((getproperty(Base, :+))((getproperty(Base, :^))((getproperty(Base, :-))(x₂, x₇), 2), (getproperty(Base, :^))((getproperty(Base, :-))(x₉, x₁₄), 2)), 1.5)))
var"##MTIIPVar#313"[17] = (getproperty(Base, :+))((getproperty(Base, :+))((getproperty(Base, :+))((getproperty(Base, :+))((getproperty(Base, :+))((getproperty(Base, :+))(0, (getproperty(Base, :/))((getproperty(Base, :*))(1, (getproperty(Base, :-))(x₁, x₃)), (getproperty(Base, :^))((getproperty(Base, :+))((getproperty(Base, :^))((getproperty(Base, :-))(x₃, x₁), 2), (getproperty(Base, :^))((getproperty(Base, :-))(x₁₀, x₈), 2)), 1.5))), (getproperty(Base, :/))((getproperty(Base, :*))(2, (getproperty(Base, :-))(x₂, x₃)), (getproperty(Base, :^))((getproperty(Base, :+))((getproperty(Base, :^))((getproperty(Base, :-))(x₃, x₂), 2), (getproperty(Base, :^))((getproperty(Base, :-))(x₁₀, x₉), 2)), 1.5))), (getproperty(Base, :/))((getproperty(Base, :*))(4, (getproperty(Base, :-))(x₄, x₃)), (getproperty(Base, :^))((getproperty(Base, :+))((getproperty(Base, :^))((getproperty(Base, :-))(x₃, x₄), 2), (getproperty(Base, :^))((getproperty(Base, :-))(x₁₀, x₁₁), 2)), 1.5))), (getproperty(Base, :/))((getproperty(Base, :*))(5, (getproperty(Base, :-))(x₅, x₃)), (getproperty(Base, :^))((getproperty(Base, :+))((getproperty(Base, :^))((getproperty(Base, :-))(x₃, x₅), 2), (getproperty(Base, :^))((getproperty(Base, :-))(x₁₀, x₁₂), 2)), 1.5))), (getproperty(Base, :/))((getproperty(Base, :*))(6, (getproperty(Base, :-))(x₆, x₃)), (getproperty(Base, :^))((getproperty(Base, :+))((getproperty(Base, :^))((getproperty(Base, :-))(x₃, x₆), 2), (getproperty(Base, :^))((getproperty(Base, :-))(x₁₀, x₁₃), 2)), 1.5))), (getproperty(Base, :/))((getproperty(Base, :*))(7, (getproperty(Base, :-))(x₇, x₃)), (getproperty(Base, :^))((getproperty(Base, :+))((getproperty(Base, :^))((getproperty(Base, :-))(x₃, x₇), 2), (getproperty(Base, :^))((getproperty(Base, :-))(x₁₀, x₁₄), 2)), 1.5)))
var"##MTIIPVar#313"[18] = (getproperty(Base, :+))((getproperty(Base, :+))((getproperty(Base, :+))((getproperty(Base, :+))((getproperty(Base, :+))((getproperty(Base, :+))(0, (getproperty(Base, :/))((getproperty(Base, :*))(1, (getproperty(Base, :-))(x₁, x₄)), (getproperty(Base, :^))((getproperty(Base, :+))((getproperty(Base, :^))((getproperty(Base, :-))(x₄, x₁), 2), (getproperty(Base, :^))((getproperty(Base, :-))(x₁₁, x₈), 2)), 1.5))), (getproperty(Base, :/))((getproperty(Base, :*))(2, (getproperty(Base, :-))(x₂, x₄)), (getproperty(Base, :^))((getproperty(Base, :+))((getproperty(Base, :^))((getproperty(Base, :-))(x₄, x₂), 2), (getproperty(Base, :^))((getproperty(Base, :-))(x₁₁, x₉), 2)), 1.5))), (getproperty(Base, :/))((getproperty(Base, :*))(3, (getproperty(Base, :-))(x₃, x₄)), (getproperty(Base, :^))((getproperty(Base, :+))((getproperty(Base, :^))((getproperty(Base, :-))(x₄, x₃), 2), (getproperty(Base, :^))((getproperty(Base, :-))(x₁₁, x₁₀), 2)), 1.5))), (getproperty(Base, :/))((getproperty(Base, :*))(5, (getproperty(Base, :-))(x₅, x₄)), (getproperty(Base, :^))((getproperty(Base, :+))((getproperty(Base, :^))((getproperty(Base, :-))(x₄, x₅), 2), (getproperty(Base, :^))((getproperty(Base, :-))(x₁₁, x₁₂), 2)), 1.5))), (getproperty(Base, :/))((getproperty(Base, :*))(6, (getproperty(Base, :-))(x₆, x₄)), (getproperty(Base, :^))((getproperty(Base, :+))((getproperty(Base, :^))((getproperty(Base, :-))(x₄, x₆), 2), (getproperty(Base, :^))((getproperty(Base, :-))(x₁₁, x₁₃), 2)), 1.5))), (getproperty(Base, :/))((getproperty(Base, :*))(7, (getproperty(Base, :-))(x₇, x₄)), (getproperty(Base, :^))((getproperty(Base, :+))((getproperty(Base, :^))((getproperty(Base, :-))(x₄, x₇), 2), (getproperty(Base, :^))((getproperty(Base, :-))(x₁₁, x₁₄), 2)), 1.5)))
var"##MTIIPVar#313"[19] = (getproperty(Base, :+))((getproperty(Base, :+))((getproperty(Base, :+))((getproperty(Base, :+))((getproperty(Base, :+))((getproperty(Base, :+))(0, (getproperty(Base, :/))((getproperty(Base, :*))(1, (getproperty(Base, :-))(x₁, x₅)), (getproperty(Base, :^))((getproperty(Base, :+))((getproperty(Base, :^))((getproperty(Base, :-))(x₅, x₁), 2), (getproperty(Base, :^))((getproperty(Base, :-))(x₁₂, x₈), 2)), 1.5))), (getproperty(Base, :/))((getproperty(Base, :*))(2, (getproperty(Base, :-))(x₂, x₅)), (getproperty(Base, :^))((getproperty(Base, :+))((getproperty(Base, :^))((getproperty(Base, :-))(x₅, x₂), 2), (getproperty(Base, :^))((getproperty(Base, :-))(x₁₂, x₉), 2)), 1.5))), (getproperty(Base, :/))((getproperty(Base, :*))(3, (getproperty(Base, :-))(x₃, x₅)), (getproperty(Base, :^))((getproperty(Base, :+))((getproperty(Base, :^))((getproperty(Base, :-))(x₅, x₃), 2), (getproperty(Base, :^))((getproperty(Base, :-))(x₁₂, x₁₀), 2)), 1.5))), (getproperty(Base, :/))((getproperty(Base, :*))(4, (getproperty(Base, :-))(x₄, x₅)), (getproperty(Base, :^))((getproperty(Base, :+))((getproperty(Base, :^))((getproperty(Base, :-))(x₅, x₄), 2), (getproperty(Base, :^))((getproperty(Base, :-))(x₁₂, x₁₁), 2)), 1.5))), (getproperty(Base, :/))((getproperty(Base, :*))(6, (getproperty(Base, :-))(x₆, x₅)), (getproperty(Base, :^))((getproperty(Base, :+))((getproperty(Base, :^))((getproperty(Base, :-))(x₅, x₆), 2), (getproperty(Base, :^))((getproperty(Base, :-))(x₁₂, x₁₃), 2)), 1.5))), (getproperty(Base, :/))((getproperty(Base, :*))(7, (getproperty(Base, :-))(x₇, x₅)), (getproperty(Base, :^))((getproperty(Base, :+))((getproperty(Base, :^))((getproperty(Base, :-))(x₅, x₇), 2), (getproperty(Base, :^))((getproperty(Base, :-))(x₁₂, x₁₄), 2)), 1.5)))
var"##MTIIPVar#313"[20] = (getproperty(Base, :+))((getproperty(Base, :+))((getproperty(Base, :+))((getproperty(Base, :+))((getproperty(Base, :+))((getproperty(Base, :+))(0, (getproperty(Base, :/))((getproperty(Base, :*))(1, (getproperty(Base, :-))(x₁, x₆)), (getproperty(Base, :^))((getproperty(Base, :+))((getproperty(Base, :^))((getproperty(Base, :-))(x₆, x₁), 2), (getproperty(Base, :^))((getproperty(Base, :-))(x₁₃, x₈), 2)), 1.5))), (getproperty(Base, :/))((getproperty(Base, :*))(2, (getproperty(Base, :-))(x₂, x₆)), (getproperty(Base, :^))((getproperty(Base, :+))((getproperty(Base, :^))((getproperty(Base, :-))(x₆, x₂), 2), (getproperty(Base, :^))((getproperty(Base, :-))(x₁₃, x₉), 2)), 1.5))), (getproperty(Base, :/))((getproperty(Base, :*))(3, (getproperty(Base, :-))(x₃, x₆)), (getproperty(Base, :^))((getproperty(Base, :+))((getproperty(Base, :^))((getproperty(Base, :-))(x₆, x₃), 2), (getproperty(Base, :^))((getproperty(Base, :-))(x₁₃, x₁₀), 2)), 1.5))), (getproperty(Base, :/))((getproperty(Base, :*))(4, (getproperty(Base, :-))(x₄, x₆)), (getproperty(Base, :^))((getproperty(Base, :+))((getproperty(Base, :^))((getproperty(Base, :-))(x₆, x₄), 2), (getproperty(Base, :^))((getproperty(Base, :-))(x₁₃, x₁₁), 2)), 1.5))), (getproperty(Base, :/))((getproperty(Base, :*))(5, (getproperty(Base, :-))(x₅, x₆)), (getproperty(Base, :^))((getproperty(Base, :+))((getproperty(Base, :^))((getproperty(Base, :-))(x₆, x₅), 2), (getproperty(Base, :^))((getproperty(Base, :-))(x₁₃, x₁₂), 2)), 1.5))), (getproperty(Base, :/))((getproperty(Base, :*))(7, (getproperty(Base, :-))(x₇, x₆)), (getproperty(Base, :^))((getproperty(Base, :+))((getproperty(Base, :^))((getproperty(Base, :-))(x₆, x₇), 2), (getproperty(Base, :^))((getproperty(Base, :-))(x₁₃, x₁₄), 2)), 1.5)))
var"##MTIIPVar#313"[21] = (getproperty(Base, :+))((getproperty(Base, :+))((getproperty(Base, :+))((getproperty(Base, :+))((getproperty(Base, :+))((getproperty(Base, :+))(0, (getproperty(Base, :/))((getproperty(Base, :*))(1, (getproperty(Base, :-))(x₁, x₇)), (getproperty(Base, :^))((getproperty(Base, :+))((getproperty(Base, :^))((getproperty(Base, :-))(x₇, x₁), 2), (getproperty(Base, :^))((getproperty(Base, :-))(x₁₄, x₈), 2)), 1.5))), (getproperty(Base, :/))((getproperty(Base, :*))(2, (getproperty(Base, :-))(x₂, x₇)), (getproperty(Base, :^))((getproperty(Base, :+))((getproperty(Base, :^))((getproperty(Base, :-))(x₇, x₂), 2), (getproperty(Base, :^))((getproperty(Base, :-))(x₁₄, x₉), 2)), 1.5))), (getproperty(Base, :/))((getproperty(Base, :*))(3, (getproperty(Base, :-))(x₃, x₇)), (getproperty(Base, :^))((getproperty(Base, :+))((getproperty(Base, :^))((getproperty(Base, :-))(x₇, x₃), 2), (getproperty(Base, :^))((getproperty(Base, :-))(x₁₄, x₁₀), 2)), 1.5))), (getproperty(Base, :/))((getproperty(Base, :*))(4, (getproperty(Base, :-))(x₄, x₇)), (getproperty(Base, :^))((getproperty(Base, :+))((getproperty(Base, :^))((getproperty(Base, :-))(x₇, x₄), 2), (getproperty(Base, :^))((getproperty(Base, :-))(x₁₄, x₁₁), 2)), 1.5))), (getproperty(Base, :/))((getproperty(Base, :*))(5, (getproperty(Base, :-))(x₅, x₇)), (getproperty(Base, :^))((getproperty(Base, :+))((getproperty(Base, :^))((getproperty(Base, :-))(x₇, x₅), 2), (getproperty(Base, :^))((getproperty(Base, :-))(x₁₄, x₁₂), 2)), 1.5))), (getproperty(Base, :/))((getproperty(Base, :*))(6, (getproperty(Base, :-))(x₆, x₇)), (getproperty(Base, :^))((getproperty(Base, :+))((getproperty(Base, :^))((getproperty(Base, :-))(x₇, x₆), 2), (getproperty(Base, :^))((getproperty(Base, :-))(x₁₄, x₁₃), 2)), 1.5)))
var"##MTIIPVar#313"[22] = (getproperty(Base, :+))((getproperty(Base, :+))((getproperty(Base, :+))((getproperty(Base, :+))((getproperty(Base, :+))((getproperty(Base, :+))(0, (getproperty(Base, :/))((getproperty(Base, :*))(2, (getproperty(Base, :-))(x₉, x₈)), (getproperty(Base, :^))((getproperty(Base, :+))((getproperty(Base, :^))((getproperty(Base, :-))(x₁, x₂), 2), (getproperty(Base, :^))((getproperty(Base, :-))(x₈, x₉), 2)), 1.5))), (getproperty(Base, :/))((getproperty(Base, :*))(3, (getproperty(Base, :-))(x₁₀, x₈)), (getproperty(Base, :^))((getproperty(Base, :+))((getproperty(Base, :^))((getproperty(Base, :-))(x₁, x₃), 2), (getproperty(Base, :^))((getproperty(Base, :-))(x₈, x₁₀), 2)), 1.5))), (getproperty(Base, :/))((getproperty(Base, :*))(4, (getproperty(Base, :-))(x₁₁, x₈)), (getproperty(Base, :^))((getproperty(Base, :+))((getproperty(Base, :^))((getproperty(Base, :-))(x₁, x₄), 2), (getproperty(Base, :^))((getproperty(Base, :-))(x₈, x₁₁), 2)), 1.5))), (getproperty(Base, :/))((getproperty(Base, :*))(5, (getproperty(Base, :-))(x₁₂, x₈)), (getproperty(Base, :^))((getproperty(Base, :+))((getproperty(Base, :^))((getproperty(Base, :-))(x₁, x₅), 2), (getproperty(Base, :^))((getproperty(Base, :-))(x₈, x₁₂), 2)), 1.5))), (getproperty(Base, :/))((getproperty(Base, :*))(6, (getproperty(Base, :-))(x₁₃, x₈)), (getproperty(Base, :^))((getproperty(Base, :+))((getproperty(Base, :^))((getproperty(Base, :-))(x₁, x₆), 2), (getproperty(Base, :^))((getproperty(Base, :-))(x₈, x₁₃), 2)), 1.5))), (getproperty(Base, :/))((getproperty(Base, :*))(7, (getproperty(Base, :-))(x₁₄, x₈)), (getproperty(Base, :^))((getproperty(Base, :+))((getproperty(Base, :^))((getproperty(Base, :-))(x₁, x₇), 2), (getproperty(Base, :^))((getproperty(Base, :-))(x₈, x₁₄), 2)), 1.5)))
var"##MTIIPVar#313"[23] = (getproperty(Base, :+))((getproperty(Base, :+))((getproperty(Base, :+))((getproperty(Base, :+))((getproperty(Base, :+))((getproperty(Base, :+))(0, (getproperty(Base, :/))((getproperty(Base, :*))(1, (getproperty(Base, :-))(x₈, x₉)), (getproperty(Base, :^))((getproperty(Base, :+))((getproperty(Base, :^))((getproperty(Base, :-))(x₂, x₁), 2), (getproperty(Base, :^))((getproperty(Base, :-))(x₉, x₈), 2)), 1.5))), (getproperty(Base, :/))((getproperty(Base, :*))(3, (getproperty(Base, :-))(x₁₀, x₉)), (getproperty(Base, :^))((getproperty(Base, :+))((getproperty(Base, :^))((getproperty(Base, :-))(x₂, x₃), 2), (getproperty(Base, :^))((getproperty(Base, :-))(x₉, x₁₀), 2)), 1.5))), (getproperty(Base, :/))((getproperty(Base, :*))(4, (getproperty(Base, :-))(x₁₁, x₉)), (getproperty(Base, :^))((getproperty(Base, :+))((getproperty(Base, :^))((getproperty(Base, :-))(x₂, x₄), 2), (getproperty(Base, :^))((getproperty(Base, :-))(x₉, x₁₁), 2)), 1.5))), (getproperty(Base, :/))((getproperty(Base, :*))(5, (getproperty(Base, :-))(x₁₂, x₉)), (getproperty(Base, :^))((getproperty(Base, :+))((getproperty(Base, :^))((getproperty(Base, :-))(x₂, x₅), 2), (getproperty(Base, :^))((getproperty(Base, :-))(x₉, x₁₂), 2)), 1.5))), (getproperty(Base, :/))((getproperty(Base, :*))(6, (getproperty(Base, :-))(x₁₃, x₉)), (getproperty(Base, :^))((getproperty(Base, :+))((getproperty(Base, :^))((getproperty(Base, :-))(x₂, x₆), 2), (getproperty(Base, :^))((getproperty(Base, :-))(x₉, x₁₃), 2)), 1.5))), (getproperty(Base, :/))((getproperty(Base, :*))(7, (getproperty(Base, :-))(x₁₄, x₉)), (getproperty(Base, :^))((getproperty(Base, :+))((getproperty(Base, :^))((getproperty(Base, :-))(x₂, x₇), 2), (getproperty(Base, :^))((getproperty(Base, :-))(x₉, x₁₄), 2)), 1.5)))
var"##MTIIPVar#313"[24] = (getproperty(Base, :+))((getproperty(Base, :+))((getproperty(Base, :+))((getproperty(Base, :+))((getproperty(Base, :+))((getproperty(Base, :+))(0, (getproperty(Base, :/))((getproperty(Base, :*))(1, (getproperty(Base, :-))(x₈, x₁₀)), (getproperty(Base, :^))((getproperty(Base, :+))((getproperty(Base, :^))((getproperty(Base, :-))(x₃, x₁), 2), (getproperty(Base, :^))((getproperty(Base, :-))(x₁₀, x₈), 2)), 1.5))), (getproperty(Base, :/))((getproperty(Base, :*))(2, (getproperty(Base, :-))(x₉, x₁₀)), (getproperty(Base, :^))((getproperty(Base, :+))((getproperty(Base, :^))((getproperty(Base, :-))(x₃, x₂), 2), (getproperty(Base, :^))((getproperty(Base, :-))(x₁₀, x₉), 2)), 1.5))), (getproperty(Base, :/))((getproperty(Base, :*))(4, (getproperty(Base, :-))(x₁₁, x₁₀)), (getproperty(Base, :^))((getproperty(Base, :+))((getproperty(Base, :^))((getproperty(Base, :-))(x₃, x₄), 2), (getproperty(Base, :^))((getproperty(Base, :-))(x₁₀, x₁₁), 2)), 1.5))), (getproperty(Base, :/))((getproperty(Base, :*))(5, (getproperty(Base, :-))(x₁₂, x₁₀)), (getproperty(Base, :^))((getproperty(Base, :+))((getproperty(Base, :^))((getproperty(Base, :-))(x₃, x₅), 2), (getproperty(Base, :^))((getproperty(Base, :-))(x₁₀, x₁₂), 2)), 1.5))), (getproperty(Base, :/))((getproperty(Base, :*))(6, (getproperty(Base, :-))(x₁₃, x₁₀)), (getproperty(Base, :^))((getproperty(Base, :+))((getproperty(Base, :^))((getproperty(Base, :-))(x₃, x₆), 2), (getproperty(Base, :^))((getproperty(Base, :-))(x₁₀, x₁₃), 2)), 1.5))), (getproperty(Base, :/))((getproperty(Base, :*))(7, (getproperty(Base, :-))(x₁₄, x₁₀)), (getproperty(Base, :^))((getproperty(Base, :+))((getproperty(Base, :^))((getproperty(Base, :-))(x₃, x₇), 2), (getproperty(Base, :^))((getproperty(Base, :-))(x₁₀, x₁₄), 2)), 1.5)))
var"##MTIIPVar#313"[25] = (getproperty(Base, :+))((getproperty(Base, :+))((getproperty(Base, :+))((getproperty(Base, :+))((getproperty(Base, :+))((getproperty(Base, :+))(0, (getproperty(Base, :/))((getproperty(Base, :*))(1, (getproperty(Base, :-))(x₈, x₁₁)), (getproperty(Base, :^))((getproperty(Base, :+))((getproperty(Base, :^))((getproperty(Base, :-))(x₄, x₁), 2), (getproperty(Base, :^))((getproperty(Base, :-))(x₁₁, x₈), 2)), 1.5))), (getproperty(Base, :/))((getproperty(Base, :*))(2, (getproperty(Base, :-))(x₉, x₁₁)), (getproperty(Base, :^))((getproperty(Base, :+))((getproperty(Base, :^))((getproperty(Base, :-))(x₄, x₂), 2), (getproperty(Base, :^))((getproperty(Base, :-))(x₁₁, x₉), 2)), 1.5))), (getproperty(Base, :/))((getproperty(Base, :*))(3, (getproperty(Base, :-))(x₁₀, x₁₁)), (getproperty(Base, :^))((getproperty(Base, :+))((getproperty(Base, :^))((getproperty(Base, :-))(x₄, x₃), 2), (getproperty(Base, :^))((getproperty(Base, :-))(x₁₁, x₁₀), 2)), 1.5))), (getproperty(Base, :/))((getproperty(Base, :*))(5, (getproperty(Base, :-))(x₁₂, x₁₁)), (getproperty(Base, :^))((getproperty(Base, :+))((getproperty(Base, :^))((getproperty(Base, :-))(x₄, x₅), 2), (getproperty(Base, :^))((getproperty(Base, :-))(x₁₁, x₁₂), 2)), 1.5))), (getproperty(Base, :/))((getproperty(Base, :*))(6, (getproperty(Base, :-))(x₁₃, x₁₁)), (getproperty(Base, :^))((getproperty(Base, :+))((getproperty(Base, :^))((getproperty(Base, :-))(x₄, x₆), 2), (getproperty(Base, :^))((getproperty(Base, :-))(x₁₁, x₁₃), 2)), 1.5))), (getproperty(Base, :/))((getproperty(Base, :*))(7, (getproperty(Base, :-))(x₁₄, x₁₁)), (getproperty(Base, :^))((getproperty(Base, :+))((getproperty(Base, :^))((getproperty(Base, :-))(x₄, x₇), 2), (getproperty(Base, :^))((getproperty(Base, :-))(x₁₁, x₁₄), 2)), 1.5)))
var"##MTIIPVar#313"[26] = (getproperty(Base, :+))((getproperty(Base, :+))((getproperty(Base, :+))((getproperty(Base, :+))((getproperty(Base, :+))((getproperty(Base, :+))(0, (getproperty(Base, :/))((getproperty(Base, :*))(1, (getproperty(Base, :-))(x₈, x₁₂)), (getproperty(Base, :^))((getproperty(Base, :+))((getproperty(Base, :^))((getproperty(Base, :-))(x₅, x₁), 2), (getproperty(Base, :^))((getproperty(Base, :-))(x₁₂, x₈), 2)), 1.5))), (getproperty(Base, :/))((getproperty(Base, :*))(2, (getproperty(Base, :-))(x₉, x₁₂)), (getproperty(Base, :^))((getproperty(Base, :+))((getproperty(Base, :^))((getproperty(Base, :-))(x₅, x₂), 2), (getproperty(Base, :^))((getproperty(Base, :-))(x₁₂, x₉), 2)), 1.5))), (getproperty(Base, :/))((getproperty(Base, :*))(3, (getproperty(Base, :-))(x₁₀, x₁₂)), (getproperty(Base, :^))((getproperty(Base, :+))((getproperty(Base, :^))((getproperty(Base, :-))(x₅, x₃), 2), (getproperty(Base, :^))((getproperty(Base, :-))(x₁₂, x₁₀), 2)), 1.5))), (getproperty(Base, :/))((getproperty(Base, :*))(4, (getproperty(Base, :-))(x₁₁, x₁₂)), (getproperty(Base, :^))((getproperty(Base, :+))((getproperty(Base, :^))((getproperty(Base, :-))(x₅, x₄), 2), (getproperty(Base, :^))((getproperty(Base, :-))(x₁₂, x₁₁), 2)), 1.5))), (getproperty(Base, :/))((getproperty(Base, :*))(6, (getproperty(Base, :-))(x₁₃, x₁₂)), (getproperty(Base, :^))((getproperty(Base, :+))((getproperty(Base, :^))((getproperty(Base, :-))(x₅, x₆), 2), (getproperty(Base, :^))((getproperty(Base, :-))(x₁₂, x₁₃), 2)), 1.5))), (getproperty(Base, :/))((getproperty(Base, :*))(7, (getproperty(Base, :-))(x₁₄, x₁₂)), (getproperty(Base, :^))((getproperty(Base, :+))((getproperty(Base, :^))((getproperty(Base, :-))(x₅, x₇), 2), (getproperty(Base, :^))((getproperty(Base, :-))(x₁₂, x₁₄), 2)), 1.5)))
var"##MTIIPVar#313"[27] = (getproperty(Base, :+))((getproperty(Base, :+))((getproperty(Base, :+))((getproperty(Base, :+))((getproperty(Base, :+))((getproperty(Base, :+))(0, (getproperty(Base, :/))((getproperty(Base, :*))(1, (getproperty(Base, :-))(x₈, x₁₃)), (getproperty(Base, :^))((getproperty(Base, :+))((getproperty(Base, :^))((getproperty(Base, :-))(x₆, x₁), 2), (getproperty(Base, :^))((getproperty(Base, :-))(x₁₃, x₈), 2)), 1.5))), (getproperty(Base, :/))((getproperty(Base, :*))(2, (getproperty(Base, :-))(x₉, x₁₃)), (getproperty(Base, :^))((getproperty(Base, :+))((getproperty(Base, :^))((getproperty(Base, :-))(x₆, x₂), 2), (getproperty(Base, :^))((getproperty(Base, :-))(x₁₃, x₉), 2)), 1.5))), (getproperty(Base, :/))((getproperty(Base, :*))(3, (getproperty(Base, :-))(x₁₀, x₁₃)), (getproperty(Base, :^))((getproperty(Base, :+))((getproperty(Base, :^))((getproperty(Base, :-))(x₆, x₃), 2), (getproperty(Base, :^))((getproperty(Base, :-))(x₁₃, x₁₀), 2)), 1.5))), (getproperty(Base, :/))((getproperty(Base, :*))(4, (getproperty(Base, :-))(x₁₁, x₁₃)), (getproperty(Base, :^))((getproperty(Base, :+))((getproperty(Base, :^))((getproperty(Base, :-))(x₆, x₄), 2), (getproperty(Base, :^))((getproperty(Base, :-))(x₁₃, x₁₁), 2)), 1.5))), (getproperty(Base, :/))((getproperty(Base, :*))(5, (getproperty(Base, :-))(x₁₂, x₁₃)), (getproperty(Base, :^))((getproperty(Base, :+))((getproperty(Base, :^))((getproperty(Base, :-))(x₆, x₅), 2), (getproperty(Base, :^))((getproperty(Base, :-))(x₁₃, x₁₂), 2)), 1.5))), (getproperty(Base, :/))((getproperty(Base, :*))(7, (getproperty(Base, :-))(x₁₄, x₁₃)), (getproperty(Base, :^))((getproperty(Base, :+))((getproperty(Base, :^))((getproperty(Base, :-))(x₆, x₇), 2), (getproperty(Base, :^))((getproperty(Base, :-))(x₁₃, x₁₄), 2)), 1.5)))
var"##MTIIPVar#313"[28] = (getproperty(Base, :+))((getproperty(Base, :+))((getproperty(Base, :+))((getproperty(Base, :+))((getproperty(Base, :+))((getproperty(Base, :+))(0, (getproperty(Base, :/))((getproperty(Base, :*))(1, (getproperty(Base, :-))(x₈, x₁₄)), (getproperty(Base, :^))((getproperty(Base, :+))((getproperty(Base, :^))((getproperty(Base, :-))(x₇, x₁), 2), (getproperty(Base, :^))((getproperty(Base, :-))(x₁₄, x₈), 2)), 1.5))), (getproperty(Base, :/))((getproperty(Base, :*))(2, (getproperty(Base, :-))(x₉, x₁₄)), (getproperty(Base, :^))((getproperty(Base, :+))((getproperty(Base, :^))((getproperty(Base, :-))(x₇, x₂), 2), (getproperty(Base, :^))((getproperty(Base, :-))(x₁₄, x₉), 2)), 1.5))), (getproperty(Base, :/))((getproperty(Base, :*))(3, (getproperty(Base, :-))(x₁₀, x₁₄)), (getproperty(Base, :^))((getproperty(Base, :+))((getproperty(Base, :^))((getproperty(Base, :-))(x₇, x₃), 2), (getproperty(Base, :^))((getproperty(Base, :-))(x₁₄, x₁₀), 2)), 1.5))), (getproperty(Base, :/))((getproperty(Base, :*))(4, (getproperty(Base, :-))(x₁₁, x₁₄)), (getproperty(Base, :^))((getproperty(Base, :+))((getproperty(Base, :^))((getproperty(Base, :-))(x₇, x₄), 2), (getproperty(Base, :^))((getproperty(Base, :-))(x₁₄, x₁₁), 2)), 1.5))), (getproperty(Base, :/))((getproperty(Base, :*))(5, (getproperty(Base, :-))(x₁₂, x₁₄)), (getproperty(Base, :^))((getproperty(Base, :+))((getproperty(Base, :^))((getproperty(Base, :-))(x₇, x₅), 2), (getproperty(Base, :^))((getproperty(Base, :-))(x₁₄, x₁₂), 2)), 1.5))), (getproperty(Base, :/))((getproperty(Base, :*))(6, (getproperty(Base, :-))(x₁₃, x₁₄)), (getproperty(Base, :^))((getproperty(Base, :+))((getproperty(Base, :^))((getproperty(Base, :-))(x₇, x₆), 2), (getproperty(Base, :^))((getproperty(Base, :-))(x₁₄, x₁₃), 2)), 1.5)))
nothing
end
# Now test against the `@taylorize` accelerated form
taylor_odef = ODEFunction{true}(taylor_f, syms = [:x₁, :x₂, :x₃, :x₄, :x₅, :x₆, :x₇, :x₈, :x₉, :x₁₀, :x₁₁, :x₁₂, :x₁₃, :x₁₄, :x₁₅, :x₁₆, :x₁₇, :x₁₈, :x₁₉, :x₂₀, :x₂₁, :x₂₂, :x₂₃, :x₂₄, :x₂₅, :x₂₆, :x₂₇, :x₂₈])
taylor_prob = ODEProblem(taylor_odef,[3.0,3.0,-1.0,-3.0,2.0,-2.0,2.0,3.0,-3.0,2.0,0,0,-4.0,4.0,0,0,0,0,0,1.75,-1.5,0,0,0,-1.25,1,0,0],(0.0,3.0))
setups = [Dict(:alg=>Tsit5())
Dict(:alg=>Vern6())
Dict(:alg=>Vern7())
Dict(:alg=>Vern9())
Dict(:alg=>TaylorMethod(12))
]
wp = WorkPrecisionSet(taylor_prob,abstols,reltols,setups;appxsol=test_sol,numruns=100,maxiters=10000,verbose=false)
plot(wp)
savefig("taylorspeed.png")
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment