MATLAB Answers

# Stopping ode45 when encountering errors and retrieving results

1 view (last 30 days)
Mohammad Farhat on 19 May 2021
Commented: Walter Roberson on 21 May 2021
%% 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?
##### 0 CommentsShowHide -1 older comments

Sign in to comment.

### Answers (2)

Alan Stevens on 19 May 2021
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
##### 4 CommentsShowHide 3 older comments
Alan Stevens on 19 May 2021
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.

Sign in to comment.

Walter Roberson on 21 May 2021
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
##### 1 CommentShowHide None
Walter Roberson on 21 May 2021
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.

Sign in to comment.

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!