感谢+分享：四自由度微分方程组matlab求解

vitihl

首先要感谢论坛里各位帮助过我的高手们，特别是siyanger关于非线性微分方程的求解的帖子， ch_j1985 在我最绝望的时候给与的帮助，以及版主eight一直以来的的严厉要求下，终于自己完成了，我的关于一个四自由度非线性方程组的求解，感谢所有帮助过我的，给我回复的所有人。现把程序贴出来，给以后需要的朋友借以参考，程序还有很多的不足之处，希望高手们继续给我提出修改的意见和建议。
不多说了，来点儿实际的！

1. function f=vdp_eq(t,y,flag,m1,m2,l,c1,c2,c3,k1,k2,k3,e1,e2,e3,e4,e11,e12,e21,e22,w,j1,j2,F0)
   
2. %y(1),y(3),y(5),y(7)为位移，y(2),y(4),y(6),y(8)为速度
   
3. f=[y(2);
   
4. (-c1*(y(2)-y(6))*(j1+e11^2*m1)-c2*(y(4)-y(8))*(j1-e11*e12*m1)-k1*(y(1)-y(5))*(j1+e11^2*m1)-k2*(y(3)-y(7))*(j1-e11*e12*m1))*2*l/((e11+e12)*m1*j1);
   
5. y(4);
   
6. (-c1*(y(2)-y(6))*(j1-e11*e12*m1)-c2*(y(4)-y(8))*(j1+e12^2*m1)-k1*(y(1)-y(5))*(j1-e11*e12*m1)-k2*(y(3)-y(7))*(j1+e12^2*m1))*2*l/((e11+e12)*m1*j1);
   
7. y(6);
   
8. (-(c1*(y(6)-y(2))+c3*y(6))*(j2+e21^2*m2)-(c2*(y(8)-y(4))+c3*y(8))*(j2-e21*e22*m2)-(k1*(y(5)-y(1))+k3*y(5))*(j2+e21^2*m2)-(k2*(y(7)-y(3))+k3*y(7))*(j2-e21*e22*m2)+(2*j2+e21*m2*(e4-e3))*F0*sin(w*t))*2*l/((e21+e22)*m2*j2);
   
9. y(8);
   
10. (-(c1*(y(6)-y(2))+c3*y(6))*(j2-e21*e22*m2)-(c2*(y(8)-y(4))+c3*y(8))*(j2+e22^2*m2)-(k1*(y(5)-y(1))+k3*y(5))*(j2-e21*e22*m2)-(k2*(y(7)-y(3))+k3*y(7))*(j2+e22^2*m2)+(2*j2+e22*m2*(e4-e3))*F0*sin(w*t))*2*l/((e21+e22)*m2*j2)];
    
11. end

复制代码

以上存为M文件

在matlab内运行

1. options=odeset('RelTol',1e-6,'AbsTol',[1e-6]);
   
2. 
   
3. %%%%Parameter make
   
4. t_end=2000; m1=3000; m2=2065; l=660; c1=5; k1=260; c2=5; k2=260; c3=30;  k3=150;
   
5. e1=0; e2=0; e3=189+e2; e4=189-e2; e11=l+e1; e12=l-e1; e21=l+e2; e22=l-e2;
   
6. w=0.143; j1=5.194e+6; j2=8.477e+6; F0=37.5;
   
7. %%%%%%%%%%%%%
   
8. 
   
9. [t,y]=ode45('vdp_eq',[0 t_end],[0 0 0 0 0 0 0 0],options,m1,m2,l,c1,c2,c3,k1,k2,k3,e1,e2,e3,e4,e11,e12,e21,e22,w,j1,j2,F0);
   
10. u1=y(:,1); u2=y(:,2); u3=y(:,3); u4=y(:,4); u5=y(:,5); u6=y(:,6); u7=y(:,7); u8=y(:,8);
    
11. figure('unit','normalized','color',[1,1,1]);
    
12. h=get(gcf); set(gcf,'Name','1','numbertitle','off');
    
13. plot(t,u1); title('位移'); xlabel('时间t');ylabel('位移x'); grid on
    
