各位木虫童鞋，我现在在调试matlab程序遇到了问题。 在关于角度theta的复杂数值积分中，调用我拥有的quadl函数时，出现了虚数。 随后基于变步长的simpson算法编制了数值积分程序，结果有虚数。 而且，我的积分程序和quadl的计算结果还有误差。 积分函数包含很多三角函数。 另外，还有e的指数函数是单身的捉迷藏。 但是，积分的过程中不会带来傅立叶变换等虚数的过程。 我不知道虚数是怎么来的。 试着画了积分函数曲线。 变动很大，为什么作为面积会出现虚数呢？ 贴上我的8个积分函数公式和各自的函数曲线。 让大神们看看，请发表意见。 论文卡在这里，一直解决不了。
%=====functiongroupoftheta====================functiongroupoftheta==========================================
functionF1=ff1(Theta ) ) ) ) ) ) ) ) )。
%======basic formula============================================================================
% toidentifywhethertherotatingaxisinsidethecircle
if rrat=0
RM=0.5*(exp ) (theta-theta0) tan ) phi ) ) rrat*exp ) ) theta-theta0) tan ) phi ) )；
r=0.5*(exp () theta-theta0) tan ) phi ) ) rrat*exp ) ) theta-theta0) tan ) phi ) )
else
RM=0.5*(exp(theta-theta0) tan ) phi ) ) rrat*exp(-) theta-theta0(tan ) phi ) )；
r=0.5*(exp () theta-theta0) tan ) phi ) )-rrat*exp(-() theta-theta0) tan ) phi ) )；
结束
a=sin(theta0)/sin (theta )-rm；
f1=2*cos(Theta ) ) ) ) r^2-a^3-0.5 ) a*RM^2) a*RM^2) r^2-0.25 ) a^3-0.5 ) a*RM^2) r )2) a .
(0.5*pi-Asin(a/r ) () 0.125*r ) 40.5*RM )2*r )2) ) )；
结束
functionF2=ff2(Theta ) ) ) ) ) ) )。
%======basic formula============================================================================
% toidentifywhethertherotatingaxisinsidethecircle
if rrat=0
RM=0.5*(exp ) (theta-theta0) tan ) phi ) ) rrat*exp ) ) theta-theta0) tan ) phi ) )；
r=0.5*(exp () theta-theta0) tan ) phi ) ) rrat*exp ) ) theta-theta0) tan ) phi )；
else
RM=0.5*(exp(theta-theta0) tan ) phi ) ) rrat*exp(-) theta-theta0(tan ) phi ) )；
r=0.5*(exp () theta-theta0) tan ) phi ) )-rrat*exp(-() theta-theta0) tan ) phi ) )；
结束
d=sin(betasthetah )/sin (betas theta ) exp ) ) thetah-theta0) tan(phi ) )-rm；
f2=2*cos(Theta ) ) ) ) r^2-d^3-0.5 ) d*RM^2) r^2) r^2-d ) 0.25 ) d^3-0.5 ) d*RM^2) r^2) .
(0.5*pi-Asin(d/r ) ) *(0.125*R^4 0.5*rm^2*R^2) )
结束
functionG1=gf1(Theta ) )。
%======basic formula=======================basic formula==============================================
% toidentifywhethertherotatingaxisinsidethecircle
if rrat=0
RM=0.5*(exp ) (theta-theta0) tan ) phi ) ) rrat*exp ) ) theta-theta0) tan ) phi ) )；
r=0.5*(exp () theta-theta0) tan ) phi ) ) rrat*exp ) ) theta-theta0) tan ) phi ) )
else
RM=0.5*(exp(theta-theta0) tan ) phi ) ) rrat*exp(-) theta-theta0(tan ) phi ) )；
r=0.5*(exp () theta-theta0) tan ) p
hi))-rrat*exp(-(theta-theta0)*tan(phi)));end
a=sin(theta0)/sin(theta)-rm;
G1=(rm^2*(R-a)+rm*(R^2-a^2)+(1/3)*(R^3-a^3))*cos(theta);
end
function G2=Gf2(theta)
%=========basic formula=============================
% TO identify whether the rotating axis is inside the circle
if rrat<=0
rm=0.5*(exp((theta-theta0)*tan(phi))-rrat*exp((theta-theta0)*tan(phi)));
R=0.5*(exp((theta-theta0)*tan(phi))+rrat*exp((theta-theta0)*tan(phi)));
else
rm=0.5*(exp((theta-theta0)*tan(phi))+rrat*exp(-(theta-theta0)*tan(phi)));
R=0.5*(exp((theta-theta0)*tan(phi))-rrat*exp(-(theta-theta0)*tan(phi)));
end
d=sin(betaS+thetah)/sin(betaS+theta)*exp((thetah-theta0)*tan(phi))-rm;
G2=(rm^2*(R-d)+rm*(R^2-d^2)+(1/3)*(R^3-d^3))*cos(theta);
end
function E1=Ef1(theta)
%=========basic formula=============================
% TO identify whether the rotating axis is inside the circle
if rrat<=0
rm=0.5*(exp((theta-theta0)*tan(phi))-rrat*exp((theta-theta0)*tan(phi)));
R=0.5*(exp((theta-theta0)*tan(phi))+rrat*exp((theta-theta0)*tan(phi)));
else
rm=0.5*(exp((theta-theta0)*tan(phi))+rrat*exp(-(theta-theta0)*tan(phi)));
R=0.5*(exp((theta-theta0)*tan(phi))-rrat*exp(-(theta-theta0)*tan(phi)));
end
a=sin(theta0)/sin(theta)-rm;
E1=(cos(theta)/(sin(theta))^3)*(R^2-a^2)^(1/2);
end
function E2=Ef2(theta)
%========= basic formula=============================
% TO identify whether the rotating axis is inside the circle
if rrat<=0
rm=0.5*(exp((theta-theta0)*tan(phi))-rrat*exp((theta-theta0)*tan(phi)));
R=0.5*(exp((theta-theta0)*tan(phi))+rrat*exp((theta-theta0)*tan(phi)));
else
rm=0.5*(exp((theta-theta0)*tan(phi))+rrat*exp(-(theta-theta0)*tan(phi)));
R=0.5*(exp((theta-theta0)*tan(phi))-rrat*exp(-(theta-theta0)*tan(phi)));
end
d=sin(betaS+thetah)/sin(betaS+theta)*exp((thetah-theta0)*tan(phi))-rm;
E2=cos(theta+betaS)/((sin(theta+betaS))^3)*(R^2-d^2)^(1/2);
end
function H1=Hf1(theta)
%===========basic formula=============================
H1=cos(theta)/(sin(theta))^3;
end
function H2=Hf2(theta)
%========= basic formula=============================
H2=cos(theta+betaS)/(sin(theta+betaS))^3;
end
==========computational parameters==================
d2r=pi/180;
theta0=20*d2r;thetah=100*d2r;rrat=0.6;brat=0.5;
betaS=45*d2r; phi=30*d2r;
A1=sin(betaS+thetah)/sin(theta0)*exp((thetah-theta0)*tan(phi));
B1=cos(betaS);
C1=sin(betaS);
thetaB=acot((A1-B1)/C1);
%===========numerical integration for W and D==============
Wd1=quadl(@Ff1,theta0,thetaB);
Wd2=quadl(@Ff2,thetaB,thetah);
Wd=Wd1+Wd2;
Wp1=quadl(@Gf1,theta0,thetaB);
Wp2=quadl(@Gf2,thetaB,thetah);
Wp=Wp1+Wp2;
Kd1=quadl(@Ef1,theta0,thetaB);
Kd2=quadl(@Ef2,thetaB,thetah);
Dd=-2*cot(phi)*(sin(theta0))^2*Kd1...
-2*cot(phi)*exp(2*(thetah-theta0)*tan(phi))*(sin(thetah+betaS))^2*Kd2;
Kp1=quadl(@Hf1,theta0,thetaB);
Kp2=quadl(@Hf2,thetaB,thetah);
Dp=-cot(phi)*(sin(theta0))^2*Kp1...
-cot(phi)*exp(2*(thetah-theta0)*tan(phi))*(sin(thetah+betaS))^2*Kp2;
disp('[Wd1,Wd2,Wp1,Wp2,Kd1,Kd2,Kp1,Kp2]=');
disp([Wd1,Wd2,Wp1,Wp2,Kd1,Kd2,Kp1,Kp2]');
八个被积函数曲线.jpg