# Differential Equation solve problem

**URL:** <https://discourse.julialang.org/t/differential-equation-solve-problem/67587>\
**Category:** Modelling & Simulations\
**Tags:** question, diffeq\
**Created:** [September 2, 2021, 3:55pm UTC](https://discourse.julialang.org/t/differential-equation-solve-problem/67587 "2021-09-02T15:55:07Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![hzgzh](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/hzgzh/32/6596_2.png) [@hzgzh](https://discourse.julialang.org/u/hzgzh)\
**Post date:** [September 2, 2021, 3:55pm UTC](https://discourse.julialang.org/t/differential-equation-solve-problem/67587/1 "2021-09-02T15:55:07Z")

</div>

I have a dae problem to solve ,but always have singualar error, could someone help?

```julia
using DifferentialEquations, Plots

function ga!(du,u,p,t)
    t2=0.02;
    τ, kp, ki, tc=p
          
    du[1] = ki * (u[5] - u[2])     
    du[2] = 1/t2 * (u[3]-u[2])
    du[3] = u[4] / tc
    du[4] = u[4] - (u[1] + kp * (u[5] - u[2]))
    du[5] = 0.0
end

u0 = [0, 0.1, 0.1, 0.0, 0.1]

p = [0.1, 1, 1, 1]

function condition(u, t, integrator)
    t - integrator.p[1]
end

function affect!(integrator)
    integrator.u[5] = 0.9
end
cb = ContinuousCallback(condition,affect!)

discrete_condition1 = (u, t, integrator) -> integrator.u[4] > 1.0 
discrete_affect1! = integrator -> integrator.u[4] = 1.0
db1 = DiscreteCallback(discrete_condition1, discrete_affect1!)

discrete_condition2 = (u, t, integrator) -> integrator.u[4] < -1.0 
discrete_affect2! = integrator -> integrator.u[4] = -1.0
db2 = DiscreteCallback(discrete_condition2, discrete_affect2!)

M = [1.0 0.0 0.0 0.0 0.0;
     0.0 1.0 0.0 0.0 0.0;
     0.0 0.0 1.0 0.0 0.0; 
     0.0 0.0 0.0 0.0 0.0;
     0.0 0.0 0.0 0.0 0.0]

f = ODEFunction((du,u,p,t)->ga!(du,u,p,t),mass_matrix = M)

tspan = (0.0, 10.0)

prob=ODEProblem(f,u0,tspan,p;callback = CallbackSet(cb, db1, db2))

sol=solve(prob,Rosenbrock23();verbose = true)

```

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [September 2, 2021, 4:10pm UTC](https://discourse.julialang.org/t/differential-equation-solve-problem/67587/2 "2021-09-02T16:10:49Z")

</div>

That doesn’t look like an Index-1 DAE? Try using the tricks in [Automated Index Reduction of DAEs · ModelingToolkit.jl](https://mtk.sciml.ai/dev/mtkitize_tutorials/modelingtoolkitize_index_reduction/)

---

<div class="post-metadata">

**Author:** ![hzgzh](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/hzgzh/32/6596_2.png) [@hzgzh](https://discourse.julialang.org/u/hzgzh)\
**Post date:** [September 3, 2021, 3:26am UTC](https://discourse.julialang.org/t/differential-equation-solve-problem/67587/3 "2021-09-03T03:26:15Z")

</div>

my system like the figure:

 ![model](https://global.discourse-cdn.com/julialang/original/3X/d/2/d2faf599b1b57d376360de7e06f714dc7d0c1e3f.png)  
u1-u5 is odefunction state, u5 is step implement by callback, the saturation of u4 is implement by discretecallback.  
as your subjust,I try the modelingtoolkitize, structural\_simplify,dae\_index\_lowering, it dont work.

```julia
using DifferentialEquations, Plots

using ModelingToolkit
using LinearAlgebra
function ga!(du,u,p,t)
    t2=0.02;
    τ, kp, ki, tc=p
    du[1] = ki * (u[4] - u[2])     
    du[2] = 1/t2 * (u[3]-u[2])

    du[3] = u[4] / tc
    du[4] = u[4] - (u[1] + kp * (u[5] - u[2]))
    du[5] = 0.0
end

u0 = [0, 0.1, 0.1, 0.0, 0.1]
p = [0.1, 1, 1, 1]

function condition(u, t, integrator)
    t - integrator.p[1]
end

function affect!(integrator)
    integrator.u[5] = 0.9
end

cb = ContinuousCallback(condition,affect!)
discrete_condition1 = (u, t, integrator) -> integrator.u[4] > 1.0 
discrete_affect1! = integrator -> integrator.u[4] = 1.0
db1 = DiscreteCallback(discrete_condition1, discrete_affect1!)

discrete_condition2 = (u, t, integrator) -> integrator.u[4] < -1.0 
discrete_affect2! = integrator -> integrator.u[4] = -1.0
db2 = DiscreteCallback(discrete_condition2, discrete_affect2!)

f = ODEFunction((du,u,p,t)->ga!(du,u,p,t); mass_matrix = Diagonal([1,1,1,0,0]))

tspan = (0.0, 10.0)

prob=ODEProblem(f,u0,tspan,p;callback = CallbackSet(cb,db1,db2))

mtk_prob = modelingtoolkitize(prob)

sys = structural_simplify(dae_index_lowering(mtk_prob))

prob = structural_simplify(dae_index_lowering(sys))

@time sol=solve(prob,Rosenbrock23();verbose = true)

plot(sol)

```

so I change the code the saturation of u4 implement in code ,it work

```julia
using DifferentialEquations, Plots
using ModelingToolkit
using LinearAlgebra
function ga!(du,u,p,t)
    t2=0.02;
    τ, kp, ki, tc=p
    po = (u[1] + kp * (u[4] - u[2]))
    po > 1 ? po = 1 : po
    po < -1 ? po = -1 : po   
    du[1] = ki * (u[4] - u[2])     
    du[2] = 1/t2 * (u[3]-u[2])
    du[3] = po / tc
end

u0 = [0, 0.1, 0.1, 0.1]
p = [0.1, 1, 1, 1]

function condition(u, t, integrator)
    t - integrator.p[1]
end

function affect!(integrator)
    integrator.u[4] = 0.9
end
cb = ContinuousCallback(condition,affect!)

f = ODEFunction((du,u,p,t)->ga!(du,u,p,t))
tspan = (0.0, 10.0)

prob=ODEProblem(f,u0,tspan,p;callback = cb)

@time sol=solve(prob,Rosenbrock23();verbose = true)

plot(sol)

```

so for common nonlinear part in ode or dae like `Blacklash,dead zone,rate limiter,saturation,relay`, which is the recommend method to implement by code it in odefunction or using callback, how about performance or recommend tutorial code for nonlinear part implement, thank @ChrisRackauckas