Stopping ode45 when encountering errors and retrieving results

%% UPDATE: Actual code now attached.

I am using an ode45 solver to solve a nonlinear ODE.

The code uses many functions, so to avoid unneccessary details, I will simply outline it:

say the differential equation is:

%%NOTE: These are not the actual expressions, nor does an error arise with
%%these. But just showing the algorithm of what I am doing.
dx = @(t,x) 2.*x.*peta(x);

and peta(x) calls various functions, say one of them is peta(x)

function [p] = peta(x)
J = eta(x);
p=J+1;
end

and say eta(x) uses fzero in one of the steps

function [spin] = eta(x)
s0=6;
G= @(x) 3.*x - 5;
spin= fzero(G,s0);
end

The error arises because of the behavior of the ODE. after some interval t, the function G in eta evaluated at the initial guess is not finite.

So I get an error in fzero, and the code is terminated. When I break fzero if G(s0) is infinite or complex, J in peta is not assigned and the integrator is also terminated.

How can I stop the integrator and retrieve the solution whenever I encounter this hurdle in fzero?

Answer 1

[Alan Stevens]

The following makes your example work. However, I suspect your real problem is more complicated. If so, you should upload it.

dx = @(t,x) 2.*x.*peta(x);

tspan = [0 1];
x0 = 1;
[t, x] = ode45(dx, tspan, x0);

plot(t,x),grid

function [p] = peta(x)
J = eta(x);
p=J+1;                 %%%%%%%%%%%%%%%
end

function [spin] = eta(x)
G= @(x) 3.*x - 5;
spin= fzero(G,x);      %%%%%%%%%%%%%%%% x rather than s0
end

[Mohammad Farhat]

Thanks Alan. But this code is not the actual code as you can tell, the expressions are nothing like this. I just wanted to show the algorithm to try to avoid these generic errors.

[Alan Stevens]

Can you upload the actual problem?

[Mohammad Farhat]

Hi @Alan . I attached the actual code. the main script is TSM.

You can see that if you run the integrator [ 0 -3.225], it works fine, and you can see the behavior of the function dropping at some t immediately to zero.

this behavior is what causes the error if you rund the integrator for a longer timespan ( say [ 0 -4]).

TSM calls IM_K22 which calls H_total.

when this behavior is reached for x, H_total using fzero gives and error because the function at the initial guess is no longer finite, and thus IM_K22 doesn't work, and thus the integrator stops.

As you might tell, I want to solve this integrator for ranges of parameters, which results in the time around which this behavior occurs to vary, so I want to fix tspan of the integrator and not manually change it every time to avoid the error.

[Alan Stevens]

The problem seems to be related to

G_orbit = @(a)  beta*sqrt(mu.*a);

in two of the files.

If a goes negative you get a complex number here which presumably messes up the remaining calculations. By replacing a by max(a,0) the program doesn't crash, but doesn't seem to end either! I don't know enough about where your equations come from to help further.

Answer 2

[Walter Roberson]

Here are some versions that you can use to get further. Not very far though.

%TSM.m

do_TSM();

function do_TSM
    
    m=0.0123;
    M=1;
    
    
    R=6371; % km
    GM=0.39860e6; GM = GM / R^3;
    Gm=0.00490e6; Gm= Gm / R^3 ;
    mu= Gm+GM;
    Gyr_to_sec= 365.25*24*3600*10^9;
    log_sigma_R= -5;
    
    
    H=2500;
    beta= M*m/(M+m);
    G_orbit = @(a)  beta*sqrt(mu.*a);
    Torque= @(a)  3/2 * Gm *m./a.^6 .* -IM_K222(a, log_sigma_R, H)*Gyr_to_sec;
    
    a0=60.142611;
    
    tend = -3.225;
    tend = -4;
    
    opts = odeset('OutputFcn', @odeplot);
    [t,jj] = ode45(@da, [ 0 tend], a0, opts);
    figure; plot(t,jj); hold on;
    
    
    function dat = da(t,a)
        if a < 0
            dat = flintmax;
        else
          dat = 2.*a./G_orbit(a).*Torque(a);
        end
    end
    
end

%H_total.m

function [spin_rate] =H_total(semi_major_axis)
    
    
    M=1;
    m= 0.0123;
    
    beta= M*m / (M+m);
    R= 6371;
    
    GM=0.39860e6; GM = GM / R^3;
    Gm=0.00490e6; Gm= Gm / R^3 ;
    mu= Gm+GM;
    C = 2/5 * M * R^2;
    
    C= C / R^2;
    
    
    h_spin = @(w) C*w;
    
    
    
    hour= 3600;
    P_rot = 23.9345*hour;
    %P_rot = 14*hour;
    Omega_E = 2*pi/P_rot;
    a0 = 384748/R;
    
    h_orbit_0= h_orbit(a0);
    G0=  h_orbit_0 + h_spin(Omega_E);
    G = @(w) h_orbit(semi_major_axis) + h_spin(w); %a in R, w in s^-1 and G in Earth mass*R^2/s
    G=@(w) G(w) - G0;
    terminator = G(Omega_E);
    %     if isinf(terminator) || ~isreal(terminator)
    %         return
    %     end
    spin_rate = fzero(G,Omega_E);
    
    function horb =  h_orbit(a)
        horb = beta*sqrt(mu*a);
    end
    
