weiniuzhu 发表于 2011-1-9 14:25

求助高手用matlab如何求变参数微分方程的解析解

本帖最后由 weiniuzhu 于 2011-1-9 19:02 编辑


1111 :下面说我的问题:,我要求的方程如上图,实际上是曲轴动力学关于扭矩的平衡方程。通过角域变换后把二阶微分方程变为一阶,
y'+P(o)*y=Q(o)=(Ap*R*(p1*f1+p2*f2+p3*f1+p4*f2)-TL)/(J+m*R^2*fsum);(1)
请问高手(1)这个一阶方程的解析解能否求出来,
      o为弧度,其他参数都是发动机常数。我想仿真在压力作用下的转速波动规律,其中pk代表第k缸的压力.由于p是变化的,我想问一下高手,此一阶方程是否存在解析解?
2222: 我查了一些帖子,有的说可以把p当做常数求解方程后再代入,不过经过我的计算,此方法不可行啊,计算结果不符合常识(常见的转速波动规律)。后来用数值方法倒是和实验结果一致,不过我还是想求出它的解析解。
3333 我的最终目的是要通过曲轴转速波动和汽缸压力来识别,m ,R ,L,Ap,等内燃机参数,关于微分方程的参数识别教程比较少,所以想把解析解求出来,通过最小二乘法就比较好做了。如果实在求不出来,只好找关于微分方程参数识别的文章了。

44444:我把我的数值解程序发出来吧,大家共同学习
%内燃机参数R=52.5e-3;%曲柄半径L=184e-3;%连杆长度Ap= 0.0071;%活塞面积J=0.5;%转动惯量m=1.957;%单个缸往复质量n1=1000;%初始转速w1=(n1*2*pi/60)^2;%初始角域角速度平方TL=0;%负载% 模拟汽缸压力数据Vp=0.002977051738358;%排量Vc=Vp/28;%扫气体积d=95e-3;%缸径
o1=-pi:17*pi/(18*59):-pi/18;%从下至点到点火提前角Vp=0.002977051738358;%排量Vc=Vp/28;%扫气体积d=95e-3;%缸径
o1=-pi:17*pi/(18*59):-pi/18;%从下至点到点火提前角R=52.5e-3;%曲柄半径;L=184e-3;%连杆长度;S=R*cos(o1)+(L^2-R^2*sin(o1).^2).^0.5;V=Vc+d^2*pi./4*(R+L-S);y1=0.8e5*(Vp./(4*V)).^1.55;%压缩冲程压力当做压缩,pv^1.55=constx0=;y0=;x1=-pi:2*pi/200:pi;p=interp1(x0,y0,x1,'cubic');p11=p(1:101);p12=p(101:201);%以第一缸点火上止点为零点,各个气缸转动一周(0-4*pi)的压力p1=;
p2=;p3=;p4=;%要求的微分方程% dy=Q-P*y
方程(1)P=8*m*R^2*(2*R^4*cos(o)^4-4*R^4*cos(o)^2+2*R^4+4*L^2*R^2*cos(o)^2-3*L^2*R^2+L^4)*cos(o)*sin(o)/(J*L^2-J*R^2+J*R^2*cos(o)^2+4*L^2*m*R^2+12*cos(o)^2*m*R^4-4*m*R^4-4*cos(o)^2*L^2*m*R^2-8*cos(o)^4*m*R^4)/(L^2-R^2+R^2*cos(o)^2);Q=(2*Ap*R*(p1*(sin(o)+1/2*R*sin(2*o)/L/(1-R^2/L^2*sin(o)^2)^(1/2))+p2*(-sin(o)+1/2*R*sin(2*o)/L/(1-R^2/L^2*sin(o)^2)^(1/2))+p3*(sin(o)+1/2*R*sin(2*o)/L/(1-R^2/L^2*sin(o)^2)^(1/2))+p4*(-sin(o)+1/2*R*sin(2*o)/L/(1-R^2/L^2*sin(o)^2)^(1/2)))-2*TL)/(J+m*R^2*(2*(sin(o)+1/2*R*sin(2*o)/L/(1-R^2/L^2*sin(o)^2)^(1/2))^2+2*(-sin(o)+1/2*R*sin(2*o)/L/(1-R^2/L^2*sin(o)^2)^(1/2))^2))然后把R LAp 等参数带入P,Q 得到dy=Q-P*y,并写出要求的微分方程M函数w22

