应《网络安全法》要求，自2017年10月1日起，未进行实名认证将不得使用互联网跟帖服务。为保障您的帐号能够正常使用，请尽快对帐号进行手机号验证，感谢您的理解与支持！

24小时热门版块排行榜    

小木虫论坛-学术科研互动平台 » 计算模拟区 » 仿真模拟 » MATLAB » 【求助】修改程序
31/1返回列表
查看: 596  |  回复: 2
只看楼主 @他人 存档 新回复提醒 (忽略) 收藏    在APP中查看

小扶摇

金虫 (正式写手)

[交流] 【求助】修改程序

目前下面的程序能得到r=kPA^ n PB^ m中的k、n、m值，但是我想把程序修改成能直接求r=k0e(-Ea/RT)PA^nPB^m，中的k0、Ea、n和m。应该如何修改程序，求大神赐教，不甚感激。
数据如下：
T=310+273.15K时
PA0 = [3.7449  3.7500  3.7548  3.7568];
PB0 = [0.0807  0.0835  0.0861  0.0871];
r  = [0.0019  0.0022  0.0024  0.0025];
T=330+273.15K时
PA0 = [3.7138  3.7215  3.7448  3.7572];
PB0 = [0.04  0.0458  0.0618  0.071];
r  = [0.0041  0.0097  0.0109  0.0111];
Ps：R是常数R=8.3145
r=kPA^nPB^m的程序
function KineticsEst3
% 动力学参数辨识: 用微分法进行反应速率分析得到速率常数k和反应级数n
% Reaction of the type -- rate = kCA^order
% order - reaction order
% rate -- reaction rate vector
% CA -- concentration vector for reactant A
% T -- vector of reaction time
% N -- number of data points
% k- reacion rate constant

clear all
clc
global negr0m PA0 PB0
PA0 = [3.7449  3.7500  3.7548  3.7568];
PB0 = [0.0807  0.0835  0.0861  0.0871];
% negr0m = [0.5  0.63  0.83  1.0  1.28  0.33  0.8  1.5  2.21  3.33];

% Problems ...
negr0m = [0.0019  0.0022  0.0024  0.0025];

% 用多变量线性回归方法估计动力学参数
y = log(negr0m);
x1 = log(PA0);
x2 = log(PB0);
y = y';
X = [ones(size(y))  x1'  x2'];
[b,bint] = regress(y,X,0.1);
k = exp(b(1))
n = b(2)
m = b(3)

% 用多变量非线性回归方法（以线性回归的结果作为非线性回归的初值）
beta0 = [k  n  m];
lb = [0 0 0];
ub = [1 3 3];
[beta,resnorm,resid,exitflag,output,lambda,jacobian] = ...
     lsqnonlin(@ObjFunc,beta0,lb,ub);
ci = nlparci(beta,resid,jacobian)

% 残差关于拟合值的残差图
negrA0c = Rate(beta,PA0,PB0);
plot(negrA0c,resid,'*')
xlabel('反应速率拟合值, torr s^-^1')
ylabel('残差R, torr s^-^1')
refline(0,0)

% 参数辨识结果
fprintf('Estimated Parameters:\n')
fprintf('\tk = %.4f ± %.4f\n',beta(1),ci(1,2)-beta(1))
fprintf('\tn = %.2f ± %.2f\n',beta(2),ci(2,2)-beta(2))
fprintf('\tm = %.2f ± %.2f\n',beta(3),ci(3,2)-beta(3))
fprintf('\tThe sum of the squares is: %.1e\n\n',resnorm)

% ------------------------------------------------------------------
function f = ObjFunc(beta)
global negr0m PA0 PB0
f = negr0m - Rate(beta,PA0,PB0);

% ------------------------------------------------------------------
function negrA = Rate(beta,PA,PB)
negrA = beta(1)*PA.^beta(2).*PB.^beta(3);   % k=beta(1),n=beta(2),m=beta(3);
回复此楼

» 猜你喜欢

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

查看全部散金贴
1楼 2010-12-10 15:32:04
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

snowman09

金虫 (小有名气)
★ ★
小扶摇(金币+10):太谢谢了，我先研究下，有问题的话还要麻烦下 2010-12-11 09:38:01
zzuwangshilei(金币+2):积极参与,辛苦了 2010-12-11 10:36:53
我把你的程序修改了一下，能运行了，但是需要增加一组实验数据，如果不懂，请给我短消息~


function KineticsEst3_1
% 动力学参数辨识: 用微分法进行反应速率分析得到速率常数k和反应级数n
% Reaction of the type -- rate = kCA^order
% order - reaction order
% rate -- reaction rate vector
% CA -- concentration vector for reactant A
% T -- vector of reaction time
% N -- number of data points
% k- reacion rate constant

