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

24小时热门版块排行榜    

小木虫论坛-学术科研互动平台 » 计算模拟区 » 仿真模拟 » MATLAB » 【求助】修改程序
31/1返回列表
查看: 595  |  回复: 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返回列表
普通表情 龙 兔 虎 猫 高级回复 (可上传附件)
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[考研] 【2026考研调剂】制药工程 284分 求相关专业调剂名额 +4 袁奂奂 2026-03-25 8/400 2026-03-25 14:32 by lbsjt
[考研] 各位老师您好：本人初试372分 +5 jj涌77 2026-03-25 6/300 2026-03-25 14:15 by mapenggao
[考研] 284求调剂 +15 Zhao anqi 2026-03-22 15/750 2026-03-25 12:51 by wht0531
[考研] 一志愿南航材料专317分求调剂 +5 炸呀炸呀炸薯条 2026-03-23 5/250 2026-03-24 16:52 by 星空星月
[考研] 321求调剂 +4 Ymlll 2026-03-24 4/200 2026-03-24 14:44 by sprinining
[基金申请] 请教下大家 2026年国家基金申请是双盲审吗？ +3 lishucheng1 2026-03-22 5/250 2026-03-24 08:22 by gltch
[考研] 335求调剂 +4 yuyu宇 2026-03-23 5/250 2026-03-23 23:49 by Txy@872106
[考研] 335分 | 材料与化工专硕 | GPA 4.07 | 有科研经历 +4 cccchenso 2026-03-23 4/200 2026-03-23 23:00 by 徐ckkk
[考研] 求调剂材料学硕080500，总分289分 5+3 @taotao 2026-03-19 21/1050 2026-03-23 10:17 by 冠c哥
[考研] 一志愿070300浙大化学358分，求调剂！ +4 酥酥鱼.. 2026-03-21 4/200 2026-03-23 08:12 by Iveryant
[考研] 求调剂一志愿海大，0703化学学硕304分，有大创项目，四级已过 +6 幸运哩哩 2026-03-22 10/500 2026-03-22 20:10 by edmund7
[考研] 一志愿中南化学（0703）总分337求调剂 +9 niko- 2026-03-19 10/500 2026-03-22 16:08 by ColorlessPI
[考研] 275求调剂 +6 shansx 2026-03-22 8/400 2026-03-22 15:27 by barlinike
[考研] 269专硕求调剂 +6 金恩贝 2026-03-21 6/300 2026-03-22 14:31 by ColorlessPI
[考研] 考研调剂 +3 呼呼？~+123456 2026-03-21 3/150 2026-03-21 20:04 by 无际的草原
[考研] 297求调剂 +3 喜欢还是不甘心 2026-03-20 3/150 2026-03-21 18:33 by 学员8dgXkO
[考研] 求助 +5 梦里的无言 2026-03-21 6/300 2026-03-21 17:51 by 学员8dgXkO
[考研] 265求调剂 +12 梁梁校校 2026-03-19 14/700 2026-03-21 13:38 by lature00
[考研] 南昌大学材料专硕311分求调剂 +6 77chaselx 2026-03-20 6/300 2026-03-21 07:24 by JourneyLucky
[考研] 321求调剂 +9 何润采123 2026-03-18 11/550 2026-03-20 23:19 by JourneyLucky
信息提示
关闭
请填处理意见
关闭