大连电子商务网站建设搜狗站长工具
小木虫的各位童鞋们，我现在调试matlab程序遇到一个问题。一个很复杂的关于角度theta的数值积分，我调用自带的quadl函数，结果出现虚数。我后来又根据变步长的simpson算法编了一个数值积分程序，结果还是有虚数，而且我的积分程序和quadl算出的结果还有一些误差。积分函数包含了很多三角函数，还有e的指数函数在里面，但是积分过程中并无傅里叶变换等可以带来虚数的过程，不知道虚数是怎么来的。我把积分函数曲线都画出来了，波动比较大，但是作为一个面积，怎么会出现虚数呢？我把我的八个积分函数表达式以及各自的函数曲线贴出来，希望大神们帮我看看，给点意见。论文卡在这里，一直得不到解决。
%========function group of theta=====================
function F1=Ff1(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;
F1=2*cos(theta)*(((R^2-a^2)^(1/2))*(0.125*R^2*a-0.25*a^3-0.5*a*rm^2+(R^2-a^2)*rm)...
+(0.5*pi-asin(a/R))*(0.125*R^4+0.5*rm^2*R^2));
end
function F2=Ff2(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;
F2=2*cos(theta)*(((R^2-d^2)^(1/2))*(0.125*R^2*d-0.25*d^3-0.5*d*rm^2+(R^2-d^2)*rm)...
+(0.5*pi-asin(d/R))*(0.125*R^4+0.5*rm^2*R^2));
end
function G1=Gf1(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;
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