为什么结果全是Nan, 而且出现Warning: Infinite or Not-a-Number value encountered.

labmat

我的程序其实是一个模型, 如果lnt的值取得很小,比如lnt=15000或者以下,这个程序就没有问题,但是我想要的是lnt=41300,这样一来就出现如题的警告,并且结果都是Nan. 好奇怪,希望大家帮我想想原因.

原程序如下:
function loss()
clear
clc

sigma0=0.1:1:200;
for k=1:length(sigma0)
k
Q3=quadl(@(t)Integ1(t,sigma0(k)),0,1);
Q4=quadl(@(t)Integ2(t,sigma0(k)),0,1);
factor(k)=Q3/Q4
end
plot(sigma0,factor)
end

function y=hanshu1(z,k3,k2,sigma,beta,lnt)
y=abs((k3/k2)*sigma*(1-cosh(beta*z)/cosh(beta*lnt/2)));
end

function y=hanshu2(z,k3,k2,sigma,beta,Rnt,lnt)
y=abs(beta*(Rnt/2)*(k3/k2)*sigma*sinh(beta*z)/cosh(beta*lnt/2));
end

function y=Integ1(t,sigma0)

lnt=15000;
Rnt=15;
Gg=440;
Eg=1030;
Ers=3;
vrs=0.3;
Rsh=17;
thickness=10;
d=20;
tauntc=0.16;
elta=2.5;
la=50000;
sigma=sigma0*sin(2*pi*t);

K=Ers*elta/((3*Ers+3*elta)-vrs*(6*Ers+6*elta));
Grs=3*elta*K*Ers/(9*K*Ers+9*K*elta-Ers*elta);
%K=-Ers*(Ers+elta)/(6*(2*Ers+elta)*(vrs-0.5));
%Grs=3*Ers*K*(Ers+elta)/(9*K*(2*Ers+elta)-Ers*(Ers+elta));

alpha=thickness/Rnt;
Eeq=2*alpha*Eg;
Geq=(1-(1-alpha)^4)*Gg;
lambda=(d+2*Rnt)^2;
k3=lambda/(lambda-pi*Rsh^2);
k2=(pi*Rnt^2/lambda)*k3+Ers/Eeq;
k1=(Ers/(2*Grs))*(Rsh-Rnt)*Rnt/lambda;
beta=sqrt(k2/(k1*lambda));
c=beta*Rnt*k3/(2*k2);
sigmas=tauntc/(c*tanh(beta*lnt/2));
sigmaf=tauntc/(c*(1-sech(beta*lnt/2))/(beta*lnt/2));
sigmas1=(k3/k2)*sigmas*(1-tanh(beta*lnt/2)/(beta*lnt/2));
sigmaf1=(k3/k2)*sigmaf*(1-tanh(beta*lnt/2)/(beta*lnt/2));
Vnt=pi*Rnt^2*lnt;
Vsh=pi*(Rsh^2-Rnt^2)*lnt;
Vrs=(d+2*Rnt)^2*la-pi*Rnt^2*lnt;
taus1=c*sigmas*(1-sech(beta*lnt/2))/(beta*lnt/2);
tauf1=tauntc;

for k=1:length(t)
if(sigma(k)<=sigmas)
leff(k)=0;
Q1=quadl(@(z)hanshu2(z,k3,k2,sigma(k),beta,Rnt,lnt),0,lnt/2);
taunt1(k)=(2/lnt)*Q1;
Q2=quadl(@(z)hanshu1(z,k3,k2,sigma(k),beta,lnt),0,lnt/2);
sigmant1(k)=(2/lnt)*Q2;
end

if((sigmas<sigma(k))&(sigma(k)<=sigmaf))
leff(k)=lnt*(sigma(k)-sigmas)/(sigmaf-sigmas);
taunt1(k)=taus1+(tauf1-taus1)*(sigma(k)-sigmas)/(sigmaf-sigmas);
sigmant1(k)=sigmas1+(sigmaf1-sigmas1)*(sigma(k)-sigmas)/(sigmaf-sigmas);
end

