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

General Usage
1

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.

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?

2

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
3

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

2 Likes
4

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
    This is done by adding;
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?

5

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.

6

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

7

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