求解2阶线性微分方程组的问题

zhangbing810822

请问解2阶线性微分方程组是不是一定要化成1阶方程组后解啊,MATLAB里有没有函数公式可套用?

[ 本帖最后由 eight 于 2007-9-15 10:23 编辑 ]

appleseed05

这是一个方法问题

解2阶线性微分方程组你可以用Newmark，welson-theta等方法
同样你也可以化成1阶方程组后用Runge-kutta

看你自己的要求，这些方法各有优缺点，你先找本结构动力学或者有限单元法看看

octopussheng

ode方法应该是目前最常用的了，如果要使用ode方法的话，就必须进行你上面说的这些步骤

qijunshuai

这位大侠，能否具体点，用Newmark法，是不是直接可以进行运算？

appleseed05

原帖由 qijunshuai 于 2007-9-16 20:06 发表

                               
登录/注册后可看大图

这位大侠，能否具体点，用Newmark法，是不是直接可以进行运算？

当然可以，结构动力学书上有现成的公式

qijunshuai

我有程序，感觉算法也不错，但输入参数后，求出的结果不对，你能否帮着看一下？难道是输入的参数不对？

qijunshuai

%计算转子系统的系数矩阵
%密度：density;
%轴单元的长度：rotor_length
%轴单元直径：rotor_dia;
%弹性模量：E
%轴段的质量矩阵：rotor_mass;
%圆盘的质量：disk_mass;
%圆盘直径：disk_dia;
%圆盘宽度：disk_width;
%圆盘的极转动惯量：disk_JP;
%圆盘的直径转动惯量： disk_JD;
%偏心距：ee1,ee2;
%轴承刚度k1x,k1y;k2x,k2y;
%轴承阻尼：c1x,cly;c2x,c2y;
%轴段单位长的质量为：ms_unit
%定义为4个节点
%其中第2个节点和第四个节点各有一圆盘。
clear
density=7850;
rotor_dia=0.01;
sl1=0.2;
sl2=0.2;
sl3=0.1;
E=2.07*10^11;
disk_dia=0.08;
disk_width=0.01;
k1x=4.7*10^8;
k2x=9.3*10^8;
c1x=2*10^4;
c2x=5*10^5;
% %偏心距
% ee1=0.0001;
% ee2=0.0003;
w=2*pi*100;
% 圆盘的质量，直径转动惯量和极转动惯量
disk_m=density*disk_width*pi*disk_dia^2/4;
disk_JD=disk_m*disk_dia^2/16;
disk_JP=2*disk_JD;
%轴段单位长质量,单位长的极转动惯量和直径转动惯量
ms_unit=density*pi*rotor_dia^2/4;
jps_unit=ms_unit*rotor_dia^2/8;
% jds_unit=
%轴段1的质量
s1_m=ms_unit*sl1;
%轴段2的质量
s2_m=ms_unit*sl2;
%轴段3的质量
s3_m=ms_unit*sl3;
% 节点2的集总质量,集中直径转动惯量、极转动惯量
node2_m=disk_m+s1_m/2+s2_m/2;
node2_JD=disk_JD+ms_unit*sl1^3/12+ms_unit*sl2/12;
node2_JP=disk_JP+jps_unit*sl1/2+jps_unit*sl2/2;
%节点4
node4_m=disk_m+s3_m/2;
node4_JD=disk_JD+ms_unit*sl3^3/12;
node4_JP=disk_JP+jps_unit*sl3/2;
%节点1
node1_m=s1_m/2;
node1_JD=ms_unit*sl1^3/12;
node1_JP=jps_unit*sl1/2;
%节点3
node3_m=s2_m/2+s3_m/2;
node3_JD=ms_unit*sl2^3/12+ms_unit*sl3/12;
node3_JP=jps_unit*sl2/2+jps_unit*sl3/2;
M1=[node1_m     0         0        0         0          0          0        0;
      0      node1_JD     0        0         0          0          0        0;
      0         0      node2_m     0         0          0          0        0;
      0         0         0      node2_JD    0          0          0        0;
      0         0         0        0       node3_m      0          0        0;
      0         0         0        0          0      node3_JD      0        0;
      0         0         0        0          0          0      node4_m     0;
      0         0         0        0          0          0         0     node4_JD
    ];
