24小时热门版块排行榜

返回列表

2011207156

金虫 (小有名气)

应助: 0 (幼儿园)
金币: 1335.7
红花: 2
帖子: 194
在线: 153.8小时
虫号: 1829583
注册: 2012-05-22
专业: 化学反应工程

[求助] matlab拟合动力学参数，非常谢谢~

function butongshike201506171
clear all;clc
format long
tspan=(0:10:390);
yexp=[0.046125761496241
0.078139373146412
0.117158327134899
0.159846402230616
0.199412869433285
0.231111460229405
0.259814231746068
0.282697037233637
0.306493661349134
0.326396561089024
0.345295684999190
0.360444423309371
0.377418750129898
0.395250807577438
0.414965633800510
0.430017989495138
0.445923212579888
0.457995753035510
0.469186259425871
0.480284344291621
0.491524402631167
0.501658351934599
0.510935126510884
0.520202070112304
0.531309455870521
0.537782255717000
0.547944130180314
0.558094213007469
0.565468701311900
0.573757577102420
0.581007002009939
0.589392949772071
0.597541754780127
0.607751809867969
0.613129267171725
0.620580135887054
0.626072239599544
0.632475341850705
0.638742438783615];
beta0=[67971 0.000839];
y0=0;
lb=[-1 -1]*1e6;
ub=[1e9 1e9];
yy=[0 yexp'];
k0=beta0;
[k,resnorm,residual,exitflag,output,lambda,jacobian] = ...
lsqnonlin(@ObjFunc,k0,lb,ub,[],tspan,y0,yexp);
fprintf('\n\n使用函数lsqnonlin()估计得到的参数值为:\n')
fprintf('\t待拟合参数 k1 = %.6f\n',k(1))
fprintf('\t待拟合参数 k2= %.6f\n',k(2))
fprintf(' 残差平方和= %.6f\n\n',resnorm)
ts=0:10: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-Square = %.6f',R2);
figure(1)
plot(ts,ys,'b',tspan,yy,'or'),legend('计算值','实验值','Location','best');
xlabel('x'),ylabel('y')

%---------------------------------------------------------
function f = ObjFunc(k,tspan,y0,yexp) % 目标函数
[t Xsim] = ode45(@KineticsEqs,tspan,y0,[],k) ;
ysim = Xsim(2:end);
size(ysim);
size(yexp);
f=ysim-yexp;
%----------------------------------------------------------

function dydt = KineticsEqs(t,y,k)
beta(1)=k(1);
beta(2)=k(2);
tspan=(0:10:390);
h=0.068031;
p=[3.472000000000000
3.519000000000000
3.548000000000000
3.589000000000000
3.634000000000000
3.676000000000000
3.710000000000000
3.740000000000000
3.765000000000000
3.789000000000000
3.811000000000000
3.832000000000000
3.847000000000000
3.865000000000000
3.885000000000000
3.906000000000000
3.921000000000000
3.938000000000000
3.952000000000000
3.964000000000000
3.977000000000000
3.988000000000000
4.000000000000000
4.010000000000000
4.020000000000000
4.032000000000000
4.039000000000000
4.050000000000000
4.061000000000000
4.069000000000000
4.078000000000000
4.087000000000000
4.095000000000000
4.105000000000000
4.115000000000000
4.122000000000000
4.129000000000000
4.135000000000000
4.142000000000000
4.150000000000000]';
D=p*h;
sp=csaps(tspan,D,1);
pp=fnder(sp);
dD=fnval(pp,tspan);
dydt= (83.38-0.0014*exp(1600*(1/277.55-1/298.15))...
*interp1(tspan,p,t,'pchip')*10000-1/beta(2)...
*interp1(tspan,D,t,'pchip')*4.229*(0.0598...
-(y*0.0150*(54.2+17.68)+(0.236204-interp1(tspan,D,t,'pchip'))...
*0.0598*4.229)/(45.07-interp1(tspan,D,t,'pchip')*4.229))...
*interp1(tspan,dD,t,'pchip')/(-45.07+interp1(tspan,D,t,'pchip')...
*4.229)/0.0150-(0.0598-(y*0.0150*(54.2+17.68)+(0.236204-...
interp1(tspan,D,t,'pchip'))*0.0598*4.229)/...
(45.07-interp1(tspan,D,t,'pchip')*4.229))*...
interp1(tspan,dD,t,'pchip')/0.0150/beta(1))/...
(1/beta(1)+1/beta(2)-1/beta(2)*interp1(tspan,D,t,'pchip')...
*71.88/(-45.07+4.229*interp1(tspan,D,t,'pchip')));
以上是求动力学参数beta(1)和beta(2)的代码，p,D,和dD是随时间变化的已知量，不知以上程序运行时interp1(tspan,D,t,'pchip')能否对应上不同时刻的D。非常谢谢。

回复此楼