matlab拟合实验心得,matlab拟合效果图(实验与拟合的比较)

function Kineticsdaifei8_6_370_shiyan_1

clear all;

clc

k0=[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0];%参数初值

lb=[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0];ub=[+inf +inf +inf +inf +inf +inf +inf +inf +inf +inf +inf +inf +inf +inf +inf +inf +inf +inf +inf];%lb、ub:参数下限和上限

x0=[1 0 0 0 0 0 0 0];

data=[0.2   0.0420 0.0853 0.0762 0.0212 0.6087 0.1522 0.0131 0.0013

0.4    0.0087 0.0217 0.0195 0.0049 0.7398 0.1843 0.0178 0.0033

0.6    0.0015 0.0047 0.0042 0.0010 0.7707 0.1912 0.0210 0.0057

0.8    0.0003 0.0012 0.0011 0.0002 0.7743 0.1915 0.0235 0.0079

];

% KineticsData1;

tspan=data(:,1);

yexp=data(:,2:end);

%使用函数fmincon()进行参数估计

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

fprintf('\n使用函数fmincon()估计得到的参数值:\n')

fprintf('\tk1=%.9f\n',k(1)),fprintf('\tk2=%.9f\n',k(2)),fprintf('\tk3=%.9f\n',k(3))

fprintf('\tk4=%.9f\n',k(4)),fprintf('\tk5=%.9f\n',k(5)),fprintf('\tk6=%.9f\n',k(6))

fprintf('\tk7=%.9f\n',k(7)),fprintf('\tk8=%.9f\n',k(8)),fprintf('\tk9=%.9f\n',k(9))

fprintf('\tk10=%.9f\n',k(10)),fprintf('\tk11=%.9f\n',k(11)),fprintf('\tk12=%.9f\n',k(12))

fprintf('\tk13=%.9f\n',k(13)),fprintf('\tk14=%.9f\n',k(14)),fprintf('\tk15=%.9f\n',k(15))

fprintf('\tk16=%.9f\n',k(16)),fprintf('\tk17=%.9f\n',k(17)),fprintf('\tk18=%.9f\n',k(18))

fprintf('\tk19=%.9f\n',k(19))

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

k_fmincon=k;

%使用函数lsqnonlin()进行参数估计

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

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

ci=nlparci(k,residual,jacobian);

fprintf('\tk1=%.9f\n',k(1)),fprintf('\tk2=%.9f\n',k(2)),fprintf('\tk3=%.9f\n',k(3))

fprintf('\tk4=%.9f\n',k(4)),fprintf('\tk5=%.9f\n',k(5)),fprintf('\tk6=%.9f\n',k(6))

fprintf('\tk7=%.9f\n',k(7)),fprintf('\tk8=%.9f\n',k(8)),fprintf('\tk9=%.9f\n',k(9))

fprintf('\tk10=%.9f\n',k(10)),fprintf('\tk11=%.9f\n',k(11)),fprintf('\tk12=%.9f\n',k(12))

fprintf('\tk13=%.9f\n',k(13)),fprintf('\tk14=%.9f\n',k(14)),fprintf('\tk15=%.9f\n',k(15))

fprintf('\tk16=%.9f\n',k(16)),fprintf('\tk17=%.9f\n',k(17)),fprintf('\tk18=%.9f\n',k(18))

fprintf('\tk19=%.9f\n',k(19))

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

%以函数fmincon()估计得到的结果为初值，使用函数lsqnonlin()进行参数估计

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()的结果为初值，使用函数lsqnonlin()估计得到的参数值:\n')

fprintf('\tk1=%.9f\n',k(1)),fprintf('\tk2=%.9f\n',k(2)),fprintf('\tk3=%.9f\n',k(3))

fprintf('\tk4=%.9f\n',k(4)),fprintf('\tk5=%.9f\n',k(5)),fprintf('\tk6=%.9f\n',k(6))

fprintf('\tk7=%.9f\n',k(7)),fprintf('\tk8=%.9f\n',k(8)),fprintf('\tk9=%.9f\n',k(9))

fprintf('\tk10=%.9f\n',k(10)),fprintf('\tk11=%.9f\n',k(11)),fprintf('\tk12=%.9f\n',k(12))

fprintf('\tk13=%.9f\n',k(13)),fprintf('\tk14=%.9f\n',k(14)),fprintf('\tk15=%.9f\n',k(15))

fprintf('\tk16=%.9f\n',k(16)),fprintf('\tk17=%.9f\n',k(17)),fprintf('\tk18=%.9f\n',k(18))

fprintf('\tk19=%.9f\n',k(19))

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

% 拟合效果图(实验与拟合的比较)

a=linspace(tspan(1),tspan(end),200);

[a b]=ode45(@KineticEqs,a,x0,[],k);

b1(:,1)=b(:,1);

figure;

plot(tspan,yexp,'o',a,b1,'b*');

hold on

% 残差关于拟合值的残差图

a=linspace(tspan(1),tspan(end),3);

[a c]=ode45(@KineticEqs,a,x0,[],k);

c1(:,1)=c(:,1);

figure;

plot(residual,'*')

xlabel('拟合(单位未知)')

ylabel('残差R(单位未知)')

refline(0,0)

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

function f=ObjFunc4Fmincon(k,x0,yexp)

tspan=[0.2 0.4 0.6 0.8];

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

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

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)+sum((y(:,5)-yexp(:,5)).^2)...

+sum((y(:,6)-yexp(:,6)).^2)+sum((y(:,7)-yexp(:,7)).^2)+sum((y(:,8)-yexp(:,8)).^2);

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

function f=ObjFunc4LNL(k,x0,yexp)

tspan=[0.2 0.4 0.6 0.8];

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

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

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

f3=y(:,3)-yexp(:,3);f4=y(:,4)-yexp(:,4);f5=y(:,5)-yexp(:,5);

f6=y(:,6)-yexp(:,6);f7=y(:,7)-yexp(:,7);f8=y(:,8)-yexp(:,8);

f=[f1;f2;f3;f4;f5;f6;f7;f8];

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

function dxdt=KineticEqs(t,x,k)

dxdt=...

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

(k(1)*x(1)-(k(5)+k(6)+k(7)+k(8)+k(9))*x(2))

(k(5)*x(2)-(k(10)+k(11)+k(12)+k(13)+k(14))*x(3))

(k(10)*x(3)-(k(15)+k(16))*x(4))

(k(2)*x(1)+k(6)*x(2)+k(11)*x(3)-(k(17)+k(18))*x(5))

(k(3)*x(1)+k(7)*x(2)+k(12)*x(3)+k(15)*x(4)+k(17)*x(5)+k(19)*x(6))

(k(8)*x(2)+k(13)*x(3)+k(16)*x(4)+k(19)*x(6))

(k(4)*x(1)+k(9)*x(2)+k(14)*x(3)+k(18)*x(5))

];