J1=[  c1x       0         0        0         0          0          0          0;
      0      node1_JP     0        0         0          0          0         0;
      0         0         0        0         0          0          0         0;
      0         0         0      node2_JP    0          0          0         0;
      0         0         0        0         c2x           0          0        0;
      0         0         0        0          0      node3_JP      0         0;
      0         0         0        0          0          0         0         0;
      0         0         0        0          0          0         0     node4_JP
    ];
EI=E*pi*rotor_dia^4/64;     %抗弯刚度

d11=12; d12=6*sl1; d13=-12;d14=d12;
d22=4*sl1^2; d23=-d12; d24=2*sl1^2;
d33=12; d34=-d12;  d44=4*sl1^2;
D_s1=[ d11 d12 d13 d14;
    d12 d22 d23 d24;
    d13 d23 d33 d34;
    d14 d24 d34 d44];
ks1=(EI/sl1^3).*D_s1;   %轴段的刚度矩阵
k11_s1=[ks1(1,1) ks1(1,2);
    ks1(2,1) ks1(2,2)];
k12_s1=[ks1(1,3) ks1(1,4);
    ks1(2,3) ks1(2,4)];
k22_s1=[ks1(3,3) ks1(3,4);
    ks1(4,3) ks1(4,4)];
%第二段轴段的刚度矩阵
d11=12; d12=6*sl2; d13=-12;d14=d12;
d22=4*sl2^2; d23=-d12; d24=2*sl2^2;
d33=12; d34=-d12;  d44=4*sl2^2;
D_s2=[ d11 d12 d13 d14;
    d12 d22 d23 d24;
    d13 d23 d33 d34;
    d14 d24 d34 d44];
ks2=(EI/sl2^3).*D_s2;  
k11_s2=[ks2(1,1) ks2(1,2);
    ks2(2,1) ks2(2,2)];
k12_s2=[ks2(1,3) ks2(1,4);
    ks2(2,3) ks2(2,4)];
k22_s2=[ks2(3,3) ks2(3,4);
    ks2(4,3) ks2(4,4)];
%第三段
d11=12; d12=6*sl3; d13=-12;d14=d12;
d22=4*sl3^2; d23=-d12; d24=2*sl3^2;
d33=12; d34=-d12;  d44=4*sl3^2;
D_s3=[ d11 d12 d13 d14;
    d12 d22 d23 d24;
    d13 d23 d33 d34;
    d14 d24 d34 d44];
ks3=(EI/sl3^3).*D_s3;  
k11_s3=[ks3(1,1) ks3(1,2);
    ks3(2,1) ks3(2,2)];
k12_s3=[ks3(1,3) ks3(1,4);
    ks3(2,3) ks3(2,4)];
k22_s3=[ks3(3,3) ks3(3,4);
    ks3(4,3) ks3(4,4)];

G1=w.*J1;
z=[0 0;
    0 0];
kx1=[k1x 0;
    0   0];
kx2=[k2x 0;
    0 0];
K2=[k11_s1+kx1      k12_s1       z            z   ;   
    k12_s1    k22_s1+k11_s2    k12_s2           z    ;   
    z        k12_s2       k22_s2+k11_s3+kx2   k12_s3 ;
    z         z             k12_s3      k22_s3  
  ];
zero=zeros(8,8);
M=[M1 zero;
    zero M1];
C=[zero G1;
    -G1 zero];
K=[K2 zero;
    zero  K2];

