ode数值积分出现复数解是怎么造成的

chunshui2003

系统运动微分方程

1. 
   
2. 
   
3. 
   
4. function  uu = equ_6_11(T,u)
   
5.   m1 = 1.44*10^5;               
   
6.   m2 = 1.0*10^4;              
   
7.   c1 = 4.3*10^5;            
   
8.   c2 = 3*10^5;            
   
9.   Ke = 5.0*10^7;         
   
10.   
    
11.   e0 = 0.8*10^-3;                                                                                
    
12.   
    
13.   Rr = 0.6;        
    
14.   Lr = 0.3;        
    
15.   delta_0 = 2.0*10^-3;  
    
16.   
    
17.   cz = 0.4*10^-3;        
    
18.   
    
19.   Rb = 0.2;         
    
20.   Lb = 0.3;                                    
    
21.   
    
22.   mu = 24*10^-3;        
    
23.   mu_0 = 4*pi*10^-7;   
    
24.   kj = 0.8;              
    
25. 
    
26. global Ij
    
27.   Fj = kj * Ij;            
    
28.   kr = 5*10^8;        
    
29.   f = 0.01;            
    
30.   
    
31.   omega = 13;            
    
32.   e = sqrt(u(1).^2 + u(3).^2);                          
    
33.   F_coeffi = Rr * Lr *pi *mu_0 *Fj^2 /delta_0.^2;
    
34.   Fx = F_coeffi * (1/2 * e + 5/8 * e.^3) * u(1)./e;
    
35.   Fy = F_coeffi * (1/2 * e + 5/8 * e.^3) * u(3)./e;      
    
36.                                                                                                                                                                            
    
37.   sigma = mu * omega * Rb * Lb * (Rb/cz)^2 * (Lb/2/Rb)^2;      
    
38.   
    
39.   A1 = u(7) + 2*u(6);
    
40.   A2 = u(5) - 2*u(8);
    
41.   
    
42.   alpha = atan(A1./A2) - pi/2* sign(A1./A2) - pi/2* sign(A1);
    
43.   
    
44. E = sqrt(u(5).^2 +u(7).^2);                    
    
45. E_minus = sqrt(1 -E.^2);
    
46. 
    
47. B1 = u(7)*cos(alpha) - u(5)*sin(alpha);
    
48. B2 = u(5)*cos(alpha) + u(7)*sin(alpha);
    
49.   
    
50.   G = 2./E_minus * ( pi/2 + atan(B1./ E_minus) );
    
51.   V = (2 + B1 * G) ./  E_minus.^2;
    
52.   S = B2 ./ ( 1 - B2.^2 );
    
53.   
    
54.   F_oil_same = -sqrt( A2.^2 + A1.^2 ) ./ E_minus.^2;
    
55.   
    
56.   f_x = F_oil_same * ( 3* u(5) *V - sin(alpha) * G - 2 *cos(alpha)*S );
    
57.   f_y = F_oil_same * ( 3* u(7) *V + cos(alpha) * G - 2 *sin(alpha)*S );   
    
58.   
    
59.   fx = sigma * f_x;
    
60.   fy = sigma * f_y;              
    
61.   
    
62. H = 1/2 * (sign(abs(e-1)) + sign(e-1));
    
63. Fx_rub = -delta_0 * (e-1)*kr./e * (u(1) - f.* u(3))*H;
    
64. Fy_rub = -delta_0 * (e-1)*kr./e * (f.* u(1) + u(3))*H;
    
65.    
    
66.   
    
67.   
    
68.   uu = zeros(8,1);
    
69.   uu(1) = u(2);
    
70.   uu(2) = -c1/(m1*omega) * u(2) - Ke/(m1*omega^2) * (u(1) - u(5)) + ...
    
71.       e0 * cos(T)/delta_0 + Fx/(m1*omega^2*delta_0)  + Fx_rub/(m1*omega^2*delta_0) ;
    
72.   uu(3) = u(4);
    
73.   uu(4) = -c1/(m1*omega) * u(4) - Ke/(m1*omega^2) * (u(3) - u(7)) + ...
    
74.       e0 * sin(T)/delta_0 + Fy/(m1*omega^2*delta_0)  + Fy_rub/(m1*omega^2*delta_0);
    
75.   uu(5) = u(6);
    
76.   uu(6) = -c2/(m2*omega) * u(6) - Ke/(2*m2*omega^2) * (u(5) - u(1)) + fx/(m2*omega^2*delta_0);
    
77.   uu(7) = u(8);
    
78.   uu(8) = -c2/(m2*omega) * u(8) - Ke/(2*m2*omega^2) * (u(7) - u(3)) + fy/(m2*omega^2*delta_0);
    
79.   
    
80.   

复制代码

求解

1. 
   
2. 
   
3. period = 2*pi;
   
4. step = period/200;
   
5. tspan = 0:step:period;
   
6. global Ij
   
7. Ij = 600;
   
8. s = [1.101267007 -0.117826981 2.501571845 0.837573562 0.360137403 -0.017557353 1.715633159 0.535182987];
   
9.   
   
10.     y0 = [s(1) s(2) s(3) s(4) s(5) s(6) s(7) s(8)];
    
11.     [t,x] = ode45('equ_6_11',tspan,y0);
    
12. 
    

复制代码

在此初值条件下得到的结果为复数，首次遇到这样的情况。如改用其他初值，如s = [0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1]则可以顺利求解。
然而，这一组初值是通过迭代获得，无法更改，其结果对后续运算也会产生影响。
请教大家，为什么会出现复数解，应采用何种方法加以解决。

jxldc.com

代码没错呀，这多试几次看看。。

chunshui2003

回复 2 # jxldc.com 的帖子

代码没有毛病。问题主要是在采用第一组初值时，得到的结果是复数的，不知道这是什么原因造成的。

www.zdscj.com|

感觉 应该是没什么问题的