|
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") |