24小时热门版块排行榜    

CyRhmU.jpeg
查看: 2225  |  回复: 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的回帖

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的回帖
查看全部 16 个回答

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的回帖

dingd

铁杆木虫 (职业作家)

【答案】应助回帖

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