首先是meshkc没有定义
其次你编这个程序最好写为m函数文件的形式
f没有定义
我只是拷了一部分代码，那个meshkc是以一个函数的格式写的，f也定义了 原帖由 无水1324 于 2007-6-29 20:31 发表
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;. ...
我不太明白你的意思，该怎么指定呢，不要笑哦，我比较笨。:@( 你这种写法通常是在m函数文件中使用的，m函数中并不直接定义，调用的时候才定义 我的函数是这样的
function f=myMulti(t,x)
%初始数据
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;
=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];
回复 #32 luoluo 的帖子
就是你就算得时刻t，取该时刻对应的km和cm 求啮合刚度和啮合阻尼的函数：function = meshkc(n)
%meshkc 计算齿轮副的啮合刚度及阻尼
%
%输入
%n 在一个啮合周期内的计算点数（缺省时计算300点）
%
%输出
%km 随时间变化的啮合刚度值（n*2）
%cm 随时间变化的啮合阻尼值(n*2)
%初始参数
z=;xi=;sigma=0.0;alpha=22.6*pi/180.0;
alpha0=20.0*pi/180.0;rb1=0.02834;
omega=;epsilon=1.6456;
pn=419.9012;
Tz=1.7544e-3;E=2.068e11;mu=0.3;
if nargin<1, n=300;
end
for i=1:2
if (z(i)>135),z (i)=135;
end
if (z(i)<25)
A(i)=11.05+0.166*(abs(z(i)-25))^0.2;
B(i)=0.32+0.05*(abs(z(i)-25))^0.35;
else
A(i)=11.05-0.166*(abs(z(i)-25))^0.2;
B(i)=0.32-0.05*(abs(z(i)-25))^0.35;
end
C(i)=(abs(z(i)-17))^0.45/38.4+1.22;
end
t01=((omega(1)+omega(2))*tan(alpha)-sqrt((z(2)+2+2*xi(2)-2*sigma)^2*...
omega(2)^2/(z(1)^2*(cos(alpha0))^2)-omega(1)^2))/omega(1)/omega(2);
nn=fix(n/epsilon);
t=0.0;
dt=Tz/(n-1);
for m=1:n
km(m,1)=t;
cm(m,1)=t;
for i=1:2
%轮齿的弯曲刚度(kw(i)(i=1,2));
lamda(i)=0.5*(z(i)+2+2*xi(i))-0.5*sqrt(1+omega(i)^2*(t+t01)^2)*...
z(i)*cos(alpha0);
kw(i)=E*(1+xi(i))^B(i)*(1+lamda(i))^C(i)/A(i);
%轮齿的接触刚度(kc(i)(i=1,2));
rho(1)=omega(1)*(t+t01)*rb1;
rho(2)=((omega(1)+omega(2))/omega(2)*tan(alpha)-omega(1)*(t+t01))*rb1;
if(i==1),j=2;else,j=1;
end
kc(i)=pi*E/((1-mu^2)*(log(pi*E*rho(i)*(rho(1)+rho(2)))-...
log(2*(1-mu^2)*rho(j)*pn)+0.814));
%单齿刚度(k(i)(i=1,2));
k(i)=kw(i)*kc(i)/(kw(i)+kc(i));
end
if(m<=nn) %单齿啮合区内啮合刚度
km(m,2)=k(1)*k(2)/(k(1)+k(2));
else %双齿啮合区内啮合刚度
km(m,2)=k(1)*k(2)/(k(1)+k(2))+km(m-nn,2);
end
t=t+dt;
gxi=0.07;
r1=55.25/2000;r2=182.75/2000;J1=3.698332e-4;J2=15.257864e-3;
cm(m,2)=2*gxi*sqrt(km(m,2)*r1^2*r2^2*J1*J2/(r1^2*J1*r2^2*J2));
end 齿轮的东西我不太懂，但是对于二阶微分方程处理成一阶的一般是这样（如下图）
然后就可以RK了 原帖由 无水1324 于 2007-6-29 21:26 发表
就是你就算得时刻t，取该时刻对应的km和cm
我明白你的意思了，我再改改看吧，谢谢了！ 原帖由 appleseed05 于 2007-6-29 21:49 发表
齿轮的东西我不太懂，但是对于二阶微分方程处理成一阶的一般是这样（如下图）
然后就可以RK了
谢谢，只不过我的刚度和阻尼是时变的，如果是常数的话我想我也能解出来吧 1.楼主说的是变系数, 并非是非线性问题.
2.不管是线性还是非线性,ode45采用的数值解法总是可以的.
3.楼主所说的离散形式,当然在离散时刻上是已知的,.当采用ode45时一般是等步长的,如果这个步长与阻尼和刚度的采样间隔一致,可以采用插值将阻尼和刚度的采样间隔之间插值,以便与ode45的步长一致. 插值方法很多,找本数值方法看一看. 原帖由 VibrationMaster 于 2007-6-30 07:12 发表
1.楼主说的是变系数, 并非是非线性问题.
2.不管是线性还是非线性,ode45采用的数值解法总是可以的.
3.楼主所说的离散形式,当然在离散时刻上是已知的,.当采用ode45时一般是等步长的,如果这个步长与阻尼和刚度的 ...
“ 采用ode45时一般是等步长的”ode45时一般是并不是等步长的。这点要注意