二阶微分方程组求解

花梨树

求解此方程组。本人解出来的全是直线，望各位高手帮助！

VibrationMaster

程序写错了? 如果短可以贴出来

花梨树

回复 2 # VibrationMaster 的帖子

M文件函数：

1. function dydt=fun3(t,y)
   
2. q1=0.032;
   
3. q2=0.032;
   
4. q3=0.048;
   
5. q4=2.353*10^-3;
   
6. q5=2.174*10^-3;
   
7. q6=5.262*10^-3;
   
8. q7=0.0029;
   
9. q8=0.0019;
   
10. q9=0.0088;
    
11. l=0.178;
    
12. v1=7000;
    
13. v2=6500;
    
14. v3=7000;
    
15. v4=7500;
    
16. u1=2.842*10^-5*(v1^2+v2^2+v3^2+v4^2);
    
17. u2=2.842*10^-5*(v4^2-v2^2)*l;
    
18. u3=2.842*10^-5*(v3^2-v1^2)*l;
    
19. u4=3.232*10^-7*(-v1^2+v2^2-v3^2+v4^2);
    
20. m=0.3;
    
21. g=9.81;
    
22. dydt=zeros(12,1);
    
23. dydt(1)=y(2);
    
24. dydt(3)=y(4);
    
25. dydt(5)=y(6);
    
26. dydt(7)=y(8);
    
27. dydt(9)=y(10);
    
28. dydt(11)=y(12);
    
29. dydt(2)=((sin(y(9))*cos(y(11))*cos(y(7))+sin(y(11))*sin(y(7)))*u1-q1*y(2))/m;
    
30. dydt(4)=((sin(y(9))*sin(y(11))*cos(y(7))-cos(y(11))*sin(y(7)))*u1-q2*y(4))/m;
    
31. dydt(6)=(cos(y(9))*cos(y(7))*u1-q3*y(6)-m*g)/m;
    
32. dydt(8)=(u2+(q5-q6)*y(10)*y(12)-q7*y(8))/q4;
    
33. dydt(10)=(u3+(q4-q6)*y(8)*y(12)-q8*y(10))/q5;
    
34. dydt(12)=(u4+(q4-q5)*y(8)*y(10)-q9*y(12))/q6;
    
35. end

复制代码

调用m文件函数：

1. [t,y]=ode45('fun4',[0:0.01:5],[0 0 0 0 0 0 0 0 0 0 0 0]);
   
2. plot(t,y(5,1),'.k',t,y(6,1),'.r',t,y(7,1),':c',t,y(8,1),':b')
   
3. title('不需要标题吧 ');
   
4. xlabel('t');
   
5. ylabel('y');

复制代码

shengshengchina

1. function s
   
2. [t,y]=ode45(@fun,[0:0.01:5],[0 0 0 0 0 0 0 0 0 0 0 0]);
   
3. size(t),size(y)
   
4. plot(t,y(:,5),'.k',t,y(:,6),'.r',t,y(:,7),':c',t,y(:,8),':b')
   
5. title('不需要标题吧');
   
6. xlabel('t');
   
7. ylabel('y');
   
8. function dydt=fun(t,y)
   
9. q1=0.032;
   
10. q2=0.032;
    
11. q3=0.048;
    
12. q4=2.353*10^-3;
    
13. q5=2.174*10^-3;
    
14. q6=5.262*10^-3;
    
15. q7=0.0029;
    
16. q8=0.0019;
    
17. q9=0.0088;
    
18. l=0.178;
    
19. v1=7000;
    
20. v2=6500;
    
21. v3=7000;
    
22. v4=7500;
    
23. u1=2.842*10^-5*(v1^2+v2^2+v3^2+v4^2);
    
24. u2=2.842*10^-5*(v4^2-v2^2)*l;
    
25. u3=2.842*10^-5*(v3^2-v1^2)*l;
    
26. u4=3.232*10^-7*(-v1^2+v2^2-v3^2+v4^2);
    
27. m=0.3;
    
28. g=9.81;
    
29. dydt=zeros(12,1);
    
30. dydt(1)=y(2);
    
31. dydt(3)=y(4);
    
32. dydt(5)=y(6);
    
33. dydt(7)=y(8);
    
34. dydt(9)=y(10);
    
35. dydt(11)=y(12);
    
36. dydt(2)=((sin(y(9))*cos(y(11))*cos(y(7))+sin(y(11))*sin(y(7)))*u1-q1*y(2))/m;
    
37. dydt(4)=((sin(y(9))*sin(y(11))*cos(y(7))-cos(y(11))*sin(y(7)))*u1-q2*y(4))/m;
    
38. dydt(6)=(cos(y(9))*cos(y(7))*u1-q3*y(6)-m*g)/m;
    
39. dydt(8)=(u2+(q5-q6)*y(10)*y(12)-q7*y(8))/q4;
    
40. dydt(10)=(u3+(q4-q6)*y(8)*y(12)-q8*y(10))/q5;
    
41. dydt(12)=(u4+(q4-q5)*y(8)*y(10)-q9*y(12))/q6;

复制代码

花梨树

回复 4 # shengshengchina 的帖子

运行结果还是不太理想啊，当改变参数V1，V2，V3，V4（可变量）.使其v1+v3=v2+v4，则输出的Z（y（3））就会上下波动，应该是一直增大或者保持水平，除非当小于某个值时才会减小。还有我看我们的主函数不一样，而你写的主函数上的 size(t),size(y)，这句什么意思、什么作用啊，还是不明白。望给与解答，谢谢你前面的帮助啊！

shengshengchina

我的size(t),size(y)，没什么用，刚开始我是想检查你的数据来着，后来发现没什么问题，但是忘了删去这语句。你说的z变量，我都不知道在哪，我也不知道你说的Z（y（3））数据是指哪个数据，另外你说的它为什么要一直增大或者保持水平，我就更不知道为什么了。

shr20604043

aaaaaaaaaaaaaa

luotao1014

{:{13}:}

boyboycom

数值解，初值很主要！