From 126dc7998314378f66adc663b8cc68b7ca1bd046 Mon Sep 17 00:00:00 2001 From: Yingbo Ma Date: Fri, 19 Aug 2022 13:12:37 -0400 Subject: [PATCH 1/5] Add up to NC5 --- src/ode/enright_pryce.jl | 469 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 469 insertions(+) create mode 100644 src/ode/enright_pryce.jl diff --git a/src/ode/enright_pryce.jl b/src/ode/enright_pryce.jl new file mode 100644 index 0000000..b198f89 --- /dev/null +++ b/src/ode/enright_pryce.jl @@ -0,0 +1,469 @@ +using ModelingToolkit +using DiffEqBase, StaticArrays, LinearAlgebra + +@variables t y(t)[1:10] +y = collect(y) +D = Differential(t) + +# Stiff +sa1sys = 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 + +sa1prob = ODEProblem{false}(sa1sys, y[1:4] .=> 1.0, (0, 20.0), dt = 1e-2) + +sa2sys = 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 + +sa2prob = ODEProblem{false}(sa2sys, y[1:9] .=> 0.0, (0, 120.0), dt = 5e-4) + +sa3sys = 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 + +sa3prob = ODEProblem{false}(sa3sys, y[1:4] .=> 1.0, (0, 20.0), dt = 1e-5) + +sa4sys = let + sa4eqs = [D(y[i]) ~ -i^5 * y[i] for i in 1:10] + ODESystem(sa4eqs, t, name = :sa4) +end + +sa4prob = ODEProblem{false}(sa4sys, y[1:10] .=> 1.0, (0, 1.0), dt = 1e-5) + +#const SA_PROBLEMS = [sa1prob, sa2prob, sa3prob, sa4prob] + +sb1sys = 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 + +sb1prob = ODEProblem{false}(sb1sys, [y[[1, 3]] .=> 1.0; y[[2, 4]] .=> 0.0;], (0, 20.0), dt = 7e-3) + + +@parameters α +sb2sys = 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 + +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 = 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 + +sc1prob = ODEProblem{false}(sc1sys, y .=> 1.0, (0, 20.0), dt = 1e-2) + + +@parameters β +sc2sys = 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 + +sc2prob = ODEProblem{false}(sc2sys, y .=> 1.0, (0, 20.0), [β => 0.1], dt = 1e-2, cse = true) +sc3prob = ODEProblem{false}(sc2sys, y .=> 1.0, (0, 20.0), [β => 1.0], dt = 1e-2, cse = true) +sc3prob = ODEProblem{false}(sc2sys, y .=> 1.0, (0, 20.0), [β => 10.0], dt = 1e-2, cse = true) +sc3prob = ODEProblem{false}(sc2sys, y .=> 1.0, (0, 20.0), [β => 20.0], dt = 1e-2, cse = true) + +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 + +sd1prob = ODEProblem{false}(sd1sys, y .=> 0.0, (0, 400.0), [β => 0.1], dt = 1.7e-2) + +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 + +sd2prob = ODEProblem{false}(sd2sys, [y[1] => 1.0; y[2:3] .=> 0.0], (0, 40.0), dt = 1e-5, cse = true) + +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 + +sd3prob = ODEProblem{false}(sd3sys, [y[1:2] .=> 1; y[3:4] .=> 0.0], (0, 20.0), dt = 2.5e-5, cse = true) + +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 + +sd4prob = ODEProblem{false}(sd4sys, [y[1:2] .=> 1; y[3] => 0.0], (0, 50.0), dt = 2.9e-4, cse = true) + +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 + +sd5prob = ODEProblem{false}(sd5sys, y[1:2] .=> 0.0, (0, 100.0), dt = 1e-4, cse = true) + +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 + +sd6prob = ODEProblem{false}(sd6sys, [y[1] => 1.0; y[2:3] .=> 0.0], (0, 1.0), dt = 3.3e-8, cse = true) + +se1sys = let + Γ = 100 + se1eqs = [D(y[1]) ~ y[2], + D(y[2]) ~ y[3], + D(y[3]) ~ y[4], + D(y[4]) ~ (y[1]^2 - sin(y[1]) - Γ^4) * y[1] + (y[2] * y[3] / (y[1]^2 + 1) - 4 * Γ^3) * y[2] + (1 - 6 * Γ^2) * y[3] + (10 * exp(-y[4]^2) - 4Γ) * y[4] + 1 + ] + + ODESystem(se1eqs, t, name = :se1) +end + +se1prob = ODEProblem{false}(se1sys, y .=> 0.0, (0, 1.0), dt = 6.8e-3, cse = true) + +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 + +se2prob = ODEProblem{false}(se2sys, [y[1] => 2.0, y[2] => 0.0], (0, 1.0), dt = 1e-3, cse = true) + +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 + +se3prob = ODEProblem{false}(se3sys, [y[1:2] .=> 1.0; y[3] => 0.0], (0, 500.0), dt = 0.02, cse = true) + +se4sys = let y = y + y = y[1:4] + U = ones(4, 4) + U[diagind(U)] .= -1 + U ./= 2 + Z = U * y + G = U' * [(Z[1]^2 - Z[2]^2) / 2, Z[1] * Z[2], Z[3]^2, Z[4]^2] + 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 + +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) + +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], + D(y[3]) ~ 7.89e-10 * y[1] - 1.1e7 * y[1] * y[3] + 1.13e3 * y[4] - 1.13e9 * y[2] * y[3], + D(y[4]) ~ 1.1e7 * y[1] * y[3] - 1.13e3 * y[4], + ] + + ODESystem(se5eqs, t, name = :se5) +end + +se5prob = ODEProblem{false}(se5sys, [y[1] => 1.76e-3; y[2:4] .=> 0.0], (0, 1000.0), dt = 1e-3, cse = true) + +sf1sys = let + k = exp(20.7 - 1500 / y[1]) + sf1eqs = [ + D(y[1]) ~ 1.3 * (y[3] - y[1]) + 10400 * k * y[2], + D(y[2]) ~ 1880 * [y[4] - y[2] * (1 + k)], + D(y[3]) ~ 1752 - 269 * y[3] + 267 * y[1], + D(y[4]) ~ 0.1 + 320 * y[2] - 321 * y[4], + ] + + ODESystem(sf1eqs, t, name = :sf1) +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 = 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 + +sf2prob = ODEProblem{false}(sf2sys, [y[1] => 1.0; y[2] => 0.0], (0, 240.0), dt = 1e-2, cse = true) + +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], + D(y[3]) ~ 1e7 * y[2] * y[1] - 1.001e4 * y[3] + 1e-3 * y[4], + D(y[4]) ~ 1e4 * y[3] - 10.001 * y[4] + 1e7 * y[2] * y[5], + D(y[5]) ~ 10 * y[4] - 1e7 * y[2] * y[5], + ] + + ODESystem(sf3eqs, t, name = :sf3) +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) + +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 + +sf4prob = ODEProblem{false}(sf4sys, [y[1] => 4.0; y[2] => 1.1; y[3] => 4.0], (0, 300.0), dt = 1e-3, cse = true) + +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], + D(y[3]) ~ 1e11 * (-9 * y[1] * y[3] + 0.001 * y[4]), + D(y[4]) ~ 1e11 * (3 * y[1] * y[2] - 0.0012 * y[4] + 9 * y[1] * y[3]), + ] + + ODESystem(sf5eqs, t, name = :sf5) +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) + +# Non-stiff +na1sys = let y = y[1] + na1eqs = D(y) ~ -y + + ODESystem(na1eqs, t, name = :na1) +end + +na1prob = ODEProblem{false}(na1sys, [y[1] => 1], (0, 20.0), cse = true) + +na2sys = let y = y[1] + na2eqs = D(y) ~ -y^3/2 + + ODESystem(na2eqs, t, name = :na2) +end + +na2prob = ODEProblem{false}(na2sys, [y[1] => 1], (0, 20.0), cse = true) + +na3sys = let y = y[1] + na3eqs = D(y) ~ y * cos(t) + + ODESystem(na3eqs, t, name = :na3) +end + +na3prob = ODEProblem{false}(na3sys, [y[1] => 1], (0, 20.0), cse = true) + +na4sys = let y = y[1] + na4eqs = D(y) ~ y/4 * (1 - y/20) + + ODESystem(na4eqs, t, name = :na4) +end + +na4prob = ODEProblem{false}(na4sys, [y[1] => 1], (0, 20.0), cse = true) + +na5sys = let y = y[1] + na5eqs = D(y) ~ (y - t) / (y + t) + + ODESystem(na5eqs, t, name = :na5) +end + +na5prob = ODEProblem{false}(na5sys, [y[1] => 4], (0, 20.0), cse = true) + +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 + +nb1prob = ODEProblem{false}(nb1sys, [y[1] => 1.0, y[2] => 3], (0, 20.0), cse = true) + +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 + +nb2prob = ODEProblem{false}(nb2sys, [y[1] => 2.0, y[2] => 0.0, y[3] => 1.0], (0, 20.0), cse = true) + +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 + +nb3prob = ODEProblem{false}(nb3sys, [y[1] => 1.0; y[2:3] .=> 0.0], (0, 20.0), cse = true) + +nb4sys = let + r = sqrt(y[1]^2 + y[2]^2) + nb4eqs = [ + D(y[1]) ~ -y[2] - (y[1] * y[3]) / r, + D(y[2]) ~ y[1] - (y[2] * y[3]) / r, + D(y[3]) ~ y[1] / r, + ] + + ODESystem(nb4eqs, t, name = :nb4) +end + +nb4prob = ODEProblem{false}(nb4sys, [y[1] => 3.0; y[2:3] .=> 0.0], (0, 20.0), cse = true) + +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 + +nb5prob = ODEProblem{false}(nb5sys, [y[1] => 0.0; y[2:3] .=> 1.0], (0, 20.0), cse = true) + +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 + +nc1prob = ODEProblem{false}(nc1sys, [y[1] => 1.0; y[2:10] .=> 0.0], (0, 20.0), cse = true) + +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 + +nc2prob = ODEProblem{false}(nc2sys, [y[1] => 1.0; y[2:10] .=> 0.0], (0, 20.0), cse = true) + +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 + +nc3prob = ODEProblem{false}(nc3sys, [y[1] => 1.0; y[2:10] .=> 0.0], (0, 20.0), cse = true) + +@variables y(t)[1:51] +y = collect(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 + +nc4prob = ODEProblem{false}(nc4sys, [y[1] => 1.0; y[2:51] .=> 0.0], (0, 20.0), cse = true) + +@variables y(t)[1:3, 1:5] +y = collect(y) +nc5sys = let + k_2 = 2.95912208286 + m_0 = 1.00000597682 + ms = [0.000954786104043 + 0.000285583733151 + 0.0000437273164546 + 0.0000517759138449 + 0.00000277777777778] + + r = [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 + +ys = [3.42947415189, 3.35386959711, 1.35494901715, + 6.64145542550, 5.97156957878, 2.18231499728, + 11.2630437207, 14.6952576794, 6.27960525067, + -30.1552268759, 1.65699966404, 1.43785752721, + -21.1238353380, 28.4465098142, 15.3882659679,] +ys′ = [-0.557160570446, 0.505696783289, 0.230578543901, + -0.415570776342, 0.365682722812, 0.169143213293, + -0.325325669158, 0.189706021964, 0.087726532278, + -0.0240476254170, -0.287659532608, -0.117219543175, + -0.176860753121, -0.216393453025, -0.0148647893090,] +y0 = y .=> reshape(ys, 3, 5) +y0′ = D.(y) .=> reshape(ys′, 3, 5) +nc5prob = ODEProblem{false}(nc5sys, [y0; y0′], (0, 20.0), cse = true) From db52c70c1fbfc7bdaef0ca66648023df708a1ab6 Mon Sep 17 00:00:00 2001 From: Yingbo Ma Date: Mon, 22 Aug 2022 15:33:15 -0400 Subject: [PATCH 2/5] Add the rest --- src/ode/enright_pryce.jl | 121 ++++++++++++++++++++++++++++++++++++++- 1 file changed, 119 insertions(+), 2 deletions(-) diff --git a/src/ode/enright_pryce.jl b/src/ode/enright_pryce.jl index b198f89..5e08604 100644 --- a/src/ode/enright_pryce.jl +++ b/src/ode/enright_pryce.jl @@ -1,4 +1,7 @@ using ModelingToolkit +using IfElse +using Symbolics +using Symbolics: unwrap using DiffEqBase, StaticArrays, LinearAlgebra @variables t y(t)[1:10] @@ -444,13 +447,13 @@ nc5sys = let 0.0000517759138449 0.00000277777777778] - r = [sum(i->y[i, j]^2, 1:3) for j in 1:5] + 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] + 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 @@ -466,4 +469,118 @@ ys′ = [-0.557160570446, 0.505696783289, 0.230578543901, -0.176860753121, -0.216393453025, -0.0148647893090,] 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) + +@variables y(t)[1:4] +y = collect(y) +@parameters ε +nd1sys = let + r = sqrt(y[1]^2 + y[2]^2)^3 + nd1eqs = [D(y[1]) ~ y[3], + D(y[2]) ~ y[4], + D(y[3]) ~ (-y[1]) / r, + D(y[4]) ~ (-y[2]) / r, + ] + ODESystem(nd1eqs, t, name = :nd1) +end + +function make_ds(nd1sys, e) + y = collect(@nonamespace nd1sys.y) + y0 = [y[1] => 1-e; y[2:3] .=> 0.0; y[4] => sqrt((1 + e) / (1 - e))] + ODEProblem{false}(nd1sys, y0, (0, 20.0), [ε => e], cse = true) +end +nd1prob = make_ds(nd1sys, 0.1) +nd2prob = make_ds(nd1sys, 0.3) +nd3prob = make_ds(nd1sys, 0.5) +nd4prob = make_ds(nd1sys, 0.7) +nd5prob = make_ds(nd1sys, 0.9) + +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 + +y0 = [y[1] => 0.6713967071418030; y[2] => 0.09540051444747446] +ne1prob = ODEProblem{false}(ne1sys, y0, (0, 20.0), cse = true) + +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 + +y0 = [y[1] => 2.0; y[2] => 0.0] +ne2prob = ODEProblem{false}(ne2sys, y0, (0, 20.0), cse = true) + +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 + +ne3prob = ODEProblem{false}(ne3sys, y[1:2] .=> 0, (0, 20.0), cse = true) + +ne4sys = let + ne4eqs = [D(y[1]) ~ y[2], + D(y[2]) ~ 0.032 - 0.4 * y[2]^2, + ] + ODESystem(ne4eqs, t, name = :ne4) +end + +ne4prob = ODEProblem{false}(ne4sys, [y[1] => 30.0, y[2] => 0.0], (0, 20.0), cse = true) + +ne5sys = let + ne5eqs = [D(y[1]) ~ y[2], + D(y[2]) ~ sqrt(1 - y[2]^2) / (25 - t),] + ODESystem(ne5eqs, t, name = :ne5) +end + +ne5prob = ODEProblem{false}(ne5sys, y[1:2] .=> 0.0, (0, 20.0), cse = true) + +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 + IfElse.ifelse(cond, 1, -1)] + ODESystem(nf1eqs, t, name = :nf1) +end + +nf1prob = ODEProblem{false}(nf1sys, y[1:2] .=> 0.0, (0, 20.0)) + +nf2sys = let + cond = term(iseven, term(floor, Int, unwrap(t), type = Int), type = Bool) + nf2eqs = [D(y[1]) ~ 55 - IfElse.ifelse(cond, 2y[1]/2, y[1]/2)] + ODESystem(nf2eqs, t, name = :nf2) +end + +nf2prob = ODEProblem{false}(nf2sys, [y[1] .=> 110.0], (0, 20.0), cse = true) + +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 + +nf3prob = ODEProblem{false}(nf3sys, y[1:2] .=> 0.0, (0, 20.0), cse = true) + +nf4sys = let + nf4eqs = [D(y[1]) ~ IfElse.ifelse(t <= 10, -2/21 - (120 * (t - 5)) / (1 + 4 * (t - 5)^2), -2y[1])] + ODESystem(nf4eqs, t, name = :nf4) +end + +nf4prob = ODEProblem{false}(nf4sys, [y[1] => 1.0], (0, 20.0), cse = true) + +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 + +nf5prob = ODEProblem{false}(nf5sys, [y[1] => 1.0], (0, 20.0), cse = true) From a15c1945226f4152c622f6c66c83300b8b45b270 Mon Sep 17 00:00:00 2001 From: Yingbo Ma Date: Mon, 22 Aug 2022 21:56:33 -0400 Subject: [PATCH 3/5] Don't use cse when it's broken Co-authored-by: Chris Elrod --- m.jl | 387 +++++++++++++++++++++++++++++++++++++++ src/ode/enright_pryce.jl | 39 ++-- 2 files changed, 407 insertions(+), 19 deletions(-) create mode 100644 m.jl diff --git a/m.jl b/m.jl new file mode 100644 index 0000000..490d0b7 --- /dev/null +++ b/m.jl @@ -0,0 +1,387 @@ +using ModelingToolkit +using LinearAlgebra +@parameters( + m₁ = 0.36, + m₂ = 0.151104, + m₃ = 0.075552, + l₁ = 0.15, + l₂ = 0.30, + J₁ = 0.002727, + J₂ = 0.0045339259, + EE = 0.20e12, + NUE= 0.30, + BB = 0.0080, + HH = 0.0080, + ρ = 7870.0, + γ = 0.0, + Ω = 150., + ) + +FACM = ρ*BB*HH*l₂ +FACK = EE*BB*HH/l₂ +FACB = BB*HH*l₂ + +MQ = zeros(Num, 4, 4) +MQ[1,1] = FACM*.5 +MQ[2,2] = FACM*.5 +MQ[3,3] = FACM*8.0 +MQ[3,4] = FACM*1.0 +MQ[4,3] = FACM*1.0 +MQ[4,4] = FACM*2.0 + +KQ = zeros(Num, 4, 4) +KQ[1,1] = FACK*π^4/24.0*(HH/l₂)^2 +KQ[2,2] = FACK*π^4*2.0/3.0*(HH/l₂)^2 +KQ[3,3] = FACK*16.0/3.0 +KQ[3,4] = -FACK*8.0/3.0 +KQ[4,3] = -FACK*8.0/3.0 +KQ[4,4] = FACK*7.0/3.0 + +BQ = zeros(Num, 4, 4) +BQ[1,3] = -FACB*16.0/π^3 +BQ[1,4] = FACB*(8.0/π^3-1.0/π) +BQ[2,4] = FACB*0.5/π +BQ[3,1] = FACB*16.0/π^3 +BQ[4,1] = -FACB*(8.0/π^3-1.0/π) +BQ[4,2] = -FACB*0.5/π + +c1 = zeros(Num, 4) +c2 = zeros(Num, 4) +c12 = zeros(Num, 4) +c21 = zeros(Num, 4) +c1[3] = FACB*2.0/3.0 +c1[4] = FACB*1.0/6.0 +c2[1] = FACB*2.0/π +c12[3] = l₂*FACB*1.0/3.0 +c12[4] = l₂*FACB*1.0/6.0 +c21[1] = l₂*FACB*1.0/π +c21[2] = -l₂*FACB*0.5/π + +DQ = zeros(Num, 4, 4) +# OPTIONAL DAMPING +# DQ[1,1] = 5.0 +# DQ[2,2] = 25.0 +# DQ[3,3] = 0.5*2.308375455264791e2 +# DQ[3,4] = -0.5*2.62688487992052e2 +# DQ[4,3] = -0.5*2.626884879920526e2 +# DQ[4,4] = 0.5*4.217421837156818e2 + +## Setup variables and initial conditions: + +# init +# Initial values: 'Close' to smooth motion, +# accelerations and multipliers consistent +# for linear stiffness term and no damping +# (ipar(1) = 0, ipar(2) = 0). +@variables( + t, + ϕ₁(t)=0.0, + ϕ₂(t)=0.0, + x₃(t)=4.500169330000000e-1, + q₁(t)=0.0, + q₂(t)=0.0, + q₃(t)=1.033398630000000e-5, + q₄(t)=1.693279690000000e-5, + dϕ₁(t)= 1.500000000000000e2, + dϕ₂(t)=-7.499576703969453e1, + dx₃(t)=-2.689386719979040e-6, + dq₁(t)= 4.448961125815990e-1, + dq₂(t)= 4.634339319238670e-3, + dq₃(t)=-1.785910760000550e-6, + dq₄(t)=-2.689386719979040e-6, + ddϕ₁(t)= 0.000000000000000e0, + ddϕ₂(t)=-1.344541576709835e-3, + ddx₃(t)=-5.062194924490193e+3, + ddq₁(t)=-6.829725665986310e-5, + ddq₂(t)= 1.813207639590617e-20, + ddq₃(t)=-4.268463266810281e+0, + ddq₄(t)= 2.098339029337557e-1, + λ₁(t) = -6.552727150584648e-8, + λ₂(t) = 3.824589509350831e+2, + λ₃(t) = -4.635908708561371e-9, + ) + + +# init1 +# Initial values: 'Close' to smooth motion, +# accelerations and multipliers consistent +# for linear stiffness term and no damping +# (ipar(1) = 0, ipar(2) = 0). +#@variables( +# t, +# ϕ₁(t) = 0.0, +# ϕ₂(t) = 0.0, +# x₃(t) = .450016933, +# q₁(t) = 0.0, +# q₂(t) = 0.0, +# q₃(t) = .103339863e-4, +# q₄(t) = .169327969e-4, +# dϕ₁(t) = .150000000e+3, +# dϕ₂(t) = -.749957670e+2, +# dx₃(t) = -.268938672e-5, +# dq₁(t) = .444896105e0, +# dq₂(t) = .463434311e-2, +# dq₃(t) = -.178591076e-5, +# dq₄(t) = -.268938672e-5, +# ddϕ₁(t) = 0.0, +# ddϕ₂(t) = -1.344541576008661e-3, +# ddx₃(t) = -5.062194923138079e+3, +# ddq₁(t) = -6.833142732779555e-5, +# ddq₂(t) = 1.449382650173157e-8, +# ddq₃(t) = -4.268463211410861e+0, +# ddq₄(t) = 2.098334687947376e-1, +# λ₁(t) = -6.397251492537153e-8, +# λ₂(t) = 3.824589508329281e+2, +# λ₃(t) = -4.376060460948886e-9, +# ) + +# init2 +# Initial values: On rigid motion trajectory, +# leading to strong oscillations. +# accelerations and multipliers consistent +# for linear stiffness term and no damping +# (ipar(1) = 0, ipar(2) = 0). +#@variables( +# t, +# ϕ₁(t) = 0.0, +# ϕ₂(t) = 0.0, +# x₃(t) = .4500e0, +# q₁(t) = 0.0, +# q₂(t) = 0.0, +# q₃(t) = 0.0, +# q₄(t) = 0.0, +# dϕ₁(t) = .150e+3, +# dϕ₂(t) = -.750e+2, +# dx₃(t) = 0.0, +# dq₁(t) = 0.0, +# dq₂(t) = 0.0, +# dq₃(t) = 0.0, +# dq₄(t) = 0.0, +# ddϕ₁(t) = 0.0, +# ddϕ₂(t) = 0.0, +# ddx₃(t) = -3.789473684210526e+3, +# ddq₁(t) = 0.0, +# ddq₂(t) = 0.0, +# ddq₃(t) = 1.924342105263158e+2, +# ddq₄(t) = 1.273026315789474e+3, +# λ₁(t) = 0.0, +# λ₂(t) = 2.863023157894737e+2, +# λ₃(t) = 0.0, +# ) + + +# Define problem equations: + +cosp1 = cos(ϕ₁) +cosp2 = cos(ϕ₂) +sinp1 = sin(ϕ₁) +sinp2 = sin(ϕ₂) +cosp12 = cos(ϕ₁-ϕ₂) +sinp12 = sin(ϕ₁-ϕ₂) + +D = Differential(t) +P = [ϕ₁, ϕ₂, x₃] +Q = [q₁, q₂, q₃, q₄] +X = [P..., Q...] +# V = D.(P) +V = [dϕ₁, dϕ₂, dx₃] +# QD = D.(Q) +QD = [dq₁, dq₂, dq₃, dq₄] +DX = [V..., QD...] +# DDX = D.(DX) +DDX = [ddϕ₁, ddϕ₂, ddx₃, ddq₁, ddq₂, ddq₃, ddq₄] + +ivs = [ + # Initial values velocity variables + D(ϕ₁) => Symbolics.getdefaultval(dϕ₁), + D(ϕ₂) => Symbolics.getdefaultval(dϕ₂), + D(x₃) => Symbolics.getdefaultval(dx₃), + D(q₁) => Symbolics.getdefaultval(q₁), + D(q₂) => Symbolics.getdefaultval(q₂), + D(q₃) => Symbolics.getdefaultval(q₃), + D(q₄) => Symbolics.getdefaultval(q₄), + # Initial values acceleration variables + D(DX[1]) => Symbolics.getdefaultval(ddϕ₁), + D(DX[2]) => Symbolics.getdefaultval(ddϕ₂), + D(DX[3]) => Symbolics.getdefaultval(ddx₃), + D(DX[4]) => Symbolics.getdefaultval(dq₁), + D(DX[5]) => Symbolics.getdefaultval(dq₂), + D(DX[6]) => Symbolics.getdefaultval(dq₃), + D(DX[7]) => Symbolics.getdefaultval(dq₄), + D(λ₁) => 0.0, + D(λ₂) => 0.0, + D(λ₃) => 0.0, + ] +for I = 1:7 + push!(ivs, D(DDX[I]) => 0.0) +end + +# Evaluate scalar products and quadratic forms. +c1TQ = dot(c1, Q) +c1TQD = dot(c1, QD) +c2TQ = dot(c2, Q) +c2TQD = dot(c2, QD) +c12TQ = dot(c12, Q) +c12TQD= dot(c12, QD) +MQQ = zeros(Num, length(Q)) +KQQ = zeros(Num, length(Q)) +DQQD = zeros(Num, length(Q)) +QTBQ = zeros(Num, length(Q)) +BQQD = zeros(Num, length(Q)) +for I = 1:length(Q) + MQQ[I] = dot(MQ[:, I], Q) + KQQ[I] = dot(KQ[:, I], Q) + DQQD[I]= dot(DQ[:, I], QD) + QTBQ[I]= dot(Q, BQ[:, I]) + BQQD[I]= dot(BQ[I, :], QD) +end +QTMQQ = dot(Q, MQQ) +QDTMQQ = dot(QD, MQQ) +QDTBQQD = dot(QD, BQQD) + +# Compute mass matrix. +AM = zeros(Num, 7, 7) # NP×NP +AM[1,1] = J₁ + m₂*l₁^2 +AM[1,2] = .5*l₁*l₂*m₂*cosp12 +AM[2,2] = J₂ +AM[1,3] = 0.0 +AM[2,3] = 0.0 +AM[3,1] = 0.0 +AM[3,2] = 0.0 +AM[3,3] = m₃ +AM[1,2] = AM[1,2] + ρ*l₁*(sinp12*c2TQ+cosp12*c1TQ) +AM[2,2] = AM[2,2] + QTMQQ + 2.0*ρ*c12TQ +for I = 1:length(Q) + AM[1,3+I] = ρ*l₁*(-sinp12*c1[I] + cosp12*c2[I]) + AM[2,3+I] = ρ*c21[I] + ρ*QTBQ[I] + AM[3,3+I] = 0.0 +end +for I = 1:length(Q) + for J = 1:I + AM[3+J,3+I] = MQ[J,I] + end +end +for I = 1:size(AM,1) + for J = I+1:size(AM,2) + AM[J,I] = AM[I,J] + end +end + +KU = 4 +KV = 0 + +# Compute constraint matrix. +QKU = KU > 0 ? Q[KU] : Num(0) +QKV = KV > 0 ? Q[KV] : Num(0) +GP = zeros(Num, 3, 7) +GP[1,1] = l₁*cosp1 +GP[1,2] = l₂*cosp2 + QKU*cosp2 - QKV*sinp2 +GP[1,3] = 0.0 +GP[2,1] = l₁*sinp1 +GP[2,2] = l₂*sinp2 + QKU*sinp2 + QKV*cosp2 +GP[2,3] = 1.0 +GP[3,1] = 1.0 +GP[3,2] = 0.0 +GP[3,3] = 0.0 +for I = 1:length(Q) + GP[1,3+I] = 0.0 + GP[2,3+I] = 0.0 + GP[3,3+I] = 0.0 +end +if KU != 0 + GP[1,3+KU] = sinp2 + GP[2,3+KU] = -cosp2 +end +if KV != 0 + GP[1,3+KV] = cosp2 + GP[2,3+KV] = sinp2 +end + +# Forces - rigid motion entries. +F = zeros(Num, 7) +F[1] = -.5*l₁*γ*(m₁+2.0*m₂)*cosp1 + + -.5*l₁*l₂*m₂*V[2]*V[2]*sinp12 +F[2] = -.5*l₂*γ*m₂*cosp2 + + .5*l₁*l₂*m₂*V[1]*V[1]*sinp12 +F[3] = 0.0 +# Superposition of flexible motion (term f^e). +F[1] += ρ*l₁*V[2]*V[2]*(-sinp12*c1TQ+cosp12*c2TQ) + + - 2.0*ρ*l₁*V[2]*(cosp12*c1TQD+sinp12*c2TQD) +F[2] += ρ*l₁*V[1]*V[1]*(sinp12*c1TQ-cosp12*c2TQ) + + - 2.0*ρ*V[2]*c12TQD - 2.0*V[2]*QDTMQQ + + - ρ*QDTBQQD - ρ*γ*(cosp2*c1TQ-sinp2*c2TQ) +# Coriolis and gravity terms flexible motion (Gamma). +for I = 1:length(Q) + F[3+I] = V[2]*V[2]*MQQ[I] + + + ρ*(V[2]*V[2]*c12[I] + + l₁*V[1]*V[1]*(cosp12*c1[I]+sinp12*c2[I]) + + 2.0*V[2]*BQQD[I]) + + - ρ*γ*(sinp2*c1[I]+cosp2*c2[I]) +end +# Stiffness + damping terms - K q - D q`. +for I = 1:length(Q) + F[3+I] = F[3+I] - KQQ[I] - DQQD[I] +end + +#if IPAR[1] == 1 +# # Nonlinear stiffness term +# FACK = 0.5*EE*BB*HH/l₂^2*π^2 +# FACB = 80.0/(π^2*9.0) +# F[4] -= FACK*(Q[1]*Q[4]-FACB*Q[2]*(-4*Q[3]+2*Q[4])) +# F[5] -= FACK*(4*Q[2]*Q[4]-FACB*Q[1]*(-4*Q[3]+2*Q[4])) +# F[6] -= FACK*4.0*FACB*Q[1]*Q[2] +# F[7] -= FACK*(0.5*Q[1]^2+2*Q[2]^2-2*FACB*Q[1]*Q[2]) +#end + +# Dynamics part II ( M*w - f + G(T)*lambda ). +DELTA = zeros(Num, 7) +for I = 1:7 + DELTA[I] = dot(AM[:,I], DDX) + + - F[I] + GP[1,I]*λ₁+GP[2,I]*λ₂+GP[3,I]*λ₃ +end + +# Acceleration level constraints. +#ALC = zeros(Num, 3) +#QDKU = KU > 0 ? QD[KU] : Num(0) +#QDKV = KV > 0 ? QD[KV] : Num(0) +#ALC = zeros(Num, 3) +#ALC[1] = -l₁*sinp1*V[1]^2 - (l₂+QKU)*sinp2*V[2]^2 + +# 2.0*V[2]*(cosp2*QDKU-sinp2*QDKV) - cosp2*V[2]^2*QKV +#ALC[2] = l₁*cosp1*V[1]^2 + (l₂+QKU)*cosp2*V[2]^2 + +# 2.0*V[2]*(sinp2*QDKU+cosp2*QDKV) - sinp2*V[2]^2*QKV +#ALC[3] = 0.0 +#for I = 1:7 +# ALC[1] += GP[1,I]*DDX[I] +# ALC[2] += GP[2,I]*DDX[I] +# ALC[3] += GP[3,I]*DDX[I] +#end + +## Velocity level constraints. +#VLC = zeros(Num, 3) +#VLC[1] = 0.0 +#VLC[2] = 0.0 +#VLC[3] = -Ω +#for I = 1:7 +# VLC[1] += GP[1,I]*DX[I] +# VLC[2] += GP[2,I]*DX[I] +# VLC[3] += GP[3,I]*DX[I] +#end + +## Position level constraints. +PLC = zeros(Num, 3) +PLC[1] = l₁*sinp1 + l₂*sinp2 + QKU*sinp2 + QKV*cosp2 +PLC[2] = x₃ - l₁*cosp1 - l₂*cosp2 + -QKU*cosp2 + QKV*sinp2 +PLC[3] = ϕ₁ - Ω*t + +eqs = Symbolics.Equation[] +for D in DELTA + push!(eqs, 0 ~ D) +end +for P in PLC + push!(eqs, 0 ~ P) +end +for I = 1:7 + push!(eqs, D(X[I]) ~ DX[I]) + push!(eqs, D(DX[I]) ~ DDX[I]) +end +eqs diff --git a/src/ode/enright_pryce.jl b/src/ode/enright_pryce.jl index 5e08604..810e711 100644 --- a/src/ode/enright_pryce.jl +++ b/src/ode/enright_pryce.jl @@ -101,10 +101,10 @@ sc2sys = let ODESystem(sc2eqs, t, name = :sc2) end -sc2prob = ODEProblem{false}(sc2sys, y .=> 1.0, (0, 20.0), [β => 0.1], dt = 1e-2, cse = true) -sc3prob = ODEProblem{false}(sc2sys, y .=> 1.0, (0, 20.0), [β => 1.0], dt = 1e-2, cse = true) -sc3prob = ODEProblem{false}(sc2sys, y .=> 1.0, (0, 20.0), [β => 10.0], dt = 1e-2, cse = true) -sc3prob = ODEProblem{false}(sc2sys, y .=> 1.0, (0, 20.0), [β => 20.0], dt = 1e-2, cse = true) +sc2prob = ODEProblem{false}(sc2sys, y .=> 1.0, (0, 20.0), [β => 0.1], dt = 1e-2) +sc3prob = ODEProblem{false}(sc2sys, y .=> 1.0, (0, 20.0), [β => 1.0], dt = 1e-2) +sc4prob = ODEProblem{false}(sc2sys, y .=> 1.0, (0, 20.0), [β => 10.0], dt = 1e-2) +sc5prob = ODEProblem{false}(sc2sys, y .=> 1.0, (0, 20.0), [β => 20.0], dt = 1e-2) sd1sys = let sd1eqs = [D(y[1]) ~ 0.2 * (y[2] - y[1]), @@ -199,7 +199,7 @@ se3sys = let ODESystem(se3eqs, t, name = :se3) end -se3prob = ODEProblem{false}(se3sys, [y[1:2] .=> 1.0; y[3] => 0.0], (0, 500.0), dt = 0.02, cse = true) +se3prob = ODEProblem{false}(se3sys, [y[1:2] .=> 1.0; y[3] => 0.0], (0, 500.0), dt = 0.02) se4sys = let y = y y = y[1:4] @@ -402,7 +402,7 @@ nc1sys = let y = y ODESystem(nc1eqs, t, name = :nc1) end -nc1prob = ODEProblem{false}(nc1sys, [y[1] => 1.0; y[2:10] .=> 0.0], (0, 20.0), cse = true) +nc1prob = ODEProblem{false}(nc1sys, [y[1] => 1.0; y[2:10] .=> 0.0], (0, 20.0)) nc2sys = let y = y n = 10 @@ -488,7 +488,7 @@ end function make_ds(nd1sys, e) y = collect(@nonamespace nd1sys.y) y0 = [y[1] => 1-e; y[2:3] .=> 0.0; y[4] => sqrt((1 + e) / (1 - e))] - ODEProblem{false}(nd1sys, y0, (0, 20.0), [ε => e], cse = true) + ODEProblem{false}(nd1sys, y0, (0, 20.0), [ε => e]) end nd1prob = make_ds(nd1sys, 0.1) nd2prob = make_ds(nd1sys, 0.3) @@ -504,7 +504,7 @@ ne1sys = let end y0 = [y[1] => 0.6713967071418030; y[2] => 0.09540051444747446] -ne1prob = ODEProblem{false}(ne1sys, y0, (0, 20.0), cse = true) +ne1prob = ODEProblem{false}(ne1sys, y0, (0, 20.0)) ne2sys = let ne2eqs = [D(y[1]) ~ y[2], @@ -523,7 +523,7 @@ ne3sys = let ODESystem(ne3eqs, t, name = :ne3) end -ne3prob = ODEProblem{false}(ne3sys, y[1:2] .=> 0, (0, 20.0), cse = true) +ne3prob = ODEProblem{false}(ne3sys, y[1:2] .=> 0, (0, 20.0)) ne4sys = let ne4eqs = [D(y[1]) ~ y[2], @@ -532,22 +532,23 @@ ne4sys = let ODESystem(ne4eqs, t, name = :ne4) end -ne4prob = ODEProblem{false}(ne4sys, [y[1] => 30.0, y[2] => 0.0], (0, 20.0), cse = true) +ne4prob = ODEProblem{false}(ne4sys, [y[1] => 30.0, y[2] => 0.0], (0, 20.0)) ne5sys = let ne5eqs = [D(y[1]) ~ y[2], - D(y[2]) ~ sqrt(1 - y[2]^2) / (25 - t),] + D(y[2]) ~ sqrt(1 + y[2]^2) / (25 - t),] ODESystem(ne5eqs, t, name = :ne5) end -ne5prob = ODEProblem{false}(ne5sys, y[1:2] .=> 0.0, (0, 20.0), cse = true) +ne5prob = ODEProblem{false}(ne5sys, y[1:2] .=> 0.0, (0, 20.0)) +@inline myifelse(x, y, z) = ifelse(x, y, z) 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 + IfElse.ifelse(cond, 1, -1)] + D(y[2]) ~ b + term(myifelse, cond, 1, -1, type=Real)] ODESystem(nf1eqs, t, name = :nf1) end @@ -555,11 +556,11 @@ nf1prob = ODEProblem{false}(nf1sys, y[1:2] .=> 0.0, (0, 20.0)) nf2sys = let cond = term(iseven, term(floor, Int, unwrap(t), type = Int), type = Bool) - nf2eqs = [D(y[1]) ~ 55 - IfElse.ifelse(cond, 2y[1]/2, y[1]/2)] + nf2eqs = [D(y[1]) ~ 55 - term(myifelse, cond, 2y[1]/2, y[1]/2, type=Real)] ODESystem(nf2eqs, t, name = :nf2) end -nf2prob = ODEProblem{false}(nf2sys, [y[1] .=> 110.0], (0, 20.0), cse = true) +nf2prob = ODEProblem{false}(nf2sys, [y[1] .=> 110.0], (0, 20.0)) nf3sys = let nf3eqs = [D(y[1]) ~ y[2], @@ -567,14 +568,14 @@ nf3sys = let ODESystem(nf3eqs, t, name = :nf3) end -nf3prob = ODEProblem{false}(nf3sys, y[1:2] .=> 0.0, (0, 20.0), cse = true) +nf3prob = ODEProblem{false}(nf3sys, y[1:2] .=> 0.0, (0, 20.0)) nf4sys = let - nf4eqs = [D(y[1]) ~ IfElse.ifelse(t <= 10, -2/21 - (120 * (t - 5)) / (1 + 4 * (t - 5)^2), -2y[1])] + 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 -nf4prob = ODEProblem{false}(nf4sys, [y[1] => 1.0], (0, 20.0), cse = true) +nf4prob = ODEProblem{false}(nf4sys, [y[1] => 1.0], (0, 20.0)) nf5sys = let c = sum(i->cbrt(i)^4, 1:19) @@ -583,4 +584,4 @@ nf5sys = let ODESystem(nf5eqs, t, name = :nf5) end -nf5prob = ODEProblem{false}(nf5sys, [y[1] => 1.0], (0, 20.0), cse = true) +nf5prob = ODEProblem{false}(nf5sys, [y[1] => 1.0], (0, 20.0)) From 08f5fbc142342aa37b49d188576dffc0b9911b5d Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Fri, 25 Oct 2024 04:47:05 +0000 Subject: [PATCH 4/5] a few updates --- lib/ODEProblemLibrary/Project.toml | 1 + .../ODEProblemLibrary/src}/enright_pryce.jl | 61 +++++++++---------- 2 files changed, 31 insertions(+), 31 deletions(-) rename {src/ode => lib/ODEProblemLibrary/src}/enright_pryce.jl (92%) diff --git a/lib/ODEProblemLibrary/Project.toml b/lib/ODEProblemLibrary/Project.toml index c341ff7..9f91d49 100644 --- a/lib/ODEProblemLibrary/Project.toml +++ b/lib/ODEProblemLibrary/Project.toml @@ -10,6 +10,7 @@ Markdown = "d6f4376e-aef5-505a-96c1-9c027394607a" ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" RuntimeGeneratedFunctions = "7e49a35a-f44a-4d26-94aa-eba1b4ca6b47" +StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" [compat] Aqua = "0.5" diff --git a/src/ode/enright_pryce.jl b/lib/ODEProblemLibrary/src/enright_pryce.jl similarity index 92% rename from src/ode/enright_pryce.jl rename to lib/ODEProblemLibrary/src/enright_pryce.jl index 810e711..79b2a66 100644 --- a/src/ode/enright_pryce.jl +++ b/lib/ODEProblemLibrary/src/enright_pryce.jl @@ -1,6 +1,5 @@ using ModelingToolkit -using IfElse -using Symbolics +using ModelingToolkit.Symbolics using Symbolics: unwrap using DiffEqBase, StaticArrays, LinearAlgebra @@ -124,7 +123,7 @@ sd2sys = let ODESystem(sd2eqs, t, name = :sd2) end -sd2prob = ODEProblem{false}(sd2sys, [y[1] => 1.0; y[2:3] .=> 0.0], (0, 40.0), dt = 1e-5, cse = true) +sd2prob = ODEProblem{false}(sd2sys, [y[1] => 1.0; y[2:3] .=> 0.0], (0, 40.0), dt = 1e-5) sd3sys = let sd3eqs = [D(y[1]) ~ y[3] - 100 * (y[1] * y[2]), @@ -135,7 +134,7 @@ sd3sys = let ODESystem(sd3eqs, t, name = :sd3) end -sd3prob = ODEProblem{false}(sd3sys, [y[1:2] .=> 1; y[3:4] .=> 0.0], (0, 20.0), dt = 2.5e-5, cse = true) +sd3prob = ODEProblem{false}(sd3sys, [y[1:2] .=> 1; y[3:4] .=> 0.0], (0, 20.0), dt = 2.5e-5) sd4sys = let sd4eqs = [D(y[1]) ~ -0.013y[1] - 1000 * (y[1] * y[3]), @@ -145,7 +144,7 @@ sd4sys = let ODESystem(sd4eqs, t, name = :sd4) end -sd4prob = ODEProblem{false}(sd4sys, [y[1:2] .=> 1; y[3] => 0.0], (0, 50.0), dt = 2.9e-4, cse = true) +sd4prob = ODEProblem{false}(sd4sys, [y[1:2] .=> 1; y[3] => 0.0], (0, 50.0), dt = 2.9e-4) sd5sys = let sd5eqs = [D(y[1]) ~ 0.01 - (1 + (y[1] + 1000) * (y[1] + 1)) * (0.01 + y[1] + y[2]), @@ -154,7 +153,7 @@ sd5sys = let ODESystem(sd5eqs, t, name = :sd5) end -sd5prob = ODEProblem{false}(sd5sys, y[1:2] .=> 0.0, (0, 100.0), dt = 1e-4, cse = true) +sd5prob = ODEProblem{false}(sd5sys, y[1:2] .=> 0.0, (0, 100.0), dt = 1e-4) sd6sys = let sd6eqs = [D(y[1]) ~ -y[1] + 10^8 * y[3] * (1 - y[1]), @@ -165,7 +164,7 @@ sd6sys = let ODESystem(sd6eqs, t, name = :sd6) end -sd6prob = ODEProblem{false}(sd6sys, [y[1] => 1.0; y[2:3] .=> 0.0], (0, 1.0), dt = 3.3e-8, cse = true) +sd6prob = ODEProblem{false}(sd6sys, [y[1] => 1.0; y[2:3] .=> 0.0], (0, 1.0), dt = 3.3e-8) se1sys = let Γ = 100 @@ -178,7 +177,7 @@ se1sys = let ODESystem(se1eqs, t, name = :se1) end -se1prob = ODEProblem{false}(se1sys, y .=> 0.0, (0, 1.0), dt = 6.8e-3, cse = true) +se1prob = ODEProblem{false}(se1sys, y .=> 0.0, (0, 1.0), dt = 6.8e-3) se2sys = let se2eqs = [D(y[1]) ~ y[2], @@ -188,7 +187,7 @@ se2sys = let ODESystem(se2eqs, t, name = :se2) end -se2prob = ODEProblem{false}(se2sys, [y[1] => 2.0, y[2] => 0.0], (0, 1.0), dt = 1e-3, cse = true) +se2prob = ODEProblem{false}(se2sys, [y[1] => 2.0, y[2] => 0.0], (0, 1.0), dt = 1e-3) se3sys = let se3eqs = [D(y[1]) ~ -(55 + y[3]) * y[1] + 65 * y[2], @@ -214,7 +213,7 @@ se4sys = let y = y ODESystem(se4eqs, t, name = :se4) 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) +se4prob = ODEProblem{false}(se4sys, [y[1] => 0.0; y[2] => -2.0; y[3:4] .=> -1.0], (0, 1000.0), dt = 1e-3) se5sys = let se5eqs = [ @@ -227,7 +226,7 @@ se5sys = let ODESystem(se5eqs, t, name = :se5) end -se5prob = ODEProblem{false}(se5sys, [y[1] => 1.76e-3; y[2:4] .=> 0.0], (0, 1000.0), dt = 1e-3, cse = true) +se5prob = ODEProblem{false}(se5sys, [y[1] => 1.76e-3; y[2:4] .=> 0.0], (0, 1000.0), dt = 1e-3) sf1sys = let k = exp(20.7 - 1500 / y[1]) @@ -241,7 +240,7 @@ sf1sys = let ODESystem(sf1eqs, t, name = :sf1) 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) +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) sf2sys = let sf2eqs = [ @@ -252,7 +251,7 @@ sf2sys = let ODESystem(sf2eqs, t, name = :sf2) end -sf2prob = ODEProblem{false}(sf2sys, [y[1] => 1.0; y[2] => 0.0], (0, 240.0), dt = 1e-2, cse = true) +sf2prob = ODEProblem{false}(sf2sys, [y[1] => 1.0; y[2] => 0.0], (0, 240.0), dt = 1e-2) sf3sys = let sf3eqs = [ @@ -266,7 +265,7 @@ sf3sys = let ODESystem(sf3eqs, t, name = :sf3) 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) +sf3prob = ODEProblem{false}(sf3sys, [y[1] => 4e-6; y[2] => 1e-6; y[3:5] .=> 0.0], (0, 100.0), dt = 1e-6) sf4sys = let sf4eqs = [ @@ -278,7 +277,7 @@ sf4sys = let ODESystem(sf4eqs, t, name = :sf4) 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) +sf4prob = ODEProblem{false}(sf4sys, [y[1] => 4.0; y[2] => 1.1; y[3] => 4.0], (0, 300.0), dt = 1e-3) sf5sys = let sf5eqs = [ @@ -291,7 +290,7 @@ sf5sys = let ODESystem(sf5eqs, t, name = :sf5) 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) +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) # Non-stiff na1sys = let y = y[1] @@ -300,7 +299,7 @@ na1sys = let y = y[1] ODESystem(na1eqs, t, name = :na1) end -na1prob = ODEProblem{false}(na1sys, [y[1] => 1], (0, 20.0), cse = true) +na1prob = ODEProblem{false}(na1sys, [y[1] => 1], (0, 20.0)) na2sys = let y = y[1] na2eqs = D(y) ~ -y^3/2 @@ -308,7 +307,7 @@ na2sys = let y = y[1] ODESystem(na2eqs, t, name = :na2) end -na2prob = ODEProblem{false}(na2sys, [y[1] => 1], (0, 20.0), cse = true) +na2prob = ODEProblem{false}(na2sys, [y[1] => 1], (0, 20.0)) na3sys = let y = y[1] na3eqs = D(y) ~ y * cos(t) @@ -316,7 +315,7 @@ na3sys = let y = y[1] ODESystem(na3eqs, t, name = :na3) end -na3prob = ODEProblem{false}(na3sys, [y[1] => 1], (0, 20.0), cse = true) +na3prob = ODEProblem{false}(na3sys, [y[1] => 1], (0, 20.0)) na4sys = let y = y[1] na4eqs = D(y) ~ y/4 * (1 - y/20) @@ -324,7 +323,7 @@ na4sys = let y = y[1] ODESystem(na4eqs, t, name = :na4) end -na4prob = ODEProblem{false}(na4sys, [y[1] => 1], (0, 20.0), cse = true) +na4prob = ODEProblem{false}(na4sys, [y[1] => 1], (0, 20.0)) na5sys = let y = y[1] na5eqs = D(y) ~ (y - t) / (y + t) @@ -332,7 +331,7 @@ na5sys = let y = y[1] ODESystem(na5eqs, t, name = :na5) end -na5prob = ODEProblem{false}(na5sys, [y[1] => 4], (0, 20.0), cse = true) +na5prob = ODEProblem{false}(na5sys, [y[1] => 4], (0, 20.0)) nb1sys = let nb1eqs = [ @@ -343,7 +342,7 @@ nb1sys = let ODESystem(nb1eqs, t, name = :nb1) end -nb1prob = ODEProblem{false}(nb1sys, [y[1] => 1.0, y[2] => 3], (0, 20.0), cse = true) +nb1prob = ODEProblem{false}(nb1sys, [y[1] => 1.0, y[2] => 3], (0, 20.0)) nb2sys = let nb2eqs = [ @@ -355,7 +354,7 @@ nb2sys = let ODESystem(nb2eqs, t, name = :nb2) end -nb2prob = ODEProblem{false}(nb2sys, [y[1] => 2.0, y[2] => 0.0, y[3] => 1.0], (0, 20.0), cse = true) +nb2prob = ODEProblem{false}(nb2sys, [y[1] => 2.0, y[2] => 0.0, y[3] => 1.0], (0, 20.0)) nb3sys = let nb3eqs = [ @@ -367,7 +366,7 @@ nb3sys = let ODESystem(nb3eqs, t, name = :nb3) end -nb3prob = ODEProblem{false}(nb3sys, [y[1] => 1.0; y[2:3] .=> 0.0], (0, 20.0), cse = true) +nb3prob = ODEProblem{false}(nb3sys, [y[1] => 1.0; y[2:3] .=> 0.0], (0, 20.0)) nb4sys = let r = sqrt(y[1]^2 + y[2]^2) @@ -380,7 +379,7 @@ nb4sys = let ODESystem(nb4eqs, t, name = :nb4) end -nb4prob = ODEProblem{false}(nb4sys, [y[1] => 3.0; y[2:3] .=> 0.0], (0, 20.0), cse = true) +nb4prob = ODEProblem{false}(nb4sys, [y[1] => 3.0; y[2:3] .=> 0.0], (0, 20.0)) nb5sys = let nb5eqs = [ @@ -392,7 +391,7 @@ nb5sys = let ODESystem(nb5eqs, t, name = :nb5) end -nb5prob = ODEProblem{false}(nb5sys, [y[1] => 0.0; y[2:3] .=> 1.0], (0, 20.0), cse = true) +nb5prob = ODEProblem{false}(nb5sys, [y[1] => 0.0; y[2:3] .=> 1.0], (0, 20.0)) nc1sys = let y = y n = 10 @@ -412,7 +411,7 @@ nc2sys = let y = y ODESystem(nc2eqs, t, name = :nc2) end -nc2prob = ODEProblem{false}(nc2sys, [y[1] => 1.0; y[2:10] .=> 0.0], (0, 20.0), cse = true) +nc2prob = ODEProblem{false}(nc2sys, [y[1] => 1.0; y[2:10] .=> 0.0], (0, 20.0)) nc3sys = let y = y n = 10 @@ -422,7 +421,7 @@ nc3sys = let y = y ODESystem(nc3eqs, t, name = :nc3) end -nc3prob = ODEProblem{false}(nc3sys, [y[1] => 1.0; y[2:10] .=> 0.0], (0, 20.0), cse = true) +nc3prob = ODEProblem{false}(nc3sys, [y[1] => 1.0; y[2:10] .=> 0.0], (0, 20.0)) @variables y(t)[1:51] y = collect(y) @@ -434,7 +433,7 @@ nc4sys = let y = y ODESystem(nc4eqs, t, name = :nc4) end -nc4prob = ODEProblem{false}(nc4sys, [y[1] => 1.0; y[2:51] .=> 0.0], (0, 20.0), cse = true) +nc4prob = ODEProblem{false}(nc4sys, [y[1] => 1.0; y[2:51] .=> 0.0], (0, 20.0)) @variables y(t)[1:3, 1:5] y = collect(y) @@ -470,7 +469,7 @@ 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) +nc5prob = ODEProblem{false}(nc5sys, [y0; y0′], (0, 20.0)) @variables y(t)[1:4] y = collect(y) @@ -514,7 +513,7 @@ ne2sys = let end y0 = [y[1] => 2.0; y[2] => 0.0] -ne2prob = ODEProblem{false}(ne2sys, y0, (0, 20.0), cse = true) +ne2prob = ODEProblem{false}(ne2sys, y0, (0, 20.0)) ne3sys = let ne3eqs = [D(y[1]) ~ y[2], From ceef2e7778410728091aecfd6ff9abfbb41aa418 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Fri, 25 Oct 2024 05:01:53 +0000 Subject: [PATCH 5/5] Get scripts running again --- lib/ODEProblemLibrary/src/enright_pryce.jl | 125 +++++++++++---------- m.jl => lib/ODEProblemLibrary/src/m.jl | 6 +- 2 files changed, 66 insertions(+), 65 deletions(-) rename m.jl => lib/ODEProblemLibrary/src/m.jl (99%) diff --git a/lib/ODEProblemLibrary/src/enright_pryce.jl b/lib/ODEProblemLibrary/src/enright_pryce.jl index 79b2a66..a031ae2 100644 --- a/lib/ODEProblemLibrary/src/enright_pryce.jl +++ b/lib/ODEProblemLibrary/src/enright_pryce.jl @@ -1,11 +1,11 @@ using ModelingToolkit +using ModelingToolkit: t_nounits as t, D_nounits as D using ModelingToolkit.Symbolics using Symbolics: unwrap using DiffEqBase, StaticArrays, LinearAlgebra -@variables t y(t)[1:10] +@variables y(t)[1:10] y = collect(y) -D = Differential(t) # Stiff sa1sys = let @@ -13,7 +13,7 @@ sa1sys = let ODESystem(sa1eqs, t, name = :sa1) end -sa1prob = ODEProblem{false}(sa1sys, y[1:4] .=> 1.0, (0, 20.0), dt = 1e-2) +sa1prob = ODEProblem{false}(structural_simplify(sa1sys), y[1:4] .=> 1.0, (0, 20.0), dt = 1e-2) sa2sys = let sa2eqs = [D(y[1]) ~ -1800 * y[1] + 900 * y[2]] @@ -24,7 +24,7 @@ sa2sys = let ODESystem(sa2eqs, t, name = :sa2) end -sa2prob = ODEProblem{false}(sa2sys, y[1:9] .=> 0.0, (0, 120.0), dt = 5e-4) +sa2prob = ODEProblem{false}(structural_simplify(sa2sys), y[1:9] .=> 0.0, (0, 120.0), dt = 5e-4) sa3sys = let sa3eqs = [ @@ -36,14 +36,14 @@ sa3sys = let ODESystem(sa3eqs, t, name = :sa3) end -sa3prob = ODEProblem{false}(sa3sys, y[1:4] .=> 1.0, (0, 20.0), dt = 1e-5) +sa3prob = ODEProblem{false}(structural_simplify(sa3sys), y[1:4] .=> 1.0, (0, 20.0), dt = 1e-5) sa4sys = let sa4eqs = [D(y[i]) ~ -i^5 * y[i] for i in 1:10] ODESystem(sa4eqs, t, name = :sa4) end -sa4prob = ODEProblem{false}(sa4sys, y[1:10] .=> 1.0, (0, 1.0), dt = 1e-5) +sa4prob = ODEProblem{false}(structural_simplify(sa4sys), y[1:10] .=> 1.0, (0, 1.0), dt = 1e-5) #const SA_PROBLEMS = [sa1prob, sa2prob, sa3prob, sa4prob] @@ -56,7 +56,7 @@ sb1sys = let ODESystem(sb1eqs, t, name = :sb1) end -sb1prob = ODEProblem{false}(sb1sys, [y[[1, 3]] .=> 1.0; y[[2, 4]] .=> 0.0;], (0, 20.0), dt = 7e-3) +sb1prob = ODEProblem{false}(structural_simplify(sb1sys), [y[[1, 3]] .=> 1.0; y[[2, 4]] .=> 0.0;], (0, 20.0), dt = 7e-3) @parameters α @@ -71,10 +71,10 @@ sb2sys = let ODESystem(sb2eqs, t, name = :sb2) end -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) +sb2prob = ODEProblem{false}(structural_simplify(sb2sys), y .=> 1.0, (0, 20.0), [α => 3], dt = 1e-2) +sb3prob = ODEProblem{false}(structural_simplify(sb2sys), y .=> 1.0, (0, 20.0), [α => 8], dt = 1e-2) +sb4prob = ODEProblem{false}(structural_simplify(sb2sys), y .=> 1.0, (0, 20.0), [α => 25], dt = 1e-2) +sb5prob = ODEProblem{false}(structural_simplify(sb2sys), y .=> 1.0, (0, 20.0), [α => 100], dt = 1e-2) sc1sys = let @@ -87,7 +87,7 @@ sc1sys = let ODESystem(sc1eqs, t, name = :sc1) end -sc1prob = ODEProblem{false}(sc1sys, y .=> 1.0, (0, 20.0), dt = 1e-2) +sc1prob = ODEProblem{false}(structural_simplify(sc1sys), y .=> 1.0, (0, 20.0), dt = 1e-2) @parameters β @@ -100,10 +100,10 @@ sc2sys = let ODESystem(sc2eqs, t, name = :sc2) end -sc2prob = ODEProblem{false}(sc2sys, y .=> 1.0, (0, 20.0), [β => 0.1], dt = 1e-2) -sc3prob = ODEProblem{false}(sc2sys, y .=> 1.0, (0, 20.0), [β => 1.0], dt = 1e-2) -sc4prob = ODEProblem{false}(sc2sys, y .=> 1.0, (0, 20.0), [β => 10.0], dt = 1e-2) -sc5prob = ODEProblem{false}(sc2sys, y .=> 1.0, (0, 20.0), [β => 20.0], dt = 1e-2) +sc2prob = ODEProblem{false}(structural_simplify(sc2sys), y .=> 1.0, (0, 20.0), [β => 0.1], dt = 1e-2) +sc3prob = ODEProblem{false}(structural_simplify(sc2sys), y .=> 1.0, (0, 20.0), [β => 1.0], dt = 1e-2) +sc4prob = ODEProblem{false}(structural_simplify(sc2sys), y .=> 1.0, (0, 20.0), [β => 10.0], dt = 1e-2) +sc5prob = ODEProblem{false}(structural_simplify(sc2sys), y .=> 1.0, (0, 20.0), [β => 20.0], dt = 1e-2) sd1sys = let sd1eqs = [D(y[1]) ~ 0.2 * (y[2] - y[1]), @@ -113,7 +113,7 @@ sd1sys = let ODESystem(sd1eqs, t, name = :sd1) end -sd1prob = ODEProblem{false}(sd1sys, y .=> 0.0, (0, 400.0), [β => 0.1], dt = 1.7e-2) +sd1prob = ODEProblem{false}(structural_simplify(sd1sys), y .=> 0.0, (0, 400.0), [β => 0.1], dt = 1.7e-2) sd2sys = let sd2eqs = [D(y[1]) ~ -0.04y[1] + 0.01 * (y[2] * y[3]), @@ -123,7 +123,7 @@ sd2sys = let ODESystem(sd2eqs, t, name = :sd2) end -sd2prob = ODEProblem{false}(sd2sys, [y[1] => 1.0; y[2:3] .=> 0.0], (0, 40.0), dt = 1e-5) +sd2prob = ODEProblem{false}(structural_simplify(sd2sys), [y[1] => 1.0; y[2:3] .=> 0.0], (0, 40.0), dt = 1e-5) sd3sys = let sd3eqs = [D(y[1]) ~ y[3] - 100 * (y[1] * y[2]), @@ -134,7 +134,7 @@ sd3sys = let ODESystem(sd3eqs, t, name = :sd3) end -sd3prob = ODEProblem{false}(sd3sys, [y[1:2] .=> 1; y[3:4] .=> 0.0], (0, 20.0), dt = 2.5e-5) +sd3prob = ODEProblem{false}(structural_simplify(sd3sys), [y[1:2] .=> 1; y[3:4] .=> 0.0], (0, 20.0), dt = 2.5e-5) sd4sys = let sd4eqs = [D(y[1]) ~ -0.013y[1] - 1000 * (y[1] * y[3]), @@ -144,7 +144,7 @@ sd4sys = let ODESystem(sd4eqs, t, name = :sd4) end -sd4prob = ODEProblem{false}(sd4sys, [y[1:2] .=> 1; y[3] => 0.0], (0, 50.0), dt = 2.9e-4) +sd4prob = ODEProblem{false}(structural_simplify(sd4sys), [y[1:2] .=> 1; y[3] => 0.0], (0, 50.0), dt = 2.9e-4) sd5sys = let sd5eqs = [D(y[1]) ~ 0.01 - (1 + (y[1] + 1000) * (y[1] + 1)) * (0.01 + y[1] + y[2]), @@ -153,7 +153,7 @@ sd5sys = let ODESystem(sd5eqs, t, name = :sd5) end -sd5prob = ODEProblem{false}(sd5sys, y[1:2] .=> 0.0, (0, 100.0), dt = 1e-4) +sd5prob = ODEProblem{false}(structural_simplify(sd5sys), y[1:2] .=> 0.0, (0, 100.0), dt = 1e-4) sd6sys = let sd6eqs = [D(y[1]) ~ -y[1] + 10^8 * y[3] * (1 - y[1]), @@ -164,7 +164,7 @@ sd6sys = let ODESystem(sd6eqs, t, name = :sd6) end -sd6prob = ODEProblem{false}(sd6sys, [y[1] => 1.0; y[2:3] .=> 0.0], (0, 1.0), dt = 3.3e-8) +sd6prob = ODEProblem{false}(structural_simplify(sd6sys), [y[1] => 1.0; y[2:3] .=> 0.0], (0, 1.0), dt = 3.3e-8) se1sys = let Γ = 100 @@ -177,7 +177,7 @@ se1sys = let ODESystem(se1eqs, t, name = :se1) end -se1prob = ODEProblem{false}(se1sys, y .=> 0.0, (0, 1.0), dt = 6.8e-3) +se1prob = ODEProblem{false}(structural_simplify(se1sys), y .=> 0.0, (0, 1.0), dt = 6.8e-3) se2sys = let se2eqs = [D(y[1]) ~ y[2], @@ -187,7 +187,7 @@ se2sys = let ODESystem(se2eqs, t, name = :se2) end -se2prob = ODEProblem{false}(se2sys, [y[1] => 2.0, y[2] => 0.0], (0, 1.0), dt = 1e-3) +se2prob = ODEProblem{false}(structural_simplify(se2sys), [y[1] => 2.0, y[2] => 0.0], (0, 1.0), dt = 1e-3) se3sys = let se3eqs = [D(y[1]) ~ -(55 + y[3]) * y[1] + 65 * y[2], @@ -198,7 +198,7 @@ se3sys = let ODESystem(se3eqs, t, name = :se3) end -se3prob = ODEProblem{false}(se3sys, [y[1:2] .=> 1.0; y[3] => 0.0], (0, 500.0), dt = 0.02) +se3prob = ODEProblem{false}(structural_simplify(se3sys), [y[1:2] .=> 1.0; y[3] => 0.0], (0, 500.0), dt = 0.02) se4sys = let y = y y = y[1:4] @@ -213,7 +213,7 @@ se4sys = let y = y ODESystem(se4eqs, t, name = :se4) end -se4prob = ODEProblem{false}(se4sys, [y[1] => 0.0; y[2] => -2.0; y[3:4] .=> -1.0], (0, 1000.0), dt = 1e-3) +se4prob = ODEProblem{false}(structural_simplify(se4sys), [y[1] => 0.0; y[2] => -2.0; y[3:4] .=> -1.0], (0, 1000.0), dt = 1e-3) se5sys = let se5eqs = [ @@ -226,7 +226,7 @@ se5sys = let ODESystem(se5eqs, t, name = :se5) end -se5prob = ODEProblem{false}(se5sys, [y[1] => 1.76e-3; y[2:4] .=> 0.0], (0, 1000.0), dt = 1e-3) +se5prob = ODEProblem{false}(structural_simplify(se5sys), [y[1] => 1.76e-3; y[2:4] .=> 0.0], (0, 1000.0), dt = 1e-3) sf1sys = let k = exp(20.7 - 1500 / y[1]) @@ -240,7 +240,7 @@ sf1sys = let ODESystem(sf1eqs, t, name = :sf1) 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) +sf1prob = ODEProblem{false}(structural_simplify(sf1sys), [y[1] => 761.0; y[2] => 0.0; y[3] => 600.0; y[4] => 0.1], (0, 1000.0), dt = 1e-4) sf2sys = let sf2eqs = [ @@ -251,7 +251,7 @@ sf2sys = let ODESystem(sf2eqs, t, name = :sf2) end -sf2prob = ODEProblem{false}(sf2sys, [y[1] => 1.0; y[2] => 0.0], (0, 240.0), dt = 1e-2) +sf2prob = ODEProblem{false}(structural_simplify(sf2sys), [y[1] => 1.0; y[2] => 0.0], (0, 240.0), dt = 1e-2) sf3sys = let sf3eqs = [ @@ -265,7 +265,7 @@ sf3sys = let ODESystem(sf3eqs, t, name = :sf3) end -sf3prob = ODEProblem{false}(sf3sys, [y[1] => 4e-6; y[2] => 1e-6; y[3:5] .=> 0.0], (0, 100.0), dt = 1e-6) +sf3prob = ODEProblem{false}(structural_simplify(sf3sys), [y[1] => 4e-6; y[2] => 1e-6; y[3:5] .=> 0.0], (0, 100.0), dt = 1e-6) sf4sys = let sf4eqs = [ @@ -277,7 +277,7 @@ sf4sys = let ODESystem(sf4eqs, t, name = :sf4) end -sf4prob = ODEProblem{false}(sf4sys, [y[1] => 4.0; y[2] => 1.1; y[3] => 4.0], (0, 300.0), dt = 1e-3) +sf4prob = ODEProblem{false}(structural_simplify(sf4sys), [y[1] => 4.0; y[2] => 1.1; y[3] => 4.0], (0, 300.0), dt = 1e-3) sf5sys = let sf5eqs = [ @@ -290,7 +290,7 @@ sf5sys = let ODESystem(sf5eqs, t, name = :sf5) 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) +sf5prob = ODEProblem{false}(structural_simplify(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) # Non-stiff na1sys = let y = y[1] @@ -299,7 +299,7 @@ na1sys = let y = y[1] ODESystem(na1eqs, t, name = :na1) end -na1prob = ODEProblem{false}(na1sys, [y[1] => 1], (0, 20.0)) +na1prob = ODEProblem{false}(structural_simplify(na1sys), [y[1] => 1], (0, 20.0)) na2sys = let y = y[1] na2eqs = D(y) ~ -y^3/2 @@ -307,7 +307,7 @@ na2sys = let y = y[1] ODESystem(na2eqs, t, name = :na2) end -na2prob = ODEProblem{false}(na2sys, [y[1] => 1], (0, 20.0)) +na2prob = ODEProblem{false}(structural_simplify(na2sys), [y[1] => 1], (0, 20.0)) na3sys = let y = y[1] na3eqs = D(y) ~ y * cos(t) @@ -315,7 +315,7 @@ na3sys = let y = y[1] ODESystem(na3eqs, t, name = :na3) end -na3prob = ODEProblem{false}(na3sys, [y[1] => 1], (0, 20.0)) +na3prob = ODEProblem{false}(structural_simplify(na3sys), [y[1] => 1], (0, 20.0)) na4sys = let y = y[1] na4eqs = D(y) ~ y/4 * (1 - y/20) @@ -323,7 +323,7 @@ na4sys = let y = y[1] ODESystem(na4eqs, t, name = :na4) end -na4prob = ODEProblem{false}(na4sys, [y[1] => 1], (0, 20.0)) +na4prob = ODEProblem{false}(structural_simplify(na4sys), [y[1] => 1], (0, 20.0)) na5sys = let y = y[1] na5eqs = D(y) ~ (y - t) / (y + t) @@ -331,7 +331,7 @@ na5sys = let y = y[1] ODESystem(na5eqs, t, name = :na5) end -na5prob = ODEProblem{false}(na5sys, [y[1] => 4], (0, 20.0)) +na5prob = ODEProblem{false}(structural_simplify(na5sys), [y[1] => 4], (0, 20.0)) nb1sys = let nb1eqs = [ @@ -342,7 +342,7 @@ nb1sys = let ODESystem(nb1eqs, t, name = :nb1) end -nb1prob = ODEProblem{false}(nb1sys, [y[1] => 1.0, y[2] => 3], (0, 20.0)) +nb1prob = ODEProblem{false}(structural_simplify(nb1sys), [y[1] => 1.0, y[2] => 3], (0, 20.0)) nb2sys = let nb2eqs = [ @@ -354,7 +354,7 @@ nb2sys = let ODESystem(nb2eqs, t, name = :nb2) end -nb2prob = ODEProblem{false}(nb2sys, [y[1] => 2.0, y[2] => 0.0, y[3] => 1.0], (0, 20.0)) +nb2prob = ODEProblem{false}(structural_simplify(nb2sys), [y[1] => 2.0, y[2] => 0.0, y[3] => 1.0], (0, 20.0)) nb3sys = let nb3eqs = [ @@ -366,7 +366,7 @@ nb3sys = let ODESystem(nb3eqs, t, name = :nb3) end -nb3prob = ODEProblem{false}(nb3sys, [y[1] => 1.0; y[2:3] .=> 0.0], (0, 20.0)) +nb3prob = ODEProblem{false}(structural_simplify(nb3sys), [y[1] => 1.0; y[2:3] .=> 0.0], (0, 20.0)) nb4sys = let r = sqrt(y[1]^2 + y[2]^2) @@ -379,7 +379,7 @@ nb4sys = let ODESystem(nb4eqs, t, name = :nb4) end -nb4prob = ODEProblem{false}(nb4sys, [y[1] => 3.0; y[2:3] .=> 0.0], (0, 20.0)) +nb4prob = ODEProblem{false}(structural_simplify(nb4sys), [y[1] => 3.0; y[2:3] .=> 0.0], (0, 20.0)) nb5sys = let nb5eqs = [ @@ -391,7 +391,7 @@ nb5sys = let ODESystem(nb5eqs, t, name = :nb5) end -nb5prob = ODEProblem{false}(nb5sys, [y[1] => 0.0; y[2:3] .=> 1.0], (0, 20.0)) +nb5prob = ODEProblem{false}(structural_simplify(nb5sys), [y[1] => 0.0; y[2:3] .=> 1.0], (0, 20.0)) nc1sys = let y = y n = 10 @@ -401,7 +401,7 @@ nc1sys = let y = y ODESystem(nc1eqs, t, name = :nc1) end -nc1prob = ODEProblem{false}(nc1sys, [y[1] => 1.0; y[2:10] .=> 0.0], (0, 20.0)) +nc1prob = ODEProblem{false}(structural_simplify(nc1sys), [y[1] => 1.0; y[2:10] .=> 0.0], (0, 20.0)) nc2sys = let y = y n = 10 @@ -411,7 +411,7 @@ nc2sys = let y = y ODESystem(nc2eqs, t, name = :nc2) end -nc2prob = ODEProblem{false}(nc2sys, [y[1] => 1.0; y[2:10] .=> 0.0], (0, 20.0)) +nc2prob = ODEProblem{false}(structural_simplify(nc2sys), [y[1] => 1.0; y[2:10] .=> 0.0], (0, 20.0)) nc3sys = let y = y n = 10 @@ -421,7 +421,7 @@ nc3sys = let y = y ODESystem(nc3eqs, t, name = :nc3) end -nc3prob = ODEProblem{false}(nc3sys, [y[1] => 1.0; y[2:10] .=> 0.0], (0, 20.0)) +nc3prob = ODEProblem{false}(structural_simplify(nc3sys), [y[1] => 1.0; y[2:10] .=> 0.0], (0, 20.0)) @variables y(t)[1:51] y = collect(y) @@ -433,7 +433,7 @@ nc4sys = let y = y ODESystem(nc4eqs, t, name = :nc4) end -nc4prob = ODEProblem{false}(nc4sys, [y[1] => 1.0; y[2:51] .=> 0.0], (0, 20.0)) +nc4prob = ODEProblem{false}(structural_simplify(nc4sys), [y[1] => 1.0; y[2:51] .=> 0.0], (0, 20.0)) @variables y(t)[1:3, 1:5] y = collect(y) @@ -487,7 +487,7 @@ end function make_ds(nd1sys, e) y = collect(@nonamespace nd1sys.y) y0 = [y[1] => 1-e; y[2:3] .=> 0.0; y[4] => sqrt((1 + e) / (1 - e))] - ODEProblem{false}(nd1sys, y0, (0, 20.0), [ε => e]) + ODEProblem{false}(structural_simplify(nd1sys), y0, (0, 20.0), [ε => e]) end nd1prob = make_ds(nd1sys, 0.1) nd2prob = make_ds(nd1sys, 0.3) @@ -503,7 +503,7 @@ ne1sys = let end y0 = [y[1] => 0.6713967071418030; y[2] => 0.09540051444747446] -ne1prob = ODEProblem{false}(ne1sys, y0, (0, 20.0)) +ne1prob = ODEProblem{false}(structural_simplify(ne1sys), y0, (0, 20.0)) ne2sys = let ne2eqs = [D(y[1]) ~ y[2], @@ -513,7 +513,7 @@ ne2sys = let end y0 = [y[1] => 2.0; y[2] => 0.0] -ne2prob = ODEProblem{false}(ne2sys, y0, (0, 20.0)) +ne2prob = ODEProblem{false}(structural_simplify(ne2sys), y0, (0, 20.0)) ne3sys = let ne3eqs = [D(y[1]) ~ y[2], @@ -522,7 +522,7 @@ ne3sys = let ODESystem(ne3eqs, t, name = :ne3) end -ne3prob = ODEProblem{false}(ne3sys, y[1:2] .=> 0, (0, 20.0)) +ne3prob = ODEProblem{false}(structural_simplify(ne3sys), y[1:2] .=> 0, (0, 20.0)) ne4sys = let ne4eqs = [D(y[1]) ~ y[2], @@ -531,7 +531,7 @@ ne4sys = let ODESystem(ne4eqs, t, name = :ne4) end -ne4prob = ODEProblem{false}(ne4sys, [y[1] => 30.0, y[2] => 0.0], (0, 20.0)) +ne4prob = ODEProblem{false}(structural_simplify(ne4sys), [y[1] => 30.0, y[2] => 0.0], (0, 20.0)) ne5sys = let ne5eqs = [D(y[1]) ~ y[2], @@ -539,27 +539,28 @@ ne5sys = let ODESystem(ne5eqs, t, name = :ne5) end -ne5prob = ODEProblem{false}(ne5sys, y[1:2] .=> 0.0, (0, 20.0)) +ne5prob = ODEProblem{false}(structural_simplify(ne5sys), y[1:2] .=> 0.0, (0, 20.0)) + +cond = term(iseven, term(floor, Int, unwrap(t), type = Int), type = Bool) -@inline myifelse(x, y, z) = ifelse(x, y, z) 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)] + D(y[2]) ~ b + ifelse(cond, 1, -1)] ODESystem(nf1eqs, t, name = :nf1) end -nf1prob = ODEProblem{false}(nf1sys, y[1:2] .=> 0.0, (0, 20.0)) +nf1prob = ODEProblem{false}(structural_simplify(nf1sys), y[1:2] .=> 0.0, (0, 20.0)) +#= nf2sys = let - cond = term(iseven, term(floor, Int, unwrap(t), type = Int), type = Bool) - nf2eqs = [D(y[1]) ~ 55 - term(myifelse, cond, 2y[1]/2, y[1]/2, type=Real)] + nf2eqs = [D(y[1]) ~ 55 - ifelse(cond, 2y[1]/2, y[1]/2)] ODESystem(nf2eqs, t, name = :nf2) end -nf2prob = ODEProblem{false}(nf2sys, [y[1] .=> 110.0], (0, 20.0)) +nf2prob = ODEProblem{false}(structural_simplify(nf2sys), [y[1] .=> 110.0], (0, 20.0)) +=# nf3sys = let nf3eqs = [D(y[1]) ~ y[2], @@ -567,14 +568,14 @@ nf3sys = let ODESystem(nf3eqs, t, name = :nf3) end -nf3prob = ODEProblem{false}(nf3sys, y[1:2] .=> 0.0, (0, 20.0)) +nf3prob = ODEProblem{false}(structural_simplify(nf3sys), y[1:2] .=> 0.0, (0, 20.0)) 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 -nf4prob = ODEProblem{false}(nf4sys, [y[1] => 1.0], (0, 20.0)) +nf4prob = ODEProblem{false}(structural_simplify(nf4sys), [y[1] => 1.0], (0, 20.0)) nf5sys = let c = sum(i->cbrt(i)^4, 1:19) @@ -583,4 +584,4 @@ nf5sys = let ODESystem(nf5eqs, t, name = :nf5) end -nf5prob = ODEProblem{false}(nf5sys, [y[1] => 1.0], (0, 20.0)) +nf5prob = ODEProblem{false}(structural_simplify(nf5sys), [y[1] => 1.0], (0, 20.0)) diff --git a/m.jl b/lib/ODEProblemLibrary/src/m.jl similarity index 99% rename from m.jl rename to lib/ODEProblemLibrary/src/m.jl index 490d0b7..0bbb55d 100644 --- a/m.jl +++ b/lib/ODEProblemLibrary/src/m.jl @@ -74,7 +74,6 @@ DQ = zeros(Num, 4, 4) # for linear stiffness term and no damping # (ipar(1) = 0, ipar(2) = 0). @variables( - t, ϕ₁(t)=0.0, ϕ₂(t)=0.0, x₃(t)=4.500169330000000e-1, @@ -179,7 +178,6 @@ sinp2 = sin(ϕ₂) cosp12 = cos(ϕ₁-ϕ₂) sinp12 = sin(ϕ₁-ϕ₂) -D = Differential(t) P = [ϕ₁, ϕ₂, x₃] Q = [q₁, q₂, q₃, q₄] X = [P..., Q...] @@ -384,4 +382,6 @@ for I = 1:7 push!(eqs, D(X[I]) ~ DX[I]) push!(eqs, D(DX[I]) ~ DDX[I]) end -eqs + +@mtkbuild sys = ODESystem(eqs, t) +someprob = ODEProblem(sys, [], (0.0,1.0)) \ No newline at end of file