24小时热门版块排行榜    

查看: 2895  |  回复: 15

abandonkkk

铁虫 (小有名气)

[求助] 请教matlab反应动力学参数估计遇到的问题,谢谢

比如做了六组实验,实验的输入变量是温度随时间的变化曲线T(t),输出是转化率X。X同T(t)的关系为
-ln(1-x)=从0积分到时间t [k0exp(Ea/RT(t)]dt;其中k0和Ea为需估计的动力学参数
请问如何解决这个问题。我看书上基本用lsqnonlin工具,但是我目前的问题是:我的六个自变量是T(t)的离散数据组,当然也可以拟合成函数,这同一般例子中的自变量是一组单纯的数不同。不知道怎么处理,谢谢!
回复此楼
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
回帖支持 ( 显示支持度最高的前 50 名 )

change0618

铁杆木虫 (著名写手)

方丈大师

【答案】应助回帖

感谢参与,应助指数 +1
CODE:
function exercise011
clear all
clc
global R
R = 8.413;
a1=[0        96.1
1        97.3
2        158.6
3        190
4        197.5
5        194.8
6        194.5
7        195.4
8        196.3
9        195.4
10        195.6
11        195.7
12        195.2
13        194.6
14        192
15        183
16        174.2
17        166.5
18        159.7
19        153.5
20        147.9];
a2=[0        96.6
1        97.8
2        98.8
3        97.9
4        97.5
5        97.9
6        171.2
7        197.9
8        196.7
9        194
10        194.7
11        195.2
12        195.9
13        195.3
14        194.9
15        195.1
16        195.5
17        192.6
18        193
19        194.1
20        195.5
21        193.6
22        192.9
23        192.2
24        193.5
25        194.8
26        192.9
27        190.5
28        183.4
29        176.1
30        169.5
31        163.6
32        158.4
33        153.6
34        149.1
35        145
36        141.2
37        137.8
38        134.5
39        131.3
40        128.4
41        125.7
42        122.9
43        120.4
44        117.9
45        115.6
46        113.4
47        111.3
48        109.3
49        107.2
50        105.3
51        103.6
52        101.8
53        100.1
54        98.5];
a3=[0        95.9
1        97.3
2        165
3        190.9
4        195.9
5        195.5
6        195.5
7        196.5
8        195.9
9        195.6
10        196.1
11        196.4
12        194.5
13        193
14        192.3
15        191.8
16        191.3
17        191
18        190.7
19        190.5
20        190.3
21        190.1
22        190.2
23        190.5
24        190.9
25        191.5
26        191
27        190.7
28        191.3
29        191.7
30        192.1
31        192.5
32        192.8
33        188.2
34        179.8
35        171.5
36        164.3
37        157.8
38        152
39        146.8
40        141.9
41        137.5
42        133.5
43        129.7
44        126.1
45        122.9
46        119.8
47        116.9
48        114.2
49        111.6
50        109.1
51        106.9
52        104.7
53        102.8
54        101.6
55        100.6];
a4=[0        96.3
1        124.7
2        184.9
3        197.3
4        194.5
5        193.7
6        194.5
7        195
8        195.4
9        194.9
10        194
11        193.3
12        192.7
13        192.2
14        191.9
15        191.6
16        191.4
17        191.1
18        191.4
19        191.5
20        191.6
21        191.8
22        191.9
23        191.9
24        191.9
25        191.9
26        191.9
27        191.9
28        191.9
29        191.9
30        191.7
31        191.7
32        191.5
33        191.6
34        191.4
35        191.3
36        191.1
37        191.1
38        191
39        190.9
40        190.7
41        190.7
42        190.4
43        183.6
44        174.9
45        166.8
46        159.3
47        152.3
48        145.6
49        139.7
50        134.3
51        129.2
52        124.6
53        120.3
54        116.3
55        112.6
56        109.2
57        105.9
58        102.8];
a5=[0        97
1        99.3
2        180.6
3        196.1
4        194.8
5        193.2
6        193.9
7        194.5
8        195
9        195.6
10        195.8
11        195.4
12        195.7
13        196.3
14        194.8
15        192.8
16        193.9
17        194.9
18        195.4
19        194
20        194.7
21        194.9
22        192.2
23        192.7
24        193.3
25        194.1
26        194.9
27        195.8
28        192.7
29        192.8
30        193.7
31        194.7
32        195.5
33        195.7
34        194.3
35        194.8
36        195.4
37        195.8
38        193.5
39        193.8
40        194.3
41        194.7
42        195.1
43        195.5
44        195.6
45        193.4
46        193.5
47        193.7
48        194.1
49        194.3
50        194.7
51        194.8
52        195
53        195.2
54        195.5
55        195.5
56        195.6
57        195.8
58        195.8
59        195.9
60        189.6
61        180.3
62        172
63        164.5
64        158
65        152
66        146.6
67        141.5
68        136.8
69        132.6
70        128.5
71        124.9
72        121.4
73        118.2
74        115.4
75        112.7
76        110.2
77        107.8
78        105.4
79        103.3
80        101.2
81        99.2];
x0= [0.232573702        0.443523301        0.522192614        0.724670708        0.942120657];
beta0=[0.5 2810];
lb = [-inf -inf];
ub = [+inf +inf];
[beta,resnorm,resid,exitflag]=lsqnonlin(@ObjFun,beta0,lb,ub,[],a1,a2,a3,a4,a5,x0)

function f = ObjFun(beta,a1,a2,a3,a4,a5,x0)
[x1,B1]=ode23s(@odefun,[a1(1,1),a1(end,1)],0,[],beta,a1);
[x2,B2]=ode23s(@odefun,[a2(1,1),a2(end,1)],0,[],beta,a2);
[x3,B3]=ode23s(@odefun,[a3(1,1),a3(end,1)],0,[],beta,a3);
[x4,B4]=ode23s(@odefun,[a4(1,1),a4(end,1)],0,[],beta,a4);
[x5,B5]=ode23s(@odefun,[a5(1,1),a5(end,1)],0,[],beta,a5);
B=[B1(end),B2(end),B3(end),B4(end),B5(end)];
xcal=1-exp(-B);
f=x0-xcal;

function f = odefun(t,x,beta,a)
global R
k0=beta(1);
Ea=beta(2);
T=Tfun(x,a);
f = k0 * exp(-Ea/R/T);

function f= Tfun(x,a)
% f=interp1(a1(:,1),a1(:,2),t,'spline');
pp=fit(a(:,1),a(:,2),'smoothingspline');
f=feval(pp,x);

15楼2012-04-24 11:02:18
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
普通回帖

dingd

铁杆木虫 (职业作家)

【答案】应助回帖

★ ★
感谢参与,应助指数 +1
zhangguangping: 金币+2, 谢谢应助! 2012-04-23 19:53:11
1stOpt应该可以解决你的问题。有数据的话放上来可以试试。
2楼2012-04-23 14:26:26
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

abandonkkk

铁虫 (小有名气)

引用回帖:
2楼: Originally posted by dingd at 2012-04-23 14:26:26:
1stOpt应该可以解决你的问题。有数据的话放上来可以试试。

你好,非常感谢,数据怎么传给你
3楼2012-04-23 15:13:46
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

abandonkkk

铁虫 (小有名气)

有5组数据
4楼2012-04-23 15:16:37
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

abandonkkk

铁虫 (小有名气)

» 本帖附件资源列表

  • 欢迎监督和反馈:小木虫仅提供交流平台,不对该内容负责。
    本内容由用户自主发布,如果其内容涉及到知识产权问题,其责任在于用户本人,如对版权有异议,请联系邮箱:xiaomuchong@tal.com
  • 附件 1 : 参数拟合.xls
  • 2012-04-23 15:19:50, 61 K
5楼2012-04-23 15:20:03
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

dingd

铁杆木虫 (职业作家)

【答案】应助回帖

就发论坛不行吗?能帮的人更多些。
6楼2012-04-23 15:20:08
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

abandonkkk

铁虫 (小有名气)

引用回帖:
6楼: Originally posted by dingd at 2012-04-23 15:20:08:
就发论坛不行吗?能帮的人更多些。

我刚才不知道怎么发,现在可以了
7楼2012-04-23 15:22:13
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

dingd

铁杆木虫 (职业作家)

【答案】应助回帖

1: 5组数据分别单独拟合得出5组k0和Ea,还是在考虑转化率的前提下一起拟合得出一组k0和Ea?转化率在公式中如何体现?
2: 公式中R为多少?最好将完整的公式图片贴一个。
8楼2012-04-23 15:38:08
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

abandonkkk

铁虫 (小有名气)

引用回帖:
8楼: Originally posted by dingd at 2012-04-23 15:38:08:
1: 5组数据分别单独拟合得出5组k0和Ea,还是在考虑转化率的前提下一起拟合得出一组k0和Ea?转化率在公式中如何体现?
2: 公式中R为多少?最好将完整的公式图片贴一个。

是用五组数据一起拟合成一组k0和Ea,
公式中的X是转化率
9楼2012-04-23 18:41:53
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

abandonkkk

铁虫 (小有名气)

zhangguangping: 编辑内容 2012-04-24 09:24
zhangguangping: , 最好引用回复或者@ID,这样对方能够及时看到。 2012-04-24 09:24:34
公式

公式

@dingd

[ Last edited by zhangguangping on 2012-4-24 at 09:24 ]
10楼2012-04-23 18:42:36
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 abandonkkk 的主题更新
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[公派出国] 售SCI一区文章,我:8.O.551.O.5.4,科目全,可伽急 +3 PbwVEMK5haaE 2026-08-25 4/200 2026-08-26 03:20 by cNXvBfCpiZOM
[考博] 售SCI一区T0P文章,我:8.O.55.1.O.5.4,科目全,可+急 +3 Qj98d49dHFb9 2026-08-25 4/200 2026-08-26 03:08 by cNXvBfCpiZOM
[考研] 售SCI一区T0P文章,我:8.O.55.1.O.5.4,科目全,可+急 +3 7K1CJE38xLG4 2026-08-25 4/200 2026-08-26 02:04 by cNXvBfCpiZOM
[论文投稿] 售SCI一区T0P文章,我:8O.55.1.O.5.4,科目齐全,可+急 +3 7K1CJE38xLG4 2026-08-25 4/200 2026-08-26 02:04 by cNXvBfCpiZOM
[基金申请] 今天务委会开完了,明天出结果吗 +18 angus9576 2026-08-25 22/1100 2026-08-26 00:41 by merchancy
[基金申请] 在坚冰还盖着北海的时候,我看到了怒放的梅花。 (金币+10) +6 ziyangfang 2026-08-25 9/450 2026-08-25 20:26 by huagongfeihu
[基金申请] 明天应该可查了!? +6 chengyan1220 2026-08-23 6/300 2026-08-25 19:45 by zfd97
[基金申请] 有没有大神帮我看看基金代码 28+4 1234567wang 2026-08-24 10/500 2026-08-25 19:15 by lfy8008
[基金申请] 某些机构,以效率低为荣,以效率低作为存在感 +9 yuleib84 2026-08-25 10/500 2026-08-25 17:14 by alexon
[基金申请] 今日不放榜?网传国自然预计 8 月 27 日可查结果 +17 医学老男孩 2026-08-20 22/1100 2026-08-25 15:36 by 医学老男孩
[基金申请] 人气不行了 +11 fansofjerry 2026-08-21 11/550 2026-08-25 11:04 by 孤独的英雄6
[教师之家] 导师吐槽:我怎么摊上了这么个极品研究生! +3 苏东坡二世 2026-08-23 3/150 2026-08-25 10:35 by shisan1313
[基金申请] 2026年的国家社科基金项目通讯评审的新规则与新动向、新挑战 +5 process2012 2026-08-23 7/350 2026-08-25 09:42 by huixian257
[基金申请] 范进中举一文的中心思想 +7 炎黄贵胄 2026-08-22 8/400 2026-08-25 08:48 by ZJTJZ
[基金申请] 只有每年这种时候来逛逛小木虫 +25 yaoyewhu2008 2026-08-20 27/1350 2026-08-25 08:01 by Equinoxhua
[基金申请] 能否退出参与的面上项目解除限项 +21 koalala 2026-08-24 24/1200 2026-08-24 19:25 by 家与远方
[教师之家] 跳槽后在研项目怎么办? +5 简单化xn 2026-08-22 10/500 2026-08-23 12:38 by 简单化xn
[基金申请] 今天放榜吗? +15 布布和一二 2026-08-19 16/800 2026-08-23 09:55 by 张春生
[基金申请] 看来今天不会放榜了? +8 chengyan1220 2026-08-21 11/550 2026-08-21 17:52 by dcqxinyang
[基金申请] 时间戳又变了 +13 wuchongjun 2026-08-20 19/950 2026-08-21 17:21 by 紫杉醇
信息提示
请填处理意见