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


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