24小时热门版块排行榜

返回列表

小扶摇

金虫 (正式写手)

应助: 0 (幼儿园)
金币: 644.6
帖子: 341
在线: 206.5小时
虫号: 918128

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

目前下面的程序能得到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);

回复此楼

» 猜你喜欢

上海工程技术大学【激光智能制造】课题组招收硕士已经有6人回复
带资进组求博导收留已经有11人回复
自荐读博已经有5人回复
求个博导看看已经有16人回复
上海工程技术大学张培磊教授团队招收博士生已经有4人回复
求助院士们，这个如何合成呀已经有4人回复
临港实验室与上科大联培博士招生1名已经有9人回复
写了一篇“相变储能技术在冷库中应用”的论文，论文内容以实验为主，投什么期刊合适？已经有6人回复
最近几年招的学生写论文不引自己组发的文章已经有11人回复
中科院杭州医学所招收博士生一名（生物分析化学、药物递送）已经有3人回复

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

查看全部散金贴

1楼 2010-12-10 15:32:04

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

snowman09

金虫 (小有名气)

应助: 0 (幼儿园)
金币: 1619.1
帖子: 148
在线: 125.1小时
虫号: 1060657

★ ★
小扶摇(金币+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的回帖

小扶摇

金虫 (正式写手)

应助: 0 (幼儿园)
金币: 644.6
帖子: 341
在线: 206.5小时
虫号: 918128

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的回帖

相关版块跳转我要订阅楼主小扶摇的主题更新

返回列表

24小时热门版块排行榜

小扶摇

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

» 猜你喜欢

» 抢金币啦！回帖就可以得到: 查看全部散金贴

snowman09

小扶摇

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

查看全部散金贴