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

24小时热门版块排行榜    

CyRhmU.jpeg
小木虫论坛-学术科研互动平台 » 计算模拟区 » 计算模拟 » 求助用matlab三次样条法算一个反应的反应速率
141/2返回列表12下一页
查看: 1892  |  回复: 13
只看楼主@他人 存档 新回复提醒 (忽略) 收藏    在APP中查看

liqianmelody

至尊木虫 (著名写手)

[求助] 求助用matlab三次样条法算一个反应的反应速率已有1人参与

本人不会用matlab这个软件，我有实验数据，请各位精通的人帮我算一下表格标黄的那列反应速率，谢谢了！数据在附件里了，PS：附件中上传时间较晚的那个是正确文件
回复此楼

» 本帖附件资源列表

  • 欢迎监督和反馈：小木虫仅提供交流平台，不对该内容负责。
    本内容由用户自主发布，如果其内容涉及到知识产权问题，其责任在于用户本人，如对版权有异议，请联系邮箱：xiaomuchong@tal.com
  • 附件 1 : 动力学实验数据.doc
    • 2014-03-26 14:46:39, 79.5 K
  • 附件 2 : 动力学实验数据.doc
    • 2014-03-26 14:48:30, 79.5 K

» 猜你喜欢

» 本主题相关价值贴推荐，对您同样有帮助:
1楼2014-03-26 14:46:44
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

月只蓝

主管区长 (职业作家)

【答案】应助回帖

★ ★ ★ ★ ★ ★ ★ ★ ★ ★
感谢参与，应助指数 +1
fegg7502: 金币+10, 3ks 2014-03-29 08:41:56
比较粗糙的办法是求数值导数，根据t CA的数据即可求出，如果需要MATLAB代码参见：http://muchong.com/bbs/viewthread.php?tid=7020495；
精确点的方法是先拟合出动力学方程中的各参数，以第一组数据（即T=325.15K这组）为例，楼主没给出反应方程式，根据三个浓度随时间的增幅，推测反应式为：
                                 A   +     B      =       C
C为生成物，初始浓度为0。忽略反应前后体系体积变化，则CA0=0.1530+0.0206=0.1736；CB0=3.4505+0.0206=3.4711；设物质A的瞬时浓度为CA，有：
                                 A        +          B            =             C
初始浓度                     CA0                CB0                        0
瞬时浓度                     CA                [CB0-(CA0-CA)]         CA0-CA

-(dCA/dt)=kCAmCBn=k0*exp(-Ea/RT)* CA^m*CB^n
得dCA/dt=-k0*exp(-Ea/RT)* CA^m* [CB0-(CA0-CA)]^n      （1）
该微分方程无法积分求得关于CA的代数式解析解，故直接对式（1）作常微分方程的拟合，求出指前因子k0、活化能Ea、两个反应级数m和n，这样把不同时间t对应的CA代入式（1）即得dCA/dt，取负得-dCA/dt。
MATLAB代码如下：
CODE:
function ode3333
clear all;clc

format long

global T
T=323.15;
tspan=[0
    30
60
90
120
180
240
300
360
420
480
540
600
];

t1=length(tspan);
yexp=[0.1530
0.1222
0.1036
0.0874
0.0619
0.0331
0.0305
0.0241
0.0185
0.0120
0.0080
0.0050
]';

t2=length(yexp);



y0=0.1736;

k0=[0.1 10000 2 2];
lb=-[1 1 1 1]*1e9;
ub=[1 1 1 1]*1e9;        

yy=[y0 yexp];

[k,resnorm,residual,exitflag,output,lambda,jacobian] = ...
    lsqnonlin(@ObjFunc,k0,lb,ub,[],tspan,y0,yexp);
ci = nlparci(k,residual,jacobian);
fprintf('\n\n使用函数lsqnonlin()估计得到的参数值为:\n')
fprintf('\t k0 = %.4f\n',k(1))
fprintf('\t Ea = %.4f\n',k(2))
fprintf('\t m = %.4f\n',k(3))
fprintf('\t n = %.4f\n',k(4))

%fprintf('\t反应级数 n = %.4f\n',k(2))
fprintf('  The sum of the squares is: %.1e\n\n',resnorm)
ts=0:1:max(tspan);

