From 4a7bd476072b9cb7b9c6423e63211ead2a093975 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Sun, 4 May 2025 00:25:33 -0400 Subject: [PATCH] Properly use symbolic constructors This was using an old form that didn't necessarily always get the order right. This PR fixes the constructors to ensure they are safe for any equation reordering, and thus fixes downstream. --- lib/ODEProblemLibrary/Project.toml | 2 +- .../src/ode_simple_nonlinear_prob.jl | 42 ++++++++----------- 2 files changed, 19 insertions(+), 25 deletions(-) diff --git a/lib/ODEProblemLibrary/Project.toml b/lib/ODEProblemLibrary/Project.toml index c341ff7..0555930 100644 --- a/lib/ODEProblemLibrary/Project.toml +++ b/lib/ODEProblemLibrary/Project.toml @@ -1,6 +1,6 @@ name = "ODEProblemLibrary" uuid = "fdc4e326-1af4-4b90-96e7-779fcce2daa5" -version = "0.1.8" +version = "0.1.9" [deps] DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e" diff --git a/lib/ODEProblemLibrary/src/ode_simple_nonlinear_prob.jl b/lib/ODEProblemLibrary/src/ode_simple_nonlinear_prob.jl index 177b4d8..4dd20c5 100644 --- a/lib/ODEProblemLibrary/src/ode_simple_nonlinear_prob.jl +++ b/lib/ODEProblemLibrary/src/ode_simple_nonlinear_prob.jl @@ -54,8 +54,7 @@ prob_ode_fitzhughnagumo = ODEProblem(fitz, [1.0; 1.0], (0.0, 1.0), eqs = [D(y) ~ μ * ((1 - x^2) * y - x), D(x) ~ y] -de = ODESystem(eqs, t; name = :van_der_pol) |> structural_simplify |> complete -van = ODEFunction(de, [y, x], [μ], jac = true, eval_module = @__MODULE__) +van = ODESystem(eqs, t; name = :van_der_pol) |> structural_simplify |> complete """ Van der Pol Equations @@ -67,11 +66,11 @@ Van der Pol Equations \\frac{dy}{dt} = μ((1-x^2)y -x) ``` -with ``μ=1.0`` and ``u_0=[0,\\sqrt{3}]`` +with ``μ=1.0`` and ``u_0=[x => \\sqrt{3}, y => 0]`` Non-stiff parameters. """ -prob_ode_vanderpol = ODEProblem(van, [0; sqrt(3)], (0.0, 1.0), 1.0) +prob_ode_vanderpol = ODEProblem(van, [y => 0, x => sqrt(3)], (0.0, 1.0), [μ => 1.0], jac=true, eval_module = @__MODULE__) """ Van der Pol Equations @@ -83,11 +82,11 @@ Van der Pol Equations \\frac{dy}{dt} = μ((1-x^2)y -x) ``` -with ``μ=10^6`` and ``u_0=[0,\\sqrt{3}]`` +with ``μ=10^6`` and ``u_0=[x => \\sqrt{3}, y => 0]`` Stiff parameters. """ -prob_ode_vanderpol_stiff = ODEProblem(van, [0; sqrt(3)], (0.0, 1.0), 1e6) +prob_ode_vanderpol_stiff = ODEProblem(van, [y => 0, x => sqrt(3)], (0.0, 1.0), [μ => 1e6], jac=true, eval_module = @__MODULE__) # ROBER @parameters k₁ k₂ k₃ @@ -96,8 +95,7 @@ prob_ode_vanderpol_stiff = ODEProblem(van, [0; sqrt(3)], (0.0, 1.0), 1e6) eqs = [D(y₁) ~ -k₁ * y₁ + k₃ * y₂ * y₃, D(y₂) ~ k₁ * y₁ - k₂ * y₂^2 - k₃ * y₂ * y₃, D(y₃) ~ k₂ * y₂^2] -de = ODESystem(eqs, t; name = :rober) |> structural_simplify |> complete -rober = ODEFunction(de, [y₁, y₂, y₃], [k₁, k₂, k₃], jac = true, eval_module = @__MODULE__) +rober = ODESystem(eqs, t; name = :rober) |> structural_simplify |> complete """ The Robertson biochemical reactions: (Stiff) @@ -118,7 +116,7 @@ Hairer Norsett Wanner Solving Ordinary Differential Equations I - Nonstiff Probl Usually solved on ``[0,1e11]`` """ -prob_ode_rober = ODEProblem(rober, [1.0; 0.0; 0.0], (0.0, 1e11), (0.04, 3e7, 1e4)) +prob_ode_rober = ODEProblem(rober, [y₁, y₂, y₃] .=> [1.0; 0.0; 0.0], (0.0, 1e11), [k₁, k₂, k₃] .=> (0.04, 3e7, 1e4), jac = true, eval_module = @__MODULE__) # Three Body const threebody_μ = big(0.012277471); @@ -177,8 +175,7 @@ prob_ode_threebody = ODEProblem(threebody, eqs = [D(y₁) ~ I₁ * y₂ * y₃, D(y₂) ~ I₂ * y₁ * y₃, D(y₃) ~ I₃ * y₁ * y₂] -de = ODESystem(eqs, t; name = :rigid_body) |> structural_simplify |> complete -rigid = ODEFunction(de, [y₁, y₂, y₃], [I₁, I₂, I₃], jac = true, eval_module = @__MODULE__) +rigid = ODESystem(eqs, t; name = :rigid_body) |> structural_simplify |> complete """ Rigid Body Equations (Non-stiff) @@ -203,7 +200,9 @@ or Hairer Norsett Wanner Solving Ordinary Differential Equations I - Nonstiff Pr Usually solved from 0 to 20. """ -prob_ode_rigidbody = ODEProblem(rigid, [1.0, 0.0, 0.9], (0.0, 20.0), (-2.0, 1.25, -0.5)) +prob_ode_rigidbody = ODEProblem(rigid, [y₁, y₂, y₃] .=> [1.0, 0.0, 0.9], (0.0, 20.0), + [I₁, I₂, I₃] .=> (-2.0, 1.25, -0.5), jac = true, + eval_module = @__MODULE__) # Pleiades Problem @@ -360,10 +359,7 @@ eqs = [D(y1) ~ -p1 * y1 + p2 * y2 + p3 * y3 + p4, p2 * y6 + p11 * y7, D(y7) ~ p10 * y6 * y8 - p12 * y7, D(y8) ~ -p10 * y6 * y8 + p12 * y7] -de = ODESystem(eqs, t; name = :hires) |> structural_simplify |> complete -hires = ODEFunction(de, [y1, y2, y3, y4, y5, y6, y7, y8], - [p1, p2, p3, p4, p5, p6, p7, p8, p9, p10, p11, p12], - jac = true) +hires = ODESystem(eqs, t; name = :hires) |> structural_simplify |> complete u0 = zeros(8) u0[1] = 1 @@ -391,10 +387,10 @@ where ``f`` is defined by Reference: [demohires.pdf](http://www.radford.edu/~thompson/vodef90web/problems/demosnodislin/Demos_Pitagora/DemoHires/demohires.pdf) Notebook: [Hires.ipynb](http://nbviewer.jupyter.org/github/JuliaDiffEq/DiffEqBenchmarks.jl/blob/master/StiffODE/Hires.ipynb) """ -prob_ode_hires = ODEProblem(hires, u0, (0.0, 321.8122), - (1.71, 0.43, 8.32, 0.0007, 8.75, - 10.03, 0.035, 1.12, 1.745, 280.0, - 0.69, 1.81)) +prob_ode_hires = ODEProblem(hires, [y1, y2, y3, y4, y5, y6, y7, y8] .=> u0, (0.0, 321.8122), + [p1, p2, p3, p4, p5, p6, p7, p8, p9, p10, p11, p12] .=> (1.71, 0.43, 8.32, 0.0007, 8.75, + 10.03, 0.035, 1.12, 1.745, 280.0, + 0.69, 1.81), jac = true, eval_module = @__MODULE__) @parameters p1 p2 p3 @variables y1(t) y2(t) y3(t) @@ -402,9 +398,7 @@ prob_ode_hires = ODEProblem(hires, u0, (0.0, 321.8122), eqs = [D(y1) ~ p1 * (y2 + y1 * (1 - p2 * y1 - y2)), D(y2) ~ (y3 - (1 + y1) * y2) / p1, D(y3) ~ p3 * (y1 - y3)] -de = ODESystem(eqs, t; name = :orego) |> structural_simplify |> complete -jac = calculate_jacobian(de) -orego = ODEFunction(de, [y1, y2, y3], [p1, p2, p3], jac = true, eval_module = @__MODULE__) +orego = ODESystem(eqs, t; name = :orego) |> structural_simplify |> complete """ Orego Problem (Stiff) @@ -424,4 +418,4 @@ where ``s=77.27``, ``w=0.161`` and ``q=8.375⋅10^{-6}``. Reference: [demoorego.pdf](http://www.radford.edu/~thompson/vodef90web/problems/demosnodislin/Demos_Pitagora/DemoOrego/demoorego.pdf) Notebook: [Orego.ipynb](http://nbviewer.jupyter.org/github/JuliaDiffEq/DiffEqBenchmarks.jl/blob/master/StiffODE/Orego.ipynb) """ -prob_ode_orego = ODEProblem(orego, [1.0, 2.0, 3.0], (0.0, 30.0), [77.27, 8.375e-6, 0.161]) +prob_ode_orego = ODEProblem(orego, [y1, y2, y3] .=> [1.0, 2.0, 3.0], (0.0, 30.0), [p1, p2, p3] .=> [77.27, 8.375e-6, 0.161], jac = true, eval_module = @__MODULE__)