24小时热门版块排行榜    

CyRhmU.jpeg
查看: 1990  |  回复: 3

kylexy3792

银虫 (小有名气)


[交流] 【求助】催化反应动力学matlab计算各基元反应的速率常数时,该如何避免较小量被忽略?

催化反应动力学matlab计算各基元反应的速率常数时,该如何避免较小量被忽略?
我按照个基元反应,依据控制速率步骤,得到相应速率表达式,从而求的各个基元反应的速率常数,但是在设定初值后,结果再次迭代最后算的结果,但是把最后的结果当初值,同时改掉其中一个变量,有些变量对最后结果影响很大,有些却几乎不影响其他结果,求高手帮忙解决!

同时求其他更好的解动力学方程组的方法!
回复此楼

» 猜你喜欢

» 本主题相关价值贴推荐,对您同样有帮助:

» 抢金币啦!回帖就可以得到:

查看全部散金贴

已阅   关注TA 给TA发消息 送TA红花 TA的回帖

kylexy3792

银虫 (小有名气)


jinkai838: 什么反应?能具体么?可以给我发文献,可以联系我 2011-03-06 08:15:20
没有人帮忙看看,解决一下么?

clc;clear;
format long
C=[18000        572.7587379        17427.24126        503.100394        63.00150801        5.429614942        1.227220877        5.34957E-05        31.80663631        28.64788272        2.468940458        0.558038702
18000        442.6621111        17557.33789        391.4982651        46.31421905        3.97596704        0.873659836        2.99887E-05        25.4952724        22.99518737        1.974082883        0.433775459
18000        395.8533613        17604.14664        353.0097887        38.82819644        3.315565494        0.699810656        2.01511E-05        23.85648271        21.55796881        1.840849281        0.388544864
18000        345.603434        17654.39657        311.2417414        31.19218479        2.643727575        0.525780251        1.48443E-05        21.89345442        19.83639301        1.681255082        0.334365283
18000        326.7498543        17673.25015        295.4076551        28.47206709        2.405087713        0.465044362        2.99887E-05        21.11229134        19.15131299        1.617746524        0.312805182
18000        268.2211638        17731.77884        245.416994        20.77111267        1.733166874        0.299890213        2.01511E-05        18.68215469        17.02008314        1.420176413        0.24573341
18000        1386.82292        16613.17708        1331.139206        101.6813129        8.431488388        1.370722461        1.48443E-05        88.51218746        80.81389331        6.701146786        1.089417668
18000        1109.037017        16890.96298        1063.169464        83.66013195        6.960106488        1.169991466        2.99887E-05        72.90022276        66.48318761        5.531070232        0.929770971
18000        557.2888809        17442.71112        532.9063777        44.3791472        3.714372299        0.661563426        2.01511E-05        38.75233847        35.2670053        2.951719351        0.525728013
18000        276.5175885        17723.48241        264.074297        22.62549864        1.899122798        0.347306091        1.48443E-05        19.78050922        17.98334006        1.509472637        0.276047995
];
P_P=C(:,4);
P_O=C(:,5);
P_D=C(:,8);
P_H=C(:,3);
P_A=C(:,7);

W_P=C(:,9);
W_O=C(:,10);
W_D=C(:,11);
W_A=C(:,12);

% W_P=(K_1.*P_P-K_a.*P_O.*P_H)./A;
% W_O=(K_1.*P_P-K_a.*P_O.*P_H-K_2.*P_O+K_b.*P_D.*P_H)./A;
% W_D=(K_2.*P_O-K_b.*P_D.*P_H)./A;

xx=[P_P P_O P_D P_H W_P W_O W_D W_A P_A];

yy=[0 0 0
    0 0 0
    0 0 0
    0 0 0
    0 0 0
    0 0 0
    0 0 0
    0 0 0
    0 0 0
    0 0 0];

