# Domain error in circuit DAE

**URL:** <https://discourse.julialang.org/t/domain-error-in-circuit-dae/73542>\
**Category:** Modelling & Simulations\
**Tags:** diffeq, differentialequation\
**Created:** [December 23, 2021, 6:55pm UTC](https://discourse.julialang.org/t/domain-error-in-circuit-dae/73542 "2021-12-23T18:55:02Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![dfcaporale](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dfcaporale/32/29216_2.png) [@dfcaporale](https://discourse.julialang.org/u/dfcaporale)\
**Post date:** [December 23, 2021, 6:55pm UTC](https://discourse.julialang.org/t/domain-error-in-circuit-dae/73542/1 "2021-12-23T18:55:02Z")

</div>

I’m trying to solve the following DAE:

\begin{cases} i\_R=C \displaystyle \frac{dV\_o}{dt} \\ i\_R= \displaystyle \frac{V\_{in}(t)-V\_o-nV\_T \ln \displaystyle (1+\displaystyle \frac{i\_R}{I\_s}) }{R} \end{cases}

The solution obtained from a simulator is:

 ![](https://global.discourse-cdn.com/julialang/original/3X/2/4/24762c09f888017834cd5c835156e365c7da4dd4.png)

Nevertheless, we want to solve the model with `DifferentialEquations.jl`, and I’ve tried with the following code:

```julia
using DifferentialEquations, Plots

f = 10E3
ω = 2*π*f
n = 1
C = 0.01E-6
VT = 25E-3
IS = 1E-6
A = 1.0
R = 10E3

function func(out, du, u, p, t)
  V0, IR = u
  dV0, _ = du

  out[1] = C*dV0 - IR
  out[2] = (A*sin(ω*t) - V0 - n*VT*log(1+IR/IS)) - R*IR
end

B = A*ω/R/(1+n*VT/R/IS)

u₀ = [0, 0.0] # V0, IR
du₀ = [0, B] # dV0, dIR

dvars = [true, false]
tspan = (0, 2E-3)
prob = DAEProblem(func, du₀, u₀, tspan, differential_vars=dvars);

```

It gives the following error:

```julia
sol = solve(prob)
DomainError with -0.33407903800223493:
log will only return a complex result if called with a complex argument. Try log(Complex(x)).

```

What do you suggest to tackle this problem?

---

<div class="post-metadata">

**Author:** ![mauro3](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mauro3/32/292_2.png) [@mauro3](https://discourse.julialang.org/u/mauro3)\
**Post date:** [December 23, 2021, 7:01pm UTC](https://discourse.julialang.org/t/domain-error-in-circuit-dae/73542/2 "2021-12-23T19:01:21Z")

</div>

The term `1+IR/IS` used in `log` goes negative which will cause the error you see. This can be either because the model is not well formulated, wrongly coded, or because the DAE solver takes too large a time step. If it is the last option, you could try `isoutofdomain` or code the objective function such that the term does not go negative.

Maybe you also need to make sure that the ICs are good for the solver to get started.

---

<div class="post-metadata">

**Author:** ![dfcaporale](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dfcaporale/32/29216_2.png) [@dfcaporale](https://discourse.julialang.org/u/dfcaporale)\
**Post date:** [December 23, 2021, 7:52pm UTC](https://discourse.julialang.org/t/domain-error-in-circuit-dae/73542/4 "2021-12-23T19:52:52Z")

</div>

@mauro3 Thank you for your answer!

I’m wondering if is it possible to use `isoutofdomain` within a `DAEProblem`?

---

<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:** [December 23, 2021, 9:17pm UTC](https://discourse.julialang.org/t/domain-error-in-circuit-dae/73542/5 "2021-12-23T21:17:32Z")

</div>

IDA will throw a warning that it’s ignoring that keyword argument. Just reduce the solver tolerance and you should be fine.

---

<div class="post-metadata">

**Author:** ![dfcaporale](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dfcaporale/32/29216_2.png) [@dfcaporale](https://discourse.julialang.org/u/dfcaporale)\
**Post date:** [December 30, 2021, 12:42am UTC](https://discourse.julialang.org/t/domain-error-in-circuit-dae/73542/6 "2021-12-30T00:42:25Z")

</div>

Thank you, @ChrisRackauckas! The `log` turns the DAE to be unstable, I guess. We are considering a workaround by slightly modifying the equation.

---

<div class="post-metadata">

**Author:** ![dfcaporale](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dfcaporale/32/29216_2.png) [@dfcaporale](https://discourse.julialang.org/u/dfcaporale)\
**Post date:** [January 22, 2022, 1:41pm UTC](https://discourse.julialang.org/t/domain-error-in-circuit-dae/73542/7 "2022-01-22T13:41:10Z")

</div>

Hi again!

I’ve been trying a workaround to my [original issue](https://discourse.julialang.org/t/domain-error-in-circuit-dae/73542), by reordering some terms in my initial model (in short, now I have an `exp` function instead of a `log`):

```julia
using DifferentialEquations, Plots
using Sundials #IDA solver

f = 10E3
ω = 2*π*f
n = 1
C = 0.01E-6
VT = 25E-3
IS = 1E-6
A = 1.0
R = 10E3
Vin(t) = A*sin(ω*t)

function func(out, du, u, p, t)
  V0, IR = u
  dV0, _ = du

  out[1] = C*dV0 - IR
  aux = (R*IR - Vin(t) + V0)/(n*VT)
  out[2] = IS*(exp(-aux) - 1) - IR
end

B = A*ω/R/(1+n*VT/R/IS)

u₀ = [0, 0.0] # V0, IR
du₀ = [0, B] # dV0, dIR

dvars = [true, false]

tspan = (0, 3.6E-3)
prob = DAEProblem(func, du₀, u₀, tspan, differential_vars=dvars);

```

The minimum `reltol` that worked is `1E-15` (maximum accuracy I can request). So next, I tried reducing the `abstol` parameter:

```julia
sol = solve(prob, reltol=1e-15, abstol=1e-17)
plot(sol, vars=(0, 1), xlab="t (s)", ylab="Vo (V)", lab="")

```

leading to the following plot that doesn’t match what I expected:

 ![plot](https://global.discourse-cdn.com/julialang/original/3X/e/f/ef3026ec73715106b2c59a395881ad7dff0843da.png)

On the other hand, when reducing `abstol` further, the solution seems to be OK, but the time horizon is shorter (for `abstol=1E-18`, the time horizon I get is 1.5E-3 s).

Do you have any further suggestions that I can try?

---

<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:** [January 25, 2022, 8:16am UTC](https://discourse.julialang.org/t/domain-error-in-circuit-dae/73542/8 "2022-01-25T08:16:25Z")

</div>

> [@dfcaporale](#):
>
> `out[2] = IS*(exp(-aux) - 1) - IR`

Why not just write this as an equality and substitute? That would improve the numerical accuracy. I would suggest using [ModelingToolkit.jl](https://github.com/SciML/ModelingToolkit.jl) for this kind of thing. Its symbolic algorithms will do exactly these types of transformations.

---

<div class="post-metadata">

**Author:** ![dfcaporale](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dfcaporale/32/29216_2.png) [@dfcaporale](https://discourse.julialang.org/u/dfcaporale)\
**Post date:** [January 25, 2022, 2:11pm UTC](https://discourse.julialang.org/t/domain-error-in-circuit-dae/73542/9 "2022-01-25T14:11:53Z")

</div>

Thanks so much, @ChrisRackauckas!

I followed the [pendulum example](https://mtk.sciml.ai/dev/mtkitize_tutorials/modelingtoolkitize_index_reduction/) from the ModelingToolkit tutorial website and got the following:

 ![Captura de pantalla 2022-01-25 11.01.29](https://global.discourse-cdn.com/julialang/original/3X/e/3/e38fad13e4aa844ed515ff9e4bbf706e72d2f8a9.png)

Looks good to me! Let me know if you want to check the code. Now I’ll explore the ModelingToolkit tool to check if I can extract some of the transformations it does to the system to solve it with tools like [ReachabilityAnalysis](https://github.com/JuliaReach/ReachabilityAnalysis.jl).

Thanks so much again!