# DifferentialEquations.jl, Error: Exponentiation yielding a complex result requires a complex argument

**URL:** <https://discourse.julialang.org/t/differentialequations-jl-error-exponentiation-yielding-a-complex-result-requires-a-complex-argument/44935>\
**Category:** General Usage\
**Tags:** package\
**Created:** [August 14, 2020, 2:32pm UTC](https://discourse.julialang.org/t/differentialequations-jl-error-exponentiation-yielding-a-complex-result-requires-a-complex-argument/44935 "2020-08-14T14:32:57Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![Gregers](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gregers/32/17339_2.png) [@Gregers](https://discourse.julialang.org/u/Gregers)\
**Post date:** [August 14, 2020, 2:32pm UTC](https://discourse.julialang.org/t/differentialequations-jl-error-exponentiation-yielding-a-complex-result-requires-a-complex-argument/44935/1 "2020-08-14T14:32:57Z")

</div>

Hey,  
I am trying to solve a boundary value problem using DifferentialEquations.jl.  
Unfortunately, I get the error:

“DomainError with -117.96192353385844:  
Exponentiation yielding a complex result requires a complex argument.  
Replace x^y with (x+0im)^y, Complex(x)^y, or similar.”

I have replicated the examples given in the tutorials, first two examples [Boundary Value Problems · DifferentialEquations.jl](https://diffeq.sciml.ai/stable/tutorials/bvp_example/).

My code is:

```julia

# Define parameter values
my = 0.04                                                          
sigma = 0.13                                                    
gamma = 5.0                                                    
rho = 0.02                                                        
r = 0.02                                                           
t_min = 3.0                                                        
t_max = 100.0                                                
tspan = (t_min,t_max)                                   

eta = (1/gamma)*rho+(1-1/gamma)*(r+0.5*my^2/(gamma*sigma^2))   
x₀ = eta^(-gamma)*(t_min^(1-gamma)/(1-gamma)) 
y₀ = eta^(-gamma)*t_min^(-gamma) value function

u₀ = [x₀, y₀]                                                     

function softhabits!(du,u,p,t)
    du[1] = u[2]                                                                                  
    du[2] = (0.5*u[2]^2*my^2/sigma^2)/((gamma/(1-gamma))*u[2]^(1-1/gamma)-rho*u[1]+r*t*u[2]) 
end 

function bc1!(residual, u, p, t)
    residual[1] = u[end][2] # u[2] equals zero at t_max
end

using OrdinaryDiffEq
using Plots

bvp3 = BVProblem(softhabits!, bc1!, u₀, tspan)
sol3 = solve(bvp3, Shooting(Vern7()))
plot(sol3)

```

This problem has a know analytical solution, which I have solved earlier same parameters as above and with this code:

```julia
function softhabits!(du,u,p,t)
    du[1] = u[2]                                                                                  
    du[2] = (0.5*u[2]^2*my^2/sigma^2)/((gamma/(1-gamma))*u[2]^(1-1/gamma)-rho*u[1]+r*t*u[2]) 
end 

prob = ODEProblem(softhabits!,u₀,tspan,maxiters=10e5)                                         
sol = solve(prob,reltol=1e-12,abstol=1e-12, saveat=1.0)                                          

```

Do some one have any clue on fixing this bug?

---

<div class="post-metadata">

**Author:** ![jonniedie](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jonniedie/32/12842_2.png) [@jonniedie](https://discourse.julialang.org/u/jonniedie)\
**Post date:** [August 14, 2020, 3:24pm UTC](https://discourse.julialang.org/t/differentialequations-jl-error-exponentiation-yielding-a-complex-result-requires-a-complex-argument/44935/2 "2020-08-14T15:24:11Z")

</div>

> [@Gregers](#):
>
> `u[2]^(1-1/gamma)`

This term will be complex when `u[2]` is negative. Rather than having exponentiation be type-unstable (which would make it slower), Julia requires the user to specify when exponentiation will lead to a complex result by passing in a `Complex` value explicitly and will throw an error otherwise. If you are expecting your `u[2]` value to be complex, you should make your initial conditions `Complex`. If you aren’t expecting it to be complex, then it shouldn’t be going negative and there might be a mistake in one of your equations.

EDIT: I didn’t notice this was a BVP when I wrote this. That could be driving your `u[2]` value negative. See the `isoutofdomain` comment below.

---

<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:** [August 14, 2020, 3:30pm UTC](https://discourse.julialang.org/t/differentialequations-jl-error-exponentiation-yielding-a-complex-result-requires-a-complex-argument/44935/3 "2020-08-14T15:30:16Z")

</div>

Have a look at `isoutofdomain` in the docs and in posts here on discourse.

---

<div class="post-metadata">

**Author:** ![Gregers](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gregers/32/17339_2.png) [@Gregers](https://discourse.julialang.org/u/Gregers)\
**Post date:** [August 17, 2020, 5:51am UTC](https://discourse.julialang.org/t/differentialequations-jl-error-exponentiation-yielding-a-complex-result-requires-a-complex-argument/44935/4 "2020-08-17T05:51:30Z")

</div>

Thanks for the reply’s jonniedie and mauro3.

I have now tried to fix the problem by:

1. Using `isoutofdomain` formualtion below

```julia
sol1 = solve(bvp1, GeneralMIRK4(), dt=0.05 , isoutofdomain=(u,p,t) -> any(y -> y < 0, u))

```

Still yields the error:  
DomainError with -59466.70272947129:  
Exponentiation yielding a complex result requires a complex argument.  
Replace x^y with (x+0im)^y, Complex(x)^y, or similar.

Is there something I forget when adding the `isoutofdomain`?

1. Using the `callback` function - stop executing when `u[2]<0`  
This is done by adding;

```julia
function condition(t,u,integrator)
  u[2] < 1e-9
end

affect!(integrator) = terminate!(integrator)

cb = ContinuousCallback(condition, affect!)

sol1 = solve(bvp1, GeneralMIRK4(), dt=0.05, callback=cb)

```

The code is colse to equal: ContinouseCallBack Example Bouning Ball  
[https://diffeq.sciml.ai/release-2.0/features/callback\_functions.html#ContinuousCallback-Examples-1](https://diffeq.sciml.ai/release-2.0/features/callback_functions.html#ContinuousCallback-Examples-1)

Still I end up with the domain error:

DomainError with -59466.70272947129:  
Exponentiation yielding a complex result requires a complex argument.  
Replace x^y with (x+0im)^y, Complex(x)^y, or similar.

Is there an explanation for how the solver accepts `u[2]<0`?

---

<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:** [August 17, 2020, 6:19am UTC](https://discourse.julialang.org/t/differentialequations-jl-error-exponentiation-yielding-a-complex-result-requires-a-complex-argument/44935/5 "2020-08-17T06:19:08Z")

</div>

You still need to make sure that the objective function does not error when it goes negative, see [Using DiffEq isoutofdomain to enforce positivity - #2 by ChrisRackauckas](https://discourse.julialang.org/t/using-diffeq-isoutofdomain-to-enforce-positivity/42443/2).

---

<div class="post-metadata">

**Author:** ![Gregers](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gregers/32/17339_2.png) [@Gregers](https://discourse.julialang.org/u/Gregers)\
**Post date:** [August 17, 2020, 9:27am UTC](https://discourse.julialang.org/t/differentialequations-jl-error-exponentiation-yielding-a-complex-result-requires-a-complex-argument/44935/6 "2020-08-17T09:27:07Z")

</div>

Thanks mauro3 for your relpy.

Is i correct that `isoutofdomain=(u,p,t) -> any(x -> x < 0, u)` ensure that all variables to be within the positive domain?

If so, is there a way to specify that `u[2] > 0 and is free u[2]`?

---

<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:** [August 17, 2020, 2:28pm UTC](https://discourse.julialang.org/t/differentialequations-jl-error-exponentiation-yielding-a-complex-result-requires-a-complex-argument/44935/7 "2020-08-17T14:28:56Z")

</div>

This works for me (albeit, I don’t know if it is right):

```julia

# Define parameter values
my = 0.04
sigma = 0.13
gamma = 5.0
rho = 0.02
r = 0.02
t_min = 3.0
t_max = 100.0
tspan = (t_min,t_max)

eta = (1/gamma)*rho+(1-1/gamma)*(r+0.5*my^2/(gamma*sigma^2))
x₀ = eta^(-gamma)*(t_min^(1-gamma)/(1-gamma))
y₀ = eta^(-gamma)*t_min^(-gamma) # value function

u₀ = [x₀, y₀]

function softhabits!(du,u,p,t)
    du[1] = u[2]
    du[2] = (0.5*u[2]^2*my^2/sigma^2)/((gamma/(1-gamma))*max(u[2],0)^(1-1/gamma)-rho*u[1]+r*t*u[2])
end

function bc1!(residual, u, p, t)
    residual[1] = u[end][2] # u[2] equals zero at t_max
end

using OrdinaryDiffEq, BoundaryValueDiffEq
using Plots

bvp3 = BVProblem(softhabits!, bc1!, u₀, tspan)
sol3 = solve(bvp3, Shooting(Vern7()), isoutofdomain=(u,p,t) -> u[2]<0)
plot(sol3)

```

Note, how I made sure that even at negative `u[2]` there is no error thrown with `max(u[2],0)`. This is necessary for some reason, presumably because `isoutofdomain` only gets called after the evaluation of the objective function (why, I don’t know).

> [@Gregers](#):
>
> If so, is there a way to specify that `u[2] > 0 and is free u[2]` ?

I don’t understand what you mean with this.