xdata = xx;
ydata = yy;
x0=[9.99E+09        1.410771372584031e+006        1.250036334940072e+010        1.550570624167701e+007        2.97E+10        8.250392906170209e+002        3.468987506591585e+002        2.790032030520837e+002        4.789732864204271e+002  9.528956235795578e+003
]
%options = optimset('TolFun',1e-3);%TolFun是函数误差
%%不知道你说的是哪个误差。如果是‘结果误差’就把上句的TolFun给改成TolX x = lsqcurvefit(fun,x0,xdata,ydata,options);
%$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
x=lsqcurvefit(@myfun310,x0',xdata,ydata)
ww10=(x(1).*x(6).*(xdata(:,1)-xdata(:,2).*xdata(:,4)./x(1).*x(2)))./(1+x(6).*xdata(:,1)+x(7).*xdata(:,2)+x(9).*xdata(:,4)+x(8).*xdata(:,3)).^2-xdata(:,5);
ww20=(x(1).*x(6).*(xdata(:,1)-xdata(:,2).*xdata(:,4)./x(1).*x(2)))./(1+x(6).*xdata(:,1)+x(7).*xdata(:,2)+x(9).*xdata(:,4)+x(8).*xdata(:,3)).^2-((x(3).*x(7).*(xdata(:,2)-xdata(:,3).*xdata(:,4)./x(3).*x(4)))./(1+x(6).*xdata(:,1)+x(7).*xdata(:,2)+x(9).*xdata(:,4)+x(8).*xdata(:,3)).^2)-xdata(:,6);
ww30=1*((x(3).*x(7).*(xdata(:,2)-xdata(:,3).*xdata(:,4)./x(3).*x(4)))./(1+x(6).*xdata(:,1)+x(7).*xdata(:,2)+x(9).*xdata(:,4)+x(8).*xdata(:,3)).^2-xdata(:,8)-xdata(:,7));
    x=lsqcurvefit(@myfun310,x0',xdata,ydata)%,options)
    ww1=(x(1).*x(6).*(xdata(:,1)-xdata(:,2).*xdata(:,4)./x(1).*x(2)))./(1+x(6).*xdata(:,1)+x(7).*xdata(:,2)+x(9).*xdata(:,4)+x(8).*xdata(:,3)).^2-xdata(:,5);
    ww2=(x(1).*x(6).*(xdata(:,1)-xdata(:,2).*xdata(:,4)./x(1).*x(2)))./(1+x(6).*xdata(:,1)+x(7).*xdata(:,2)+x(9).*xdata(:,4)+x(8).*xdata(:,3)).^2-((x(3).*x(7).*(xdata(:,2)-xdata(:,3).*xdata(:,4)./x(3).*x(4)))./(1+x(6).*xdata(:,1)+x(7).*xdata(:,2)+x(9).*xdata(:,4)+x(8).*xdata(:,3)).^2)-xdata(:,6);
    ww3=1*((x(3).*x(7).*(xdata(:,2)-xdata(:,3).*xdata(:,4)./x(3).*x(4)))./(1+x(6).*xdata(:,1)+x(7).*xdata(:,2)+x(9).*xdata(:,4)+x(8).*xdata(:,3)).^2-xdata(:,8)-xdata(:,7));
    if abs(ww1-ww10)>10^(-6) |abs(ww2-ww20)>10^(-6)| abs(ww3-ww30)>10^(-6)
        ww10=ww1;
        ww20=ww2;
        ww30=ww3;
        x=lsqcurvefit(@myfun310,x0',xdata,ydata)%,options)
    end
    K_1=x(1),K_a=x(2),K_2=x(3),K_b=x(4),K=x(5),K_P=x(6),K_O=x(7),K_D=x(8),K_H=x(9),K_A=x(10)
    KK_1=x(1)/x(2)
    KK_2=x(3)/x(4)
    x0=[x(1) x(2) x(3) x(4) x(5) x(6) x(7) x(8) x(9) x(10)]
%$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$%$$$$$$$$$$
方程为:
function F=myfun310(x,xdata)
format long
F=[(x(1).*x(6).*(xdata(:,1)-xdata(:,2).*xdata(:,4)./x(1).*x(2)))./(1+x(6).*xdata(:,1)+x(7).*xdata(:,2)+x(9).*xdata(:,4)+x(8).*xdata(:,3)).^2-xdata(:,5),(x(1).*x(6).*(xdata(:,1)-xdata(:,2).*xdata(:,4)./x(1).*x(2)))./(1+x(6).*xdata(:,1)+x(7).*xdata(:,2)+x(9).*xdata(:,4)+x(8).*xdata(:,3)).^2-((x(3).*x(7).*(xdata(:,2)-xdata(:,3).*xdata(:,4)./x(3).*x(4)))./(1+x(6).*xdata(:,1)+x(7).*xdata(:,2)+x(9).*xdata(:,4)+x(8).*xdata(:,3)).^2)-xdata(:,6),1*((x(3).*x(7).*(xdata(:,2)-xdata(:,3).*xdata(:,4)./x(3).*x(4)))./(1+x(6).*xdata(:,1)+x(7).*xdata(:,2)+x(9).*xdata(:,4)+x(8).*xdata(:,3)).^2-xdata(:,8)-xdata(:,7))]
有没有好心人帮忙看看啊?

[ Last edited by kylexy3792 on 2011-3-3 at 20:02 ]
2楼2011-03-03 20:00:42
已阅   关注TA 给TA发消息 送TA红花 TA的回帖

kylexy3792

银虫 (小有名气)


没人帮忙解决!
3楼2011-03-05 22:40:16
已阅   关注TA 给TA发消息 送TA红花 TA的回帖

zt_chem

木虫 (正式写手)


★ ★
kylexy3792(金币+5): 谢谢,我去试试看,如果可以计算出结果,再给你追加金币! 2011-03-06 14:02:49
jinkai838(金币+2): 不错!可以具体发帖 2011-03-06 14:42:48
试着用lsqnonlin优化效果是否能好些?
4楼2011-03-06 08:37:35
已阅   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 kylexy3792 的主题更新
信息提示
请填处理意见