%%------------------------------------  算法           
u=zeros(16,1);                 
v=zeros(16,1);
x(:,1)=u;                                     %位移                       
x1(:,1)=v;                                    %速度
e1=0.0005 ;                                  %偏心距
e2=0.0003;
delt=0.5;
alpha=0.25;
dt=2*pi/(1024*8);                            %积分步长
%%%%计算积分常数
b0=1/(alpha*dt^2);
b1=delt/(alpha*dt);

b2=1/(alpha*dt);
b3=(1/(2*alpha))-1;
b4=(delt/alpha)-1;
b5=(dt/2)*((delt/alpha)-2);
b6=dt*(1-delt);
b7=delt*dt;     
t_max=4;                 
ii=1;
t(1)=0;                 
q=zeros(16,1);
while t(ii)<t_max                        
     fux1(ii)=disk_m*e1*w*w*cos(w*t(ii));       %不平衡力
     fuy1(ii)=disk_m*e1*w*w*sin(w*t(ii));       %不平衡力
%      fux2(ii)=disk_m*e2*w*w*cos(w*t(ii));       %不平衡力
%      fuy2(ii)=disk_m*e2*w*w*sin(w*t(ii));       %不平衡力
     Q(:,ii)=zeros(16,1);
     Q1(:,ii)=zeros(16,1);
     Q1(3,ii)=fux1(ii);
     Q1(11,ii)=fuy1(ii);
%      Q1(7,ii)=fux2(ii);
%      Q1(15,ii)=fuy2(ii);
     Q(:,ii)=Q1(:,ii);
     Ccc=C;
     
     K1=K+b0*M+b1*Ccc ;        %公式
     
     if ii==1
        x2(:,ii)=inv(M)*(Q(:,ii)-K*x(:,ii)-Ccc*x1(:,ii));
        q(:,ii)=Q(:,ii)+M*(b0*x(:,ii)+b2*x1(:,ii)+b3*x2(:,ii))+Ccc*(b1*x(:,ii)+b4*x1(:,ii)+b5*x2(:,ii));
        x(:,ii+1)=K1\q(:,ii);
        x2(:,ii+1)=b0*(x(:,ii+1)-x(:,ii))-b2*x1(:,ii)-b3*x2(:,ii);
        x1(:,ii+1)=x1(:,ii)+b6*x2(:,ii)+b7*x2(:,ii+1);
      else
        q(:,ii+1)=Q(:,ii)+M*(b0*x(:,ii)+b2*x1(:,ii)+b3*x2(:,ii))+Ccc*(b1*x(:,ii)+b4*x1(:,ii)+b5*x2(:,ii));
        x(:,ii+1)=K1\q(:,ii+1);                                      %位移
        x2(:,ii+1)=b0*(x(:,ii+1)-x(:,ii))-b2*x1(:,ii)-b3*x2(:,ii);   %加速度
        x1(:,ii+1)=x1(:,ii)+b6*x2(:,ii)+b7*x2(:,ii+1);               %速度   
      end
       ii=ii+1;
       t(ii)=t(ii-1)+dt;
end

[ 本帖最后由 qijunshuai 于 2007-9-17 21:13 编辑 ]

qijunshuai

哪位高人帮着运行下。我不知道问题出在哪儿了！！多谢！！

花如月

自己根据出错提示，进行修改吧。置顶帖子里有常见的程序错误和解决办法

qijunshuai

运行没有错误，但结果不大对！！

咕噜噜

你怎么知道结果不对，结果贴出来看看，说明一下为什么不对

qijunshuai

原帖由 咕噜噜 于 2007-9-18 15:23 发表

                               
登录/注册后可看大图

你怎么知道结果不对，结果贴出来看看，说明一下为什么不对

值都成了无穷大了！求得值越来越大！！

[ 本帖最后由 eight 于 2007-9-19 10:07 编辑 ]

octopussheng

这种情况两种原因：
1。你的方程是刚性的，用ode15试试，ode45已经不适用了
2。修改你的系统参数，直到结果比较接近理论值～