ode中是否应该选择每一次计算的最后一行值作为下一次计算初值

chunshui2003

选择最后值作为下次初值

直接采用初值初值

上面两幅图分别为采用每一次计算的最后一行值作为下一次计算初值，以及使用最开始初值得到的分岔计算结果，有很大差异。
应该采用哪一种方式是更加合理的选择呢？

gghhjj

这个个人认为应该有两种可能
（1）系统本身就是对初值极为敏感的系统，不同的初值就可能导致系统解相差很大，这个你可以在多取几个初值算一下看看是否属于这种情况；
（2）可能庞加莱截面选取不一致产生的结果，你可以把计算稳定后的相图做出来比较一下看看是否属于这种情况；
（3）计算程序上存在问题；
目前个人只能想到这三种可能性

chunshui2003

回复 2 # gghhjj 的帖子

学长你好，按照你的建议我又重新算了一下：
1 采用默认精度reltol=1e-3,abstol=1e-6。初值分别采用0.001、0.005、0.01和0.05。对于y0保持不变，记作“直接结果”；对于y0在每一次循环后变化，记作“间接结果”。当采用直接结果时，系统在y0=0.001、0.01和0.05时规律基本相同，只是幅值上有差异，而y0=0.005与其他三者均不同，规律类似于间接结果；当采用间接结果时，y0的取值对分岔结果没有实质影响，规律几乎相同甚至幅值也相似。
2 这次我把控制参数的变化范围扩大到原来的3倍多（虽然实际上不可能出现），然后分别绘制其直接和间接结果的分岔图以及cz=0.0040时的轨迹图和poincare映射图。

直接分岔图.bmp

间接分岔图.bmp

cz=0.00400_直接poincare.bmp

cz=0.00400_直接轨迹.bmp

cz=0.00400_间接轨迹.bmp

cz=0.00400_间接poincare.bmp

chunshui2003

程序如下，和在本版块http://forum.vibunion.com/thread-98188-1-1.html中的类似，区别在于这里采用内嵌函数，麻烦学长给看一下：

1. 
   
2. 
   
3. function  [time,y] = qiujie
   
4. tic
   
5. period = 2*pi;               %%%因为对时间t采用了无量纲化，cos(oemga*t)变成了cos(T)，所以周期变为2*pi/1，暂时这样
   
6. tspan = (0:period/100:1000*period);
   
7. y0 = [0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001];
   
8. 
   
9. e0 = 0.6*10^-3; Lb = 100*10^-3;         %Situ_H_imp
   
10.       
    
11. cz = 0.0002:0.0001:0.0200;
    
12. 
    
13. options = odeset('Reltol',1e-3,'Abstol',1e-6);
    
14. 
    
15. row = length(cz);
    
16. column = length(80100:100:100001);
    
17. 
    
18. U = zeros(row,column);
    
19. 
    
20. for i = 1:length(cz)
    
21.     disp(cz(i))
    
22. [time,y] = ode15s(@equ_elastic_without_UMP,tspan,y0,options,e0,cz(i),Lb);
    
23. y0 = y(end,:);
    
24. 
    
25. U(i,:) = y(80100:100:end,1)';
    
26. 
    
27. plot(y(80100:end,1),y(80100:end,3));
    
28. saveas(gcf,sprintf('Lb=1.00,e0=0.0006,without_UMP,cz=%6.5f_traj_Situ_H_imp.bmp',cz(i)),'bmp');             % Situ_H
    
29. 
    
30. 
    
31. 
    
32. plot(y(80100:100:end,1),y(80100:100:end,3),'r.','markersize',5);
    
33. saveas(gcf,sprintf('Lb=1.00,e0=0.0006,without_UMP,cz=%6.5f_poincare_Situ_H_imp.bmp',cz(i)),'bmp');             % Situ_H
    
34. 
    
35. 
    
36. end
    
37. 
    
38. plot(cz,U,'r.','markersize',5);
    
39. saveas(gcf,'分岔图.bmp','bmp');
    
40. toc
    
41. 
    
42. 
    
43. 
    
44. 
    
45. %%%%%%%%         运动微分方程          %%%%%%%%%
    
46. 
    
47. function  uu = equ_elastic_without_UMP(T,u,e0,cz,Lb)
    
48. 
    
49. m1 = 60;               %转子质量 kg
    
50.   m2 = 25;               %轴颈质量 kg
    
51.   c1 = 4000;             %转子阻尼 N.s/m
    
52.   c2 = 1200;             %轴承阻尼 N.s/m
    
53.   Ke = 6.2*10^6;          %转轴刚度 N/m
    
54. %   e0 = 0.6*10^-3;        %转子质量偏心 m
    
55. %   Rr = 60*10^-3;         %转子半径 m
    
56. %   Lr = 150*10^-3;        %转子长 m
    
57.   delta_0 = 4.5*10^-3;   %均匀气隙大小 m
    
