多自由度微分方程的求解

octopussheng

个人觉得四自由度的近似解（如用多尺度、摄动法）本身就是一个大难题了

luoluo

我明白大家的意思了，谢谢各位，我先自己求解一下看看，遇到问题再和大家讨论

咕噜噜

难度当然不小啊，非线性尤其是强非线性多自由度求解大于2自由度都不容易

drbethune

:@o 如果只有离散值的话,即使作出了数值解,这个解的精度也很难判别吧

咕噜噜

离散值？刚度和阻尼不是时间的函数吗

drbethune

当然是时间的函数,但也许没有现成的表达式,只有少数几个时刻的测量结果(或者是其他的原因)

咕噜噜

看楼主所给方程的情况不像离散状况

appleseed05

原帖由 luoluo 于 2007-6-28 18:41 发表
我建立了一个四自由度的振动方程:
    [M]{x''}+[C(t)]{x'}+[K(t)]{x}={F}
其中的刚度和阻尼都是时变的，对于这类方程，该如何求解呢？
用matlab自带的ode45可以求解吗？
哪位朋友有类似方程的matlab源码啊 ...

可以的，把它变成一阶的状态方程，找本转子动力学的书看看，上面有讲

[ 本帖最后由 mjhzhjg 于 2007-7-8 22:20 编辑 ]

luoluo

还是解不出来，我的方程化为一阶方程组如下：
% x1'=x2
% x2'=(-c1y*x2-k1y*x1-km*(x1+r1*x5-x3+r2*x7)-cm*(x2+r1*x6-x4+r2*x8))/m1
% x3'=x4
% x4'=(-c2y*x4-k2y*x3+km*(x1+r1*x5-x3+r2*x7)+cm*(x2+r1*x6-x4+r2*x8))/m2
% x5'=x6
% x6'=(-r1*(km*(x1+r1*x5-x3+r2*x7)+cm*(x2+r1*x6-x4+r2*x8))-M1)/I1
% x7'=x8
% x8'=(r2*(km*(x1+r1*x5-x3+r2*x7)+cm*(x2+r1*x6-x4+r2*x8))-M2)/I2

luoluo

其中的km和cm通过一个函数得到，
[km,cm]=meshkc(n)；算出来的km和cm都是维数为300的数组，其余的参数都是常数

luoluo

我写的将原方程组转化为一阶方程的程序如下：
m1=0.96;I1=4.3659e-4;m2=2.88;I2=8.3602e-3;  %质量和转动惯量
z1=19;z2=48; %主动轮和从动轮的齿数
k1y=6.56e7;k2y=6.56e7;c1y=1.8e5;c2y=1.8e5;   %支承刚度和阻尼
%km=1.3e6;cm=3.62;  %啮合刚度和阻尼
n=300;
[km,cm]=meshkc(n); % 一个啮合周期内的km和cm值
f1=30;   %输入轴的转动频率
T1=1/f1; %输入轴旋转一周的时间
Tz=T1/z1; %一个啮合周期
nsum = n*z1; %总的采样点数
r1=0.02834;r2=0.07160;  %大、小齿轮基圆半径
L=0.016;    %齿宽
Cr=1.6456;  %重合度
M1=11.9;M2=48.8;    %输入输出扭矩

%将原微分方程组化为一阶方程组，并写成myMulti.m
f=[x(2);(-c1y*x(2)-k1y*x(1)-km*(x(1)+r1*x(5)-x(3)+r2*x(7))-cm*(x(2)+r1*x(6)-x(4)+r2*x(8)))/m1;...
    x(4);(-c2y*x(4)-k2y*x(3)+km*(x(1)+r1*x(5)-x(3)+r2*x(7))+cm*(x(2)+r1*x(6)-x(4)+r2*x(8)))/m2;...
    x(6);(-r1*(km*(x(1)+r1*x(5)-x(3)+r2*x(7))+cm*(x(2)+r1*x(6)-x(4)+r2*x(8)))-M1)/I1;...
    x(8);(r2*(km*(x(1)+r1*x(5)-x(3)+r2*x(7))+cm*(x(2)+r1*x(6)-x(4)+r2*x(8)))-M2)/I2];
运行的时候报错误：Error using ==> vertcat
CAT arguments dimensions are not consistent.

Error in ==> myMulti at 19
f=[x(2);(-c1y*x(2)-k1y*x(1)-km*(x(1)+r1*x(5)-x(3)+r2*x(7))-cm*(x(2)+r1*x(6)-x(4)+r2*x(8)))，应该还是km和cm处理的不对吧，可是该怎么处理呢？

无水1324

km和cm应该在不同的时刻选取一个值

luoluo

原帖由 无水1324 于 2007-6-29 18:48 发表
km和cm应该在不同的时刻选取一个值

km和cm就是不同时刻取得值啊，它们都是和时间有关系的一些离散的值，我的km和cm都是n*2的矩阵，在我的程序中n取了300，km和cm的第一列代表的是时间，第二列代表的是km和cm的值

无水1324

f=[x(2);(-c1y*x(2)-k1y*x(1)-km*(x(1)+r1*x(5)-x(3)+r2*x(7))-cm*(x(2)+r1*x(6)-x(4)+r2*x(8)))/m1;...
    x(4);(-c2y*x(4)-k2y*x(3)+km*(x(1)+r1*x(5)-x(3)+r2*x(7))+cm*(x(2)+r1*x(6)-x(4)+r2*x(8)))/m2;...
    x(6);(-r1*(km*(x(1)+r1*x(5)-x(3)+r2*x(7))+cm*(x(2)+r1*x(6)-x(4)+r2*x(8)))-M1)/I1;...
    x(8);(r2*(km*(x(1)+r1*x(5)-x(3)+r2*x(7))+cm*(x(2)+r1*x(6)-x(4)+r2*x(8)))-M2)/I2];

那么你在f中必须指定

咕噜噜

首先是meshkc没有定义
其次你编这个程序最好写为m函数文件的形式
f没有定义