求助关于使用ode45求解裂纹转子微分方程时的出错

bryant

微分方程中一个系数是开闭函数，就是需要调用这个函数
主函数：
global e xi Omega omega k;%定义全局变量
e=0.1;xi=0.01;Omega=0.5;omega=1;k=0.1;%对参数赋值
[t,u]=ode45('ode3',[0:0.1:50],[1 0 0 0]);
子函数：
function du=ode3(t,u)%解微分方程组函数
global e xi Omega omega k
alpha=pi/6;
beta=pi/6;
theta=omega*t+beta;
ff=fcrack(alpha,theta);
du=zeros(4,1);%u1=x,u2=x',u3=y,u4=y'
du(1)=u(2);
du(2)=ff*k/Omega^2*sin(2*omega*t)/2*u(3)-2*xi/Omega*u(2)-(1/Omega^2-ff*k/Omega^2*cos(omega*t)^2)*u(1)+1/Omega^2+e*cos(omega*t);
du(3)=u(4);
du(4)=ff*k/Omega^2*sin(2*omega*t)/2*u(1)-2*xi/Omega*u(4)-(1/Omega^2-ff*k/Omega^2*sin(omega*t)^2)*u(3)+e*sin(omega*t);

function f=fcrack(alpha,theta)%开闭函数
theta2=mod(theta,2*pi);%令theta2的范围在0-2pi之内
if -pi/2+alpha<=theta2&theta2<=pi/2-alpha
    f=1;
elseif pi/2-alpha<=theta2&theta2<=pi/2+alpha
    f=1/2+1/2*cos((theta2-pi/2+alpha)*pi/(2*alpha));
elseif pi/2+alpha<=theta2&theta2<=3*pi/2-alpha;
    f=0;
elseif 3*pi/2-alpha<=theta2&theta2<=3*pi/2+alpha;
    f=1/2+1/2*cos((theta2-3*pi/2-alpha)*pi/(2*alpha));
end

运行时错误信息：One or more output arguments not assigned during call to 'ode3 (fcrack)'.
请问程序应该怎么修改才正确

[ 本帖最后由 eight 于 2008-5-8 22:10 编辑 ]

ch_j1985

把红色的“;”去掉试试看

if -pi/2+alpha<=theta2&theta2<=pi/2-alpha
    f=1;
elseif pi/2-alpha<=theta2&theta2<=pi/2+alpha
    f=1/2+1/2*cos((theta2-pi/2+alpha)*pi/(2*alpha));
elseif pi/2+alpha<=theta2&theta2<=3*pi/2-alpha;
    f=0;
elseif 3*pi/2-alpha<=theta2&theta2<=3*pi/2+alpha;
    f=1/2+1/2*cos((theta2-3*pi/2-alpha)*pi/(2*alpha));
end

bryant

去掉了，还是那个问题。谢谢关注

[ 本帖最后由 eight 于 2008-5-8 22:00 编辑 ]

ch_j1985

运行之后出现如下错误：

Error in ==> fcrack at 2
theta2=mod(theta,2*pi);%令theta2的范围在0-2pi之内
??? Output argument "f" (and maybe others) not assigned during call to "E:\fcrack.m
(fcrack)".
Error in ==> ode3 at 6
ff=fcrack(alpha,theta);
Error in ==> ode45 at 324
    f(:,2) = feval(odeFcn,t+hA(1),y+f*hB(:,1),odeArgs{:});

应该是开闭函数有问题，是不是应该把问题的有关知识介绍一下？

bryant

开闭函数是转子运转过程中由于裂纹的开闭而对系统刚度的影响，与时间有关系。

[ 本帖最后由 bryant 于 2008-5-8 16:24 编辑 ]

eight

请在置顶帖中找答案，已经不是第一次提醒的了，这个帖子 Undefined function or variable 'e'.出错求助 就叫你看置顶帖，我们整理的成果不是用来摆设的

[ 本帖最后由 eight 于 2008-5-8 22:13 编辑 ]

haitao0318fly

:@o 已经解决了吗？？？？？

快乐

请问，你在建模时用的系数omige是如何取值的?是否用到了求解局部柔度？
本人也在建立裂纹转子的模型，可是在一步求局部柔度的大的二次积分就被难住了，还请各位朋友帮忙解决一下，在此先谢谢了！！！

mjtjiang

回复 快乐 的帖子

可能现在你已经做出来了吧，不过还是说一下我做的。我是用quad2dggen双重积分的，matlab里可能没有这个函数，可以在网上下载NIT工具箱。

玉林

这个开闭函数有问题吧。theta2的取值范围是0～2*pi,但是开闭函数的范围是-pi/3~5*pi/3.