Instability with DifferentialEquations.jl


#1

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);
            
            household_premium = (a_ - a)/f[2] + delta - delta_ + risk_premium;
                   
            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;
    rk = r + risk_premium;
    r_k = r + household_premium;
    
    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

#2

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