2
$\begingroup$

I'm trying to solve a differential algebraic equation in Julia's modelingTookKit.JL, where the vector field has the form f(X) = 0.

I found an example of a DAE in the below link modelingToolkit.JL DAE](https://benchmarks.sciml.ai/html/DAE/ROBERDAE.html) with follwing code. However, when I replace say

D(y₁) ~ -k₁*y₁ + k₃*y₂*y₃

by

0 ~ -k₁*y₁ + k₃*y₂*y₃ - D(y₁)

I get an error. If someone knows of or could provide an example where the the equation of the vector field is only implicitly defined i.e. f(X) = 0, it would be much appreciated.

using OrdinaryDiffEq, DiffEqDevTools, Sundials, ModelingToolkit, ODEInterfaceDiffEq,
      Plots, DASSL, DASKR
using LinearAlgebra

@variables t y₁(t)=1.0 y₂(t)=0.0 y₃(t)=0.0
@parameters k₁=0.04 k₂=3e7 k₃=1e4
D = Differential(t)
eqs = [
  D(y₁) ~ -k₁*y₁ + k₃*y₂*y₃
  D(y₂) ~  k₁*y₁ - k₃*y₂*y₃ - k₂*y₂^2
  0 ~  y₁ + y₂ + y₃ - 1
]
@named sys = ODESystem(eqs)
simpsys = structural_simplify(sys)
mmprob = ODEProblem(sys,[],(0.0,1e5))
daeprob = DAEProblem(sys,[D(y₁)=>-0.04,
                              D(y₂)=>0.04,
                              D(y₃)=>0.0],[],(0.0,1e5))
odaeprob = ODAEProblem(simpsys,[],(0.0,1e5))

ref_sol = solve(daeprob,IDA(),abstol=1/10^14,reltol=1/10^14);
ode_ref_sol = solve(odaeprob,CVODE_BDF(),abstol=1/10^14,reltol=1/10^14);

probs = [mmprob,daeprob,odaeprob]
refs = [ref_sol,ref_sol,ode_ref_sol];
$\endgroup$

1 Answer 1

3
$\begingroup$
0 ~ -k₁*y₁ + k₃*y₂*y₃ - D(y₁)

is handled in this PR. This feature will be included in the next release of ModelingToolkit.jl. Also, note that you should use

simpsys = structural_simplify(sys)
mmprob = ODEProblem(simpsys,[],(0.0,1e5))

The structural_simplify function doesn't mutate the system.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.