Add up to NC5

Author: YingboMa

src/ode/enright_pryce.jl

@@ -1,10 +1,11 @@
 using ModelingToolkit
-using DiffEqBase, StaticArrays
+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)
@@ -162,3 +163,307 @@ sd6sys = let
 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)