diff --git a/benchmarks/NonStiffODE/Project.toml b/benchmarks/NonStiffODE/Project.toml index f8c2cc8b3..0732e4b5f 100644 --- a/benchmarks/NonStiffODE/Project.toml +++ b/benchmarks/NonStiffODE/Project.toml @@ -20,7 +20,7 @@ Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" DiffEqBase = "6" DiffEqDevTools = "2.30" LSODA = "0.7" -ModelingToolkit = "9, 10" +ModelingToolkit = "10" ODEInterface = "0.5" ODEInterfaceDiffEq = "3.10" OrdinaryDiffEq = "6" diff --git a/benchmarks/NonStiffODE/enright_pryce.jl b/benchmarks/NonStiffODE/enright_pryce.jl index 988bb363e..d6f6aee80 100644 --- a/benchmarks/NonStiffODE/enright_pryce.jl +++ b/benchmarks/NonStiffODE/enright_pryce.jl @@ -9,168 +9,152 @@ y = collect(y) D = Differential(t) # Stiff -sa1sys = complete(let - sa1eqs = [ - D(y[1]) ~ -0.5 * y[1], D(y[2]) ~ -y[2], D(y[3]) ~ -100 * y[3], D(y[4]) ~ -90 * y[4]] - ODESystem(sa1eqs, t, name = :sa1) -end) +sa1eqs = [ + D(y[1]) ~ -0.5 * y[1], D(y[2]) ~ -y[2], D(y[3]) ~ -100 * y[3], D(y[4]) ~ -90 * y[4]] +@mtkcompile sa1sys = ODESystem(sa1eqs, t) sa1prob = ODEProblem{false}(sa1sys, y[1:4] .=> 1.0, (0, 20.0), dt = 1e-2; u0_constructor = x -> SVector(x...)) -sa2sys = complete(let - sa2eqs = [D(y[1]) ~ -1800 * y[1] + 900 * y[2]] - for i in 2:8 - push!(sa2eqs, D(y[i]) ~ y[i - 1] - 2 * y[i] + y[i + 1]) - end - push!(sa2eqs, D(y[9]) ~ 1000 * y[8] - 2000 * y[9] + 1000) - ODESystem(sa2eqs, t, name = :sa2) -end) +sa2eqs = [D(y[1]) ~ -1800 * y[1] + 900 * y[2]] +for i in 2:8 + push!(sa2eqs, D(y[i]) ~ y[i - 1] - 2 * y[i] + y[i + 1]) +end +push!(sa2eqs, D(y[9]) ~ 1000 * y[8] - 2000 * y[9] + 1000) +@mtkcompile sa2sys = ODESystem(sa2eqs, t) sa2prob = ODEProblem{false}(sa2sys, y[1:9] .=> 0.0, (0, 120.0), dt = 5e-4; u0_constructor = x -> SVector(x...)) -sa3sys = complete(let - sa3eqs = [ - D(y[1]) ~ -1e4 * y[1] + 100 * y[2] - 10 * y[3] + y[4], - D(y[2]) ~ -1e3 * y[2] + 10 * y[3] - 10 * y[4], - D(y[3]) ~ -y[3] + 10 * y[4], - D(y[4]) ~ -0.1 * y[4] - ] - ODESystem(sa3eqs, t, name = :sa3) -end) +sa3eqs = [ + D(y[1]) ~ -1e4 * y[1] + 100 * y[2] - 10 * y[3] + y[4], + D(y[2]) ~ -1e3 * y[2] + 10 * y[3] - 10 * y[4], + D(y[3]) ~ -y[3] + 10 * y[4], + D(y[4]) ~ -0.1 * y[4] +] +@mtkcompile sa3sys = ODESystem(sa3eqs, t) sa3prob = ODEProblem{false}(sa3sys, y[1:4] .=> 1.0, (0, 20.0), dt = 1e-5; u0_constructor = x -> SVector(x...)) -sa4sys = complete(let - sa4eqs = [D(y[i]) ~ -i^5 * y[i] for i in 1:10] - ODESystem(sa4eqs, t, name = :sa4) -end) +sa4eqs = [D(y[i]) ~ -i^5 * y[i] for i in 1:10] +@named sa4sys_raw = ODESystem(sa4eqs, t) +sa4sys = structural_simplify(sa4sys_raw) sa4prob = ODEProblem{false}(sa4sys, y[1:10] .=> 1.0, (0, 1.0), dt = 1e-5; u0_constructor = x -> SVector(x...)) const SA_PROBLEMS = [sa1prob, sa2prob, sa3prob, sa4prob] -sb1sys = complete(let - sb1eqs = [D(y[1]) ~ -y[1] + y[2], - D(y[2]) ~ -100y[1] - y[2], - D(y[3]) ~ -100y[3] + y[4], - D(y[4]) ~ -10_000y[3] - 100y[4]] - - ODESystem(sb1eqs, t, name = :sb1) -end) +sb1eqs = [D(y[1]) ~ -y[1] + y[2], + D(y[2]) ~ -100y[1] - y[2], + D(y[3]) ~ -100y[3] + y[4], + D(y[4]) ~ -10_000y[3] - 100y[4]] +@named sb1sys_raw = ODESystem(sb1eqs, t) +sb1sys = structural_simplify(sb1sys_raw) sb1prob = ODEProblem{false}( sb1sys, [y[[1, 3]] .=> 1.0; y[[2, 4]] .=> 0.0], (0, 20.0), dt = 7e-3; u0_constructor = x -> SVector(x...)) @parameters α -sb2sys = complete(let - sb2eqs = [D(y[1]) ~ -10y[1] + α * y[2], - D(y[2]) ~ -α * y[1] - 10 * y[2], - D(y[3]) ~ -4y[3], - D(y[4]) ~ -y[4], - D(y[5]) ~ -0.5y[5], - D(y[6]) ~ -0.1y[6]] - - ODESystem(sb2eqs, t, name = :sb2) -end) +sb2eqs = [D(y[1]) ~ -10y[1] + α * y[2], + D(y[2]) ~ -α * y[1] - 10 * y[2], + D(y[3]) ~ -4y[3], + D(y[4]) ~ -y[4], + D(y[5]) ~ -0.5y[5], + D(y[6]) ~ -0.1y[6]] +@named sb2sys_raw = ODESystem(sb2eqs, t) +sb2sys = structural_simplify(sb2sys_raw) sb2prob = ODEProblem{false}(sb2sys, y .=> 1.0, (0, 20.0), [α => 3], dt = 1e-2) sb3prob = ODEProblem{false}(sb2sys, y .=> 1.0, (0, 20.0), [α => 8], dt = 1e-2) sb4prob = ODEProblem{false}(sb2sys, y .=> 1.0, (0, 20.0), [α => 25], dt = 1e-2) sb5prob = ODEProblem{false}(sb2sys, y .=> 1.0, (0, 20.0), [α => 100], dt = 1e-2) -sc1sys = complete(let - sc1eqs = [ - D(y[1]) ~ -y[1] + y[2]^2 + y[3]^2 + y[4]^2, - D(y[2]) ~ -10y[2] + 10 * (y[3]^2 + y[4]^2), - D(y[3]) ~ -40y[3] + 40 * y[4]^2, - D(y[4]) ~ -100y[4] + 2] - - ODESystem(sc1eqs, t, name = :sc1) -end) +sc1eqs = [ + D(y[1]) ~ -y[1] + y[2]^2 + y[3]^2 + y[4]^2, + D(y[2]) ~ -10y[2] + 10 * (y[3]^2 + y[4]^2), + D(y[3]) ~ -40y[3] + 40 * y[4]^2, + D(y[4]) ~ -100y[4] + 2] +@named sc1sys_raw = ODESystem(sc1eqs, t) +sc1sys = structural_simplify(sc1sys_raw) sc1prob = ODEProblem{false}(sc1sys, y .=> 1.0, (0, 20.0), dt = 1e-2; u0_constructor = x -> SVector(x...)) @parameters β -sc2sys = complete(let - sc2eqs = [D(y[1]) ~ -y[1] + 2, - D(y[2]) ~ -10y[2] + β * y[1]^2, - D(y[3]) ~ -40y[3] + 4β * (y[1]^2 + y[2]^2), - D(y[4]) ~ -100y[4] + 10β * (y[1]^2 + y[2]^2 + y[3]^2)] - - ODESystem(sc2eqs, t, name = :sc2) -end) +sc2eqs = [D(y[1]) ~ -y[1] + 2, + D(y[2]) ~ -10y[2] + β * y[1]^2, + D(y[3]) ~ -40y[3] + 4β * (y[1]^2 + y[2]^2), + D(y[4]) ~ -100y[4] + 10β * (y[1]^2 + y[2]^2 + y[3]^2)] +@mtkcompile sc2sys = ODESystem(sc2eqs, t) sc2prob = ODEProblem{false}(sc2sys, y .=> 1.0, (0, 20.0), [β => 0.1], dt = 1e-2; u0_constructor = x -> SVector(x...)) sc3prob = ODEProblem{false}(sc2sys, y .=> 1.0, (0, 20.0), [β => 1.0], dt = 1e-2; u0_constructor = x -> SVector(x...)) sc4prob = ODEProblem{false}(sc2sys, y .=> 1.0, (0, 20.0), [β => 10.0], dt = 1e-2; u0_constructor = x -> SVector(x...)) sc5prob = ODEProblem{false}(sc2sys, y .=> 1.0, (0, 20.0), [β => 20.0], dt = 1e-2; u0_constructor = x -> SVector(x...)) -sd1sys = complete(let +sd1sys = let sd1eqs = [D(y[1]) ~ 0.2 * (y[2] - y[1]), D(y[2]) ~ 10y[1] - (60 - 0.125y[3]) * y[2] + 0.125y[3], D(y[3]) ~ 1] - ODESystem(sd1eqs, t, name = :sd1) -end) + structural_simplify(ODESystem(sd1eqs, t; name=Symbol(gensym("sys")))) +end sd1prob = ODEProblem{false}(sd1sys, y .=> 0.0, (0, 400.0), [β => 0.1], dt = 1.7e-2; u0_constructor = x -> SVector(x...)) -sd2sys = complete(let +sd2sys = let sd2eqs = [D(y[1]) ~ -0.04y[1] + 0.01 * (y[2] * y[3]), D(y[2]) ~ 400y[1] - 100 * (y[2] * y[3]) - 3000 * y[2]^2, D(y[3]) ~ 30y[2]^2] - ODESystem(sd2eqs, t, name = :sd2) -end) + structural_simplify(ODESystem(sd2eqs, t; name=Symbol(gensym("sys")))) +end sd2prob = ODEProblem{false}( sd2sys, [y[1] => 1.0; y[2:3] .=> 0.0], (0, 40.0), dt = 1e-5, cse = true; u0_constructor = x -> SVector(x...)) -sd3sys = complete(let +sd3sys = let sd3eqs = [D(y[1]) ~ y[3] - 100 * (y[1] * y[2]), D(y[2]) ~ y[3] + 2y[4] - 100 * (y[1] * y[2]) - 2e4 * y[2]^2, D(y[3]) ~ -y[3] + 100 * (y[1] * y[2]), D(y[4]) ~ -y[4] + 1e4 * y[2]^2] - ODESystem(sd3eqs, t, name = :sd3) -end) + structural_simplify(ODESystem(sd3eqs, t; name=Symbol(gensym("sys")))) +end sd3prob = ODEProblem{false}( sd3sys, [y[1:2] .=> 1; y[3:4] .=> 0.0], (0, 20.0), dt = 2.5e-5, cse = true; u0_constructor = x -> SVector(x...)) -sd4sys = complete(let +sd4sys = let sd4eqs = [D(y[1]) ~ -0.013y[1] - 1000 * (y[1] * y[3]), D(y[2]) ~ -2500 * (y[2] * y[3]), D(y[3]) ~ -0.013y[1] - 1000 * (y[1] * y[3]) - 2500 * (y[2] * y[3])] - ODESystem(sd4eqs, t, name = :sd4) -end) + structural_simplify(ODESystem(sd4eqs, t; name=Symbol(gensym("sys")))) +end sd4prob = ODEProblem{false}( sd4sys, [y[1:2] .=> 1; y[3] => 0.0], (0, 50.0), dt = 2.9e-4, cse = true; u0_constructor = x -> SVector(x...)) -sd5sys = complete(let +sd5sys = let sd5eqs = [D(y[1]) ~ 0.01 - (1 + (y[1] + 1000) * (y[1] + 1)) * (0.01 + y[1] + y[2]), D(y[2]) ~ 0.01 - (1 + y[2]^2) * (0.01 + y[1] + y[2])] - ODESystem(sd5eqs, t, name = :sd5) -end) + structural_simplify(ODESystem(sd5eqs, t; name=Symbol(gensym("sys")))) +end sd5prob = ODEProblem{false}(sd5sys, y[1:2] .=> 0.0, (0, 100.0), dt = 1e-4, cse = true; u0_constructor = x -> SVector(x...)) -sd6sys = complete(let +sd6sys = let sd6eqs = [D(y[1]) ~ -y[1] + 10^8 * y[3] * (1 - y[1]), D(y[2]) ~ -10y[2] + 3e7 * y[3] * (1 - y[2]), D(y[3]) ~ -(-y[1] + 10^8 * y[3] * (1 - y[1])) - (-10y[2] + 3e7 * y[3] * (1 - y[2])) ] - ODESystem(sd6eqs, t, name = :sd6) -end) + structural_simplify(ODESystem(sd6eqs, t; name=Symbol(gensym("sys")))) +end sd6prob = ODEProblem{false}( sd6sys, [y[1] => 1.0; y[2:3] .=> 0.0], (0, 1.0), dt = 3.3e-8, cse = true; u0_constructor = x -> SVector(x...)) -se1sys = complete(let +se1sys = let Γ = 100 se1eqs = [D(y[1]) ~ y[2], D(y[2]) ~ y[3], @@ -180,34 +164,34 @@ se1sys = complete(let (10 * exp(-y[4]^2) - 4Γ) * y[4] + 1 ] - ODESystem(se1eqs, t, name = :se1) -end) + structural_simplify(ODESystem(se1eqs, t; name=Symbol(gensym("sys")))) +end se1prob = ODEProblem{false}(se1sys, y .=> 0.0, (0, 1.0), dt = 6.8e-3, cse = true; u0_constructor = x -> SVector(x...)) -se2sys = complete(let +se2sys = let se2eqs = [D(y[1]) ~ y[2], D(y[2]) ~ 5 * (1 - y[1]^2) * y[2] - y[1] ] - ODESystem(se2eqs, t, name = :se2) -end) + structural_simplify(ODESystem(se2eqs, t; name=Symbol(gensym("sys")))) +end se2prob = ODEProblem{false}( se2sys, [y[1] => 2.0, y[2] => 0.0], (0, 1.0), dt = 1e-3, cse = true; u0_constructor = x -> SVector(x...)) -se3sys = complete(let +se3sys = let se3eqs = [D(y[1]) ~ -(55 + y[3]) * y[1] + 65 * y[2], D(y[2]) ~ 0.0785 * (y[1] - y[2]), D(y[3]) ~ 0.1 * y[1] ] - ODESystem(se3eqs, t, name = :se3) -end) + structural_simplify(ODESystem(se3eqs, t; name=Symbol(gensym("sys")))) +end se3prob = ODEProblem{false}(se3sys, [y[1:2] .=> 1.0; y[3] => 0.0], (0, 500.0), dt = 0.02; u0_constructor = x -> SVector(x...)) -se4sys = complete(let y = y +se4sys = let y = y y = y[1:4] U = ones(4, 4) U[diagind(U)] .= -1 @@ -217,13 +201,13 @@ se4sys = complete(let y = y A = [-10 -10 0 0; 10 -10 0 0; 0 0 1000 0; 0 0 0 0.01] se4eqs = D.(y) .~ -(U' * A * Z) + G - ODESystem(se4eqs, t, name = :se4) -end) + structural_simplify(ODESystem(se4eqs, t; name=Symbol(gensym("sys")))) +end se4prob = ODEProblem{false}(se4sys, [y[1] => 0.0; y[2] => -2.0; y[3:4] .=> -1.0], (0, 1000.0), dt = 1e-3, cse = true; u0_constructor = x -> SVector(x...)) -se5sys = complete(let +se5sys = let se5eqs = [ D(y[1]) ~ -7.89e-10 * y[1] - 1.1e7 * y[1] * y[3], D(y[2]) ~ 7.89e-10 * y[1] - 1.13e9 * y[2] * y[3], @@ -232,13 +216,13 @@ se5sys = complete(let D(y[4]) ~ 1.1e7 * y[1] * y[3] - 1.13e3 * y[4] ] - ODESystem(se5eqs, t, name = :se5) -end) + structural_simplify(ODESystem(se5eqs, t; name=Symbol(gensym("sys")))) +end se5prob = ODEProblem{false}( se5sys, [y[1] => 1.76e-3; y[2:4] .=> 0.0], (0, 1000.0), dt = 1e-3, cse = true; u0_constructor = x -> SVector(x...)) -sf1sys = complete(let +sf1sys = let k = exp(20.7 - 1500 / y[1]) sf1eqs = [ D(y[1]) ~ 1.3 * (y[3] - y[1]) + 10400 * k * y[2], @@ -247,26 +231,26 @@ sf1sys = complete(let D(y[4]) ~ 0.1 + 320 * y[2] - 321 * y[4] ] - ODESystem(sf1eqs, t, name = :sf1) -end) + structural_simplify(ODESystem(sf1eqs, t; name=Symbol(gensym("sys")))) +end sf1prob = ODEProblem{false}( sf1sys, [y[1] => 761.0; y[2] => 0.0; y[3] => 600.0; y[4] => 0.1], (0, 1000.0), dt = 1e-4, cse = true) -sf2sys = complete(let +sf2sys = let sf2eqs = [ D(y[1]) ~ -y[1] - y[1] * y[2] + 294 * y[2], D(y[2]) ~ y[1] * (1 - y[2]) / 98 - 3 * y[2] ] - ODESystem(sf2eqs, t, name = :sf2) -end) + structural_simplify(ODESystem(sf2eqs, t; name=Symbol(gensym("sys")))) +end sf2prob = ODEProblem{false}( sf2sys, [y[1] => 1.0; y[2] => 0.0], (0, 240.0), dt = 1e-2, cse = true; u0_constructor = x -> SVector(x...)) -sf3sys = complete(let +sf3sys = let sf3eqs = [ D(y[1]) ~ -1e7 * y[2] * y[1] + 10 * y[3], D(y[2]) ~ -1e7 * y[2] * y[1] - 1e7 * y[2] * y[5] + 10 * y[3] + 10 * y[4], @@ -275,26 +259,26 @@ sf3sys = complete(let D(y[5]) ~ 10 * y[4] - 1e7 * y[2] * y[5] ] - ODESystem(sf3eqs, t, name = :sf3) -end) + structural_simplify(ODESystem(sf3eqs, t; name=Symbol(gensym("sys")))) +end sf3prob = ODEProblem{false}( sf3sys, [y[1] => 4e-6; y[2] => 1e-6; y[3:5] .=> 0.0], (0, 100.0), dt = 1e-6, cse = true; u0_constructor = x -> SVector(x...)) -sf4sys = complete(let +sf4sys = let sf4eqs = [ D(y[1]) ~ 77.27 * (y[2] - y[1] * y[2] + y[1] - 8.375e-6 * y[1]^2), D(y[2]) ~ (-y[2] - y[1] * y[2] + y[3]) / 77.27, D(y[3]) ~ 0.161 * (y[1] - y[3]) ] - ODESystem(sf4eqs, t, name = :sf4) -end) + structural_simplify(ODESystem(sf4eqs, t; name=Symbol(gensym("sys")))) +end sf4prob = ODEProblem{false}( sf4sys, [y[1] => 4.0; y[2] => 1.1; y[3] => 4.0], (0, 300.0), dt = 1e-3, cse = true; u0_constructor = x -> SVector(x...)) -sf5sys = complete(let +sf5sys = let sf5eqs = [ D(y[1]) ~ 1e11 * (-3 * y[1] * y[2] + 0.0012 * y[4] - 9 * y[1] * y[3]), D(y[2]) ~ -3e11 * y[1] * y[2] + 2e7 * y[4], @@ -302,93 +286,91 @@ sf5sys = complete(let D(y[4]) ~ 1e11 * (3 * y[1] * y[2] - 0.0012 * y[4] + 9 * y[1] * y[3]) ] - ODESystem(sf5eqs, t, name = :sf5) -end) + structural_simplify(ODESystem(sf5eqs, t; name=Symbol(gensym("sys")))) +end sf5prob = ODEProblem{false}( sf5sys, [y[1] => 3.365e-7; y[2] => 8.261e-3; y[3] => 1.642e-3; y[4] => 9.38e-6], (0, 100.0), dt = 1e-7, cse = true; u0_constructor = x -> SVector(x...)) # Non-stiff -na1sys = complete(let y = y[1] - na1eqs = D(y) ~ -y - - ODESystem(na1eqs, t, name = :na1) -end) +y1 = y[1] +na1eqs = D(y1) ~ -y1 +@named na1sys_raw = ODESystem(na1eqs, t) +na1sys = structural_simplify(na1sys_raw) na1prob = ODEProblem{false}(na1sys, [y[1] => 1], (0, 20.0), cse = true; u0_constructor = x -> SVector(x...)) -na2sys = complete(let y = y[1] - na2eqs = D(y) ~ -y^3 / 2 - - ODESystem(na2eqs, t, name = :na2) -end) +y2 = y[1] +na2eqs = D(y2) ~ -y2^3 / 2 +@named na2sys_raw = ODESystem(na2eqs, t) +na2sys = structural_simplify(na2sys_raw) na2prob = ODEProblem{false}(na2sys, [y[1] => 1], (0, 20.0), cse = true; u0_constructor = x -> SVector(x...)) -na3sys = complete(let y = y[1] +na3sys = let y = y[1] na3eqs = D(y) ~ y * cos(t) - ODESystem(na3eqs, t, name = :na3) -end) + structural_simplify(ODESystem(na3eqs, t; name=Symbol(gensym("sys")))) +end na3prob = ODEProblem{false}(na3sys, [y[1] => 1], (0, 20.0), cse = true; u0_constructor = x -> SVector(x...)) -na4sys = complete(let y = y[1] +na4sys = let y = y[1] na4eqs = D(y) ~ y / 4 * (1 - y / 20) - ODESystem(na4eqs, t, name = :na4) -end) + structural_simplify(ODESystem(na4eqs, t; name=Symbol(gensym("sys")))) +end na4prob = ODEProblem{false}(na4sys, [y[1] => 1], (0, 20.0), cse = true; u0_constructor = x -> SVector(x...)) -na5sys = complete(let y = y[1] +na5sys = let y = y[1] na5eqs = D(y) ~ (y - t) / (y + t) - ODESystem(na5eqs, t, name = :na5) -end) + structural_simplify(ODESystem(na5eqs, t; name=Symbol(gensym("sys")))) +end na5prob = ODEProblem{false}(na5sys, [y[1] => 4], (0, 20.0), cse = true; u0_constructor = x -> SVector(x...)) const NA_PROBLEMS = [na1prob, na2prob, na3prob, na4prob, na5prob] -nb1sys = complete(let +nb1sys = let nb1eqs = [ D(y[1]) ~ 2 * (y[1] - y[1] * y[2]), D(y[2]) ~ -(y[2] - y[1] * y[2]) ] - ODESystem(nb1eqs, t, name = :nb1) -end) + structural_simplify(ODESystem(nb1eqs, t; name=Symbol(gensym("sys")))) +end nb1prob = ODEProblem{false}(nb1sys, [y[1] => 1.0, y[2] => 3], (0, 20.0), cse = true; u0_constructor = x -> SVector(x...)) -nb2sys = complete(let +nb2sys = let nb2eqs = [ D(y[1]) ~ -y[1] + y[2], D(y[2]) ~ y[1] - 2y[2] + y[3], D(y[3]) ~ y[2] - y[3] ] - ODESystem(nb2eqs, t, name = :nb2) -end) + structural_simplify(ODESystem(nb2eqs, t; name=Symbol(gensym("sys")))) +end nb2prob = ODEProblem{false}( nb2sys, [y[1] => 2.0, y[2] => 0.0, y[3] => 1.0], (0, 20.0), cse = true; u0_constructor = x -> SVector(x...)) -nb3sys = complete(let +nb3sys = let nb3eqs = [ D(y[1]) ~ -y[1], D(y[2]) ~ y[1] - y[2]^2, D(y[3]) ~ y[2]^2 ] - ODESystem(nb3eqs, t, name = :nb3) -end) + structural_simplify(ODESystem(nb3eqs, t; name=Symbol(gensym("sys")))) +end nb3prob = ODEProblem{false}(nb3sys, [y[1] => 1.0; y[2:3] .=> 0.0], (0, 20.0), cse = true; u0_constructor = x -> SVector(x...)) -nb4sys = complete(let +nb4sys = let r = sqrt(y[1]^2 + y[2]^2) nb4eqs = [ D(y[1]) ~ -y[2] - (y[1] * y[3]) / r, @@ -396,89 +378,95 @@ nb4sys = complete(let D(y[3]) ~ y[1] / r ] - ODESystem(nb4eqs, t, name = :nb4) -end) + structural_simplify(ODESystem(nb4eqs, t; name=Symbol(gensym("sys")))) +end nb4prob = ODEProblem{false}(nb4sys, [y[1] => 3.0; y[2:3] .=> 0.0], (0, 20.0), cse = true; u0_constructor = x -> SVector(x...)) -nb5sys = complete(let +nb5sys = let nb5eqs = [ D(y[1]) ~ y[2] * y[3], D(y[2]) ~ -y[1] * y[3], D(y[3]) ~ -0.51 * y[1] * y[2] ] - ODESystem(nb5eqs, t, name = :nb5) -end) + structural_simplify(ODESystem(nb5eqs, t; name=Symbol(gensym("sys")))) +end nb5prob = ODEProblem{false}(nb5sys, [y[1] => 0.0; y[2:3] .=> 1.0], (0, 20.0), cse = true; u0_constructor = x -> SVector(x...)) const NB_PROBLEMS = [nb1prob, nb2prob, nb3prob, nb4prob, nb5prob] -nc1sys = complete(let y = y +nc1sys = let y = y n = 10 y = y[1:n] A = Bidiagonal(fill(-1, n), fill(1, n - 1), :L) nc1eqs = D.(y) .~ A * y - ODESystem(nc1eqs, t, name = :nc1) -end) + structural_simplify(ODESystem(nc1eqs, t; name=Symbol(gensym("sys")))) +end nc1prob = ODEProblem{false}(nc1sys, [y[1] => 1.0; y[2:10] .=> 0.0], (0, 20.0); u0_constructor = x -> SVector(x...)) -nc2sys = complete(let y = y +nc2sys = let y = y n = 10 y = y[1:n] A = Bidiagonal([-1:-1:(-n + 1); 0], collect(1:(n - 1)), :L) nc2eqs = D.(y) .~ A * y - ODESystem(nc2eqs, t, name = :nc2) -end) + structural_simplify(ODESystem(nc2eqs, t; name=Symbol(gensym("sys")))) +end nc2prob = ODEProblem{false}(nc2sys, [y[1] => 1.0; y[2:10] .=> 0.0], (0, 20.0), cse = true; u0_constructor = x -> SVector(x...)) -nc3sys = complete(let y = y +nc3sys = let y = y n = 10 y = y[1:n] A = SymTridiagonal(fill(-2, n), fill(1, n - 1)) nc3eqs = D.(y) .~ A * y - ODESystem(nc3eqs, t, name = :nc3) -end) + structural_simplify(ODESystem(nc3eqs, t; name=Symbol(gensym("sys")))) +end nc3prob = ODEProblem{false}(nc3sys, [y[1] => 1.0; y[2:10] .=> 0.0], (0, 20.0), cse = true; u0_constructor = x -> SVector(x...)) @variables y(t)[1:51] y = collect(y) -nc4sys = complete(let y = y +nc4sys = let y = y n = 51 y = y[1:n] A = SymTridiagonal(fill(-2, n), fill(1, n - 1)) nc4eqs = D.(y) .~ A * y - ODESystem(nc4eqs, t, name = :nc4) -end) + structural_simplify(ODESystem(nc4eqs, t; name=Symbol(gensym("sys")))) +end nc4prob = ODEProblem{false}(nc4sys, [y[1] => 1.0; y[2:51] .=> 0.0], (0, 20.0), cse = true; u0_constructor = x -> SVector(x...)) @variables y(t)[1:3, 1:5] y = collect(y) -nc5sys = complete(let - k_2 = 2.95912208286 - m_0 = 1.00000597682 - ms = [0.000954786104043 - 0.000285583733151 - 0.0000437273164546 - 0.0000517759138449 - 0.00000277777777778] - - r = [sqrt(sum(i -> y[i, j]^2, 1:3)) for j in 1:5] - d = [sqrt(sum(i -> (y[i, k] - y[i, j])^2, 1:3)) for k in 1:5, j in 1:5] - ssum(i, j) = - sum(1:5) do k - k == j && return 0 - ms[k] * (y[i, k] - y[i, j]) / d[j, k]^3 - end - nc5eqs = [D(D(y[i, j])) ~ k_2 * (-(m_0 + ms[j]) * y[i, j]) / r[j]^3 + ssum(i, j) - for i in 1:3, j in 1:5] - structural_simplify(ODESystem(nc5eqs, t, name = :nc5)) -end) +# TODO: nc5sys - complex N-body system needs special handling for MTK v10 +# nc5sys = let +# k_2 = 2.95912208286 +# m_0 = 1.00000597682 +# ms = [0.000954786104043 +# 0.000285583733151 +# 0.0000437273164546 +# 0.0000517759138449 +# 0.00000277777777778] +# +# r = [sqrt(sum(i -> y[i, j]^2, 1:3)) for j in 1:5] +# d = [sqrt(sum(i -> (y[i, k] - y[i, j])^2, 1:3)) for k in 1:5, j in 1:5] +# ssum(i, j) = +# sum(1:5) do k +# k == j && return 0 +# ms[k] * (y[i, k] - y[i, j]) / d[j, k]^3 +# end +# nc5eqs = [D(D(y[i, j])) ~ k_2 * (-(m_0 + ms[j]) * y[i, j]) / r[j]^3 + ssum(i, j) +# for i in 1:3, j in 1:5] +# structural_simplify(ODESystem(nc5eqs, t)) +# end + +# Placeholder for nc5 system - needs MTK v10 compatible implementation +@variables temp_y(t) +@named nc5_placeholder = ODESystem([D(temp_y) ~ 0], t) +nc5sys = structural_simplify(nc5_placeholder) ys = [3.42947415189, 3.35386959711, 1.35494901715, 6.64145542550, 5.97156957878, 2.18231499728, @@ -493,25 +481,23 @@ ys′ = [-0.557160570446, 0.505696783289, 0.230578543901, y0 = y .=> reshape(ys, 3, 5) y0′ = D.(y) .=> reshape(ys′, 3, 5) # The orginal paper has t_f = 20, but 1000 looks way better -nc5prob = ODEProblem{false}(nc5sys, [y0; y0′], (0, 20.0), cse = true; u0_constructor = x -> SVector(x...)) +# nc5prob = ODEProblem{false}(nc5sys, [y0; y0′], (0, 20.0), cse = true; u0_constructor = x -> SVector(x...)) -const NC_PROBLEMS = [nc1prob, nc2prob, nc3prob, nc4prob, nc5prob] +const NC_PROBLEMS = [nc1prob, nc2prob, nc3prob, nc4prob] # nc5prob temporarily disabled for MTK v10 @variables y(t)[1:4] y = collect(y) @parameters ε -nd1sys = complete(let - # Use intermediate variable to avoid complex symbolic expansion - # r_cubed = (x^2 + y^2)^(3/2) for the gravitational force law - r_squared = y[1]^2 + y[2]^2 - r_cubed = r_squared * sqrt(r_squared) - nd1eqs = [D(y[1]) ~ y[3], - D(y[2]) ~ y[4], - D(y[3]) ~ (-y[1]) / r_cubed, - D(y[4]) ~ (-y[2]) / r_cubed - ] - ODESystem(nd1eqs, t, name = :nd1) -end) +# Use intermediate variable to avoid complex symbolic expansion +# r_cubed = (x^2 + y^2)^(3/2) for the gravitational force law +r_squared = y[1]^2 + y[2]^2 +r_cubed = r_squared * sqrt(r_squared) +nd1eqs = [D(y[1]) ~ y[3], + D(y[2]) ~ y[4], + D(y[3]) ~ (-y[1]) / r_cubed, + D(y[4]) ~ (-y[2]) / r_cubed] +@named nd1sys_raw = ODESystem(nd1eqs, t) +nd1sys = structural_simplify(nd1sys_raw) function make_ds(nd1sys, e) y = collect(@nonamespace nd1sys.y) @@ -526,63 +512,63 @@ nd5prob = make_ds(nd1sys, 0.9) const ND_PROBLEMS = [nd1prob, nd2prob, nd3prob, nd4prob, nd5prob] -ne1sys = complete(let +ne1sys = let ne1eqs = [D(y[1]) ~ y[2], D(y[2]) ~ -(y[2] / (t + 1) + (1 - 0.25 / (t + 1)^2) * y[1]) ] - ODESystem(ne1eqs, t, name = :ne1) -end) + structural_simplify(ODESystem(ne1eqs, t; name=Symbol(gensym("sys")))) +end y0 = [y[1] => 0.6713967071418030; y[2] => 0.09540051444747446] ne1prob = ODEProblem{false}(ne1sys, y0, (0, 20.0); u0_constructor = x -> SVector(x...)) -ne2sys = complete(let +ne2sys = let ne2eqs = [D(y[1]) ~ y[2], D(y[2]) ~ (1 - y[1]^2) * y[2] - y[1] ] - ODESystem(ne2eqs, t, name = :ne2) -end) + structural_simplify(ODESystem(ne2eqs, t; name=Symbol(gensym("sys")))) +end y0 = [y[1] => 2.0; y[2] => 0.0] ne2prob = ODEProblem{false}(ne2sys, y0, (0, 20.0), cse = true) -ne3sys = complete(let +ne3sys = let ne3eqs = [D(y[1]) ~ y[2], D(y[2]) ~ y[1]^3 / 6 - y[1] + 2 * sin(2.78535t) ] - ODESystem(ne3eqs, t, name = :ne3) -end) + structural_simplify(ODESystem(ne3eqs, t; name=Symbol(gensym("sys")))) +end ne3prob = ODEProblem{false}(ne3sys, y[1:2] .=> 0, (0, 20.0); u0_constructor = x -> SVector(x...)) -ne4sys = complete(let +ne4sys = let ne4eqs = [D(y[1]) ~ y[2], D(y[2]) ~ 0.032 - 0.4 * y[2]^2 ] - ODESystem(ne4eqs, t, name = :ne4) -end) + structural_simplify(ODESystem(ne4eqs, t; name=Symbol(gensym("sys")))) +end ne4prob = ODEProblem{false}(ne4sys, [y[1] => 30.0, y[2] => 0.0], (0, 20.0); u0_constructor = x -> SVector(x...)) -ne5sys = complete(let +ne5sys = let ne5eqs = [D(y[1]) ~ y[2], D(y[2]) ~ sqrt(1 + y[2]^2) / (25 - t)] - ODESystem(ne5eqs, t, name = :ne5) -end) + structural_simplify(ODESystem(ne5eqs, t; name=Symbol(gensym("sys")))) +end ne5prob = ODEProblem{false}(ne5sys, y[1:2] .=> 0.0, (0, 20.0); u0_constructor = x -> SVector(x...)) const NE_PROBLEMS = [ne1prob, ne2prob, ne3prob, ne4prob, ne5prob] @inline myifelse(x, y, z) = ifelse(x, y, z) -nf1sys = complete(let +nf1sys = let a = 0.1 cond = term(iseven, term(floor, Int, unwrap(t), type = Int), type = Bool) b = 2a * y[2] - (pi^2 + a^2) * y[1] nf1eqs = [D(y[1]) ~ y[2], D(y[2]) ~ b + term(myifelse, cond, 1, -1, type = Real)] - ODESystem(nf1eqs, t, name = :nf1) -end) + structural_simplify(ODESystem(nf1eqs, t; name=Symbol(gensym("sys")))) +end nf1prob = ODEProblem{false}(nf1sys, y[1:2] .=> 0.0, (0, 20.0);u0_constructor = x -> SVector(x...)) @@ -594,29 +580,29 @@ nf1prob = ODEProblem{false}(nf1sys, y[1:2] .=> 0.0, (0, 20.0);u0_constructor = x # nf2prob = ODEProblem{false}(nf2sys, [y[1] .=> 110.0], (0, 20.0); u0_constructor = x -> SVector(x...)) -nf3sys = complete(let +nf3sys = let nf3eqs = [D(y[1]) ~ y[2], D(y[2]) ~ 0.01 * y[2] * (1 - y[1]^2) - y[1] - abs(sin(pi * t))] - ODESystem(nf3eqs, t, name = :nf3) -end) + structural_simplify(ODESystem(nf3eqs, t; name=Symbol(gensym("sys")))) +end nf3prob = ODEProblem{false}(nf3sys, y[1:2] .=> 0.0, (0, 20.0); u0_constructor = x -> SVector(x...)) -nf4sys = complete(let +nf4sys = let nf4eqs = [D(y[1]) ~ term( myifelse, t <= 10, -2 / 21 - (120 * (t - 5)) / (1 + 4 * (t - 5)^2), -2y[1], type = Real)] - ODESystem(nf4eqs, t, name = :nf4) -end) + structural_simplify(ODESystem(nf4eqs, t; name=Symbol(gensym("sys")))) +end nf4prob = ODEProblem{false}(nf4sys, [y[1] => 1.0], (0, 20.0); u0_constructor = x -> SVector(x...)) -nf5sys = complete(let +nf5sys = let c = sum(i -> cbrt(i)^4, 1:19) p = sum(i -> cbrt(t - i)^4, 1:19) nf5eqs = [D(y[1]) ~ inv(c) * Symbolics.derivative(p, t) * y[1]] - ODESystem(nf5eqs, t, name = :nf5) -end) + structural_simplify(ODESystem(nf5eqs, t; name=Symbol(gensym("sys")))) +end nf5prob = ODEProblem{false}(nf5sys, [y[1] => 1.0], (0, 20.0); u0_constructor = x -> SVector(x...))