| 查看: 527 | 回复: 2 | |||
| 当前只显示满足指定条件的回帖,点击这里查看本话题的所有回帖 | |||
[交流]
【求助】修改程序
|
|||
|
目前下面的程序能得到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); |
» 猜你喜欢
孩子确诊有中度注意力缺陷
已经有13人回复
三甲基碘化亚砜的氧化反应
已经有4人回复
请问下大家为什么这个铃木偶联几乎不反应呢
已经有5人回复
请问有评职称,把科研教学业绩算分排序的高校吗
已经有5人回复
2025冷门绝学什么时候出结果
已经有3人回复
天津工业大学郑柳春团队欢迎化学化工、高分子化学或有机合成方向的博士生和硕士生加入
已经有4人回复
康复大学泰山学者周祺惠团队招收博士研究生
已经有6人回复
AI论文写作工具:是科研加速器还是学术作弊器?
已经有3人回复
论文投稿,期刊推荐
已经有4人回复
硕士和导师闹得不愉快
已经有13人回复
» 抢金币啦!回帖就可以得到:
硫化物全固态电池的产业化破局:手套箱如何实现全线稳定制造
+1/86
加拿大/英属哥伦比亚大学曹彦凯课题组招收全奖博士/博后 [机器学习/优化/控制方向]
+1/86
招贤纳博(已结束)
+1/75
Call for papers,征稿
+1/70
南昌大学药学博士招生
+1/46
燕山大学2026级考核制博士招生——电容去离子方向(超级电容器衍生应用)
+1/36
捷克布拉格查理大学(QS260)招收第一性原理计算方向博士生
+1/34
双一流大学湘潭大学“化工过程模拟与强化”国家地方联合工程研究中心招收各类博士生
+1/34
浙大城市学院向铮教授与浙江大学联合招聘博后
+1/27
2026年度智能交通课题组诚招理工科背景博士
+1/27
武汉大学医学部周翌丹课题组科研助理、博士后招聘
+1/8
长江大学武汉校区诚招新能源博士(工程热物理、电气、油气)-2025
+1/8
伦敦大学学院(University College London)机械工程系博士招生/CSC联培招生
+1/7
澳门科技大学2026年数学博士招生—杨钧翔助理教授计算物理与数学课题组
+1/6
香港城市大学 Microsystems and Semiconductor Technology Lab 博士生招生启事
+1/5
诚招博士后及研究人员
+1/5
沙特法赫德国王石油与矿产大学(KFUPM)膜分离课题组招生
+1/2
澳洲皇家墨尔本理工RMIT招收网安方向CSC PhD
+1/2
中国科学技术大学 精准智能化学重点实验室 武建昌课题组招聘博士,博士后
+1/2
南京航空航天大学核科学与技术方向招收博士生
+1/2
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 ] |
3楼2010-12-15 13:30:40
★ ★
小扶摇(金币+10):太谢谢了,我先研究下,有问题的话还要麻烦下 2010-12-11 09:38:01
zzuwangshilei(金币+2):积极参与,辛苦了 2010-12-11 10:36:53
小扶摇(金币+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 |
2楼2010-12-11 09:28:21














回复此楼