11
2- using DiffEqDevTools, Sundials, ODEInterfaceDiffEq,
3- Plots, DASSL, DASKR
4- using ModelingToolkit, OrdinaryDiffEq
2+ using DiffEqDevTools, ODEInterfaceDiffEq, Plots
3+ using ModelingToolkit, OrdinaryDiffEq, Symbolics
54using ModelingToolkit: t_nounits as t, D_nounits as D
65using LinearAlgebra
76
@@ -71,16 +70,11 @@ u0 = [y₁ => 0.0
7170 y₇ => 6.0
7271 y₈ => 0.0 ]
7372
74- @mtkbuild sys = ODESystem (eqs, t)
73+ @mtkcompile sys = System (eqs, t)
7574tspan = (0.0 , 0.2 )
7675mtkprob = ODEProblem (sys, u0, tspan)
7776ref_sol = solve (mtkprob, Rodas5P (), abstol = 1e-10 , reltol = 1e-10 )
7877
79- du = mtkprob. f (mtkprob. u0, mtkprob. p, 0.0 )
80- du0 = D .(unknowns (sys)) .=> du
81- daeprob = DAEProblem (sys, du0, [], tspan)
82- dae_ref_sol = solve (daeprob, IDA (), abstol = 1 / 10 ^ 7 , reltol = 1 / 10 ^ 7 )
83-
8478function transamp (du, u, p, t)
8579 y₁, y₂, y₃, y₄, y₅, y₆, y₇, y₈ = u
8680 Uₑ = 0.1 sin (200 π * t)
@@ -116,7 +110,7 @@ function transamp(du, u, p, t)
116110 nothing
117111end
118112
119- dirMassMatrix = Float64 .(ModelingToolkit . unwrap .(substitute .(
113+ dirMassMatrix = Float64 .(Symbolics . value .(substitute .(
120114 [- C₁ C₁ 0 0 0 0 0 0
121115 C₁ - C₁ 0 0 0 0 0 0
122116 0 0 - C₂ 0 0 0 0 0
@@ -130,8 +124,8 @@ mmf = ODEFunction(transamp, mass_matrix = dirMassMatrix)
130124mmprob = ODEProblem (mmf, [0.0 , 3.0 , 3.0 , 6.0 , 3.0 , 3.0 , 6.0 , 0.0 ], tspan)
131125mm_refsol = solve (mmprob, Rodas5 (), reltol = 1e-12 , abstol = 1e-12 )
132126
133- probs = [mtkprob, daeprob, mmprob]
134- refs = [ref_sol, ref_sol, mm_refsol];
127+ probs = [mtkprob, mmprob]
128+ refs = [ref_sol, mm_refsol];
135129
136130
137131plot (ref_sol, idxs = [y₁, y₂, y₃, y₄, y₅, y₆, y₇, y₈])
@@ -145,10 +139,8 @@ reltols = 1.0 ./ 10.0 .^ (1:4);
145139setups = [Dict (:prob_choice => 1 , :alg => Rodas4 ()),
146140 Dict (:prob_choice => 1 , :alg => FBDF ()),
147141 Dict (:prob_choice => 1 , :alg => QNDF ()),
148- Dict (:prob_choice => 1 , :alg => radau ()),
142+ Dict (:prob_choice => 2 , :alg => radau ()),
149143 Dict (:prob_choice => 1 , :alg => RadauIIA5 ()),
150- Dict (:prob_choice => 2 , :alg => DFBDF ()),
151- Dict (:prob_choice => 2 , :alg => IDA ())
152144]
153145
154146wp = WorkPrecisionSet (probs, abstols, reltols, setups;
@@ -160,13 +152,10 @@ abstols = 1.0 ./ 10.0 .^ (6:8)
160152reltols = 1.0 ./ 10.0 .^ (2 : 4 );
161153setups = [Dict (:prob_choice => 1 , :alg => Rosenbrock23 ()),
162154 Dict (:prob_choice => 1 , :alg => Rodas4 ()),
163- Dict (:prob_choice => 2 , :alg => IDA ()),
164- Dict (:prob_choice => 3 , :alg => Rodas5P ()),
165- Dict (:prob_choice => 3 , :alg => Rodas4 ()),
166- Dict (:prob_choice => 3 , :alg => rodas ()),
167- Dict (:prob_choice => 3 , :alg => FBDF ()),
168- Dict (:prob_choice => 2 , :alg => IDA ()),
169- Dict (:prob_choice => 2 , :alg => DASKR. daskr ())
155+ Dict (:prob_choice => 2 , :alg => Rodas5P ()),
156+ Dict (:prob_choice => 2 , :alg => Rodas4 ()),
157+ Dict (:prob_choice => 2 , :alg => rodas ()),
158+ Dict (:prob_choice => 2 , :alg => FBDF ()),
170159]
171160wp = WorkPrecisionSet (probs, abstols, reltols, setups;
172161 save_everystep = false , appxsol = refs, maxiters = Int (1e5 ), numruns = 10 )
@@ -179,10 +168,8 @@ setups = [Dict(:prob_choice => 1, :alg=>Rosenbrock23()),
179168 Dict (:prob_choice => 1 , :alg => Rodas4 ()),
180169 Dict (:prob_choice => 1 , :alg => FBDF ()),
181170 Dict (:prob_choice => 1 , :alg => QNDF ()),
182- Dict (:prob_choice => 1 , :alg => radau ()),
171+ Dict (:prob_choice => 2 , :alg => radau ()),
183172 Dict (:prob_choice => 1 , :alg => RadauIIA5 ()),
184- Dict (:prob_choice => 2 , :alg => DFBDF ()),
185- Dict (:prob_choice => 2 , :alg => IDA ())
186173]
187174wp = WorkPrecisionSet (probs, abstols, reltols, setups; error_estimate = :l2 ,
188175 save_everystep = false , appxsol = refs, maxiters = Int (1e5 ), numruns = 10 )
@@ -193,13 +180,10 @@ abstols = 1.0 ./ 10.0 .^ (6:8)
193180reltols = 1.0 ./ 10.0 .^ (2 : 4 );
194181setups = [Dict (:prob_choice => 1 , :alg => Rosenbrock23 ()),
195182 Dict (:prob_choice => 1 , :alg => Rodas4 ()),
196- Dict (:prob_choice => 2 , :alg => IDA ()),
197- Dict (:prob_choice => 3 , :alg => Rodas5P ()),
198- Dict (:prob_choice => 3 , :alg => Rodas4 ()),
199- Dict (:prob_choice => 3 , :alg => rodas ()),
200- Dict (:prob_choice => 3 , :alg => FBDF ()),
201- Dict (:prob_choice => 2 , :alg => IDA ()),
202- Dict (:prob_choice => 2 , :alg => DASKR. daskr ())
183+ Dict (:prob_choice => 2 , :alg => Rodas5P ()),
184+ Dict (:prob_choice => 2 , :alg => Rodas4 ()),
185+ Dict (:prob_choice => 2 , :alg => rodas ()),
186+ Dict (:prob_choice => 2 , :alg => FBDF ()),
203187]
204188wp = WorkPrecisionSet (probs, abstols, reltols, setups; error_estimate = :l2 ,
205189 save_everystep = false , appxsol = refs, maxiters = Int (1e5 ), numruns = 10 )
@@ -210,16 +194,13 @@ abstols = 1.0 ./ 10.0 .^ (7:12)
210194reltols = 1.0 ./ 10.0 .^ (4 : 9 )
211195
212196setups = [Dict (:prob_choice => 1 , :alg => Rodas5P ()),
213- Dict (:prob_choice => 3 , :alg => Rodas5P ()),
197+ Dict (:prob_choice => 2 , :alg => Rodas5P ()),
214198 Dict (:prob_choice => 1 , :alg => Rodas4 ()),
215- Dict (:prob_choice => 3 , :alg => Rodas4 ()),
199+ Dict (:prob_choice => 2 , :alg => Rodas4 ()),
216200 Dict (:prob_choice => 1 , :alg => FBDF ()),
217201 Dict (:prob_choice => 1 , :alg => QNDF ()),
218- Dict (:prob_choice => 1 , :alg => radau ()),
202+ Dict (:prob_choice => 2 , :alg => radau ()),
219203 Dict (:prob_choice => 1 , :alg => RadauIIA5 ()),
220- Dict (:prob_choice => 2 , :alg => DFBDF ()),
221- Dict (:prob_choice => 2 , :alg => IDA ()),
222- Dict (:prob_choice => 2 , :alg => DASKR. daskr ())
223204]
224205
225206wp = WorkPrecisionSet (probs, abstols, reltols, setups;
0 commit comments