System of equation, domain error: initial value setting problem?

Specific Domains Numerics
1

Good afternoon,

I have been blocked for a very long time on this so any help would really help me as I am completely stuck there. I am new to Julia so please pardon my errors.

I have a system of equations (very simplified below) and keep getting “domain error”. It works for certain initial values, or for some equations. I noticed that whenever the value of equations evaluated at the initial values are below zero, it seems to generate a domain error.

For instance, when i change equation F[1] of the sytem from F[1]= x[2] * D , it works i.e. I do not get domain error (- the algorithm does not converge but as this is a simplified version of my true system, I am not concerned about that for now).

But when I change it to F[1]= x[2] * D + 2 , I get a domain error (so this is even before looking into convergence for a solution).

Would anyone know why? Understanding what is going on here would help me solve this error and the future ones. Many thanks for your help!

using QuantEcon
using NLsolve
using ForwardDiff

Re=5.0
Rl=2.0
E=4.0
D=0.1

function foc2(θlow::AbstractFloat,
                D::AbstractFloat) 
    
v(c) = 4 * c
c21(θ) = 1/ (1 - θ)   
val2(θ) = θ * D * v(c21(θ))
val22 = quadgk(val2, 0, θlow)[1]

return val22 
    
end 

 function f!(F ,x) 
            F[1]= x[2] * D
            F[2]= foc2(x[9], D)
            F[3]= D - x[3] - x[1] 
            F[4]= E - x[4] - D 
            F[5]= x[2] * D - x[5] * ((Re-Rl)/Rl) 
            F[6]= x[4] - x[6] - x[5]
            F[7]= x[8] -  x[5] / x[1]  
            F[8]= x[9] - x[3]/ (D * x[2]) 
            F[9]= x[2] * D * x[7] + x[2] * D * (1 - x[7]) * (x[8]/Re) - x[3] - x[1] * x[8] 
       
end 
        
        
    initial_x = [D; D; 0.0; E/2.0; 0.01; E/2.0; 0.3; 0.1; 0.1 ] 
    initial_F = similar(initial_x)
        
    df = OnceDifferentiable(f!, initial_x, initial_F)
       
    results1_bk = mcpsolve(df, [0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01 , 0.0], 
                        [D, D, D, E, E, D, 0.9, 0.99, 0.9], initial_x )
2

Worked for me in v0.6.4 when I put F[1]= x[2] * D + 2.0. No “domain error”.

Please confirm you have checked for divide by zero etc and anything else you’ve tried. I always specify the type of arguments in functions as it sometimes helps. I also habitually write 0.0 for a float and 0 for an integer. If nothing else it makes code less ambiguous. I know Julia does the conversions so it’s largely unnecessary.

1 Like
3

Thanks a lot for your help! I tried F[1]= x[2] * D - 2.0 (os changing to .0 as you suggested and i still get a domain error. I am using Julia 0.6.3 . I am trying tu update to v0.6.4 to see if this is an error in my version of Julia maybe (but it would be surprising, right?)

Yes i had to correct for domain error coming indeed from dividing by 0. I also had another issue: some integrals were not defined for values of unknow near zero that I (probably artificially) solved by imposing the unknowns to be above 0.01 instead of zero.

Thanks again.

4

I’m also able to reproduce in 0.7 for F[1]= x[2] * D - 2 (previously you wrote F[1]= x[2] * D + 2, which did not produce the error). I had to add using QuadGK to get it to work on 0.7

I think the reason it’s not working is that this line:

c21(θ) = 1/ (1 - θ)

produces NaN when θ = 1, which occurs if θlow > 1, which seems to happen for this particular system of equations. This integral does diverge, so it is perhaps not too surprising you’re seeing this behavior.

1 Like