| 查看: 1184 | 回复: 3 | ||
[求助]
动力学参数拟合
|
|
原料是：甲醇、CO、氧气。主反应：甲醇、CO、氧气合成DMC，副反应：CO和氧气合成CO2. 要干的活：拟合出生成DMC和CO2的动力学参数（指前因子、活化能、级数） 动力学方程为dy1/dw=f(y1、y2)和dy2/dw=f(y1、y2),具体形式见附图。 整理后的实验数据：N0(原料初始的摩尔数)、P（压力）、t(温度)、x1(甲醇初始的摩尔分数)、x2(CO初始的摩尔分数)、x3(O2初始的摩尔分数)、y1(DMC产物的摩尔分数)、y2(O2产物的摩尔分数)、W=0.375(催化剂的质量)。 关键难度在：主要是两个常微分方程、数据多不知如何进行运算。下面是自己编的程序，出现问题并且解决不上来，希望大神指教 function YYHJIA clear all;clc global t T P N0 X1 X2 X3 ; t=[140 140 140.2 143.4 150 150.1 150.2 151.2 153.7 159.2 160 161.3 140.2 141.1 149.7 150.3 150.8 159.7 161.2 138 150.3 160 160.2 162.3]; %温度 T=t+273.15; P=[0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.3 0.3 0.3 0.3 0.3]; %压力 N0=[0.3091 0.3260 0.2796 0.2856 0.2801 0.3260 0.3091 0.2853 0.3190 0.2853 0.3962 0.3190 0.3085 0.2793 0.2804 0.2856 0.3083 0.3085 0.2788 0.2368 0.26 0.2597 0.3207 0.2344]; %原料初始摩尔 x1=[0.624 0.5 0.5310 0.5710 0.53 0.5 0.624 0.572 0.511 0.572 0.524 0.511 0.625 0.531 0.529 0.571 0.625 0.625 0.532 0.501 0.57 0.571 0.555 0.506]; %原料甲醇占的摩尔分数 x2=[0.25 0.334 0.349 0.287 0.348 0.334 0.25 0.286 0.345 0.286 0.352 0.345 0.25 0.349 0.348 0.286 0.25 0.25 0.35 0.334 0.287 0.288 0.334 0.337]; %原料CO摩尔分数 x3=[0.127 0.166 0.121 0.142 0.122 0.166 0.127 0.142 0.144 0.142 0.124 0.144 0.125 0.12 0.123 0.143 0.124 0.125 0.118 0.165 0.142 0.141 0.111 0.157]; %原料O2的摩尔分数 Y1=[0.015 0.016 0.015 0.017 0.02 0.021 0.021 0.021 0.021 0.025 0.023 0.025 0.018 0.018 0.025 0.026 0.026 0.03 0.03 0.023 0.028 0.036 0.034 0.038]; %出口DMC的摩尔数 Y2=[0.018 0.02 0.018 0.033 0.04 0.041 0.042 0.042 0.045 0.054 0.053 0.054 0.031 0.032 0.05 0.051 0.05 0.059 0.058 0.043 0.051 0.063 0.061 0.066]; %出口CO的摩尔数 k0=[8e4 2e4 -0.5 0.5 1 4e5 3e4 0.5 0.5 ];%(指前因子、活化能、级数)初始量 lb=[0 0 -3 -3 -3 0 0 -3 -3 ];%下限 ub=[+inf +inf +inf +inf +inf +inf +inf +inf +inf];%上限 %使用lsqnonlin函数进行参数估计 [k,resnorm,residual,exitflag,output,lambda,jacobian]=lsqnonlin(@ObjFunc,k0,lb,ub,[],T,P,N0,X1,X2,X3,Y1,Y2); ci=nlparci(k,residual,jacobian); function f=ObjFunc(k,T,P,N0,X1,X2,X3,Y1,Y2) wspan=[0 0.375]; Y0=[Y1' Y2'] [m,n]=size(Y0); Ycal = zeros(m,n); for i=1:m [w yy]=ode45(@Equations,wspan,Y0(i,  ,[],k);Ycal(i, =yy(end,;end f=Ycal-Y0; function dydw=Equations(w,yy,k) %如何表达呢？这么多数字 r1=k(1)*exp*(-k(2)./T).*P.^(k(3)+k(4)+k(5))*... (yy(1).*(1.5*x1-2)+0.5*yy(2).*x1).^(k(3))*... (yy(1).*(1.5*x2-1)+yy(2).*(0.5*x2-1)+x2).^(k(4))*... (yy(1).*(1.5*x3-0.5)+yy(2).*(0.3*x3-0.5)+x3).^k(5); r2=k(6)*exp*(-k(7)./T).*P.^(k(8)+k(9))*... (yy(1).*(1.5*x2-1)+yy(2).*(0.5*x2-1)+x2).^k(8)*... (yy(1).*(1.5*x3-0.5)+yy(2).*(0.3*x3-0.5)+x3).^k(9); dy1dw=(1+1.5*yy(1)+0.5*yy(2)).*(r1.*(1.5*yy(1)+1)+0.5*yy(1).*r2)./N0; dy2dw=(1+1.5*yy(1)+0.5*yy(2)).*(r2.*(1+0.5*yy(2))+1.5*r1.*yy(2))./N0; dydw=[dy1dw;dy2dw]; 出现的问题如图：  出现的问题.jpg  方程形式.jpg  方程具体形式2.jpg |
回复此楼» 猜你喜欢
建议基金发布提前给出明确的时间点
已经有14人回复
-
范进中举一文的中心思想
已经有5人回复
-
2026国自然函评费到账
已经有17人回复
-
让我中一个面上吧！
已经有15人回复
-
2026年的国家社科基金项目通讯评审的新规则与新动向、新挑战
已经有5人回复
-
跳槽后在研项目怎么办？
已经有10人回复
-
什么时候开奖？
已经有10人回复
-
今天放榜吗？
已经有16人回复
-
93BebMhtakh前后11位开头都是大写
已经有8人回复
-
只有每年这种时候来逛逛小木虫
已经有26人回复
» 本主题相关价值贴推荐，对您同样有帮助:
- 这样的动力学方程能够有1stopt拟合出参数来么
已经有4人回复
- 使用MATLAB拟合动力学参数问题
已经有18人回复
- 动力学方程组参数求解问题
已经有9人回复
- 使用matlab最优化方法拟合获得多个动力学参数中的问题
已经有4人回复
- 求大神帮忙matlab拟合函数求参数
已经有9人回复
- 求助1stopt拟合动力学参数
已经有4人回复
- 拜求！matlab拟合动力学参数
已经有14人回复
- 求助动力学方程？？？
已经有3人回复
- 求助用matlab拟合动力学方程
已经有13人回复
- 动力学拟合求助
已经有4人回复
- 反应动力学的拟合有什么意义？
已经有12人回复
- 通过实验数据拟合，求解公式中的参数
已经有14人回复
- 二级动力学 拟合
已经有9人回复
- 怎样用一组参数同时拟合两个曲线--matlab
已经有5人回复
- 动力学方程参数估计方法
已经有14人回复
- 动力学参数拟合
已经有26人回复
- matlab 拟合反应动力学参数结果很差。大家帮忙看一下
已经有14人回复
- 为啥我的一级动力学方程拟合出来的数据R是负数呢？
已经有16人回复
- matlab拟合拟均相动力学参数
已经有6人回复
- 【求助】用matlab最优化方法进行参数拟合
已经有17人回复
- 【求助】拟合动力学方程求助
已经有13人回复
信彼南山
木虫 (著名写手)
- 应助: 33 (小学生)
- 金币: 4142.9
- 散金: 1221
- 红花: 16
- 帖子: 1178
- 在线: 233.5小时
- 虫号: 1133529
- 注册: 2010-10-27
- 专业: 导航、制导与传感技术
2楼2015-03-05 01:28:52
3楼2015-03-05 08:35:46
|
function YYHJIA clear all;clc global T P N0 x1 x2 x3; t=[140 140 140.2 143.4 150 150.1 150.2 151.2 153.7 159.2 160 161.3 140.2 141.1 149.7 150.3 150.8 159.7 161.2 138 150.3 160 160.2 162.3]; %温度 T=t+273.15; P=[0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.3 0.3 0.3 0.3 0.3]; %压力 N0=[0.3091 0.3260 0.2796 0.2856 0.2801 0.3260 0.3091 0.2853 0.3190 0.2853 0.3962 0.3190 0.3085 0.2793 0.2804 0.2856 0.3083 0.3085 0.2788 0.2368 0.26 0.2597 0.3207 0.2344]; %原料初始摩尔 x1=[0.624 0.5 0.5310 0.5710 0.53 0.5 0.624 0.572 0.511 0.572 0.524 0.511 0.625 0.531 0.529 0.571 0.625 0.625 0.532 0.501 0.57 0.571 0.555 0.506]; %原料甲醇占的摩尔分数 x2=[0.25 0.334 0.349 0.287 0.348 0.334 0.25 0.286 0.345 0.286 0.352 0.345 0.25 0.349 0.348 0.286 0.25 0.25 0.35 0.334 0.287 0.288 0.334 0.337]; %原料CO摩尔分数 x3=[0.127 0.166 0.121 0.142 0.122 0.166 0.127 0.142 0.144 0.142 0.124 0.144 0.125 0.12 0.123 0.143 0.124 0.125 0.118 0.165 0.142 0.141 0.111 0.157]; %原料O2的摩尔分数 Y1=[0.015 0.016 0.015 0.017 0.02 0.021 0.021 0.021 0.021 0.025 0.023 0.025 0.018 0.018 0.025 0.026 0.026 0.03 0.03 0.023 0.028 0.036 0.034 0.038]; %出口DMC的摩尔数 Y2=[0.018 0.02 0.018 0.033 0.04 0.041 0.042 0.042 0.045 0.054 0.053 0.054 0.031 0.032 0.05 0.051 0.05 0.059 0.058 0.043 0.051 0.063 0.061 0.066]; %出口CO的摩尔数 k0=[8e4 2e4 -0.5 0.5 1 4e5 3e4 0.5 0.5 ];%(指前因子、活化能、级数)初始量 lb=[0 0 -3 -3 -3 0 0 -3 -3 ];%下限 ub=[+inf +inf +inf +inf +inf +inf +inf +inf +inf];%上限 %使用lsqnonlin函数进行参数估计 [k,resnorm,residual,exitflag,output,lambda,jacobian]=lsqnonlin(@ObjFunc,k0,lb,ub,[],T,P,N0,x1,x2,x3,Y1,Y2); ci=nlparci(k,residual,jacobian); function f=ObjFunc(k,T,P,N0,x1,x2,x3,Y1,Y2) wspan=[0 0.375]; Y0=[Y1' Y2']; X=[N0' P' T' x1' x2' x3' Y1' Y2'] [m,n]=size(Y0); Ycal = zeros(m,n); for i=1:m [w yy]=ode45(@Equations,wspan,Y0(i, ,[],k);Ycal(i, =yy(end,;end f=Ycal-Y0; function dydw=Equations(w,yy,k,N0,P,T,x1,x2,x3) global T P N0 x1 x2 x3; %如何表达呢？这么多数字 r1=k(1)*exp(-k(2)./T).*P.^(k(3)+k(4)+k(5))*... (yy(1).*(1.5*x1-2)+0.5*yy(2).*x1).^(k(3))*... (yy(1).*(1.5*x2-1)+yy(2).*(0.5*x2-1)+x2).^(k(4))*... (yy(1).*(1.5*x3-0.5)+yy(2).*(0.3*x3-0.5)+x3).^k(5); r2=k(6)*exp(-k(7)./T).*P.^(k(8)+k(9))*... (yy(1).*(1.5*x2-1)+yy(2).*(0.5*x2-1)+x2).^k(8)*... (yy(1).*(1.5*x3-0.5)+yy(2).*(0.3*x3-0.5)+x3).^k(9); dy1dw=(1+1.5*yy(1)+0.5*yy(2)).*(r1.*(1.5*yy(1)+1)+0.5*yy(1).*r2)./N0; dy2dw=(1+1.5*yy(1)+0.5*yy(2)).*(r2.*(1+0.5*yy(2))+1.5*r1.*yy(2))./N0; dydw=[dy1dw;dy2dw]; 图像为出现的错误结果  图1.jpg |
4楼2015-03-05 09:36:48
40