感谢+分享：四自由度微分方程组matlab求解
首先要感谢论坛里各位帮助过我的高手们，特别是siyanger关于非线性微分方程的求解的帖子， ch_j1985 在我最绝望的时候给与的帮助，以及版主eight一直以来的的严厉要求下，终于自己完成了，我的关于一个四自由度非线性方程组的求解，感谢所有帮助过我的，给我回复的所有人。现把程序贴出来，给以后需要的朋友借以参考，程序还有很多的不足之处，希望高手们继续给我提出修改的意见和建议。不多说了，来点儿实际的！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)
%y(1),y(3),y(5),y(7)为位移，y(2),y(4),y(6),y(8)为速度
f=[y(2);
(-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);
(-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);
y(6);
(-(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);
y(8);
(-(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)];
end以上存为M文件
在matlab内运行options=odeset('RelTol',1e-6,'AbsTol',);
%%%%Parameter make
t_end=2000; 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;
%%%%%%%%%%%%%
=ode45('vdp_eq',,,options,m1,m2,l,c1,c2,c3,k1,k2,k3,e1,e2,e3,e4,e11,e12,e21,e22,w,j1,j2,F0);
u1=y(:,1); u2=y(:,2); u3=y(:,3); u4=y(:,4); u5=y(:,5); u6=y(:,6); u7=y(:,7); u8=y(:,8);
figure('unit','normalized','color',);
h=get(gcf); set(gcf,'Name','1','numbertitle','off');
plot(t,u1); title('位移'); xlabel('时间t');ylabel('位移x'); grid on
figure('unit','normalized','color',);
h=get(gcf); set(gcf,'Name','2','numbertitle','off');
plot(u1,u2); title('图2'); xlabel('位移x');ylabel('速度dx'); grid on
figure('unit','normalized','color',);
h=get(gcf); set(gcf,'Name','3','numbertitle','off');
plot(t,u3); title('位移'); xlabel('时间t');ylabel('位移x'); grid on
figure('unit','normalized','color',);
h=get(gcf); set(gcf,'Name','4','numbertitle','off');
plot(u3,u4); title('图4'); xlabel('位移x');ylabel('速度dx'); grid on
figure('unit','normalized','color',);
h=get(gcf); set(gcf,'Name','5','numbertitle','off');
plot(t,u5); title('位移'); xlabel('时间t');ylabel('位移x'); grid on
figure('unit','normalized','color',);
h=get(gcf); set(gcf,'Name','6','numbertitle','off');
plot(u5,u6); title('图6'); xlabel('位移x');ylabel('速度dx'); grid on
figure('unit','normalized','color',);
h=get(gcf); set(gcf,'Name','7','numbertitle','off');
plot(t,u7); title('位移'); xlabel('时间t');ylabel('位移x'); grid on
figure('unit','normalized','color',);
h=get(gcf); set(gcf,'Name','8','numbertitle','off');
plot(u7,u8); title('图8'); xlabel('位移x');ylabel('速度dx');
grid on这个输出结果是我需要用的，时间设定的比较长，会运行比较长的时间。
[ 本帖最后由 ChaChing 于 2009-3-8 12:00 编辑 ] 我想问下，你的那些初值的数值都是怎么取的？如：
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 编辑 ]
关于参数地选择
因为这是一个实体模型的方程组，所以所有参数地选择是按照工程实际的参数来定的，就是说，实体模型里是多少，数学模型里同样也选用多少 怎么得到下面的式子？(-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 于 2009-3-8 07:40 编辑 ]
回复 5楼 beehappy 的帖子
呃。。。好像确实没有非线性项。 确实是一个很好的资源。非线性越来越热点了真心收藏！谢谢！！ 多谢楼主分享 请问如果要加入白噪声的路面激励该怎么办呢？ 同上，求白噪声路面激励函数.....
可以这样定义吗？
fun1是路面激励函数，编程函数
将F0在计算程序里赋值为F0=fun1，可以这样吗？
我也没发现非线性项.....小弟初来乍到，求指教，诸位大侠见笑了
感谢 学习下 vdp_eq中的方程组时怎么化出来的啊 每天来这里的目的就是为了学习，期待着有一天也能为这个论坛做一点贡献。
页:
[1]