声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 3226|回复: 7

[编程技巧] 求助高手用matlab如何求变参数微分方程的解析解

[复制链接]
发表于 2011-1-9 14:25 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

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

D9E_SLUY~)37K4HZ{HEMZDJ.jpg
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=const
x0=[o1, 0, pi/18, 1, pi/2, 2, pi ];
y0=[y1, 2.3e6, 3.6e6, 1.5e6, 0.6e6, 0.3e6, 0.04e6];
x1=-pi:2*pi/200:pi;
p=interp1(x0,y0,x1,'cubic');
p11=p(1:101);p12=p(101:201);
%以第一缸点火上止点为零点,各个气缸转动一周(0-4*pi)的压力
p1=[p12(1:100),zeros(1,200),p11(1:100)];
p2=[p11(1:100),p12(1:100),zeros(1,200)];
p3=[zeros(1,100),p11(1:100),p12(1:100),zeros(1,100)];
p4=[zeros(1,200),p11(1:100), p12(1:100)];
%要求的微分方程
% 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);
空间代码
clear
load pdata.txt;
p1=pdata(:,1); p2=pdata (:,2); p3=pdata(:,3); p4=pdata(:,4);
y0=(1000*2*pi/60)^2;yy=y0
tspan=0:4*pi/399:4*pi;
for i=

1:length(tspan)-1
[t,y]=ode45(@w22,[tspan(i),tspan(i+1)],y0,[],p1(i),p2(i),p3(i),p4(i));
y0=y(end,: );
yy=[yy;y0];
end
w=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(负载)为零(即扭矩和不平衡),所以求得到的转速为不断上升
)}@U$L0@N73$}BTU0HA6L4T.jpg
我的课题是关于汽缸压力识别的,有兴趣的同学可加我QQ(三九三七三77三三一零,共同讨论。
关于曲轴瞬时转速可参考文献
1:利用转速波动信号在线识别内燃机气缸压力的研究(刘世元)
2:基于非线性动力学模型的柴油机瞬时转速仿真和缸内压力重构研究







file:///C:/DOCUME%7E1/ADMINI%7E1/LOCALS%7E1/Temp/msohtml1/02/clip_image002.gif
clip_image002.jpg
D9E_SLUY~)37K4HZ{HEMZDJ.jpg
回复
分享到:

使用道具 举报

发表于 2011-1-9 18:23 | 显示全部楼层
这个问题好像挺专业的,等高手啊!
 楼主| 发表于 2011-1-10 18:30 | 显示全部楼层
我已经想到一些思路,不过并不是最好的啊:可以先把汽缸压力模拟成角度的函数,这样整个式子只含有一个变量o,不过汽缸压力模型并不知道,所以还是有困难的。
发表于 2011-1-10 22:24 | 显示全部楼层
不过汽缸压力模型并不知道?
这个好像是专业问题,不是MATLAB自身的问题!
 楼主| 发表于 2011-1-11 00:29 | 显示全部楼层
呵呵,这个要做变参量的微分方程参数识别,其实也是matlab编程的问题哦,
发表于 2011-1-11 20:55 | 显示全部楼层
看完这两帖, LZ就了解版主的意思:@)

建议提问的网友分清 编程问题 和 专业问题
http://forum.vibunion.com/forum/ ... p;extra=&page=1
提问的智慧!!!!(发帖前请认真阅读)
http://forum.vibunion.com/forum/viewthread.php?tid=21991
发表于 2011-4-28 09:51 | 显示全部楼层
我现在做的问题恨这个相似,似乎不好解决!
 楼主| 发表于 2011-5-7 16:26 | 显示全部楼层
回复 7 # 龙马 的帖子

我本来是想先求解微分方程,然后识别方程的系数,后来尝试直接识别微分方程的系数,差不多做出来了

评分

1

查看全部评分

您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-11-29 08:33 , Processed in 0.084729 second(s), 23 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表