齿轮（变系数）微分方程该如何求解？已给自编程序，请指教

xiaowu

齿轮振动微分方程

图（1）为齿轮动力学振动微分方程

其中x为要求的位移响应，Zeta为齿轮副相对阻尼系数，Zeta =0.07，m为等效质量m=3；Fn为法向载荷Fn=100N。

k（t）为齿轮副的时变啮合刚度

对k(t)我做了关于时间t的傅里叶变换：

k=750527161.016+1.2466e+008*sin(2*pi*f0*t-0.86045)+8.3314e+007*sin(4*pi*f0*t-0.2292)+4.9789e+007*sin(6*pi*f0*t+0.41062)+1.7784e+007*sin(8*pi*f0*t+1.0182)+7.633e+006*sin(10*pi*f0*t-1.2305);

这个方程是变系数微分方程，其中就k（t）一个时变系数，k（t)该如何表示？是否可把k（t）作为t的函数直接代入方程？

下面是我编的M文件和求解程序，各位大虾帮我看看处理的对不对？谢谢

M文件

function Yd= fun7( t,x,k,zeta,m,Fn)
% FUN7.M： 时变微分方程例子
Fn=100;
zeta=0.07;
m=3;
n1=3000;
z1=45;
We=2*pi*z1*n1/60;
Tz=0.00044444;
f0=1/Tz;
k=750527161.016+1.2466e+008*sin(2*pi*f0*t-0.86045)+8.3314e+007*sin(4*pi*f0*t-0.2292)+4.9789e+007*sin(6*pi*f0*t+0.41062)+1.7784e+007*sin(8*pi*f0*t+1.0182)+7.633e+006*sin(10*pi*f0*t-1.2305);

dx=zeros(2,1);
dx(1)=x(2);
Yd=[x(2);Fn/m-2*zeta*sqrt(k/m)*x(2)-k/m*x(1)];

求解程序：

clear
m=3;
Tz=0.00044444;
f0=1/Tz;
tspan=linspace(0,0.00044444,120);
Y0=[1.5040e-007;0];
[t,XX]=ode45('fun7',tspan,Y0);
X=XX(:,1);
Y=XX(:,2);
plot(t,X),xlabel('经历时间/(s)'),ylabel('位移x'),figure;
plot(t,Y),xlabel('经历时间/(s)'),ylabel('速度v'),figure;
plot(X,Y)

greatchina

回复 1 # xiaowu 的帖子

你好，请问你的这个问题解决了吗？我的问题和你的一样，也是需要求解一个非线性方程。

happy

楼主的程序对不对我不知道，主要是看方程有没有写错之类的
但是调用格式是对的，系数可以通过求解的函数传递进去

帮助文件中的第三个例子就是类似的一个问题

zhouyang664

同意楼上，matlab--help ode45下的example3就是带参数传递的微分方程，建议楼主先学习一下，

xiaowu

谢谢各位的解答，我得出了结果了，程序没错，得出的曲线和实际有差异是因为一个参数写错了其中m=0.061658就对了。谢谢各位的热心帮助。

xiaowu

回复 2 # greatchina 的帖子

我得到结果了，程序没错，一个参数取错了。

greatchina

回复 6 # xiaowu 的帖子

你好，能不能给我传个程序呀，我好好学习一下。我是做链传动的，和你做的差不多。谢谢
great-china123@163.com

happy

一楼不是已经贴程序了吗？

hiv5

我在ODE45看帮助怎么只有两个例子呢？

ChaChing

回复 9 # hiv5 的帖子

版本不同吧!?
我的R2009a裡头有

fawcgzmg

1楼的方程从何而来？

xiaowu

回复 11 # fawcgzmg 的帖子

一个论文中看到的，需要可以发给你。

fawcgzmg

回复 12 # xiaowu 的帖子

请发到fawcgzmg@126.com

fawcgzmg

xiaowu 发表于 2010-10-25 11:13

                               
登录/注册后可看大图

图（1）为齿轮动力学振动微分方程其中x为要求的位移响应，Zeta为齿轮副相对阻尼系数，Zeta =0.07， ...

这个方程没考虑误差啊

maple1520

回复 12 # xiaowu 的帖子

请发给我，谢谢，maple1520@163.com