12楼: Originally posted by dingd at 2012-04-23 20:56:39:
这是公式：
43/99/291104_1335185785_534.jpg

12楼: Originally posted by dingd at 2012-04-23 20:56:39:
这是公式：
43/99/291104_1335185785_534.jpg

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);