求助用matlab三次样条法算一个反应的反应速率 - 计算模拟

版块导航: 正在加载中...

登录注册

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

24小时热门版块排行榜

返回列表

justin534

新虫 (初入文坛)

应助: 0 (幼儿园)
金币: 175.5
帖子: 22
在线: 24小时
虫号: 3005267
注册: 2014-02-28
性别: GG
专业: 化学反应工程

???????:

10?: Originally posted by ????? at 2014-03-29 15:11:44
?????????????????????????????????
?????????????dingd ??1stopt???????????????????...

???????,???!????????}??!??????????????????W?,?F????????3??,??6??k1 k2 k3 k4 k5 k6!
R1=k(1)*DEC*PA-(k(1)/k(2))*EPC*EA
R2=k(3)*EPC*PA-(k(3)/k(4))*DPC*EA
R3=k(5)*EPC.^2-(k(5)/k(6))*DPC*DEC
???????????????a!???????RUN,??????????!

function ODE_parafit2
%-------data
yexp=[5.14  12.76  45.41  21.52  3.31;
5.36  13.20  46.94  20.90  3.17;
5.51  13.47  48.68  20.01  3.08;
5.78  14.03  48.16  19.34  3.04;
6.07  14.33  47.80  18.65  2.94]';
x0=[93.7 6.3 0 0 0];
k0 = [1 1 1 1 1 1];
lb = [0 0 0 0 0 0];
ub = [1 1 1 1 1 1]*1e5;
tspan = [0 38.46 50 62.5 100 166.666 250 500];
%-----------------
[k,fval,flag] =fmincon(@ObjFuncFmincon,k0,[],[],[],[],lb,ub,[],[],tspan,x0,yexp);
fprintf('\n\n??ú???fmincon()?????????????:\n')
fprintf('\tk1 = %.16f \n',k(1))
fprintf('\tk2 = %.16f \n',k(2))
fprintf('\tk3 = %.16f \n',k(3))
fprintf('\tk4 = %.16f \n',k(4))
fprintf('\tk5 = %.16f \n',k(5))
fprintf('\tk6 = %16f \n',k(6))
fprintf('The sum of squares is: %.16f \n\n',fval)

[k,resnorm,residual,exitflag,output,lambda,jacobian] =lsqnonlin(@ObjFunc,k0,lb,ub,[],tspan,x0,yexp);
ci = nlparci(k,residual,jacobian);
fprintf('\n\n??ú???lsqnonlin()?????????????:\n')
fprintf('\tk1 = %.16f \n',k(1))
fprintf('\tk2 = %.16f \n',k(2))
fprintf('\tk3 = %.16f \n',k(3))
fprintf('\tk4 = %.16f \n',k(4))
fprintf('\tk5 = %.16f \n',k(5))
fprintf('\tk6 = %16f \n',k(6))
fprintf('The sum of squares is: %.16f \n\n',resnorm)