end

The change to IM_K222.m was just reformatting it to make it more clear where the nested functions are:

%IM_K222.m

function [Imk22] = IM_K222(semi_major_axis, log_sigma_R, H)
    
    
    log_sigmaR_low=log_sigma_R;
    log_sigmaR_high=log_sigma_R;
    log_sigmaR = log_sigmaR_low:log_sigmaR_high;
    N_sigmaR=length(log_sigmaR);
    sigmaR = 10.^(log_sigmaR);
    
    
    N_max=30;
    
    R= 6378e3; %m
    %m % motoyama: 2600/ 3900 / 5200
    g= 9.81; %m/s^2
    rho= 1022; %kg/m^3
    N = 1e-3; %s^-1
    cs = 1545; %m s^-1
    M_planet=5.9724e24; %kg
    M_moon=	7.346e22; %kg
    hour= 3600;
    P_rot = 23.9345*hour;
    Omega_E = 2*pi/P_rot;
    N2 = N^2;
    cs2 = cs^2;
    M_ocean = 4*pi*(R^2)*H*rho;
    G	= 6.67384e-11;
    
    %the code should take radius in meters
    a_EM= semi_major_axis*R; %m
    n_orb = sqrt(G*(M_planet+M_moon)/(a_EM^3));
    Omega= H_total(a_EM/1e3/(R/1e3));
    
    sigma= 2*(Omega - n_orb);
    Nomega= length(sigma);
    
    
    
    
    m=2;
    NSmin=30;
    NS=max(NSmin, 2*N_max);
    Lambda= zeros(1,N_max);
    Avector= zeros(N_max,N_max);
    Bvector= zeros(N_max,N_max);
    
    C2n2_tab= zeros(Nomega, N_max);
    Lambda_tab= zeros(Nomega, N_max);
    Qn=zeros(Nomega,N_max);%HERE
    Qtot=zeros(Nomega,N_sigmaR);
    
    
    for js=1:N_sigmaR
        nu= 2*Omega./(sigma - 1i*sigmaR(js));
        for kk=1:Nomega
            nuk=nu(kk);
            [Lambda,Avector,Bvector]= Hough(m,nuk,NS,N_max);
            for ik=1:N_max
                C2n2_tab(kk,ik) = Avector(1,ik)*Bvector(ik,1);
            end
            Lambda_tab(kk,:) = Lambda;
            
            for n=1:N_max
                %Qn(kk,n) = QLove_uni(sigma(kk),R,H,g,N2,cs2,sigmaf(js),Lambda(n),1,1);
                Qn(kk,n) = QLove_uni_neutre(sigma(kk),R,H,g,cs2,sigmaR(js),Lambda(n));
            end
            for n=1:N_max
                Qtot(kk,js) = Qtot(kk,js) + Qn(kk,n)*C2n2_tab(kk,n);
            end
            
        end
        ImQtot = imag(Qtot);
        Imk22 = (G*M_ocean)/(5*R).*ImQtot;
        k22  = (G*M_ocean)/(5*R).*Qtot;
        Q= abs( k22./Imk22);
        absImk22 = abs(Imk22);
        %   plot(omega, log10(Q)); hold on;
        
        
    end
    
    
    function [Lambda,Avector,Bvector]= Hough(m,nuk,NS,Nfunc)
        [Lambda_brut,Avector_brut]  = H_even(m,nuk,NS);
        Avector=Avector_brut(1:Nfunc,1:Nfunc);
        Lambda=Lambda_brut(1:Nfunc);
        Bvector_brut=inv(Avector_brut);
        Bvector=Bvector_brut(1:Nfunc,1:Nfunc);
    end
    function [Lambda,Avector] = H_even(m,nuk,NS)
        absm = abs(m);
        S=zeros(NS,NS);
        indices = 1:NS;
        k = absm+2*(indices-1);
        diagg = M_r(k,m,nuk);
        bdiagg = L_r(k(1:NS-1),m,nuk);
        
        S(1,1) = diagg(1);
        
        for j=2: NS
            S(j,j-1) = bdiagg(j-1);
            S(j-1,j) = bdiagg(j-1);
            S(j,j) = diagg(j);
        end
        
        [Avector,X] = eig(S);
        X=diag(X);
        Lambda = 1./(X*nuk.^2);
    end
    function [Mr] = M_r(r,s,nuk)
        rr = 2*r;
        snu = s*nuk;
        rp1 = r+1;
        rp2 = r+2;
        r2 = r.^2;
        
        M1 = -(snu-r.*rp1)./((r.*rp1*nuk).^2);
        M2 = (r+s+1).*(r-s+1).*rp2.^2 ./ (rp1.^2 .* (rr+3).*(rr+1).*(snu-rp1.*rp2));
        M3 = (r-1).^2 .*(r2-s^2)./(r2.*(4.*r2-1).*(snu-r.*(r-1)));
        Mr=M1+M2+M3;
    end
    function [Lr] = L_r(r,s,nuk)
        
        rr = 2*r;
        snu = s*nuk;
        
        num = sqrt( (r+s+1).*(r+s+2).*(r-s+1).*(r-s+2));
        den = (rr+3).*sqrt((rr+1).*(rr+5)).*(snu-(r+1).*(r+2));
        
        Lr= num./den;
    end
    function [Qtot]= QLove_uni(sigma, R,H,g,N2,cs2,sigmaf,Lambdan,op_Qxi,op_Qrho)
        H2=H^2;
        R2=R^2;
        sigmati = sigma - 1i*sigmaf;
        ssti = sigma*sigmati;
        sg2=2*g/R;
        ss2=cs2/R2;
        
        sigmagn2 = g*Lambdan/(2*R);
        Ksph = -2*H/R;
        Ksph=0;
        
        epsn = R2*ssti/(Lambdan*cs2);
        etan = 1 - epsn;
        gHcs2 = g*H/cs2;
        N2Hg = N2*H/g;
        
        tau  = gHcs2 + N2Hg;
        delta = (tau + Ksph)/2;
        gam = delta - tau;
        
        
        kn2 = (Lambdan*H2)/R2 * (N2/ssti - 1)*etan +(gHcs2 + Ksph)*N2Hg - delta^2;
        kn = sqrt(kn2);
        coskn = cos(kn);
        sinkn = sin(kn);
        
        An = (gHcs2 - N2Hg + Ksph)/2;
        Bn = (gHcs2*(2*ssti/N2 - 1) - N2Hg + Ksph)/2;
        Cn = An + H*(N2 - ssti)/g;
        den = kn*(Cn-An)*coskn - (An*Cn +kn2)*sinkn;
        Andel = An - delta;
        Cndel = Cn - delta;
        expd= exp(delta);
        expt = exp(tau);
        expg = exp(gam);
        gam2 = gam^2;
        
        
        Qtot = (1+1i)*0;
        Qxi = Qtot;
        Qrho = Qtot;
        if (op_Qxi==1)
            Qxi0 = - Lambdan/(R2*ssti*(kn2+delta^2));
            Dn = kn*(An-delta)*(An-Cn)*expd;
            En = (Cn-delta)*(An^2+kn2);
            Qxi = Qxi0*(An - delta + (Dn + En*sinkn)/den);
        end
        
        
        if (op_Qrho==1)
            gk2 = gam2 + kn2;
            Bndel = Bn - delta;
            c1 = gam*(Bn-Cn) + Bn*Cn + kn2;
            c2 = gam*(Bn-An)+An*Bn+kn2;
            c3 = gam*(An*Bn+ kn2) + kn2*(An-Bn);
            c4 = gam*(Bn*Cn+kn2) + kn2*(Cn-Bn);
            
            Gn = kn*( Andel*expd*c1 + Cndel*c2/expg );
            Hn = -kn*( expt*Andel*c1 + Cndel*c2 );
            Kn = Cndel*c3 - expt*Andel*c4;
            Qrho0 = N2*H*Lambdan/(g*R2*ssti*(kn2+delta^2));
            Qrho = - Qrho0*((expt-1)*Bndel/tau + (Gn + Hn*coskn + Kn*sinkn)/(den*gk2));
        end
        Qtot = Qxi + Qrho;
    end
    function [Qtot] = QLove_uni_neutre(sigma,R,H,g,cs2,sigmaf,Lambdan)
        R2=R^2;
        kh2 = Lambdan / R2;
        epsgs = (g*H)/(2*cs2);
        Delta = sqrt( (H*g*kh2/(1-epsgs)^2 - (sigmaf/2)^2)*1i/1i);
        sigma1 =  - Delta + 1i*(sigmaf/2);
        sigma2 =   Delta + 1i*sigmaf/2;
        
        Qtot = -Lambdan/( R2) / ((sigma-sigma1)*(sigma-sigma2));
        
    end
    
end

[Walter Roberson]

That is the current code. If you run it, you will see that the values drop down steeply in the plot and then the program will look like it hangs. It does not actually hang: it just switches into a step size on the order of 1E-10 trying to find coefficients to make progress.

The Runge-Kutta process involves evaluating at several places that are carefully constructed with respect to each other in order to match derivatives in endpoints of derivatives of projection models. The value of the function at the probe locations becomes part of the information to determine where the function is probed. In areas where the function is changing rapidly, the probe locations can end up being relatively far away. You can end up probing at locations that yield large positive values, and also probing at locations that yield large negative values, and those might mathematically mostly cancel out, leaving you hardly moved. Balancing a pin on-end.

What can you do? Well, if you set the options to have a larger minimum step size such as 1e-9, then it would give up quickly complaining about a likely singularity... which would be better than hunting for hours near the singularity.

Your code takes eigenvalues of a array with minimum dimension 30, so I am not going to try to do a symbolic analysis to determine theoretical singularities.