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/lib/ODEProblemLibrary/src/enright_pryce.jl b/lib/ODEProblemLibrary/src/enright_pryce.jl new file mode 100644 index 0000000..a031ae2 --- /dev/null +++ b/lib/ODEProblemLibrary/src/enright_pryce.jl @@ -0,0 +1,587 @@ +using ModelingToolkit +using ModelingToolkit: t_nounits as t, D_nounits as D +using ModelingToolkit.Symbolics +using Symbolics: unwrap +using DiffEqBase, StaticArrays, LinearAlgebra + +@variables y(t)[1:10] +y = collect(y) + +# 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}(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]] + 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}(structural_simplify(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}(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}(structural_simplify(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}(structural_simplify(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}(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 + 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}(structural_simplify(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}(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]), + 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}(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]), + 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}(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]), + 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}(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]), + 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}(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]), + D(y[2]) ~ 0.01 - (1 + y[2]^2) * (0.01 + y[1] + y[2])] + + ODESystem(sd5eqs, t, name = :sd5) +end + +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]), + 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}(structural_simplify(sd6sys), [y[1] => 1.0; y[2:3] .=> 0.0], (0, 1.0), dt = 3.3e-8) + +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}(structural_simplify(se1sys), y .=> 0.0, (0, 1.0), dt = 6.8e-3) + +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}(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], + 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}(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] + 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}(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 = [ + 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}(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]) + 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}(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 = [ + 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}(structural_simplify(sf2sys), [y[1] => 1.0; y[2] => 0.0], (0, 240.0), dt = 1e-2) + +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}(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 = [ + 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}(structural_simplify(sf4sys), [y[1] => 4.0; y[2] => 1.1; y[3] => 4.0], (0, 300.0), dt = 1e-3) + +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}(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] + na1eqs = D(y) ~ -y + + ODESystem(na1eqs, t, name = :na1) +end + +na1prob = ODEProblem{false}(structural_simplify(na1sys), [y[1] => 1], (0, 20.0)) + +na2sys = let y = y[1] + na2eqs = D(y) ~ -y^3/2 + + ODESystem(na2eqs, t, name = :na2) +end + +na2prob = ODEProblem{false}(structural_simplify(na2sys), [y[1] => 1], (0, 20.0)) + +na3sys = let y = y[1] + na3eqs = D(y) ~ y * cos(t) + + ODESystem(na3eqs, t, name = :na3) +end + +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) + + ODESystem(na4eqs, t, name = :na4) +end + +na4prob = ODEProblem{false}(structural_simplify(na4sys), [y[1] => 1], (0, 20.0)) + +na5sys = let y = y[1] + na5eqs = D(y) ~ (y - t) / (y + t) + + ODESystem(na5eqs, t, name = :na5) +end + +na5prob = ODEProblem{false}(structural_simplify(na5sys), [y[1] => 4], (0, 20.0)) + +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}(structural_simplify(nb1sys), [y[1] => 1.0, y[2] => 3], (0, 20.0)) + +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}(structural_simplify(nb2sys), [y[1] => 2.0, y[2] => 0.0, y[3] => 1.0], (0, 20.0)) + +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}(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) + 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}(structural_simplify(nb4sys), [y[1] => 3.0; y[2:3] .=> 0.0], (0, 20.0)) + +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}(structural_simplify(nb5sys), [y[1] => 0.0; y[2:3] .=> 1.0], (0, 20.0)) + +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}(structural_simplify(nc1sys), [y[1] => 1.0; y[2:10] .=> 0.0], (0, 20.0)) + +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}(structural_simplify(nc2sys), [y[1] => 1.0; y[2:10] .=> 0.0], (0, 20.0)) + +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}(structural_simplify(nc3sys), [y[1] => 1.0; y[2:10] .=> 0.0], (0, 20.0)) + +@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}(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) +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, 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) +# The orginal paper has t_f = 20, but 1000 looks way better +nc5prob = ODEProblem{false}(nc5sys, [y0; y0′], (0, 20.0)) + +@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}(structural_simplify(nd1sys), y0, (0, 20.0), [ε => e]) +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}(structural_simplify(ne1sys), y0, (0, 20.0)) + +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}(structural_simplify(ne2sys), y0, (0, 20.0)) + +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}(structural_simplify(ne3sys), y[1:2] .=> 0, (0, 20.0)) + +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}(structural_simplify(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),] + ODESystem(ne5eqs, t, name = :ne5) +end + +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) + +nf1sys = let + a = 0.1 + b = 2a * y[2] - (pi^2 + a^2) * y[1] + nf1eqs = [D(y[1]) ~ y[2], + D(y[2]) ~ b + ifelse(cond, 1, -1)] + ODESystem(nf1eqs, t, name = :nf1) +end + +nf1prob = ODEProblem{false}(structural_simplify(nf1sys), y[1:2] .=> 0.0, (0, 20.0)) + +#= +nf2sys = let + nf2eqs = [D(y[1]) ~ 55 - ifelse(cond, 2y[1]/2, y[1]/2)] + ODESystem(nf2eqs, t, name = :nf2) +end + +nf2prob = ODEProblem{false}(structural_simplify(nf2sys), [y[1] .=> 110.0], (0, 20.0)) +=# + +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}(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}(structural_simplify(nf4sys), [y[1] => 1.0], (0, 20.0)) + +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}(structural_simplify(nf5sys), [y[1] => 1.0], (0, 20.0)) diff --git a/lib/ODEProblemLibrary/src/m.jl b/lib/ODEProblemLibrary/src/m.jl new file mode 100644 index 0000000..0bbb55d --- /dev/null +++ b/lib/ODEProblemLibrary/src/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)=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(ϕ₁-ϕ₂) + +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 + +@mtkbuild sys = ODESystem(eqs, t) +someprob = ODEProblem(sys, [], (0.0,1.0)) \ No newline at end of file