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;

% ¶¯Á¦Ñ§Êý¾Ý
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;

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 ...

4Â¥: Originally posted by liqianmelody at 2014-03-26 21:21:56
ллÄã°ïÎÒ½âÊÍÄÇô¶à£¬ÄúÂ¥ÉÏ×ö³öÁ˶¯Á¦Ñ§·½³ÌµÄ½á¹ûºÍÎÒ֮ǰÕÒÈË×öµÄ²»Ì«Ò»Ñù£¬Äú×ö³öµÄ½á¹û²»Ì«¸úÎÒµÄʵ¼ÊÒ»Ñù£¬Ö®Ç°ÊÇͨ¹ý1stopt×ö³öµÄ£¬Ã»ÓÐ×ö³ö·´Ó¦ËÙÂÊ£¬ËùÒÔÕâ´Î¾ÍÏë×ö³öÕâ¸ö£¬ÎÒÊÇ×ö·ÖÀëµÄ£¬¶ÔmatlabÒ»µã ...

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       ...

6Â¥: Originally posted by justin534 at 2014-03-28 12:34:46
²»ºÃÒâ˼°æ´ó,ÎÒÐÞ¸ÄÁËÒ»ÏÂÄãµÄ¾Ž´a,µ«ÊDz¶àÖ»¸Ä”µ“þ¶øÒÑ,µ«²»Öª•þºÍ²»ÄÜrun,¿É·ñŽÍÎÒ¿´Ò»ÏÂ,ÖxÖx!


function ode3333
clear all;clc

format long

tspan=;

t1=length(tspan);
yexp=';

t2=len ...

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

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 ...

9?: Originally posted by justin534 at 2014-03-29 12:28:35
?x?x???,????kf2????????x???????,?????????N??????kf2?????,?M??Y??????e??,???????????????????},??????lb ub??^?g?O???????r??,?M??Y??????!...