%%  plot figure
tspan = [0 38.46 50 62.5 100 166.666 250 500];
for i=1:5
[t x] = ode45(@KineticsEqs,tspan,x0,[],k);
fy(i,

=x(end,

;
end
Ypre=fy(

;
Yexp=yexp(

;
plot(Ypre,'or',Yexp,'b.')
legend('prediction','experiment')
%----------fmincon fun-------------------
function f = ObjFuncFmincon(k,x0,yexp)
tspan = [0 38.46 50 62.5 100 166.666 250 500];
for i=1:5
[t x] = ode45(@KineticsEqs,tspan,x0,[],k);
fy(i)=sum((yexp(i,

-x(end,

).^2);
end
f = sum(fy) ;
% -------------------lsqnonlin fun-------------------
function f = ObjFunc(k,x0,yexp)
tspan = [0 38.46 50 62.5 100 166.666];
for i=1:5
[t x] = ode45(@KineticsEqs,tspan,x0,[],k);
fy(i,

=yexp(i,

-x(end,

;
end
f =fy(

;
%-----------odefun--------------------------------
function dYdt = KineticsEqs(t,y,k)
k1=k(1);k2=k(2);k3=k(3);k4=k(4);k5=k(5);k6=k(6);
DEC=y(1);PA=y(2);DPC=y(3);EA=y(4);EPC=y(5);
R1=k(1)*DEC*PA-(k(1)/k(2))*EPC*EA;
R2=k(3)*EPC*PA-(k(3)/k(4))*DPC*EA;
R3=k(5)*EPC.^2-(k(5)/k(6))*DPC*DEC;
dy1dt=-R1+R3;
dy2dt=-R1-R2;
dy3dt=R2+R3;
dy4dt=R1+R2;
dy5dt=R1-R2-R3;
dYdt=[dy1dt;dy2dt;dy3dt;dy4dt;dy5dt];

赞一下

回复此楼

11楼2014-04-24 14:50:42

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

justin534

新虫 (初入文坛)

应助: 0 (幼儿园)
金币: 175.5
帖子: 22
在线: 24小时
虫号: 3005267
注册: 2014-02-28
性别: GG
专业: 化学反应工程

引用回帖:

10楼: Originally posted by 月只蓝 at 2014-03-29 15:11:44
可能是方程的问题，也可能是实验数据的问题。
你可以请实习版主dingd 用1stopt软件做一下，来验证一下。...

sorry,剛剛的編碼有誤,這個是OK的編碼!還有data是隨意打的,所以只是請你幫我看可不可以RUN!

function ODE_parafit2
%-------data
yexp=[5.14  12.76  45.41  21.52  3.31 1 1;
5.36  13.20  46.94  20.90  3.17 1 1;
5.51  13.47  48.68  20.01  3.08 1 1;
5.78  14.03  48.16  19.34  3.04 1 1;
6.07  14.33  47.80  18.65  2.94 1 1]';
x0=[93.7 6.3 0 0 0];
k0 = [1 1 1 1 1 1];
lb = [0 0 0 0 0 0];
ub = [1 1 1 1 1 1]*1e5;
tspan = [0 38.46 50 62.5 100 166.666 250 500];
%-----------------
[k,fval,flag] =fmincon(@ObjFuncFmincon,k0,[],[],[],[],lb,ub,[],[],tspan,x0,yexp);
fprintf('\n\n使用函数fmincon()估计得到的参数值为:\n')
fprintf('\tk1 = %.16f \n',k(1))
fprintf('\tk2 = %.16f \n',k(2))
fprintf('\tk3 = %.16f \n',k(3))
fprintf('\tk4 = %.16f \n',k(4))
fprintf('\tk5 = %.16f \n',k(5))
fprintf('\tk6 = %16f \n',k(6))
fprintf('The sum of squares is: %.16f \n\n',fval)

[k,resnorm,residual,exitflag,output,lambda,jacobian] =lsqnonlin(@ObjFunc,k0,lb,ub,[],tspan,x0,yexp);
ci = nlparci(k,residual,jacobian);
fprintf('\n\n使用函数lsqnonlin()估计得到的参数值为:\n')
fprintf('\tk1 = %.16f \n',k(1))
fprintf('\tk2 = %.16f \n',k(2))
fprintf('\tk3 = %.16f \n',k(3))
fprintf('\tk4 = %.16f \n',k(4))
fprintf('\tk5 = %.16f \n',k(5))
fprintf('\tk6 = %16f \n',k(6))
fprintf('The sum of squares is: %.16f \n\n',resnorm)

%%  plot figure
tspan = [0 38.46 50 62.5 100 166.666 250 500];
for i=1:8
[t x] = ode45(@KineticsEqs,tspan,x0,[],k);
fy(i,

=x(end,

;
end
Ypre=fy(

;
Yexp=yexp(

;
plot(Ypre,'or',Yexp,'b.')
legend('prediction','experiment')
%----------fmincon fun-------------------
function f = ObjFuncFmincon(k,x0,yexp)
tspan = [0 38.46 50 62.5 100 166.666 250 500];
for i=1:8
[t x] = ode45(@KineticsEqs,tspan,x0,[],k);
fy(i)=sum((yexp(i,

-x(end,

).^2);
end
f = sum(fy) ;
% -------------------lsqnonlin fun-------------------
function f = ObjFunc(k,x0,yexp)
tspan = [0 38.46 50 62.5 100 166.666];
for i=1:8
[t x] = ode45(@KineticsEqs,tspan,x0,[],k);
fy(i,

=yexp(i,

-x(end,

;
end
f =fy(

;
%-----------odefun--------------------------------

function dydt = KineticsEqs(t,y,k)
k1=k(1);k2=k(2);k3=k(3);k4=k(4);k5=k(5);k6=k(6);
DEC=y(1);PA=y(2);DPC=y(3);EA=y(4);EPC=y(5);
R1=k(1)*DEC*PA-(k(1)/k(2))*EPC*EA;
R2=k(3)*EPC*PA-(k(3)/k(4))*DPC*EA;
R3=k(5)*EPC.^2-(k(5)/k(6))*DPC*DEC;
dy1dt=-R1+R3;
dy2dt=-R1-R2;
dy3dt=R2+R3;
dy4dt=R1+R2;
dy5dt=R1-R2-R3;
dydt=[dy1dt;dy2dt;dy3dt;dy4dt;dy5dt];

赞一下

回复此楼

12楼2014-04-24 15:11:47

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

WULONGLOVE

铁虫 (初入文坛)

应助: 0 (幼儿园)
金币: 111
帖子: 3
在线: 10.9小时
虫号: 2352853
注册: 2013-03-16
专业: 建筑学

引用回帖:

2楼: Originally posted by 月只蓝 at 2014-03-26 20:00:43
比较粗糙的办法是求数值导数，根据t CA的数据即可求出，如果需要MATLAB代码参见：http://muchong.com/bbs/viewthread.php?tid=7020495；
精确点的方法是先拟合出动力学方程中的各参数，以第一组数据（即T=325.15K这组 ...

大神， lsqnonlin(@ObjFunc,k0,lb,ub,[],tspan,y0,yexp); 中的y0代表啥意思啊

赞一下

回复此楼

13楼2014-12-03 17:15:12

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