if(sigma(k)>sigmaf)
leff(k)=lnt;
taunt1(k)=tauntc;
sigmant1(k)=sigmaf1;
end
end
sigmars1=(sigma*(d+2*Rnt)^2-pi*Rnt^2*sigmant1)/((d+2*Rnt)^2-pi*Rsh^2);
epsilongrs1=sigmars1/Ers;
epsilongnt1=sigmant1/Eeq;
dWnt=tauntc*2*pi*Rnt*leff.^2.*(epsilongrs1-epsilongnt1);
dWeq=Vsh*taunt1.^2/(6*Grs);
y=dWnt;
end

function y=Integ2(t,sigma0)
lnt=15000;
Rnt=15;
Gg=440;
Eg=1030;
Ers=3;
vrs=0.3;
Rsh=17;
thickness=10;
d=20;
tauntc=0.16;
elta=2.5;
la=50000;
sigma=sigma0*sin(2*pi*t);

K=Ers*elta/((3*Ers+3*elta)-vrs*(6*Ers+6*elta));
Grs=3*elta*K*Ers/(9*K*Ers+9*K*elta-Ers*elta);
% K2=-Ers*(Ers+elta)/(6*(2*Ers+elta)*(vrs-0.5));
% Grs2=3*Ers*K2*(Ers+elta)/(9*K2*(2*Ers+elta)-Ers*(Ers+elta));

alpha=thickness/Rnt;
Eeq=2*alpha*Eg;
Geq=(1-(1-alpha)^4)*Gg;
lambda=(d+2*Rnt)^2;
k3=lambda/(lambda-pi*Rsh^2);
k2=(pi*Rnt^2/lambda)*k3+Ers/Eeq;
k1=(Ers/(2*Grs))*(Rsh-Rnt)*Rnt/lambda;
beta=sqrt(k2/(k1*lambda));
c=beta*Rnt*k3/(2*k2);
sigmas=tauntc/(c*tanh(beta*lnt/2));
sigmaf=tauntc/(c*(1-sech(beta*lnt/2))/(beta*lnt/2));
sigmas1=(k3/k2)*sigmas*(1-tanh(beta*lnt/2)/(beta*lnt/2));
sigmaf1=(k3/k2)*sigmaf*(1-tanh(beta*lnt/2)/(beta*lnt/2));
Vnt=pi*Rnt^2*lnt;
Vsh=pi*(Rsh^2-Rnt^2)*lnt;
Vrs=(d+2*Rnt)^2*la-pi*Rnt^2*lnt;
taus1=c*sigmas*(1-sech(beta*lnt/2))/(beta*lnt/2);
tauf1=tauntc;

for k=1:length(t)
if(sigma(k)<=sigmas)
leff(k)=0;
Q1=quadl(@(z)hanshu2(z,k3,k2,sigma(k),beta,Rnt,lnt),0,lnt/2);
taunt1(k)=(2/lnt)*Q1;
Q2=quadl(@(z)hanshu1(z,k3,k2,sigma(k),beta,lnt),0,lnt/2);
sigmant1(k)=(2/lnt)*Q2;
end

if((sigmas<sigma(k))&(sigma(k)<=sigmaf))
leff(k)=lnt*(sigma(k)-sigmas)/(sigmaf-sigmas);
taunt1(k)=taus1+(tauf1-taus1)*(sigma(k)-sigmas)/(sigmaf-sigmas);
sigmant1(k)=sigmas1+(sigmaf1-sigmas1)*(sigma(k)-sigmas)/(sigmaf-sigmas);
end

if(sigma(k)>sigmaf)
leff(k)=lnt;
taunt1(k)=tauntc;
sigmant1(k)=sigmaf1;
end
end
sigmars1=(sigma*(d+2*Rnt)^2-pi*Rnt^2*sigmant1)/((d+2*Rnt)^2-pi*Rsh^2);
epsilongrs1=sigmars1/Ers;
epsilongnt1=sigmant1/Eeq;
W=(sigmant1.^2/(2*Eeq)+taunt1.^2/(2*Geq))*Vnt+Vsh*taunt1.^2/(6*Grs)+Vrs*sigmars1.^2/(2*Ers);
y=W;

end