# 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[2])

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

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

# VOLATILITY COMPUTATION
sigma_eta_eta = sigma*(psi - eta)/amplification;  # sigma_eta * eta
sigma_q = sigma_eta_eta*df[2]/f[2];
sigma_theta = sigma_eta_eta*df[1]/f[1];
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[2] - Phi + delta - sigma*sigma_q + risk_premium;
mu_eta_eta = -(psi - eta)*(sigma + sigma_q)*(sigma + sigma_q + sigma_theta) + eta*(a - iota)/f[2] + eta*(1.0 - psi)*(delta_ - delta);
qpp = 2.0*(mu_q*f[2] - df[2]*mu_eta_eta)/sigma_eta_eta^2;
thetapp = 2.0*((rho - r)*f[1] - df[1]*mu_eta_eta)/sigma_eta_eta^2;

leverage = psi/eta;

ddf[1] = thetapp
ddf[2] = qpp

end

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

function condition(t,u,integrator)
u[2] - qmax
u[3]
u[4]
end

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

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

for iter in 1:50

df0[2] = ( 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[2];  # increase q'(0)
else        # if q(eta) reached qmax or theta'(0) reached 0
QR = df0[2];  # 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.