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

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 https://diffeq.sciml.ai/stable/tutorials/bvp_example/.

My code is:

``````
# 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:

``````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?

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.

3 Likes

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

2 Likes

Thanks for the reply’s jonniedie and mauro3.

I have now tried to fix the problem by:

1. Using `isoutofdomain` formualtion below
``````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`
``````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

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`?

You still need to make sure that the objective function does not error when it goes negative, see Using DiffEq isoutofdomain to enforce positivity.

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]`?

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

``````
# 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).

I don’t understand what you mean with this.

2 Likes