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

24小时热门版块排行榜    

小木虫论坛-学术科研互动平台 » 计算模拟区 » 仿真模拟 » MATLAB » 【求助】修改程序
31/1返回列表
查看: 586  |  回复: 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返回列表
普通表情 龙 兔 虎 猫 高级回复 (可上传附件)
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[考研] 考研复试调剂，过国家线的同学都可报名 +4 黑！在干嘛 2026-02-28 5/250 2026-03-02 16:48 by chock1337
[考研] 清华大学 材料与化工 353分求调剂 +3 awaystay 2026-03-02 4/200 2026-03-02 16:35 by chuocheng
[考研] 295求调剂。一志愿报考郑州大学化学工艺学硕，总分295分 +5 yl1 2026-03-02 5/250 2026-03-02 15:24 by sucesssucess
[考研] 303求调剂 +5 今夏不夏 2026-03-01 5/250 2026-03-02 15:01 by 向上的胖东
[考研] 290分材料工程085601求调剂 数二英一 +3 llx0610 2026-03-02 3/150 2026-03-02 14:15 by yc258
[考研] 化工专硕348，一志愿985求调剂 +6 弗格个 2026-02-28 9/450 2026-03-02 14:09 by liyongv
[考研] 化工270求调剂 +8 什么名字qwq 2026-03-02 8/400 2026-03-02 13:03 by houyaoxu
[考研] 一志愿山东大学材料与化工325求调剂 +3 半截的诗0927 2026-03-02 3/150 2026-03-02 12:58 by houyaoxu
[基金申请] 成果系统访问量大，请15分钟后再尝试。由此给您造成的不便，敬请谅解。 +5 xhuama 2026-03-02 5/250 2026-03-02 12:34 by stidwellNK
[考研] 265分求调剂不调专业和学校有行学上就 +6 礼堂丁真258 2026-02-28 9/450 2026-03-02 12:04 by 52hz~~
[考研] 268求调剂 +4 简单点0 2026-03-02 5/250 2026-03-02 11:54 by ms629
[考研] 求调剂 +3 熬夜的猫头鹰 2026-03-02 3/150 2026-03-02 11:45 by 刘兵
[考研] 材料类求调剂 +11 wana_kiko 2026-02-28 14/700 2026-03-02 08:46 by 聪明的大松鼠
[基金申请] 本子写完了，给DS兄弟看了，得了92分 +3 Doma 2026-03-01 7/350 2026-03-02 00:00 by jnzsy
[硕博家园] 博士自荐 +7 科研狗111 2026-02-26 11/550 2026-03-01 22:24 by 哲平L
[考研] 298求调剂 +6 axyz3 2026-02-28 6/300 2026-03-01 19:00 by 18137688336
[考研] 291分工科求调剂 +9 science饿饿 2026-03-01 10/500 2026-03-01 18:55 by 18137688336
[考研] 285求调剂 +8 满头大汗的学生 2026-02-28 8/400 2026-03-01 16:47 by caszguilin
[考研] 298求调剂 +9 人间唯你是清欢 2026-02-28 12/600 2026-03-01 14:23 by Ducount.Y
[论文投稿] 求助coordination chemistry reviews 的写作模板 10+3 ljplijiapeng 2026-02-27 4/200 2026-03-01 09:07 by babero
信息提示
关闭
请填处理意见
关闭