function dy=w22(o,y,p1,p2,p3,p4)dy=(1719005963368809/2305843009213693952*p1*(sin(o)+105/736*sin(2*o)/(1-2933139525741949/36028797018963968*sin(o)^2)^(1/2))+1719005963368809/2305843009213693952*p2*(-sin(o)+105/736*sin(2*o)/(1-2933139525741949/36028797018963968*sin(o)^2)^(1/2))+1719005963368809/2305843009213693952*p3*(sin(o)+105/736*sin(2*o)/(1-2933139525741949/36028797018963968*sin(o)^2)^(1/2))+1719005963368809/2305843009213693952*p4*(-sin(o)+105/736*sin(2*o)/(1-2933139525741949/36028797018963968*sin(o)^2)^(1/2)))/(1/2+6218836978571121/576460752303423488*(sin(o)+105/736*sin(2*o)/(1-2933139525741949/36028797018963968*sin(o)^2)^(1/2))^2+6218836978571121/576460752303423488*(-sin(o)+105/736*sin(2*o)/(1-2933139525741949/36028797018963968*sin(o)^2)^(1/2))^2)-y*(55775835223465209299070696186807/85070591730234615865843651857942052864*cos(o)^4+629338605529600761106491833156601/42535295865117307932921825928971026432*cos(o)^2+3235856465456370238136197978809399/85070591730234615865843651857942052864)*cos(o)*sin(o)/(2393779504991444819/147573952589676412928+3809512810466140650883377/4611686018427387904000000000*cos(o)^2-17552045488319130019/147573952589676412928000*cos(o)^4)/(8963892640724197/288230376151711744+6355479794145243/2305843009213693952*cos(o)^2);
空间代码
clearload pdata.txt;p1=pdata(:,1); p2=pdata (:,2); p3=pdata(:,3); p4=pdata(:,4);
y0=(1000*2*pi/60)^2;yy=y0tspan=0:4*pi/399:4*pi;
for i=
1:length(tspan)-1
=ode45(@w22,,y0,[],p1(i),p2(i),p3(i),p4(i));
y0=y(end,: );
yy=;endw=yy.^0.5;%求得的yy为角域角速度的平方
n=w*60/(2*pi);%求转速file:///C:/DOCUME%7E1/ADMINI%7E1/LOCALS%7E1/Temp/%29%7D@U$L0@N73$%7DBTU0HA6L4T.jpg
plot(tspan,yy)由于我设置TL(负载)为零(即扭矩和不平衡),所以求得到的转速为不断上升我的课题是关于汽缸压力识别的,有兴趣的同学可加我QQ(三九三七三77三三一零,共同讨论。
关于曲轴瞬时转速可参考文献
1:利用转速波动信号在线识别内燃机气缸压力的研究(刘世元)
2:基于非线性动力学模型的柴油机瞬时转速仿真和缸内压力重构研究






file:///C:/DOCUME%7E1/ADMINI%7E1/LOCALS%7E1/Temp/msohtml1/02/clip_image002.gif

zhouyang664 发表于 2011-1-9 18:23

这个问题好像挺专业的,等高手啊!

weiniuzhu 发表于 2011-1-10 18:30

我已经想到一些思路,不过并不是最好的啊:可以先把汽缸压力模拟成角度的函数,这样整个式子只含有一个变量o,不过汽缸压力模型并不知道,所以还是有困难的。

zhouyang664 发表于 2011-1-10 22:24

不过汽缸压力模型并不知道?
这个好像是专业问题,不是MATLAB自身的问题!

weiniuzhu 发表于 2011-1-11 00:29

呵呵,这个要做变参量的微分方程参数识别,其实也是matlab编程的问题哦,

ChaChing 发表于 2011-1-11 20:55

看完这两帖, LZ就了解版主的意思:@)

建议提问的网友分清 编程问题 和 专业问题
http://forum.vibunion.com/forum/viewthread.php?tid=36746&extra=&page=1
提问的智慧!!!!(发帖前请认真阅读)
http://forum.vibunion.com/forum/viewthread.php?tid=21991

龙马 发表于 2011-4-28 09:51

我现在做的问题恨这个相似,似乎不好解决!

weiniuzhu 发表于 2011-5-7 16:26

回复 7 # 龙马 的帖子

我本来是想先求解微分方程,然后识别方程的系数,后来尝试直接识别微分方程的系数,差不多做出来了
页: [1]
查看完整版本: 求助高手用matlab如何求变参数微分方程的解析解