clear all
clc
global negr0m PA0 PB0 RT
PA0 = [3.7449  3.7500  3.7548  3.7568 3.7588 ];
PB0 = [0.0807  0.0835  0.0861  0.0871 0.0891];
RT=8.3*(310+273)*ones(1,5);
% negr0m = [0.5  0.63  0.83  1.0  1.28  0.33  0.8  1.5  2.21  3.33];

% Problems ...
negr0m = [0.0019  0.0022  0.0024  0.0025 0.0028];

% 用多变量线性回归方法估计动力学参数
y = log(negr0m);
x1 = log(PA0);
x2 = log(PB0);
x3 =log(8.3*(310+273))*ones(1,5);
y = y';
X = [ones(size(y))  x1'  x2' x3'];
[b,bint] = regress(y,X,0.1);
k0 = exp(b(1))
n = b(2)
m = b(3)
Ea = b(4)

% 用多变量非线性回归方法（以线性回归的结果作为非线性回归的初值）
beta0 = [k0  n  m Ea];
lb = [0 0 0 -10];
ub = [1 3 3 10];
[beta,resnorm,resid,exitflag,output,lambda,jacobian] = ...
     lsqnonlin(@ObjFunc,beta0,lb,ub);
ci = nlparci(beta,resid,jacobian)

% 残差关于拟合值的残差图
negrA0c = Rate(beta,PA0,PB0,RT);
plot(negrA0c,resid,'*')
xlabel('反应速率拟合值, torr s^-^1')
ylabel('残差R, torr s^-^1')
refline(0,0)

% 参数辨识结果
fprintf('Estimated Parameters:\n')
fprintf('\tk = %.4f ± %.4f\n',beta(1),ci(1,2)-beta(1))
fprintf('\tn = %.2f ± %.2f\n',beta(2),ci(2,2)-beta(2))
fprintf('\tm = %.2f ± %.2f\n',beta(3),ci(3,2)-beta(3))
fprintf('\tEa = %.2f ± %.2f\n',beta(4),ci(4,2)-beta(4))
fprintf('\tThe sum of the squares is: %.1e\n\n',resnorm)

%------------------------------------------------------------------
function f = ObjFunc(beta)
global negr0m PA0 PB0 RT
f = negr0m - Rate(beta,PA0,PB0,RT);

% ------------------------------------------------------------------
function negrA = Rate(beta,PA,PB,RT)
negrA = beta(1)*PA.^beta(2).*PB.^beta(3).*exp(beta(4)./RT);   % k=beta(1),n=beta(2),m=beta(3);3.75883.7588
赞一下(1人) 回复此楼
2楼2010-12-11 09:28:21
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

小扶摇

金虫 (正式写手)
robert2020:建议使用“引用回复该帖”，不然对方收不到你的信息。 2010-12-16 09:15:35
clear all
clc
global negr0m PA0 PB0 RT
PA0 = [3.7377 3.7449  3.7500  3.7548  3.7568  3.5690 3.5859 3.5970 3.6046 3.6119 3.6161];
PB0 = [0.0769 0.0807  0.0835  0.0861  0.0871  0.0389 0.0468 0.0520 0.0554 0.0581 0.0601];
RT1=8.3145*(310+273.15)*ones(1,5);
RT2=8.3145*(330+273.15)*ones(1,6);
RT=[RT1 RT2];                                       
% negr0m = [0.5  0.63  0.83  1.0  1.28  0.33  0.8  1.5  2.21  3.33];

% Problems ...
negr0m = [0.0023 0.0024  0.0024  0.0025  0.0025  0.0040 0.0042  0.0044  0.0046 0.0048 0.0049];
% 用多变量线性回归方法估计动力学参数
y = log(negr0m);
x1 = log(PA0);
x2 = log(PB0);
x3 =log(RT);
y = y';
X = [ones(size(y))  x1'  x2' x3'];
[b,bint] = regress(y,X,0.1);
k0 = exp(b(1))
n = b(2)
m = b(3)
Ea = b(4)

% 用多变量非线性回归方法（以线性回归的结果作为非线性回归的初值）
beta0 = [k0  n  m Ea];
lb = [0 0 0 -10];
ub = [1 3 3 10];
[beta,resnorm,resid,exitflag,output,lambda,jacobian] = ...
     lsqnonlin(@ObjFunc,beta0,lb,ub);
ci = nlparci(beta,resid,jacobian)

% 残差关于拟合值的残差图
negrA0c = Rate(beta,PA0,PB0,RT);
plot(negrA0c,resid,'*')
xlabel('反应速率拟合值, torr s^-^1')
ylabel('残差R, torr s^-^1')
refline(0,0)

% 参数辨识结果
fprintf('Estimated Parameters:\n')
fprintf('\tk = %.4f ± %.4f\n',beta(1),ci(1,2)-beta(1))
fprintf('\tn = %.2f ± %.2f\n',beta(2),ci(2,2)-beta(2))
fprintf('\tm = %.2f ± %.2f\n',beta(3),ci(3,2)-beta(3))
fprintf('\tEa = %.2f ± %.2f\n',beta(4),ci(4,2)-beta(4))
fprintf('\tThe sum of the squares is: %.1e\n\n',resnorm)

%------------------------------------------------------------------
function f = ObjFunc(beta)
global negr0m PA0 PB0 RT
f = negr0m - Rate(beta,PA0,PB0,RT);

% ------------------------------------------------------------------
function negrA = Rate(beta,PA,PB,RT)
negrA = beta(1)*PA.^beta(2).*PB.^beta(3).*exp(-beta(4)./RT);   % k=beta(1),n=beta(2),m=beta(3);Ea=beta(4);





为什么最终的结果和
lb = [0 0 0 -10];
ub = [1 3 3 10];
有很大关系，该如何解决

[ Last edited by 小扶摇 on 2010-12-15 at 13:35 ]
赞一下(1人) 回复此楼
3楼2010-12-15 13:30:40
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 小扶摇 的主题更新
31/1返回列表
普通表情 龙 兔 虎 猫 高级回复 (可上传附件)
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[考研] 求b区院校调剂 +3 周56 2026-03-24 4/200 2026-03-25 13:17 by 周56
[考研] 求调剂 一志愿 本科 北科大 化学 343 +4 13831862839 2026-03-24 5/250 2026-03-25 09:47 by 无际的草原
[考研] 0856材料专硕353求调剂 +6 NIFFFfff 2026-03-20 6/300 2026-03-24 21:02 by hello七七
[考研] 070300化学求调剂 +9 苑豆豆 2026-03-20 9/450 2026-03-24 17:15 by licg0208
[材料工程] 一志愿C9材料与化工专业总分300求调剂 +4 曼111 2026-03-24 5/250 2026-03-24 15:44 by 星空星月
[考研] 307求调剂 +3 余意卿 2026-03-21 6/300 2026-03-24 15:03 by 余意卿
[考研] 291求调剂 +8 hhhhxn.. 2026-03-23 8/400 2026-03-23 23:15 by peike
[考研] 一志愿武理材料工程348求调剂 +6 ￣^￣゜汗 2026-03-19 9/450 2026-03-23 19:53 by pswait
[考研] 316求调剂 +7 梁茜雯 2026-03-19 7/350 2026-03-23 16:21 by lingjue
[考研] 接收2026硕士调剂(学硕+专硕) +4 allen-yin 2026-03-23 6/300 2026-03-23 15:04 by 汪！？！
[考研] 323求调剂 +6 洼小桶 2026-03-18 6/300 2026-03-23 00:29 by king123！
[考研] 石河子大学（211、双一流）硕博研究生长期招生公告 +3 李子目 2026-03-22 3/150 2026-03-22 21:01 by 怎么释怀
[考研] 一志愿西安交通大学材料工程专业 282分求调剂 +11 枫桥ZL 2026-03-18 13/650 2026-03-22 20:26 by edmund7
[考研] 寻找调剂 +4 倔强芒? 2026-03-21 4/200 2026-03-22 16:14 by 木托莫露露
[考研] 一志愿中南化学（0703）总分337求调剂 +9 niko- 2026-03-19 10/500 2026-03-22 16:08 by ColorlessPI
[考研] 324求调剂 +6 lucky呀呀呀鸭 2026-03-20 6/300 2026-03-22 16:01 by ColorlessPI
[考研] 275求调剂 +6 shansx 2026-03-22 8/400 2026-03-22 15:27 by barlinike
[考研] 336求调剂 +5 rmc8866 2026-03-21 5/250 2026-03-21 17:24 by 学员8dgXkO
[考研] 290求调剂 +7 ^O^乜 2026-03-19 7/350 2026-03-20 21:43 by JourneyLucky
[考研] 本科郑州大学物理学院，一志愿华科070200学硕，346求调剂 +4 我不是一根葱 2026-03-18 4/200 2026-03-19 09:11 by 浮云166
信息提示
关闭
请填处理意见
关闭