From f7a365b561e25652aab0f75fff28cbc96c61fe63 Mon Sep 17 00:00:00 2001 From: Wen Wei Tseng Date: Thu, 28 Dec 2023 15:25:38 +0800 Subject: [PATCH] fix --- Manifest.toml | 175 ++++++++++++++++------------------------ docs/figs/ch1-07.jl | 65 ++++++++------- docs/figs/ch1-09.jl | 94 ++++++++++++--------- docs/figs/ch4-01.jl | 8 +- docs/figs/ch4-15.jl | 6 +- docs/figs/ch7-07.jl | 71 ++++++++++------ docs/intro-03-diffeq.jl | 146 +++++++++++++++++---------------- 7 files changed, 293 insertions(+), 272 deletions(-) diff --git a/Manifest.toml b/Manifest.toml index 55b79524..d60d7f2e 100644 --- a/Manifest.toml +++ b/Manifest.toml @@ -20,25 +20,6 @@ git-tree-sha1 = "faa260e4cb5aba097a73fab382dd4b5819d8ec8c" uuid = "1520ce14-60c1-5f80-bbc7-55ef81b5835c" version = "0.4.4" -[[deps.Accessors]] -deps = ["CompositionsBase", "ConstructionBase", "Dates", "InverseFunctions", "LinearAlgebra", "MacroTools", "Test"] -git-tree-sha1 = "a7055b939deae2455aa8a67491e034f735dd08d3" -uuid = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697" -version = "0.1.33" - - [deps.Accessors.extensions] - AccessorsAxisKeysExt = "AxisKeys" - AccessorsIntervalSetsExt = "IntervalSets" - AccessorsStaticArraysExt = "StaticArrays" - AccessorsStructArraysExt = "StructArrays" - - [deps.Accessors.weakdeps] - AxisKeys = "94b1ba4f-4ee9-5380-92f1-94cde586c3c5" - IntervalSets = "8197267c-284f-5f27-9208-e0e47529a953" - Requires = "ae029012-a4dd-5104-9daa-d747884805df" - StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" - StructArrays = "09ab397b-f2b6-538f-b94a-2f83cf4a842a" - [[deps.Adapt]] deps = ["LinearAlgebra", "Requires"] git-tree-sha1 = "cde29ddf7e5726c9fb511f340244ea3481267608" @@ -288,15 +269,6 @@ git-tree-sha1 = "02d2316b7ffceff992f3096ae48c7829a8aa0638" uuid = "b152e2b5-7a66-4b01-a709-34e65c35f657" version = "0.1.3" -[[deps.CompositionsBase]] -git-tree-sha1 = "802bb88cd69dfd1509f6670416bd4434015693ad" -uuid = "a33af91c-f02d-484b-be07-31d278c5ca2b" -version = "0.1.2" -weakdeps = ["InverseFunctions"] - - [deps.CompositionsBase.extensions] - CompositionsBaseInverseFunctionsExt = "InverseFunctions" - [[deps.ConcreteStructs]] git-tree-sha1 = "f749037478283d372048690eb3b5f92a79432b34" uuid = "2569d6c7-a4a2-43d3-a901-331e8e4be471" @@ -363,9 +335,9 @@ uuid = "ade2ca70-3891-5945-98fb-dc099432e06a" [[deps.DelayDiffEq]] deps = ["ArrayInterface", "DataStructures", "DiffEqBase", "LinearAlgebra", "Logging", "OrdinaryDiffEq", "Printf", "RecursiveArrayTools", "Reexport", "SciMLBase", "SimpleNonlinearSolve", "SimpleUnPack"] -git-tree-sha1 = "e40378efd2af7658d0a0579aa9e15b17137368f4" +git-tree-sha1 = "9534eac1eb01f9fbc9319799cabfdb6ec9c56e42" uuid = "bcd4f6db-9728-5f36-b5f7-82caef46ccdb" -version = "5.44.0" +version = "5.45.0" [[deps.DelimitedFiles]] deps = ["Mmap"] @@ -375,9 +347,9 @@ version = "1.9.1" [[deps.DiffEqBase]] deps = ["ArrayInterface", "DataStructures", "DocStringExtensions", "EnumX", "EnzymeCore", "FastBroadcast", "ForwardDiff", "FunctionWrappers", "FunctionWrappersWrappers", "LinearAlgebra", "Logging", "Markdown", "MuladdMacro", "Parameters", "PreallocationTools", "PrecompileTools", "Printf", "RecursiveArrayTools", "Reexport", "SciMLBase", "SciMLOperators", "Setfield", "SparseArrays", "Static", "StaticArraysCore", "Statistics", "Tricks", "TruncatedStacktraces"] -git-tree-sha1 = "8775b80752e9656000ab3800adad8ee22c9cb8f6" +git-tree-sha1 = "f45adda33e3eb24ae1428dfa9b57ded8248ced99" uuid = "2b5f629d-d688-5b77-993f-72d75c75574e" -version = "6.145.0" +version = "6.145.4" [deps.DiffEqBase.extensions] DiffEqBaseChainRulesCoreExt = "ChainRulesCore" @@ -405,9 +377,9 @@ version = "6.145.0" [[deps.DiffEqCallbacks]] deps = ["DataStructures", "DiffEqBase", "ForwardDiff", "Functors", "LinearAlgebra", "Markdown", "NLsolve", "Parameters", "RecipesBase", "RecursiveArrayTools", "SciMLBase", "StaticArraysCore"] -git-tree-sha1 = "d0b94b3694d55e7eedeee918e7daee9e3b873399" +git-tree-sha1 = "cf334da651a6e42c50e1477d6ab978f1b8be3057" uuid = "459566f4-90b8-5000-8ac3-15dfb0a30def" -version = "2.35.0" +version = "2.36.1" weakdeps = ["OrdinaryDiffEq", "Sundials"] [[deps.DiffEqNoiseProcess]] @@ -436,9 +408,9 @@ version = "1.15.1" [[deps.DifferentialEquations]] deps = ["BoundaryValueDiffEq", "DelayDiffEq", "DiffEqBase", "DiffEqCallbacks", "DiffEqNoiseProcess", "JumpProcesses", "LinearAlgebra", "LinearSolve", "NonlinearSolve", "OrdinaryDiffEq", "Random", "RecursiveArrayTools", "Reexport", "SciMLBase", "SteadyStateDiffEq", "StochasticDiffEq", "Sundials"] -git-tree-sha1 = "19a5b6314715139ddefea4108a105bb9b90dc4fb" +git-tree-sha1 = "8864b6a953eeba7890d23258aca468d90ca73fd6" uuid = "0c46a032-eb83-5123-abaf-570d42b7fbaa" -version = "7.11.0" +version = "7.12.0" [[deps.Distances]] deps = ["LinearAlgebra", "Statistics", "StatsAPI"] @@ -824,12 +796,6 @@ weakdeps = ["Statistics"] [deps.IntervalSets.extensions] IntervalSetsStatisticsExt = "Statistics" -[[deps.InverseFunctions]] -deps = ["Test"] -git-tree-sha1 = "68772f49f54b479fa88ace904f6127f0a3bb2e46" -uuid = "3587e190-3f89-42d0-90ee-14403ec27112" -version = "0.1.12" - [[deps.InvertedIndices]] git-tree-sha1 = "0dc7b50b8d436461be01300fd8cd45aa0274b038" uuid = "41ab1584-1d38-5bbf-9106-f11c6c58b48f" @@ -951,12 +917,6 @@ version = "0.16.1" DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" SymEngine = "123dc426-2d89-5057-bbad-38513e3affd8" -[[deps.LatticeRules]] -deps = ["Random"] -git-tree-sha1 = "7f5b02258a3ca0221a6a9710b0a0a2e8fb4957fe" -uuid = "73f95e8e-ec14-4e6a-8b18-0d2e271c4e55" -version = "0.0.1" - [[deps.LayoutPointers]] deps = ["ArrayInterface", "LinearAlgebra", "ManualMemory", "SIMDTypes", "Static", "StaticArrayInterface"] git-tree-sha1 = "62edfee3211981241b57ff1cedf4d74d79519277" @@ -1177,6 +1137,12 @@ git-tree-sha1 = "78f6e33434939b0ac9ba1df81e6d005ee85a7396" uuid = "a3b82374-2e81-5b9e-98ce-41277c0e4c87" version = "2.1.0" +[[deps.MaybeInplace]] +deps = ["ArrayInterface", "LinearAlgebra", "MacroTools", "SparseArrays"] +git-tree-sha1 = "a85c6a98c9e5a2a7046bc1bb89f28a3241e1de4d" +uuid = "bb5d69b7-63fc-4a16-80bd-7e42200c7bdb" +version = "0.1.1" + [[deps.MbedTLS]] deps = ["Dates", "MbedTLS_jll", "MozillaCACerts_jll", "NetworkOptions", "Random", "Sockets"] git-tree-sha1 = "c067a280ddc25f196b5e7df3877c6b226d390aaf" @@ -1204,9 +1170,9 @@ uuid = "a63ad114-7e13-5084-954f-fe012c677804" [[deps.ModelingToolkit]] deps = ["AbstractTrees", "ArrayInterface", "Combinatorics", "Compat", "ConstructionBase", "DataStructures", "DiffEqBase", "DiffEqCallbacks", "DiffRules", "Distributed", "Distributions", "DocStringExtensions", "DomainSets", "ForwardDiff", "FunctionWrappersWrappers", "Graphs", "IfElse", "InteractiveUtils", "JuliaFormatter", "JumpProcesses", "LabelledArrays", "Latexify", "Libdl", "LinearAlgebra", "MLStyle", "MacroTools", "NaNMath", "OrdinaryDiffEq", "PrecompileTools", "RecursiveArrayTools", "Reexport", "RuntimeGeneratedFunctions", "SciMLBase", "Serialization", "Setfield", "SimpleNonlinearSolve", "SparseArrays", "SpecialFunctions", "StaticArrays", "SymbolicIndexingInterface", "SymbolicUtils", "Symbolics", "URIs", "UnPack", "Unitful"] -git-tree-sha1 = "2f872327288342d1e3edc538484b4f317e3cf59c" +git-tree-sha1 = "944c2d44b08b08035569043b27c6cd9a2d2476fc" uuid = "961ee093-0014-501f-94e3-6117800e7a78" -version = "8.73.2" +version = "8.74.1" [deps.ModelingToolkit.extensions] MTKBifurcationKitExt = "BifurcationKit" @@ -1260,20 +1226,34 @@ uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908" version = "1.2.0" [[deps.NonlinearSolve]] -deps = ["ADTypes", "ArrayInterface", "ConcreteStructs", "DiffEqBase", "EnumX", "FastBroadcast", "FiniteDiff", "ForwardDiff", "LineSearches", "LinearAlgebra", "LinearSolve", "PrecompileTools", "RecursiveArrayTools", "Reexport", "SciMLBase", "SciMLOperators", "SimpleNonlinearSolve", "SparseArrays", "SparseDiffTools", "StaticArraysCore", "UnPack"] -git-tree-sha1 = "6166ccd8f79c93c636ca61ab4cd18f555932563d" +deps = ["ADTypes", "ArrayInterface", "ConcreteStructs", "DiffEqBase", "EnumX", "FastBroadcast", "FastClosures", "FiniteDiff", "ForwardDiff", "LazyArrays", "LineSearches", "LinearAlgebra", "LinearSolve", "MaybeInplace", "PrecompileTools", "Printf", "RecursiveArrayTools", "Reexport", "SciMLBase", "SciMLOperators", "SimpleNonlinearSolve", "SparseArrays", "SparseDiffTools", "StaticArrays", "UnPack"] +git-tree-sha1 = "72b036b728461272ae1b1c3f7096cb4c319d8793" uuid = "8913a72c-1f9b-4ce2-8d82-65094dcecaec" -version = "2.8.2" +version = "3.4.0" [deps.NonlinearSolve.extensions] NonlinearSolveBandedMatricesExt = "BandedMatrices" NonlinearSolveFastLevenbergMarquardtExt = "FastLevenbergMarquardt" + NonlinearSolveFixedPointAccelerationExt = "FixedPointAcceleration" NonlinearSolveLeastSquaresOptimExt = "LeastSquaresOptim" + NonlinearSolveMINPACKExt = "MINPACK" + NonlinearSolveNLsolveExt = "NLsolve" + NonlinearSolveSIAMFANLEquationsExt = "SIAMFANLEquations" + NonlinearSolveSpeedMappingExt = "SpeedMapping" + NonlinearSolveSymbolicsExt = "Symbolics" + NonlinearSolveZygoteExt = "Zygote" [deps.NonlinearSolve.weakdeps] BandedMatrices = "aae01518-5342-5314-be14-df237901396f" FastLevenbergMarquardt = "7a0df574-e128-4d35-8cbd-3d84502bf7ce" + FixedPointAcceleration = "817d07cb-a79a-5c30-9a31-890123675176" LeastSquaresOptim = "0fc2ff8b-aaa3-5acd-a817-1944a5e08891" + MINPACK = "4854310b-de5a-5eb6-a2a5-c1dee2bd17f9" + NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56" + SIAMFANLEquations = "084e46ad-d928-497d-ad5e-07fa361a48c4" + SpeedMapping = "f1835b91-879b-4a3f-a438-e4baacf14412" + Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" + Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f" [[deps.OffsetArrays]] git-tree-sha1 = "6a731f2b5c03157418a20c12195eb4b74c8f8621" @@ -1336,10 +1316,10 @@ uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d" version = "1.6.3" [[deps.OrdinaryDiffEq]] -deps = ["ADTypes", "Adapt", "ArrayInterface", "DataStructures", "DiffEqBase", "DocStringExtensions", "ExponentialUtilities", "FastBroadcast", "FastClosures", "FiniteDiff", "ForwardDiff", "FunctionWrappersWrappers", "IfElse", "InteractiveUtils", "LineSearches", "LinearAlgebra", "LinearSolve", "Logging", "LoopVectorization", "MacroTools", "MuladdMacro", "NLsolve", "NonlinearSolve", "Polyester", "PreallocationTools", "PrecompileTools", "Preferences", "RecursiveArrayTools", "Reexport", "SciMLBase", "SciMLNLSolve", "SciMLOperators", "SimpleNonlinearSolve", "SimpleUnPack", "SparseArrays", "SparseDiffTools", "StaticArrayInterface", "StaticArrays", "TruncatedStacktraces"] -git-tree-sha1 = "5f9e7ce227d0e447c3803cc05ef5d8f75f84b9ea" +deps = ["ADTypes", "Adapt", "ArrayInterface", "DataStructures", "DiffEqBase", "DocStringExtensions", "ExponentialUtilities", "FastBroadcast", "FastClosures", "FillArrays", "FiniteDiff", "ForwardDiff", "FunctionWrappersWrappers", "IfElse", "InteractiveUtils", "LineSearches", "LinearAlgebra", "LinearSolve", "Logging", "LoopVectorization", "MacroTools", "MuladdMacro", "NonlinearSolve", "Polyester", "PreallocationTools", "PrecompileTools", "Preferences", "RecursiveArrayTools", "Reexport", "SciMLBase", "SciMLOperators", "SimpleNonlinearSolve", "SimpleUnPack", "SparseArrays", "SparseDiffTools", "StaticArrayInterface", "StaticArrays", "TruncatedStacktraces"] +git-tree-sha1 = "7857772d839bf0773f5db9e456149c832ba5facb" uuid = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" -version = "6.59.3" +version = "6.65.0" [[deps.PCRE2_jll]] deps = ["Artifacts", "Libdl"] @@ -1500,16 +1480,6 @@ git-tree-sha1 = "9ebcd48c498668c7fa0e97a9cae873fbee7bfee1" uuid = "1fd47b50-473d-5c70-9696-f719f8f3bcdc" version = "2.9.1" -[[deps.QuasiMonteCarlo]] -deps = ["Accessors", "ConcreteStructs", "LatticeRules", "LinearAlgebra", "Primes", "Random", "Requires", "Sobol", "StatsBase"] -git-tree-sha1 = "cc086f8485bce77b6187141e1413c3b55f9a4341" -uuid = "8a4e6c94-4038-4cdc-81c3-7e6ffdb2a71b" -version = "0.3.3" -weakdeps = ["Distributions"] - - [deps.QuasiMonteCarlo.extensions] - QuasiMonteCarloDistributionsExt = "Distributions" - [[deps.REPL]] deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"] uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb" @@ -1559,18 +1529,20 @@ uuid = "01d81517-befc-4cb6-b9ec-a95719d0359c" version = "0.6.12" [[deps.RecursiveArrayTools]] -deps = ["Adapt", "ArrayInterface", "DocStringExtensions", "GPUArraysCore", "IteratorInterfaceExtensions", "LinearAlgebra", "RecipesBase", "Requires", "StaticArraysCore", "Statistics", "SymbolicIndexingInterface", "Tables"] -git-tree-sha1 = "d7087c013e8a496ff396bae843b1e16d9a30ede8" +deps = ["Adapt", "ArrayInterface", "DocStringExtensions", "GPUArraysCore", "IteratorInterfaceExtensions", "LinearAlgebra", "RecipesBase", "Requires", "SparseArrays", "StaticArraysCore", "Statistics", "SymbolicIndexingInterface", "Tables"] +git-tree-sha1 = "17bc8c517dd0ced987f7c6c89baabf798b74d9f0" uuid = "731186ca-8d62-57ce-b412-fbd966d074cd" -version = "2.38.10" +version = "3.3.0" [deps.RecursiveArrayTools.extensions] + RecursiveArrayToolsFastBroadcastExt = "FastBroadcast" RecursiveArrayToolsMeasurementsExt = "Measurements" RecursiveArrayToolsMonteCarloMeasurementsExt = "MonteCarloMeasurements" RecursiveArrayToolsTrackerExt = "Tracker" RecursiveArrayToolsZygoteExt = "Zygote" [deps.RecursiveArrayTools.weakdeps] + FastBroadcast = "7034ab61-46d4-4ed7-9d0f-46aef9175898" Measurements = "eff96d63-e80a-5855-80a2-b1b0885c5ab7" MonteCarloMeasurements = "0987c9cc-fe09-11e8-30f0-b96dd679fdca" Tracker = "9f7883ad-71c0-57eb-9f7f-b5c9e6d3789c" @@ -1645,10 +1617,10 @@ uuid = "476501e8-09a2-5ece-8869-fb82de89a1fa" version = "0.6.42" [[deps.SciMLBase]] -deps = ["ADTypes", "ArrayInterface", "CommonSolve", "ConstructionBase", "Distributed", "DocStringExtensions", "EnumX", "FillArrays", "FunctionWrappersWrappers", "IteratorInterfaceExtensions", "LinearAlgebra", "Logging", "Markdown", "PrecompileTools", "Preferences", "Printf", "QuasiMonteCarlo", "RecipesBase", "RecursiveArrayTools", "Reexport", "RuntimeGeneratedFunctions", "SciMLOperators", "StaticArraysCore", "Statistics", "SymbolicIndexingInterface", "Tables", "TruncatedStacktraces"] -git-tree-sha1 = "32ea825941f7b58a6f48268f4b76971ae8eb9eec" +deps = ["ADTypes", "ArrayInterface", "CommonSolve", "ConstructionBase", "Distributed", "DocStringExtensions", "EnumX", "FillArrays", "FunctionWrappersWrappers", "IteratorInterfaceExtensions", "LinearAlgebra", "Logging", "Markdown", "PrecompileTools", "Preferences", "Printf", "RecipesBase", "RecursiveArrayTools", "Reexport", "RuntimeGeneratedFunctions", "SciMLOperators", "StaticArraysCore", "Statistics", "SymbolicIndexingInterface", "Tables", "TruncatedStacktraces"] +git-tree-sha1 = "f46e4dfdf24b07db0c4476a0fe730f9bef5e136a" uuid = "0bca4576-84f4-4d90-8ffe-ffa030f20462" -version = "2.10.0" +version = "2.12.1" [deps.SciMLBase.extensions] SciMLBaseChainRulesCoreExt = "ChainRulesCore" @@ -1667,12 +1639,6 @@ version = "2.10.0" RCall = "6f49c342-dc21-5d91-9882-a32aef131414" Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f" -[[deps.SciMLNLSolve]] -deps = ["DiffEqBase", "LineSearches", "NLsolve", "Reexport", "SciMLBase"] -git-tree-sha1 = "765b788339abd7d983618c09cfc0192e2b6b15fd" -uuid = "e9a6253c-8580-4d32-9898-8661bb511710" -version = "0.1.9" - [[deps.SciMLOperators]] deps = ["ArrayInterface", "DocStringExtensions", "Lazy", "LinearAlgebra", "Setfield", "SparseArrays", "StaticArraysCore", "Tricks"] git-tree-sha1 = "51ae235ff058a64815e0a2c34b1db7578a06813d" @@ -1716,16 +1682,16 @@ uuid = "777ac1f9-54b0-4bf8-805c-2214025038e7" version = "1.1.0" [[deps.SimpleNonlinearSolve]] -deps = ["ArrayInterface", "DiffEqBase", "FiniteDiff", "ForwardDiff", "LinearAlgebra", "PrecompileTools", "Reexport", "SciMLBase", "StaticArraysCore"] -git-tree-sha1 = "69b1a53374dd14d7c165d98cb646aeb5f36f8d07" +deps = ["ADTypes", "ArrayInterface", "ConcreteStructs", "DiffEqBase", "FiniteDiff", "ForwardDiff", "LinearAlgebra", "MaybeInplace", "PrecompileTools", "Reexport", "SciMLBase", "StaticArraysCore"] +git-tree-sha1 = "8d672bd91dc432fb286b6d4bcf1a5dc417e932a3" uuid = "727e6d20-b764-4bd8-a329-72de5adea6c7" -version = "0.1.25" +version = "1.2.0" [deps.SimpleNonlinearSolve.extensions] - SimpleNonlinearSolveNNlibExt = "NNlib" + SimpleNonlinearSolvePolyesterForwardDiffExt = "PolyesterForwardDiff" [deps.SimpleNonlinearSolve.weakdeps] - NNlib = "872c559c-99b0-510c-b3b7-b6c96a88d5cd" + PolyesterForwardDiff = "98d1487c-24ca-40b6-b7ab-df2af84e126b" [[deps.SimpleTraits]] deps = ["InteractiveUtils", "MacroTools"] @@ -1738,20 +1704,14 @@ git-tree-sha1 = "58e6353e72cde29b90a69527e56df1b5c3d8c437" uuid = "ce78b400-467f-4804-87d8-8f486da07d0a" version = "1.1.0" -[[deps.Sobol]] -deps = ["DelimitedFiles", "Random"] -git-tree-sha1 = "5a74ac22a9daef23705f010f72c81d6925b19df8" -uuid = "ed01d8cd-4d21-5b2a-85b4-cc3bdc58bad4" -version = "1.5.0" - [[deps.Sockets]] uuid = "6462fe0b-24de-5631-8697-dd941f90decc" [[deps.SortingAlgorithms]] deps = ["DataStructures"] -git-tree-sha1 = "5165dfb9fd131cf0c6957a3a7605dede376e7b63" +git-tree-sha1 = "66e0a8e672a0bdfca2c3f5937efb8538b9ddc085" uuid = "a2af1166-a08f-5f64-846c-94a0d3cef48c" -version = "1.2.0" +version = "1.2.1" [[deps.SparseArrays]] deps = ["Libdl", "LinearAlgebra", "Random", "Serialization", "SuiteSparse_jll"] @@ -1845,17 +1805,20 @@ deps = ["HypergeometricFunctions", "IrrationalConstants", "LogExpFunctions", "Re git-tree-sha1 = "f625d686d5a88bcd2b15cd81f18f98186fdc0c9a" uuid = "4c63d2b9-4356-54db-8cca-17b64c39e42c" version = "1.3.0" -weakdeps = ["ChainRulesCore", "InverseFunctions"] [deps.StatsFuns.extensions] StatsFunsChainRulesCoreExt = "ChainRulesCore" StatsFunsInverseFunctionsExt = "InverseFunctions" + [deps.StatsFuns.weakdeps] + ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" + InverseFunctions = "3587e190-3f89-42d0-90ee-14403ec27112" + [[deps.SteadyStateDiffEq]] -deps = ["DiffEqBase", "DiffEqCallbacks", "LinearAlgebra", "NLsolve", "Reexport", "SciMLBase"] -git-tree-sha1 = "2ca69f4be3294e4cd987d83d6019037d420d9fc1" +deps = ["ConcreteStructs", "DiffEqBase", "DiffEqCallbacks", "LinearAlgebra", "Reexport", "SciMLBase"] +git-tree-sha1 = "a735fd5053724cf4de31c81b4e2cc429db844be5" uuid = "9672c7b4-1e72-59bd-8a11-6ac3964bc41f" -version = "1.16.1" +version = "2.0.1" [[deps.StochasticDiffEq]] deps = ["Adapt", "ArrayInterface", "DataStructures", "DiffEqBase", "DiffEqNoiseProcess", "DocStringExtensions", "FiniteDiff", "ForwardDiff", "JumpProcesses", "LevyArea", "LinearAlgebra", "Logging", "MuladdMacro", "NLsolve", "OrdinaryDiffEq", "Random", "RandomNumbers", "RecursiveArrayTools", "Reexport", "SciMLBase", "SciMLOperators", "SparseArrays", "SparseDiffTools", "StaticArrays", "UnPack"] @@ -1897,22 +1860,21 @@ uuid = "fb77eaff-e24c-56d4-86b1-d163f2edb164" version = "5.2.2+0" [[deps.SymbolicIndexingInterface]] -deps = ["DocStringExtensions"] -git-tree-sha1 = "f8ab052bfcbdb9b48fad2c80c873aa0d0344dfe5" +git-tree-sha1 = "be414bfd80c2c91197823890c66ef4b74f5bf5fe" uuid = "2efcf032-c050-4f8e-a9bb-153293bab1f5" -version = "0.2.2" +version = "0.3.1" [[deps.SymbolicUtils]] -deps = ["AbstractTrees", "Bijections", "ChainRulesCore", "Combinatorics", "ConstructionBase", "DataStructures", "DocStringExtensions", "DynamicPolynomials", "IfElse", "LabelledArrays", "LinearAlgebra", "MultivariatePolynomials", "NaNMath", "Setfield", "SparseArrays", "SpecialFunctions", "StaticArrays", "TimerOutputs", "Unityper"] -git-tree-sha1 = "2f3fa844bcd33e40d8c29de5ee8dded7a0a70422" +deps = ["AbstractTrees", "Bijections", "ChainRulesCore", "Combinatorics", "ConstructionBase", "DataStructures", "DocStringExtensions", "DynamicPolynomials", "IfElse", "LabelledArrays", "LinearAlgebra", "MultivariatePolynomials", "NaNMath", "Setfield", "SparseArrays", "SpecialFunctions", "StaticArrays", "SymbolicIndexingInterface", "TimerOutputs", "Unityper"] +git-tree-sha1 = "849b1dfb1680a9e9f2c6023f79a49b694fb6d0da" uuid = "d1185830-fcd6-423d-90d6-eec64667417b" -version = "1.4.0" +version = "1.5.0" [[deps.Symbolics]] -deps = ["ArrayInterface", "Bijections", "ConstructionBase", "DataStructures", "DiffRules", "Distributions", "DocStringExtensions", "DomainSets", "DynamicPolynomials", "Groebner", "IfElse", "LaTeXStrings", "LambertW", "Latexify", "Libdl", "LinearAlgebra", "LogExpFunctions", "MacroTools", "Markdown", "NaNMath", "PrecompileTools", "RecipesBase", "RecursiveArrayTools", "Reexport", "Requires", "RuntimeGeneratedFunctions", "SciMLBase", "Setfield", "SparseArrays", "SpecialFunctions", "StaticArrays", "SymbolicUtils"] -git-tree-sha1 = "28f55dcd865e4a97f262fc476306cee14e8d4651" +deps = ["ArrayInterface", "Bijections", "ConstructionBase", "DataStructures", "DiffRules", "Distributions", "DocStringExtensions", "DomainSets", "DynamicPolynomials", "Groebner", "IfElse", "LaTeXStrings", "LambertW", "Latexify", "Libdl", "LinearAlgebra", "LogExpFunctions", "MacroTools", "Markdown", "NaNMath", "PrecompileTools", "RecipesBase", "Reexport", "Requires", "RuntimeGeneratedFunctions", "SciMLBase", "Setfield", "SparseArrays", "SpecialFunctions", "StaticArrays", "SymbolicIndexingInterface", "SymbolicUtils"] +git-tree-sha1 = "80e67b76699b1ec130ac60946b3d06017f58141f" uuid = "0c5d862f-8b57-4792-8d23-62f2024744c7" -version = "5.11.0" +version = "5.13.0" [deps.Symbolics.extensions] SymbolicsSymPyExt = "SymPy" @@ -2023,12 +1985,15 @@ deps = ["Dates", "LinearAlgebra", "Random"] git-tree-sha1 = "3c793be6df9dd77a0cf49d80984ef9ff996948fa" uuid = "1986cc42-f94f-5a68-af5c-568840ba703d" version = "1.19.0" -weakdeps = ["ConstructionBase", "InverseFunctions"] [deps.Unitful.extensions] ConstructionBaseUnitfulExt = "ConstructionBase" InverseFunctionsUnitfulExt = "InverseFunctions" + [deps.Unitful.weakdeps] + ConstructionBase = "187b0558-2788-49d3-abe0-74a17ed4e7c9" + InverseFunctions = "3587e190-3f89-42d0-90ee-14403ec27112" + [[deps.UnitfulLatexify]] deps = ["LaTeXStrings", "Latexify", "Unitful"] git-tree-sha1 = "e2d817cc500e960fdbafcf988ac8436ba3208bfd" diff --git a/docs/figs/ch1-07.jl b/docs/figs/ch1-07.jl index ace10843..dfd37d19 100644 --- a/docs/figs/ch1-07.jl +++ b/docs/figs/ch1-07.jl @@ -7,8 +7,7 @@ For Figures 1.7, 7.13, 7.14, 7.15 ===# using DifferentialEquations -using LabelledArrays -using SimpleUnPack +using ModelingToolkit using Plots Plots.default(linewidth=2) @@ -20,37 +19,43 @@ hil(x, k) = x / (x + k) hil(x, k, n) = hil(x^n, k^n) exprel(x) = x / expm1(x) -# Collins switch model -function collins!(D, u, p, t) - @unpack a1, a2, β, γ, i1, i2 = p - @unpack s1, s2 = u - D.s1 = a1 * hil(1 + i2, s2, β) - s1 - D.s2 = a2 * hil(1 + i1, s1, γ) - s2 - return nothing +# Define the problem +function build_collins(;name) + @parameters begin + a1=3.0 + a2=2.5 + β=4.0 + γ=4.0 + end + + @variables begin + t + s1(t)=0.075 + s2(t)=2.5 + i1(t) + i2(t) + end + + D = Differential(t) + + eqs = [ + D(s1) ~ a1 * hil(1 + i2, s2, β) - s1 + D(s2) ~ a2 * hil(1 + i1, s1, γ) - s2 + i2 ~ 10 * (10 < t) * (t < 20) + i1 ~ 10 * (30 < t) * (t < 40) + ] + sys = ODESystem(eqs; name) + return structural_simplify(sys) end -# events -on_10!(integrator) = integrator.p.i2 = 10. -on_20!(integrator) = integrator.p.i2 = 0. -on_30!(integrator) = integrator.p.i1 = 10. -on_40!(integrator) = integrator.p.i1 = 0. +#--- +@named sys = build_collins() -events = CallbackSet( - PresetTimeCallback(10., on_10!), - PresetTimeCallback(20., on_20!), - PresetTimeCallback(30., on_30!), - PresetTimeCallback(40., on_40!), -) - -# Prepare -ps = LVector(a1=3.0, a2=2.5, β=4.0, γ=4.0, i1=0.0, i2=0.0) -u0 = LVector(s1=0.075, s2=2.5) -tend = 50.0 - -# solve and visualize -prob = ODEProblem(collins!, u0, tend, ps) -sol = solve(prob, callback=events) +# Solve the problem +tspan = (0.0, 50.0) +prob = ODEProblem(sys, [], tspan, []) +sol = solve(prob, tstops=10:10:40) +# Visual fig = plot(sol, legend=:right, xlabel = "Time", ylabel="Concentration", title="Fig 1.7") - fig |> PNG diff --git a/docs/figs/ch1-09.jl b/docs/figs/ch1-09.jl index 4f683e09..f9add091 100644 --- a/docs/figs/ch1-09.jl +++ b/docs/figs/ch1-09.jl @@ -6,8 +6,8 @@ Hodgkin-Huxley model ===# using DifferentialEquations -using LabelledArrays -using SimpleUnPack +using ModelingToolkit +using IfElse using Plots Plots.default(linewidth=2) @@ -23,49 +23,71 @@ exprel(x) = x / expm1(x) _istim(t) = ifelse(20 <= t <= 21, -6.6, 0.0) + ifelse(60 <= t <= 61, -6.9, 0.0) # HH Neuron model -function hh!(D, u, p, t) - @unpack G_N_BAR, E_N, G_K_BAR, E_K, G_LEAK, E_LEAK, C_M = p - @unpack v, m, h, n = u - mα = exprel(-0.10 * (v + 35)) - mβ = 4.0 * exp(-(v + 60) / 18.0) - hα = 0.07 * exp(- ( v + 60) / 20) - hβ = 1 / (exp(-(v+30)/10) + 1) - nα = 0.1 * exprel(-0.1 * (v+50)) - nβ = 0.125 * exp( -(v+60) / 80) - iNa = G_N_BAR * (v - E_N) * (m^3) * h - iK = G_K_BAR * (v - E_K) * (n^4) - iLeak = G_LEAK * (v - E_LEAK) - iStim = _istim(t) - D.v = -(iNa + iK + iLeak + iStim) / C_M - D.m = -(mα + mβ) * m + mα - D.h = -(hα + hβ) * h + hα - D.n = -(nα + nβ) * n + nα - return nothing +function build_hh(;name) + @parameters begin + E_N = 55.0 ## Reversal potential of Na (mV) + E_K = -72.0 ## Reversal potential of K (mV) + E_LEAK = -49.0 ## Reversal potential of leaky channels (mV) + G_N_BAR = 120.0 ## Max. Na channel conductance (mS/cm^2) + G_K_BAR = 36.0 ## Max. K channel conductance (mS/cm^2) + G_LEAK = 0.30 ## Max. leak channel conductance (mS/cm^2) + C_M = 1.0 ## membrane capacitance (uF/cm^2)) + end + + @variables begin + t + mα(t) + mβ(t) + hα(t) + hβ(t) + nα(t) + nβ(t) + iNa(t) + iK(t) + iLeak(t) + iStim(t) + v(t) = -59.8977 + m(t) = 0.0536 + h(t) = 0.5925 + n(t) = 0.3192 + end + + D = Differential(t) + + eqs = [ + mα ~ exprel(-0.10 * (v + 35)), + mβ ~ 4.0 * exp(-(v + 60) / 18.0), + hα ~ 0.07 * exp(- ( v + 60) / 20), + hβ ~ 1 / (exp(-(v+30)/10) + 1), + nα ~ 0.1 * exprel(-0.1 * (v+50)), + nβ ~ 0.125 * exp( -(v+60) / 80), + iNa ~ G_N_BAR * (v - E_N) * (m^3) * h, + iK ~ G_K_BAR * (v - E_K) * (n^4), + iLeak ~ G_LEAK * (v - E_LEAK), + iStim ~ -6.6 * (20 < t) * (t < 21) + (-6.9) * (60 < t) * (t < 61), + D(v) ~ -(iNa + iK + iLeak + iStim) / C_M, + D(m) ~ -(mα + mβ) * m + mα, + D(h) ~ -(hα + hβ) * h + hα, + D(n) ~ -(nα + nβ) * n + nα, + ] + + sys = ODESystem(eqs; name) + return structural_simplify(sys) end #--- -ps = ( - E_N = 55.0, ## Reversal potential of Na (mV) - E_K = -72.0, ## Reversal potential of K (mV) - E_LEAK = -49.0, ## Reversal potential of leaky channels (mV) - G_N_BAR = 120.0, ## Max. Na channel conductance (mS/cm^2) - G_K_BAR = 36.0, ## Max. K channel conductance (mS/cm^2) - G_LEAK = 0.30, ## Max. leak channel conductance (mS/cm^2) - C_M = 1.0 ## membrane capacitance (uF/cm^2)) -) -u0 = LVector(v=-59.8977, m=0.0536, h=0.5925, n=0.3192) tend = 100.0 - -#--- -prob = ODEProblem(hh!, u0, tend, ps) +@named sys = build_hh() +prob = ODEProblem(sys, [], tend, []) #--- sol = solve(prob, tstops=[20., 60.]) #--- -p1 = plot(sol, idxs=[:v], ylabel="Membrane potential (mV)", xlabel="", legend=false, title="Fig 1.9") -p2 = plot(sol, idxs = [:m, :h, :n], xlabel="") -p3 = plot(_istim, sol.t, xlabel="Time (ms)", ylabel="Current", labels="Stimulation current") +@unpack v, m, h, n, iStim = sys +p1 = plot(sol, idxs=[v], ylabel="Membrane potential (mV)", xlabel="", legend=false, title="Fig 1.9") +p2 = plot(sol, idxs = [m, h, n], xlabel="") +p3 = plot(sol, idxs = [iStim], xlabel="Time (ms)", ylabel="Current", labels="Stimulation current") fig = plot(p1, p2, p3, layout=(3, 1), size=(600, 900), leftmargin=5*Plots.mm) fig |> PNG diff --git a/docs/figs/ch4-01.jl b/docs/figs/ch4-01.jl index 983758e5..f87a2e79 100644 --- a/docs/figs/ch4-01.jl +++ b/docs/figs/ch4-01.jl @@ -60,7 +60,7 @@ fig |> PNG # ## Fig. 4.2 B (Phase plot) -fig = plot( sols[1], idxs=(:a, :b), +fig = plot( sols[1], idxs=(1, 2), xlabel="[A]", ylabel="[B]", aspect_ratio=:equal, legend=nothing, title="Fig. 4.2 B (Phase plot)", @@ -86,7 +86,7 @@ fig |> PNG fig = plot(title="Fig. 4.3 B (Phase plot)") for sol in sols - plot!(fig, sol, idxs=(:a, :b), legend=nothing) + plot!(fig, sol, idxs=(1, 2), legend=nothing) end plot!(fig, xlabel="[A]", ylabel="[B]", xlims=(0., 2.), ylims=(0., 2.), aspect_ratio=:equal) @@ -120,7 +120,7 @@ fig = plot(title="Fig. 4.4 A (Phase plot with vector field)") quiver!(fig, xx, yy, quiver=∂F44, line=(:lightgrey), arrow=(:closed), aspect_ratio=:equal) for sol in sols - plot!(fig, sol, idxs=(:a, :b), linealpha=0.7, legend=nothing) + plot!(fig, sol, idxs=(1, 2), linealpha=0.7, legend=nothing) end plot!(fig, size=(600, 600), xlims=(rxy[1], rxy[end]), ylims=(rxy[1], rxy[end]), xlabel="[A]", ylabel="[B]") @@ -133,7 +133,7 @@ fig = plot(title="Fig. 4.5 A (Phase plot with nullclines)") ## Phase plots for sol in sols - plot!(fig, sol, idxs=(:a, :b), linealpha=0.7, lab=nothing) + plot!(fig, sol, idxs=(1, 2), linealpha=0.7, lab=nothing) end ## nullclines diff --git a/docs/figs/ch4-15.jl b/docs/figs/ch4-15.jl index a1ab4ae0..ef01c490 100644 --- a/docs/figs/ch4-15.jl +++ b/docs/figs/ch4-15.jl @@ -83,7 +83,7 @@ plot!(fig, Float64[], Float64[], line=(:black), label="A nullcline") contour!(fig, 0:0.01:4, 0:0.01:4, ∂B, levels=[0], cbar=false, line=(:black, :dash)) plot!(fig, Float64[], Float64[], line=(:black, :dash), label="B nullcline") for sol in sols - plot!(fig, sol, idxs=(:a, :b), label=nothing) + plot!(fig, sol, idxs=(1, 2), label=nothing) end plot!(fig, xlim=(0, 4), ylim=(0, 4), legend=:topright, size=(600, 600), xlabel="[A]", ylabel="[B]") @@ -140,7 +140,7 @@ plot!(fig, Float64[], Float64[], line=(:black), label="A nullcline") contour!(fig, 0:0.01:4, 0:0.01:4, ∂B, levels=[0], cbar=false, line=(:black, :dash)) plot!(fig, Float64[], Float64[], line=(:black, :dash), label="B nullcline") for sol in sols - plot!(fig, sol, idxs=(:a, :b), label=nothing) + plot!(fig, sol, idxs=(1, 2), label=nothing) end plot!(fig, xlim=(0, 4), ylim=(0, 4), legend=:topright, size=(600, 600), xlabel="[A]", ylabel="[B]") @@ -151,7 +151,7 @@ fig |> PNG sol = solve(ODEProblem(model415!, LVector(a=2.0, b=1.5), 10.0, ps2)) fig = plot(title="Fig 4.17") -plot!(fig, sol, idxs=(:a, :b), label=nothing, arrow=:closed) +plot!(fig, sol, idxs=(1, 2), label=nothing, arrow=:closed) quiver!(fig, xx, yy, quiver=∂F416, line=(:lightgrey), arrow=(:closed), aspect_ratio=:equal) contour!(fig, 1:0.01:3, 1:0.01:3, ∂A, levels=[0], cbar=false, line=(:black)) plot!(fig, identity, 0, 0, line=(:black), label="A nullcline") diff --git a/docs/figs/ch7-07.jl b/docs/figs/ch7-07.jl index bdb5a862..6ae25927 100644 --- a/docs/figs/ch7-07.jl +++ b/docs/figs/ch7-07.jl @@ -15,11 +15,11 @@ PNG(fig) = display("image/png", fig) #--- rn = @reaction_network begin - a1 / (1+RToverK1*(K2/(K2+L))^4), 0 --> M + a1 / (1 + RToverK1 * (K2 / (K2 + L))^4), 0 --> M δM, M --> 0 (c1 * M, δY), 0 <--> Y mm(Le, kL * Y, KML), 0 --> L - mm(L, 2 * kg * (Y/4), KMg), L ⇒ 0 + mm(L, 2 * kg * (Y / 4), KMg), L => 0 δL, L --> 0 end @@ -30,20 +30,20 @@ setdefaults!(rn, [ :δY => 0.03, :δL => 0.02, :a1 => 0.29, - :K2 => 2.92*1e6, + :K2 => 2.92 * 1e6, :RToverK1 => 213.2, :c1 => 18.8, - :kL => 6*1e4, + :kL => 6 * 1e4, :KML => 680, - :kg => 3.6*1e3, - :KMg => 7*1e5, + :kg => 3.6 * 1e3, + :KMg => 7 * 1e5, :Le => 0.0, :M => 0.01, :Y => 0.1, :L => 0.0 ]) -osys = convert(ODESystem, rn; remove_conserved = true) +osys = convert(ODESystem, rn; remove_conserved=true) equations(osys) @@ -52,14 +52,25 @@ equations(osys) @unpack Le = osys idx = findfirst(isequal(Le), parameters(osys)) -cb1 = PresetTimeCallback([500.], i -> begin i.p[idx] = 50.0; set_proposed_dt!(i, 0.01) end) -cb2 = PresetTimeCallback([1000.], i -> begin i.p[idx] = 100.0; set_proposed_dt!(i, 0.01) end) -cb3 = PresetTimeCallback([1500.], i -> begin i.p[idx] = 150.0; set_proposed_dt!(i, 0.01) end) -cb4 = PresetTimeCallback([2000.], i -> begin i.p[idx] = 0.0; set_proposed_dt!(i, 0.01) end) -cbs = CallbackSet(cb1, cb2, cb3, cb4) - -prob = ODEProblem(osys, [], (0., 2500.)) -sol = solve(prob, callback=cbs) +cb1 = PresetTimeCallback([500.0], i -> begin + i.p[idx] = 50.0 + set_proposed_dt!(i, 0.01) +end) +cb2 = PresetTimeCallback([1000.0], i -> begin + i.p[idx] = 100.0 + set_proposed_dt!(i, 0.01) +end) +cb3 = PresetTimeCallback([1500.0], i -> begin + i.p[idx] = 150.0 + set_proposed_dt!(i, 0.01) +end) +cb4 = PresetTimeCallback([2000.0], i -> begin + i.p[idx] = 0.0 + set_proposed_dt!(i, 0.01) +end) + +prob = ODEProblem(osys, [], (0.0, 2500.0)) +sol = solve(prob, callback=CallbackSet(cb1, cb2, cb3, cb4)) @unpack M, Y, L = osys fig = plot(sol, idxs=[Y], xlabel="Time (min)", title="Fig 7.7 (A)", label="β-galactosidase monomer") @@ -84,7 +95,7 @@ fig |> PNG # Compare the original model and the modified model rn_mod = @reaction_network begin - a1 / (1+RToverK1*(K2/(K2+L))^4), 0 --> M + a1 / (1 + RToverK1 * (K2 / (K2 + L))^4), 0 --> M δM, M --> 0 (c1 * M, δY), 0 <--> Y mm(Le, kL * 4 * Enz, KML), 0 --> L @@ -98,13 +109,13 @@ setdefaults!(rn_mod, [ :δY => 0.03, :δL => 0.02, :a1 => 0.29, - :K2 => 2.92*1e6, + :K2 => 2.92 * 1e6, :RToverK1 => 213.2, :c1 => 18.8, - :kL => 6*1e4, + :kL => 6 * 1e4, :KML => 680, - :kg => 3.6*1e3, - :KMg => 7*1e5, + :kg => 3.6 * 1e3, + :KMg => 7 * 1e5, :Le => 0.0, :M => 0.01, :Y => 0.1, @@ -112,7 +123,7 @@ setdefaults!(rn_mod, [ :Enz => 40.0 ]) -osys_mod = convert(ODESystem, rn_mod; remove_conserved = true) +osys_mod = convert(ODESystem, rn_mod; remove_conserved=true) equations(osys_mod) #--- @@ -129,19 +140,27 @@ idx_mod = findfirst(isequal(Le), parameters(prob_mod.f.sys)) lerange = range(0, 100, 101) eprob = EnsembleProblem(prob; - prob_func=(prob, i, repeat) -> begin prob.p[idx] = lerange[i]; prob end, - output_func=(sol, i) -> (sol[Y]/4, false) + prob_func=(prob, i, repeat) -> begin + prob.p[idx] = lerange[i] + prob + end, + output_func=(sol, i) -> (sol[Y] / 4, false) ) eprob_mod = EnsembleProblem(prob_mod; - prob_func=(prob, i, repeat) -> begin prob.p[idx_mod] = lerange[i]; prob end, - output_func=(sol, i) -> (sol[Y]/4, false) + prob_func=(prob, i, repeat) -> begin + prob.p[idx_mod] = lerange[i] + prob + end, + output_func=(sol, i) -> (sol[Y] / 4, false) ) sim = solve(eprob, DynamicSS(Rodas5()); trajectories=length(lerange)) sim_mod = solve(eprob_mod, DynamicSS(Rodas5()); trajectories=length(lerange)) -fig = plot(lerange, [sim sim_mod], label=["Original" "Modified"], +fig = plot(lerange, sim.u, label="Original") +fig = plot!(fig, lerange, sim_mod.u, label="Modified") +fig = plot!(fig, xlabel="External lactose concentration (μM)", ylabel="β-galactosidase", title="Fig 7.7 (B)" diff --git a/docs/intro-03-diffeq.jl b/docs/intro-03-diffeq.jl index 2f44d772..965d88f8 100644 --- a/docs/intro-03-diffeq.jl +++ b/docs/intro-03-diffeq.jl @@ -2,11 +2,12 @@ # Solving differential equations in Julia -## Standard procedures +**Standard procedures** - Define a model function representing the right-hand-side (RHS) of the sysstem. - Out-of-place form: `f(u, p, t)` where `u` is the state variable(s), `p` is the parameter(s), and `t` is the independent variable (usually time). The output is the right hand side (RHS) of the differential equation system. - In-place form: `f!(du, u, p, t)`, where the output is saved to `du`. The rest is the same as the out of place form. The in-place form has potential performance benefits since it allocates less than the out-of-place (`f(u, p, t)`) counterpart. + - Using ModelingToolkit.jl : define equations and build an ODE system. - Initial conditions (`u0`) for the state variable(s). - (Optional) define parameter(s) `p`. - Define a problem (e.g. `ODEProblem`) using the modeling function (`f`), initial conditions (`u0`), simulation time span (`tspan == (tstart, tend)`), and parameter(s) `p`. @@ -16,17 +17,6 @@ Documentation: -===# - -using DifferentialEquations -using Plots -Plots.default(linewidth=2) - -# PNG output in Literate.jl -PNG(fig) = display("image/png", fig) - -#=== - ### Exponential decay model The concentration of a decaying nuclear isotope could be described as an exponential decay: @@ -43,8 +33,14 @@ $$ ===# -# Model function, in the out-of-place form `f(u, p, t)` +using DifferentialEquations +using Plots +Plots.default(linewidth=2) +# PNG output in Literate.jl +PNG(fig) = display("image/png", fig) + +# Model function, in the out-of-place form `f(u, p, t)` expdecay(u, p, t) = p * u p = -1.0 ## Rate of exponential decay @@ -68,7 +64,7 @@ sol(1.0) # Time points sol.t -# Solutions at respective time points +# Solutions at corresponding time points sol.u #=== @@ -126,6 +122,38 @@ fig = plot(sol, label=["S" "I" "R"], legend=:right) fig |> PNG +#=== + +## Using ModelingToolkit.jl (recommended) + +[ModelingToolkit.jl](https://mtk.sciml.ai/dev/) is a high-level package for symbolic-numeric modeling and simulation in the Julia ecosystem. +===# + +using ModelingToolkit +using DifferentialEquations +using Plots +Plots.default(linewidth=2) + +# ### Exponential decay model + +@parameters λ ## Decaying rate constant +@variables t C(t) ## Time and concentration +D = Differential(t) ## Differential operator + +# Define an ODE with equations +@named expdecaySys = ODESystem([D(C) ~ -λ*C]) + +#--- +u0 = [C => 1.0] +p = [λ => 1.0] +tspan = (0.0, 2.0) + +prob = ODEProblem(expdecaySys, u0, tspan, p) +sol = solve(prob) + +fig = plot(sol) +fig |> PNG + #=== ### Lorenz system @@ -138,90 +166,70 @@ $$ \frac{dz}{dt} &= xy - \beta z \end{align} $$ - -In this example, we will use [LabelledArrays.jl](https://github.com/SciML/LabelledArrays.jl) to get DSL-like syntax. ===# - -using LabelledArrays using DifferentialEquations +using ModelingToolkit using Plots Plots.default(linewidth=2) -#--- -function lorenz!(du,u,p,t) - du.x = p.σ*(u.y-u.x) - du.y = u.x*(p.ρ-u.z) - u.y - du.z = u.x*u.y - p.β*u.z +function build_lorentz(; name) + @parameters begin + σ=10.0 + ρ=28.0 + β=8/3 + end + + @variables begin + t ## Independent variable (time) + x(t)=1.0 ## Independent variable (time) + y(t)=0.0 ## Independent variable (time) + z(t)=0.0 ## Independent variable (time) + end + + D = Differential(t) + + eqs = [ + D(x) ~ σ * ( y -x ), + D(y) ~ x * (ρ - z) - y, + D(z) ~ x * y - β * z + ] + + sys = ODESystem(eqs; name) + return sys end #--- -u0 = LVector(x=1.0, y=0.0, z=0.0) -p = LVector(σ=10.0, ρ=28.0, β=8/3) tspan = (0.0, 100.0) -prob = ODEProblem(lorenz!,u0,tspan,p) +@named sys = build_lorentz() +prob = ODEProblem(sys,[],tspan, []) sol = solve(prob) # x-y-z time-series fig = plot(sol) - fig |> PNG -# `idxs=(1, 2, 3)` makes a phase plot with 1st, 2nd, and the 3rd state variable. With `LabelledArrays`, you can use symbols instead of index numbers. - -fig = plot(sol, idxs=(:x, :y, :z)) +#--- +fig = plot(sol, idxs=[sys.y]) fig |> PNG -# The zeroth variable in `idxs` is the independent variable (usually time). The below command plots the time series of the second state variable (`y`). - -fig = plot(sol, idxs=(0, 2)) +# `idxs=(sys.x, sys.y, sys.z)` makes a phase plot with 1st, 2nd, and the 3rd state variable. +fig = plot(sol, idxs=(sys.x, sys.y, sys.z), label=false, size=(800,800)) fig |> PNG # ## Saving simulation results - using DataFrames using CSV df = DataFrame(sol) CSV.write("lorenz.csv", df) -#=== - -## Using ModelingToolkit.jl (advanced) - -[ModelingToolkit.jl](https://mtk.sciml.ai/dev/) is a high-level package for symbolic-numeric modeling and simulation ni the Julia DiffEq ecosystem. - -===# +# ### SIR model using DifferentialEquations using ModelingToolkit using Plots Plots.default(linewidth=2) -# ### Exponential decay model - -@parameters λ ## Decaying rate constant -@variables t C(t) ## Time and concentration -D = Differential(t) ## Differential operator - -# Define ODE equation(s) -eqs = [D(C) ~ -λ*C] - -# Define an ODE system -@named expdecaySys = ODESystem(eqs) - -#--- -u0 = [C => 1.0] -p = [λ => 1.0] -tspan = (0.0, 2.0) - -prob = ODEProblem(expdecaySys, u0, tspan, p) -sol = solve(prob) - -fig = plot(sol) -fig |> PNG - -# ### SIR model - @parameters β γ @variables t s(t) i(t) r(t) D = Differential(t) @@ -278,8 +286,10 @@ sir_rn = @reaction_network begin γ, I --> R end -p = [:β => 1.0, :γ => 0.3] -u0 = [:S => 0.99, :I => 0.01, :R => 0.00] +@unpack β, γ, S, I, R = sir_rn + +p = [β => 1.0, γ => 0.3] +u0 = [S => 0.99, I => 0.01, R => 0.00] tspan = (0., 20.) prob = ODEProblem(sir_rn, u0, tspan, p)