# Instability with DifferentialEquations.jl

I am trying to solve a second order ODE equation by using SecondOrderODEProblem as problem type. Consistently, I have encountered WARNING: Instability detected. Aborting.

I have tried various solver and still, the problem persists. What are the errors in my code? The same problem is being solved in Matlab without any error.

``````function odefnct(eta,f,df,ddf)

(Phi, iota) = investment(f)

# calculate psi
psi_L = eta
psi_R = min(f/df + eta, 1.0)

for n in 1:50
psi = (psi_L + psi_R)/2.0;
amplification = 1.0 - df/f*(psi - eta);

# VOLATILITY COMPUTATION
sigma_eta_eta = sigma*(psi - eta)/amplification;  # sigma_eta * eta
sigma_q = sigma_eta_eta*df/f;
sigma_theta = sigma_eta_eta*df/f;
risk_premium = - sigma_theta*(sigma + sigma_q);

if household_premium > 0.0 # households want to hold more
psi_R = psi;
else
psi_L = psi;
end
end

mu_q = r - (a - iota)/f - Phi + delta - sigma*sigma_q + risk_premium;
mu_eta_eta = -(psi - eta)*(sigma + sigma_q)*(sigma + sigma_q + sigma_theta) + eta*(a - iota)/f + eta*(1.0 - psi)*(delta_ - delta);
qpp = 2.0*(mu_q*f - df*mu_eta_eta)/sigma_eta_eta^2;
thetapp = 2.0*((rho - r)*f - df*mu_eta_eta)/sigma_eta_eta^2;

leverage = psi/eta;

ddf = thetapp
ddf = qpp

end

etaspan = (0.0,1.0)
f0 = [1.0 , q_]
df0 = [-1e+10, 1e+8]

function condition(t,u,integrator)
u - qmax
u
u
end

affect!(integrator) = terminate!(integrator)
cb = ContinuousCallback(condition,affect!)

QL = 0.0 ; QR = 1e+15
global sol

for iter in 1:50

df0 = ( QL +  QR)/ 2;  # this is q'(0)
prob = SecondOrderODEProblem(odefnct, f0, df0, etaspan)
sol = solve(prob, callback = cb)
u1,u2,u3,u4 = sol.u[end]

if abs(u4) < 1e-8  # if q'(eta) has reached zero, we
QL = df0;  # increase q'(0)
else        # if q(eta) reached qmax or theta'(0) reached 0
QR = df0;  # we reduce q'(0)
end

end
``````

What’s `q_`? Also, `investment` isn’t defined.

What method are you using to solve it in MATLAB? Given those initial velocities, is it stiff? Since you didn’t declare a solving method or what type of equation it is, it’s defaulting to an explicit non-stiff solver (like `ode45`). If you’re using something like `ode15s` in MATLAB, you should change the algorithm you’re using in Julia to be something similar (or just pass `alg_hints=[:stiff]`)

1 Like

#some constants
const Ip1=12.1299/27.2
const Ip2=0.78
const I=7.0e13
const omg1=0.057
#const omg2=0.0
#const ε=0.0
#const CEP=0.0

const T1=2.0π/omg1
#const T2=2.0
π/omg2
const E0=sqrt(I/(3.51e16))

I am so sorry that i have made a mistake when I use MATLAB to solve the equation. when I use the same initial condition, I found that ode45() ode15s()in Matlab still cant solve it. When I try to find where the problem is, it seems that when the time t arrive a certain value t0, the step next will go wrong. However, the sol before t0 is Ok,i have calculated du in time t0, there is nothing wrong with them, i think, if du don’t have problems, the iteration is stable ,and i just can’t understand why the solver throw a warning here.

1 Like

This doesn’t seem to be related to the original question, so please make a different post. The issue is certainly in the model somewhere: look for asymptotic properties that could lead to the issue.