Instability with DifferentialEquations.jl

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

1 Like
3

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

4

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
5

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.