1 view (last 30 days)

Show older comments

%% 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?

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

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.

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

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.

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

Start Hunting!