CODE:

function KineticsEst5 % 動(dòng)力學(xué)ODE方程模型的參數(shù)估計(jì)

%

%

%

% The variables y here are y(1)=xB, y(2)=xoNB, y(3)=xmNB,y(4)=xpNB,y(5)=xDNB .

clear all

clc

k0 = [5 5 5 5 5]; % 參數(shù)初值

lb = [0 0 0 0 0]; % 參數(shù)下限

ub = [+inf +inf +inf +inf +inf]; % 參數(shù)上限

x0 = [0 0 0 0 0 0];

KineticsData;

yexp = ExpData(:,2:6); % yexp: 實(shí)驗(yàn)數(shù)據(jù)[x1 x4 x5 x6]

% 使用函數(shù)fmincon()進(jìn)行參數(shù)估計(jì)

[k,fval,flag] = fmincon(@ObjFunc4Fmincon,k0,[],[],[],[],lb,ub,[],[],x0,yexp);

fprintf('\n使用函數(shù)fmincon()估計(jì)得到的參數(shù)值為:\n')

fprintf('\tk1 = %.4f\n',k(1))

fprintf('\tk2 = %.4f\n',k(2))

fprintf('\tk3 = %.4f\n',k(3))

fprintf('\tk4 = %.4f\n',k(4))

fprintf('\tk5 = %.4f\n',k(5))

fprintf(' The sum of the squares is: %.1e\n\n',fval)

k_fmincon = k;

% 使用函數(shù)lsqnonlin()進(jìn)行參數(shù)估計(jì)

[k,resnorm,residual,exitflag,output,lambda,jacobian] = ...

lsqnonlin(@ObjFunc4LNL,k0,lb,ub,[],x0,yexp);

ci = nlparci(k,residual,jacobian);

fprintf('\n\n使用函數(shù)lsqnonlin()估計(jì)得到的參數(shù)值為:\n')

Output

% 以函數(shù)fmincon()估計(jì)得到的結(jié)果為初值，使用函數(shù)lsqnonlin()進(jìn)行參數(shù)估計(jì)

k0 = k_fmincon;

[k,resnorm,residual,exitflag,output,lambda,jacobian] = ...

lsqnonlin(@ObjFunc4LNL,k0,lb,ub,[],x0,yexp);

ci = nlparci(k,residual,jacobian);

fprintf('\n\n以fmincon()的結(jié)果為初值，使用函數(shù)lsqnonlin()估計(jì)得到的參數(shù)值為:\n')

Output

% ------------------------------------------------------------------

function f = ObjFunc4Fmincon(k,x0,yexp)

tspan = [0.00 : 5.03 : 15.09];

[t x] = ode45(@KineticEqs,tspan,x0,[],k);

y(:,1) = x(:,1);

y(:,2:4) = x(:,4:6);

f = sum((y(:,1)-yexp(:,1)).^2) + sum((y(:,2)-yexp(:,2)).^2) ...

+ sum((y(:,3)-yexp(:,3)).^2) + sum((y(:,4)-yexp(:,4)).^2);

% ------------------------------------------------------------------

function f = ObjFunc4LNL(k,x0,yexp)

tspan = [0.00 : 5.03 : 15.09];

[t x] = ode45(@KineticEqs,tspan,x0,[],k);

y(:,1) = x(:,1);

y(:,2:4) = x(:,4:6);

f1 = y(:,1) - yexp(:,1);

f2 = y(:,2) - yexp(:,2);

f3 = y(:,3) - yexp(:,3);

f4 = y(:,4) - yexp(:,4);

f = [f1; f2; f3; f4];

% ------------------------------------------------------------------

function dxdt = KineticEqs(t,x,k)

M=1.003;

dxdt = ...

[ ( (k(1)+k(2)+k(3))*(1-x(1))*M*(1-x(2)-x(3)-x(4)-2*x(5)) )

( M*(1-x(2)-x(3)-x(4)-2*x(5))*(k(1)*(1-x(1))-k(4)*x(2)) )

( k(2)*(1-x(1))*M*(1-x(2)-x(3)-x(4)-2*x(5)) )

( M*(1-x(2)-x(3)-x(4)-2*x(5))*(k(3)*(1-x(1))-k(5)*x(4)) )

( M*(1-x(2)-x(3)-x(4)-2*x(5))*(k(4)*x(2)+k(5)*x(4)) )

];

創(chuàng)作挑戰(zhàn)賽新人創(chuàng)作獎(jiǎng)勵(lì)來(lái)咯，堅(jiān)持創(chuàng)作打卡瓜分現(xiàn)金大獎(jiǎng)

總結(jié)

以上是生活随笔為你收集整理的matlab数值很小出错,求大神帮忙解决一下，用MATLAB求解动力学数据总是出错~ - 计算模拟 - 小木虫 - 学术 科研 互动社区...的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問(wèn)題。

如果覺(jué)得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。