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

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(θ))

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 )

``````

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

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.

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