[ts ys]=ode45(@KineticsEqs,ts,y0,[],k);
[ttt XXsim] = ode45(@KineticsEqs,tspan,y0,[],k);
y=XXsim(2:end);
xexp=yexp;
R2=1-sum((xexp'-y).^2)./sum((xexp'-mean(y)).^2);
fprintf('\n\t相关系数之平方R^2 = %.6f',R2);
k0=k(1);
Ea=k(2);
m=k(3);
n=k(4);
y=yexp';
dydt_a=-k0*exp(-Ea/8.314/T)*y.^m.*(3.2975+y).^n;
[-dydt_a]


figure(1)
plot(ts,ys,'b',tspan,yy,'or'),legend('计算值','实验值','Location','best'),


%---------------------------------------------------------
function f = ObjFunc(k,tspan,y0,yexp)           % 目标函数
[t Xsim] = ode45(@KineticsEqs,tspan,y0,[],k) ;
Xsim;
ysim = Xsim(2:end);
f=ysim-yexp';
%----------------------------------------------------------

function dydt = KineticsEqs(t,y,k)
global T
k0=k(1);
Ea=k(2);
m=k(3);
n=k(4);
dydt=-k0*exp(-Ea/8.314/T)*y.^m.*(3.2975+y).^n;

计算结果显示：
k0 = 923.7021
         Ea = 9079.9827
         m = 1.1387
         n = -6.7529
相关系数之平方R^2 = 0.994517

-dCA/dt：

1.0e-003 *

   0.865245002514304
   0.711666673454071
   0.611807848360504
   0.520628349512815
   0.369919263546430
   0.192217346733282
   0.176046594470231
   0.136397181433893
   0.102089651326710
   0.063193844805022
   0.040152229316972
   0.023656194149690
拟合结果见附图1。
其他T下面的数据，修改代码中T、tspan、yexp对应的数据即可。

如果楼主嫌麻烦，直接用前面说的求数值导数的方法也可以。
求助用matlab三次样条法算一个反应的反应速率
附图1.jpg
赞一下 回复此楼
MATLAB、MS小问题、普通问题请发帖求助！时间精力有限，恕不接受无偿私信求助。
2楼2014-03-26 20:00:43
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

月只蓝

主管区长 (职业作家)

【答案】应助回帖

★ ★
fegg7502: 金币+2, 应助指数+1, 3ks 2014-03-29 08:42:05
注意到你要求用三次样条法，三次B样条法求数值导数的例子参见下列链接的第12楼
http://muchong.com/bbs/viewthread.php?tid=6744599&fpage=1&target=self&page=2
涉及代码为”
CODE:
% 动力学数据
t = [0  20  40  60  120  180  300];
CAm = [10  8  6  5  3  2  1];

% 用最小二乘样条拟合法计算微分dCA/dt--使用不经过实验点的B样条插值函数
knots = 3;
K = 3;                  % 三次B样条
sp = spap2(knots,K,t,CAm);
pp = fnder(sp);         % 计算B样条函数的导函数
dCAdt = fnval(pp,t)    % 计算t处的导函数值
rAm = dCAdt;

“
该代码转自原链接第12楼。
赞一下 回复此楼
MATLAB、MS小问题、普通问题请发帖求助！时间精力有限，恕不接受无偿私信求助。
3楼2014-03-26 20:19:21
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

liqianmelody

至尊木虫 (著名写手)
★
fegg7502: 金币+1, 鼓励交流 2014-03-29 08:42:18
引用回帖:
3楼: Originally posted by 月只蓝 at 2014-03-26 20:19:21
注意到你要求用三次样条法，三次B样条法求数值导数的例子参见下列链接的第12楼
http://muchong.com/bbs/viewthread.php?tid=6744599&fpage=1&target=self&page=2
涉及代码为”
% 动力学数据
t = ;
C ...

谢谢你帮我解释那么多，您楼上做出了动力学方程的结果和我之前找人做的不太一样，您做出的结果不太跟我的实际一样，之前是通过1stopt做出的，没有做出反应速率，所以这次就想做出这个，我是做分离的，对matlab一点也不懂，所以插入什么代码的也不会，就想请人帮我做出结果，我看您对此比较精通，如果您可以帮我做出，感激涕零，顺便可以追加金币，谢谢了！！！
赞一下 回复此楼
4楼2014-03-26 21:21:56
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

月只蓝

主管区长 (职业作家)

【答案】应助回帖

★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ...
liqianmelody: 金币+50, ★★★★★最佳答案 2014-03-27 11:15:20
fegg7502: 金币+2, 应助指数+1, 3ks 2014-03-29 08:42:30
引用回帖:
4楼: Originally posted by liqianmelody at 2014-03-26 21:21:56
谢谢你帮我解释那么多，您楼上做出了动力学方程的结果和我之前找人做的不太一样，您做出的结果不太跟我的实际一样，之前是通过1stopt做出的，没有做出反应速率，所以这次就想做出这个，我是做分离的，对matlab一点 ...