14. figure('unit','normalized','color',[1,1,1]);
    
15. h=get(gcf); set(gcf,'Name','2','numbertitle','off');
    
16. plot(u1,u2); title('图2'); xlabel('位移x');ylabel('速度dx'); grid on
    
17. figure('unit','normalized','color',[1,1,1]);
    
18. h=get(gcf); set(gcf,'Name','3','numbertitle','off');
    
19. plot(t,u3); title('位移'); xlabel('时间t');ylabel('位移x'); grid on
    
20. figure('unit','normalized','color',[1,1,1]);
    
21. h=get(gcf); set(gcf,'Name','4','numbertitle','off');
    
22. plot(u3,u4); title('图4'); xlabel('位移x');ylabel('速度dx'); grid on
    
23. figure('unit','normalized','color',[1,1,1]);
    
24. h=get(gcf); set(gcf,'Name','5','numbertitle','off');
    
25. plot(t,u5); title('位移'); xlabel('时间t');ylabel('位移x'); grid on
    
26. figure('unit','normalized','color',[1,1,1]);
    
27. h=get(gcf); set(gcf,'Name','6','numbertitle','off');
    
28. plot(u5,u6); title('图6'); xlabel('位移x');ylabel('速度dx'); grid on
    
29. figure('unit','normalized','color',[1,1,1]);
    
30. h=get(gcf); set(gcf,'Name','7','numbertitle','off');
    
31. plot(t,u7); title('位移'); xlabel('时间t');ylabel('位移x'); grid on
    
32. figure('unit','normalized','color',[1,1,1]);
    
33. h=get(gcf); set(gcf,'Name','8','numbertitle','off');
    
34. plot(u7,u8); title('图8'); xlabel('位移x');ylabel('速度dx');
    
35. grid on

复制代码

这个输出结果是我需要用的，时间设定的比较长，会运行比较长的时间。

[ 本帖最后由 ChaChing 于 2009-3-8 12:00 编辑 ]

cgx11209

我想问下，你的那些初值的数值都是怎么取的？如：
m1=3000; m2=2065; l=660; c1=5; k1=260;
c2=5; k2=260; c3=30;  k3=150;
e1=0; e2=0; e3=189+e2; e4=189-e2;
e11=l+e1; e12=l-e1; e21=l+e2; e22=l-e2;
w=0.143; j1=5.194e+6; j2=8.477e+6; F0=37.5;

[ 本帖最后由 ChaChing 于 2009-3-8 11:52 编辑 ]

vitihl

因为这是一个实体模型的方程组，所以所有参数地选择是按照工程实际的参数来定的，就是说，实体模型里是多少，数学模型里同样也选用多少

yuanli

怎么得到下面的式子？
(-c1*(y(2)-y(6))*(j1+e11^2*m1)-c2*(y(4)-y(8))*(j1-e11*e12*m1)-k1*(y(1)-y(5))*(j1+e11^2*m1)-k2*(y(3)-y(7))*(j1-e11*e12*m1))*2*l/((e11+e12)*m1*j1);
y(4);

beehappy

受益匪浅，非常谢谢论坛里的高手！想问问版主或者其他大侠，这个四自由度非线性方程组中的非线性指的是方程右边的正弦函数吗？非线性体现在哪里呢？

[ 本帖最后由 beehappy 于 2009-3-8 07:40 编辑 ]

sogooda

呃。。。好像确实没有非线性项。

甜炒栗子

确实是一个很好的资源。非线性越来越热点了
真心收藏！谢谢！！

05071109

多谢楼主分享

jinghy88

请问如果要加入白噪声的路面激励该怎么办呢？

812268685

同上，求白噪声路面激励函数.....

可以这样定义吗？
fun1是路面激励函数，编程函数
将F0在计算程序里赋值为F0=fun1，可以这样吗？

我也没发现非线性项.....小弟初来乍到，求指教，诸位大侠见笑了

WF1987

感谢   学习下  

WF1987

vdp_eq  中的方程组时怎么化出来的啊  

yap988

每天来这里的目的就是为了学习，期待着有一天也能为这个论坛做一点贡献。