matlab结果中有虚数,Matlab数值积分结果是虚数

最新推荐文章于 2024-01-26 09:26:04 发布

聂瓦

最新推荐文章于 2024-01-26 09:26:04 发布

阅读量4.5k

点赞数 1

文章标签： matlab结果中有虚数

博主在使用MATLAB进行复杂角度theta的数值积分时遇到了结果为虚数的难题，尝试了quadl函数和自编的 Simpson 算法，但结果仍有虚数且存在误差。积分函数涉及三角函数和指数函数，但未进行傅里叶变换等可能引入虚数的操作。博主分享了积分函数代码并寻求帮助，希望能解决这个阻碍论文进展的问题。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

小木虫的各位童鞋们，我现在调试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]');

t-5642162-1-pid-2

八个被积函数曲线.jpg