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

24小时热门版块排行榜    

小木虫论坛-学术科研互动平台 » 计算模拟区 » 仿真模拟 » MATLAB » 【求助】修改程序
51/1返回列表
查看: 542  |  回复: 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的回帖

小扶摇

金虫 (正式写手)
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的回帖
查看全部 3 个回答

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的回帖
查看全部 3 个回答
普通表情 龙 兔 虎 猫 高级回复 (可上传附件)
信息提示
关闭
请填处理意见
关闭