24小时热门版块排行榜

返回列表

当前只显示满足指定条件的回帖，点击这里查看本话题的所有回帖

lxyy

银虫 (小有名气)

应助: 0 (幼儿园)
金币: 237.1
散金: 480
红花: 9
帖子: 182
在线: 210.5小时
虫号: 1131417
注册: 2010-10-25
专业: 化学反应工程

[求助] matlab拟合方程参数时初值的选择

A+B=C  A+W=R  R+B=C+W
先考虑的三个反应都是可逆反应，有以下几个式子：编写程序如下
%30℃
clear all
clc
global keq1 keq2 t0 c0

keq1 = 3.94;       % 平衡常数k1
keq3 = 1.08; % 平衡常数k3

t0 = [5,10, 15, 30, 45, 60,90,120,160,200,240,300,360];
c0 = [5.939601,5.203973,0,0.022484,0.078956]       % A B C R W 组成的初始浓度
ci=[3.742508811 4.213566482 0.990406938 0.252498405 0.090415512;
3.045716838 3.387013322 1.816960097 0.215546741 0.081072022;
2.504754523 2.901991595 2.301981824 0.193427102 0.099462604;
1.65098702 2.036686653 3.167286766 0.170236582 0.095144044;
1.271633833 1.655460403 3.548513016 0.124730781 0.10334349;
1.129887499 1.630815394 3.573158025 0.129798695 0.097764543;
1.01239208 1.460748979 3.74322444 0.144649391 0.073614958;
1.007417603 1.429180292 3.774793127 0.123477124 0.094249307;
0.990005723 1.499216939 3.704756481 0.135524067 0.076072022;
0.970073726 1.493966669 3.71000675 0.130791851 0.076387812;
0.96903024 1.429705429 3.77426799 0.119510017 0.101202216;
0.980873385 1.457988124 3.745985296 0.124189274 0.097022161;
0.98477262 1.453402965 3.750570455 0.1200936 0.0762608];

% cA cC cR 对应时间数值

k0=[1 1 1]; % k1+, k2, k3+ 初始向值
lb = [0 0 0];
ub = [1000  1000 500];  % 上下限

% 使用函数lsqnonlin()进行参数估计
[k,resnorm,residual,exitflag] = lsqnonlin(@ObjFunc,k0,lb,ub,[],ci);

k1plus=k(1);
k1minus= k1plus/keq1;
k2=k(2);
k3plus=k(3);
k3minus= k3plus/keq3;

% ------------------------------------------------------------------
function f = ObjFunc(k,ci) % 目标函数
global t0  c0

[t,c_cal] = ode45(@Euqations,t0,c0,[],k);
f= sum((c_cal-ci)^2);

% ------------------------------------------------------------------
function dcdt = Euqations(t,c,k)
global keq1 keq2

k1plus=k(1);
k1minus= k1plus/keq1;
k2=k(2);
k3plus=k(3);
k3minus= k3plus/keq3;

cA=c(:,1);
cB=c(:,2);
cC=c(:,3);
cR=c(:,4);
cW=c(:,5);

dcAdt = k1plus*cA*cB-k1minus*cC+k2*cA*cW;
dcCdt = -k1plus*cA*cB+k1minus*cC+k3plus*cB*cR-k3minus*cC*cW;
dcRdt = k2*cA*cW-k3plus*cB*cR+k3minus*cC*cW;
dcdt = [dcAdt; dcCdt;dcRdt];
但是k的初值一直影响结果，程序要怎么改一下初值对结果影响不大（这个程序有问题，）

回复此楼

» 猜你喜欢

第一性原理计算方向2026级博士申请 PRB*1，四级484 已经有1人回复
求助VISSIM破解版软件已经有0人回复
物理学I论文润色/翻译怎么收费? 已经有263人回复
求2026年在台湾举行的物理和材料领域国际学术会议信息已经有0人回复
求国际会议网站已经有1人回复
求取一些关于纳米材料和纳米技术相关的英文PPT。已经有0人回复
【复旦大学】二维材料方向招收2026年博士研究生1名已经有0人回复
北京纳米能源与系统研究所王中林院士/曹南颖研究员课题组2026级硕/博/博后招生已经有10人回复
荷兰Utrecht University超快太赫兹光谱王海教授课题招收2026 CSC博士生已经有15人回复
反铁磁体中的磁性切换：两种不同的机制已成功可视化已经有0人回复

1楼2012-03-28 10:24:19

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

fielding

木虫 (正式写手)

应助: 5 (幼儿园)
金币: 3544.5
散金: 135
红花: 4
帖子: 672
在线: 219.3小时
虫号: 822444
注册: 2009-08-05
性别: GG
专业: 能源化工

【答案】应助回帖

★ ★ ★ ★ ★ ★ ★ ★ ★ ★
感谢参与，应助指数 +1
lxyy: 金币+10, ★★★很有帮助 2012-03-30 12:30:33

function KineticsEst
clear all
clc

x0 = [4.639975 5.289842 0 0.039974 0.063806];       % A B C R W 组成的初始浓
k0=[1 1 1]; % k1, k2, k3 初始向值
lb = [0 0 0];
ub = [1  +inf +inf];  % 上下限
KineticsData3;
yexp = Kinetics(:,2:6);
% 使用函数lsqnonlin()进行参数估计
[k,resnorm,residual,exitflag,output,lambda,jacobian]=...
lsqnonlin(@ObjFunc,k0,lb,ub,[],x0,yexp);

fprintf('\n\n使用函数lsqnonlin()估计得到的参数值为:\n')
fprintf('\tk1 = %.4f ± %.4f\n',k(1))
fprintf('\tk2 = %.4f ± %.4f\n',k(2))
fprintf('\tk3 = %.4f ± %.4f\n',k(3))
% ------------------------------------------------------------------
function f = ObjFunc(k,x0,yexp) % 目标函数
tspan = [0;5;10;15;30;46;60;90;120;160;200;240;300;360];
[t x]  = ode45(@Euqations,tspan,x0,[],k);
y(:,1) = x(:,1);
y(:,2:5) = x(:,2:5);
f1 = y(:,1) - yexp(:,1);
f2 = y(:,2) - yexp(:,2);
f3 = y(:,3) - yexp(:,3);
f4 = y(:,4) - yexp(:,4);
f5 = y(:,5) - yexp(:,5);
f = [f1; f2; f3; f4;f5];
% ------------------------------------------------------------------
function dcdt = Euqations(t,x,k)

r1 = k(1)*x(1)*x(2)-k(1)*x(3)/3.94;
r2 = k(2)*x(1)*x(5);
r3 = k(3)*x(2)*x(4)-k(3)*x(3)*x(5)/1.08;
dcAdt = -r1-r2;
dcBdt = -r1-r3;
dcCdt = r1+r3;
dcRdt = r2-r3;
dcWdt = -r2+r3;
dcdt = [dcAdt;dcBdt; dcCdt;dcRdt;dcWdt;];