求解2阶线性微分方程组的问题
请问解2阶线性微分方程组是不是一定要化成1阶方程组后解啊,MATLAB里有没有函数公式可套用?[ 本帖最后由 eight 于 2007-9-15 10:23 编辑 ] 这是一个方法问题
解2阶线性微分方程组你可以用Newmark，welson-theta等方法
同样你也可以化成1阶方程组后用Runge-kutta
看你自己的要求，这些方法各有优缺点，你先找本结构动力学或者有限单元法看看 ode方法应该是目前最常用的了，如果要使用ode方法的话，就必须进行你上面说的这些步骤
回复 #2 appleseed05 的帖子
这位大侠，能否具体点，用Newmark法，是不是直接可以进行运算？ 原帖由 qijunshuai 于 2007-9-16 20:06 发表 http://www.chinavib.com/forum/images/common/back.gif这位大侠，能否具体点，用Newmark法，是不是直接可以进行运算？
当然可以，结构动力学书上有现成的公式
回复 #5 appleseed05 的帖子
我有程序，感觉算法也不错，但输入参数后，求出的结果不对，你能否帮着看一下？难道是输入的参数不对？回复 #4 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;
zeroK2];
%%------------------------------------算法
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 编辑 ]
回复 #7 qijunshuai 的帖子
哪位高人帮着运行下。我不知道问题出在哪儿了！！多谢！！回复 #8 qijunshuai 的帖子
自己根据出错提示，进行修改吧。置顶帖子里有常见的程序错误和解决办法回复 #9 花如月 的帖子
运行没有错误，但结果不大对！！ 你怎么知道结果不对，结果贴出来看看，说明一下为什么不对 原帖由 咕噜噜 于 2007-9-18 15:23 发表 http://www.chinavib.com/forum/images/common/back.gif你怎么知道结果不对，结果贴出来看看，说明一下为什么不对
值都成了无穷大了！求得值越来越大！！
[ 本帖最后由 eight 于 2007-9-19 10:07 编辑 ]
回复 #12 qijunshuai 的帖子
这种情况两种原因：1。你的方程是刚性的，用ode15试试，ode45已经不适用了
2。修改你的系统参数，直到结果比较接近理论值～
页:
[1]