大神求教这个警告是什么原因造成的，如何修改程序？

Rotor2014

大神求教这个警告是什么原因造成的，如何修改程序？
用ode45求解系统方程
for j=1:L
j
[t ,Y]=ode45(@weifenfangcheng,[0:(2*pi/200):(288*2*pi),[0.01 0.001 0.01 0.001 0.01 0.001 0.01 0.001])；
.....
当j=2时出现警告：“在t=1.007596e-2处失败。在时间t处，若不将步长降至允许的最小值（2.775558e-17）一下，积分公差要求无法满足。”
in ode45 at 308
下标索引必须为正整数类型或逻辑类型。

另外，我是在求解方程里面加了一项函数，在不添加前方程可以运行得到结果，由于多考虑其他因素添加了一项后就出现了错误。请问出现这种错误的原因是什么啊，怎么解决啊？

Rotor2014

有没有知道，指点一下

ChaChing

好奇问下, 报错讯息是中文?

另提醒下
求助完整格式：出错代码和出错提示

Rotor2014

ChaChing 发表于 2014-12-29 10:07
好奇问下, 报错讯息是中文?

另提醒下

用的是最新的中文版matlab2014b，报错提示是中文的。

Rotor2014

ChaChing 发表于 2014-12-29 10:07
好奇问下, 报错讯息是中文?

另提醒下

程序如下，主要问题就是我在程序里面未添加Fa这一项之前是可以进行计算得到图，添加Fa之后就出现警告，终止运行。

1. function fangchengfenchatu
   
2. tic
   
3. N=256;
   
4. M=200;
   
5. RelTol = 1e-6;                %相对误差Relative tolerance
   
6. AbsTol = 1e-6;                %绝对误差Absolute tolerance
   
7. options = odeset('RelTol',RelTol,'AbsTol',ones(1,8)*AbsTol);  %需要改变的参数‘ones(1,4)’
   
8. ommiga=200:10:2500;
   
9. L=length(ommiga);
   
10. for j=1:L
    
11. j
    
12.   [t,Y] = ode45( @fangcheng001,[0:(2*pi/M):(288*2*pi)],[0.01 0.001 0.01 0.001 0.01 0.001 0.01 0.001],options,ommiga(j));
    
13.   y1=Y(:,1);
    
14. y11=y1((end-N*M):end);
    
15. 
    
16. for i=1:N
    
17.   
    
18.    Point1(i,j)=y11(1+(i-1)*M);
    
19. end
    
20.   
    
21. end
    
22. 
    
23. ommiga = ommiga';
    
24. 
    
25. plot(ommiga,Point1,'.k','MarkerSize',1)
    
26. 
    
27. toc

复制代码

1. function dy=fangcheng001(t,y,ommiga)  
   
2. %ommiga=450;%转速
   
3. m1=4;
   
4. m2=32.1;
   
5. c1=1050;
   
6. c2=2100;
   
7. c=0.00011;
   
8. R=0.025;
   
9. L=0.012;
   
10. e=0.00005;
    
11. miu=0.018;
    
12. k=2.5*10^7;
    
13. s=miu*ommiga*R*L*(R/c)^2*(L/2/R)^2;
    
14. M=s/(m1*ommiga^2*c);
    
15. %求解Fa;
    
16. xi=1.2e-3;
    
17. Rt=0.5;
    
18. Rb=0.37;
    
19. B1=35/360*2*pi;
    
20. B2=40/360*2*pi;
    
21. xii=0.83;
    
22. V=30;
    
23. p0=11.8;
    
24. C=V^2*sin(B1)*p0*(cos(B1)+xii*B2);
    
25. A1=(Rt^2-Rb^2)^2*pi*C*Rt/(Rt^2-Rb^2+2*Rt*xi)^2;
    
26. A3=3*(Rt^2-Rb^2)^2*pi*C*Rt^3/(Rt^2-Rb^2+2*Rt*xi)^4;
    
27. E=sqrt(y(5)^2+y(7)^2)/xi;
    
28. Fa=A1*E+A3^3*E^3;
    
29. %无量纲方程
    
30. 
    
31. dy=[ y(2);
    
32.     -c1/(m1*ommiga)*y(2)-k/(m1*ommiga^2)*(y(1)-y(5))+fx(y(1),y(3),y(2),y(4))*M;
    
33.     y(4);
    
34.     -c1/(m1*ommiga)*y(4)-k/(m1*ommiga^2)*(y(3)-y(7))-9.8/(c*ommiga^2)+fy(y(1),y(3),y(2),y(4))*M;
    
35.     y(6);
    
36.     -c2/(m2*ommiga)*y(6)-2*k/(m2*ommiga^2)*(y(5)-y(1))+e*cos(t)/c+Fa*cos(t)/(c*m2*ommiga^2);
    
37.     y(8);
    
38.     -c2/(m2*ommiga)*y(8)-2*k/(m2*ommiga^2)*(y(7)-y(3))+e*sin(t)/c-9.8/(c*ommiga^2)+Fa*sin(t)/(c*m2*ommiga^2);];
    
39. end

复制代码

chybeyond

Rotor2014 发表于 2014-12-29 11:07
程序如下，主要问题就是我在程序里面未添加Fa这一项之前是可以进行计算得到图，添加Fa之后就出现警告，终 ...

fx,fy表示什么

mrli1990

新人报到，向前辈致敬！

Rotor2014

chybeyond 发表于 2014-12-30 09:01
fx,fy表示什么

忘记写了  不好意思

1. %油膜力程序
   
2. function out=fx(y1,y3,y2,y4)
   
3. alpha=atan((y3+2*y2)/(y1-2*y4))-pi/2*sign((y3+2*y2)/(y1-2*y4))-pi/2*sign(y3+2*y2);
   
4. G=2*(pi/2+atan((y3*cos(alpha)-y1*sin(alpha))/(sqrt(1-y1^2-y3^2))))/(sqrt(1-y1^2-y3^2));
   
5. V=(2+(y3*cos(alpha)-y1*sin(alpha))*G)/(1-y1^2-y3^2);
   
6. S=(y1*cos(alpha)+y3*sin(alpha))/(1-(y1*cos(alpha)+y3*sin(alpha))^2);
   
7. K=sqrt((y1-2*y4)^2+(y3+2*y2)^2)/(1-y1^2-y3^2);
   
8. out=K*(3*y1*V-sin(alpha)*G-2*cos(alpha)*S);
   
9. end
   
10. function out=fy(y1,y3,y2,y4)
    
11. alpha=atan((y3+2*y2)/(y1-2*y4))-pi/2*sign((y3+2*y2)/(y1-2*y4))-pi/2*sign(y3+2*y2);
    
12. G=2*(pi/2+atan((y3*cos(alpha)-y1*sin(alpha))/(sqrt(1-y1^2-y3^2))))/(sqrt(1-y1^2-y3^2));
    
13. V=(2+(y3*cos(alpha)-y1*sin(alpha))*G)/(1-y1^2-y3^2);
    
14. S=(y1*cos(alpha)+y3*sin(alpha))/(1-(y1*cos(alpha)+y3*sin(alpha))^2);
    
15. K=sqrt((y1-2*y4)^2+(y3+2*y2)^2)/(1-y1^2-y3^2);
    
16. out=K*(3*y3*V+cos(alpha)*G-2*sin(alpha)*S);
    
17. end

复制代码

chybeyond

y1=Y(:,1);  
y11=y1((end-N*M):end);
y1只有一个值，所以y11=y1((end-N*M):end)相当于y11=y1((1-256*200):1)=y1(-51199:1)，
y1下标索引只能是正值或者逻辑数值

Rotor2014

chybeyond 发表于 2014-12-30 12:45
y1=Y(:,1);  
y11=y1((end-N*M):end);
y1只有一个值，所以y11=y1((end-N*M):end)相当于y11=y1((1-256*200 ...

非常感谢您的指导，我没有添加Fa之前是可以得到结果的。还有就是提示“当j=2时出现警告：“在t=1.007596e-2处失败。在时间t处，若不将步长降至允许的最小值（2.775558e-17）一下，积分公差要求无法满足。”应该怎么处理啊？

想变成蝴蝶的鱼

请问楼主解决了没有？