58. %   cz = 0.2*10^-3;        %轴颈间隙 m
    
59.   Rb = 50*10^-3;         %轴承半径 m
    
60. %   Lb = 20*10^-3;        %轴承长 m
    
61.   mu = 18*10^-3;         %绝对润滑油粘度 无单位
    
62. 
    
63.   
    
64.   omega = 13;            %大轴旋转角速度
    
65.   
    
66. %   global cz
    
67.   
    
68.   
    
69.   sigma = mu * omega * Rb * Lb * (Rb/cz)^2 * (Lb/2/Rb)^2;        %Sommerfeld修正系数
    
70.   
    
71.   A1 = u(7) + 2*u(6);
    
72.   A2 = u(5) - 2*u(8);
    
73.   
    
74.   alpha = atan(A1./A2) - pi/2* sign(A1./A2) - pi/2* sign(A1);
    
75.   
    
76. E = sqrt(u(5).^2 +u(7).^2);                     %轴颈轴心轨迹
    
77. E_minus = sqrt(1 -E.^2);
    
78. 
    
79. B1 = u(7)*cos(alpha) - u(5)*sin(alpha);
    
80. B2 = u(5)*cos(alpha) + u(7)*sin(alpha);
    
81.   
    
82.   G = 2./E_minus * ( pi/2 + atan(B1./ E_minus) );
    
83.   V = (2 + B1 * G) ./  E_minus.^2;
    
84.   S = B2 ./ ( 1 - B2.^2 );
    
85.   
    
86.   F_oil_same = -sqrt( A2.^2 + A1.^2 ) ./ E_minus.^2;
    
87.   
    
88.   f_x = F_oil_same * ( 3* u(5) *V - sin(alpha) * G - 2 *cos(alpha)*S );
    
89.   f_y = F_oil_same * ( 3* u(7) *V + cos(alpha) * G - 2 *sin(alpha)*S );     %油膜力的无量纲表达式
    
90.   
    
91.   fx = sigma * f_x;
    
92.   fy = sigma * f_y;              %非线性油膜力
    
93.   
    
94.   uu = zeros(8,1);
    
95.   uu(1) = u(2);
    
96.   uu(2) = -c1/(m1*omega) * u(2) - Ke/(m1*omega^2) * u(1) + Ke * cz/(m1*omega^2*delta_0)*u(5) + e0 * cos(T)/delta_0 ;
    
97.   uu(3) = u(4);
    
98.   uu(4) = -c1/(m1*omega) * u(4) - Ke/(m1*omega^2) * u(3) + Ke * cz/(m1*omega^2*delta_0)*u(7) + e0 * sin(T)/delta_0 ;
    
99.   uu(5) = u(6);
    
100.   uu(6) = -c2/(m2*omega) * u(6) + Ke * delta_0/(2*m2*omega^2*cz) * u(1) - Ke/(2*m2*omega^2) * u(5) + fx/(m2*omega^2 *cz);
     
101.   uu(7) = u(8);
     
102.   uu(8) = -c2/(m2*omega) * u(8) + Ke * delta_0/(2*m2*omega^2*cz) * u(3) - Ke/(2*m2*omega^2) * u(7) + fy/(m2*omega^2 *cz);
     
103. 
     

复制代码

octopussheng

可看看龙格库塔法的计算过程，肯定要用到前一步的计算值的。

chunshui2003

回复 5 # octopussheng 的帖子

谢谢学长的回答，我想了想确实应该是这样。毕竟最后的结果是稳态的。
另外，再请教一下，根据上面的最后的轨迹图和poincare图能否判断此时系统处于混沌状态。因为看着映射图像，但轨迹图又好像是拟周期运动。

octopussheng

不像是混沌运动的庞加莱映射。
也不像拟周期，因为庞加莱映射不封闭

chunshui2003

回复 7 # octopussheng 的帖子

学长你好，这个东西我已经弄了好久了。但正如你所说，映射图和轨迹图总也对不上，我作图的部分应该是没有问题吧

1. 
   
2. tspan = (0:period/100:1000*period);
   
3. plot(y(80100:end,1),y(80100:3))          % 轴心轨迹图
   
4. plot(y(80100:100:end,1),y(80100:100:end,3))        % poincare映射图
   

复制代码

一个周期划分100个点，1000个周期，选取稳态后的积分点进行绘制。我现在有点糊涂了，难道是选取点出了问题吗？
请学长帮忙分析一下原因。

octopussheng

plot(y(80100:end,1),y(80100:3)) 似有问题

是否应该是plot(y(80100:end,1),y(80100:end,3))

chunshui2003

回复 9 # octopussheng 的帖子

这里是笔误。忘记end了，y(1)是x1,y(3)是y1。另外，学长，希望你看一下我发的另一个帖子，帮忙判断一下哪种分岔程序更为合理。
http://forum.vibunion.com/thread-98401-1-1.html

octopussheng

绘图部分应是没有问题的。