你说的这个帖子么：
http://muchong.com/bbs/viewthread.php?tid=6727491
实习版主dingd已经把各个参数求出来了：
k0        870504.566962617
ea        53899.8216594799
m        1.02066293612806
n        1.08210287262186
再把CA代回-k0*exp(-Ea/(R*T1))*CA^m*CB^n不就得到想要的-dCA/dt了么？
赞一下 回复此楼
MATLAB、MS小问题、普通问题请发帖求助！时间精力有限，恕不接受无偿私信求助。
5楼2014-03-26 21:29:35
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

justin534

新虫 (初入文坛)
★
fegg7502: 金币+1, 鼓励交流 2014-03-29 08:42:40
引用回帖:
5楼: Originally posted by 月只蓝 at 2014-03-26 21:29:35
你说的这个帖子么：
http://muchong.com/bbs/viewthread.php?tid=6727491
实习版主dingd已经把各个参数求出来了：
k0        870504.566962617
ea        53899.8216594799
m        1.02066293612806
n       ...

不好意思版大,我修改了一下你的編碼,但是差不多只改數據而已,但不知會和不能run,可否幫我看一下,謝謝!


function ode3333
clear all;clc

format long

tspan=[50
62.5
80
100
125
166.666
250
333.3333
500
];

t1=length(tspan);
yexp=[1.414
1.1369
0.91816
0.81575
0.7424811
0.6323514
0.6139
0.6109
0.60789
]';

t2=length(yexp);



y0=1.9880715;

k0=[0.1 1];
lb=[0 0 ];
ub=[1 1 ]*1e9;        

yy=[y0 yexp];

[k,resnorm,residual,exitflag,output,lambda,jacobian] = ...
    lsqnonlin(@ObjFunc,k0,lb,ub,[],tspan,y0,yexp);
ci = nlparci(k,residual,jacobian);
fprintf('\n\n使用函?lsqnonlin()估?得到的??值?:\n')
fprintf('\t kf1 = %.4f\n',k(1))
fprintf('\t kf2 = %.4f\n',k(2))
fprintf('  The sum of the squares is: %.1e\n\n',resnorm)
ts=0:1:max(tspan);

[ts ys]=ode45(@KineticsEqs,ts,y0,[],k);
[ttt XXsim] = ode45(@KineticsEqs,tspan,y0,[],k);
y=XXsim(2:end);
xexp=yexp;
R2=1-sum((xexp'-y).^2)./sum((xexp'-mean(y)).^2);
fprintf('\n\t相?系?之平方R^2 = %.6f',R2);
kf1=k(1);
kf2=k(2);
y=yexp';
dydt_a=-(k(1)*y*((y+3.97614).^2)-(k(1)/k(2)).*(y-1.9880715).*((y-1.9880715).^2));
[-dydt_a]


figure(1)
plot(ts,ys,'b',tspan,yy,'or'),legend('?算值','??值','Location','best'),


%---------------------------------------------------------
function f = ObjFunc(k,tspan,y0,yexp)           % 目?函?
[t Xsim] = ode45(@KineticsEqs,tspan,y0,[],k) ;
Xsim;
ysim = Xsim(2:end);
f=ysim-yexp';
%----------------------------------------------------------

function dydt = KineticsEqs(t,y,k)
kf1=k(1);
kf2=k(2);
dydt=-(kf1*y*((y+3.97614).^2)-(kf1/kf2)*(y-1.9880715)*((y-1.9880715).^2));
赞一下 回复此楼
6楼2014-03-28 12:34:46
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

月只蓝

主管区长 (职业作家)
★
fegg7502: 金币+1, 3ks 2014-03-29 08:42:52
引用回帖:
6楼: Originally posted by justin534 at 2014-03-28 12:34:46
不好意思版大,我修改了一下你的編碼,但是差不多只改數據而已,但不知會和不能run,可否幫我看一下,謝謝!


function ode3333
clear all;clc

format long

tspan=;

t1=length(tspan);
yexp=';

t2=len ...

t=0时，y=1.9880715么？
以下代码按上述设定编写：
CODE:
function ode3333
clear all;clc

format long

tspan=[0
    50
62.5
80
100
125
166.666
250
333.3333
500
];

t1=length(tspan);
yexp=[1.414
1.1369
0.91816
0.81575
0.7424811
0.6323514
0.6139
0.6109
0.60789
]';

t2=length(yexp);



y0=1.9880715;

k0=[1e-6 -0.1];
lb=[1 1]*-1e9;
ub=[1 0];        

yy=[y0 yexp];

[k,resnorm,residual,exitflag,output,lambda,jacobian] = ...
    lsqnonlin(@ObjFunc,k0,lb,ub,[],tspan,y0,yexp);
ci = nlparci(k,residual,jacobian);
fprintf('\n\n使用函数lsqnonlin()估计得到:\n')
fprintf('\t kf1 = %.9f\n',k(1))
fprintf('\t kf2 = %.9f\n',k(2))
fprintf('  The sum of the squares is: %.1e\n\n',resnorm)
ts=0:1:max(tspan);

[ts ys]=ode45(@KineticsEqs,ts,y0,[],k);
[ttt XXsim] = ode45(@KineticsEqs,tspan,y0,[],k);
y=XXsim(2:end);
xexp=yexp;
R2=1-sum((xexp'-y).^2)./sum((xexp'-mean(y)).^2);
fprintf('\n\t相关系数之平方R^2 = %.6f',R2);
kf1=k(1);
kf2=k(2);
y=yexp';
%dydt_a=-(k(1)*y*((y+3.97614).^2)-(k(1)/k(2)).*(y-1.9880715).*((y-1.9880715).^2));
%[-dydt_a]


figure(1)
plot(ts,ys,'b',tspan,yy,'or'),legend('计算值','实验值','Location','best'),


%---------------------------------------------------------
function f = ObjFunc(k,tspan,y0,yexp)           % 目?函?
[t Xsim] = ode45(@KineticsEqs,tspan,y0,[],k) ;
Xsim;
ysim = Xsim(2:end);
f=ysim-yexp';
%----------------------------------------------------------

function dydt = KineticsEqs(t,y,k)
kf1=k(1);
kf2=k(2);
dydt=-(kf1*y*((y+3.97614).^2)-(kf1/kf2)*(y-1.9880715)*((y-1.9880715).^2));

计算结果：
使用函数lsqnonlin()估计得到:
         kf1 = 0.000315091
         kf2 = -0.232971187
  The sum of the squares is: 4.9e-002


        相关系数之平方R^2 = 0.923278>>
求助用matlab三次样条法算一个反应的反应速率-1
附图1.png
赞一下 回复此楼
MATLAB、MS小问题、普通问题请发帖求助！时间精力有限，恕不接受无偿私信求助。
7楼2014-03-28 16:31:03
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

pdl9527

专家顾问 (小有名气)
★
fegg7502: 金币+1, 鼓励交流 2014-03-29 08:42:59
如果其他温度求出来的参数差别很大怎么办，能不能求得不同温度的参数然后求平均值做呢？
赞一下 回复此楼
8楼2014-03-28 21:42:17
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

justin534

新虫 (初入文坛)
引用回帖:
7楼: Originally posted by 月只蓝 at 2014-03-28 16:31:03
t=0时，y=1.9880715么？
以下代码按上述设定编写：


function ode3333
clear all;clc

format long

tspan=;

t1=length(tspan);
yexp=';

t2=length(yexp);



y0=1.9880715;

k0=;
lb=*-1 ...

謝謝版大,可是kf2值物理意義上是正值,是不是有什麼方法將kf2趨於正值,擬合結果又不錯的,又或者�?驗數據本身就有問題,因為我將lb ub的區間設為正值的時候,擬合結果就很差!
赞一下 回复此楼
9楼2014-03-29 12:28:35
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

月只蓝

主管区长 (职业作家)
引用回帖:
9楼: Originally posted by justin534 at 2014-03-29 12:28:35
謝謝版大,可是kf2值物理意義上是正值,是不是有什麼方法將kf2趨於正值,擬合結果又不錯的,又或者�?驗數據本身就有問題,因為我將lb ub的區間設為正值的時候,擬合結果就很差!...

可能是方程的问题，也可能是实验数据的问题。
你可以请实习版主dingd 用1stopt软件做一下，来验证一下。
赞一下 回复此楼
MATLAB、MS小问题、普通问题请发帖求助！时间精力有限，恕不接受无偿私信求助。
10楼2014-03-29 15:11:44
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 liqianmelody 的主题更新
141/2返回列表12下一页
普通表情 龙 兔 虎 猫
信息提示
关闭
请填处理意见
关闭