sqnonlin和lsqcurvefit 微分法和积分法来估算微分动力学方程的参数

1、工程应用数学作业（lsqnonlin和lsqcurvefit）材研1301目的：学会使用lsqnonlin 和lsqcurvefit两种不同的最优化函数，并比较用微分法和积分法来估算微分动力学方程的参数。过程：一般的动力学方程如下：，分别用微分法和积分法进行反应速率分析得到速率常数k和反应级数n。1、lsqnonlin微分法function KineticsEst1_Diff % 动力学参数辨识: 用微分法进行反应速率分析得到速率常数k和反应级数nclear allclc % 动力学数据t = 0 20 40 60 120 180 300;CAm = 10 8 6 5 3 2 1; % 用最小

2、二乘样条拟合法计算微分dCA/dt-使用不经过实验点的B样条插值函数knots = 3;K = 3; % 三次B样条sp = spap2(knots,K,t,CAm);pp = fnder(sp); % 计算B样条函数的导函数dCAdt = fnval(pp,t) % 计算t处的导函数值 % 绘制浓度拟合曲线ti = linspace(t(1),t(end),200);CAi = fnval(sp,ti);plot(t,CAm,'ro',ti,CAi,'b-')xlabel('t')ylabel('C_A')legend('

3、;实验值','B样条拟合') % 非线性拟合beta0 = 0.0053 1.39;beta,resnorm,residual,exitflag,output,lambda,jacobian = . lsqnonlin(OptObjFunc,beta0,rAm,CAm); ci = nlparci(beta,residual,jacobian); % 参数辨识结果fprintf('Estimated Parameters:n')fprintf('tk = %.4f ± %.4fn',beta(1),ci(1,2)-beta(1)

4、fprintf('tn = %.2f ± %.2fn',beta(2),ci(2,2)-beta(2)fprintf(' The sum of the squares is: %.1enn',sum(residual.2) % 绘制反应速率拟合曲线figureplot(t,rAm,'ro',t,Rate(CAm,beta),'b*')xlabel('t')ylabel('dC_Adt')legend('Experiment','Kinetic Model')

5、 % -function f = OptObjFunc(beta,rAm,CAm)rAc = Rate(CAm,beta);f = rAc - rAm; % -function rA = Rate(CA,beta)rA = -beta(1)*CA.beta(2); % -rA = -dCA/dt = k*CAn, 其中k=beta(1), n=beta(2)输入>> KineticsEst1_Diff输出dCAdt = -0.1258 -0.0977 -0.0696 -0.0502 -0.0179 -0.0132 -0.0039Local minimum possible.lsqn

6、onlin stopped because the final change in the sum of squares relative to its initial value is less than the default value of the function tolerance.<stopping criteria details>Estimated Parameters:k = 0.0055 ± 0.0025n = 1.37 ± 0.22 The sum of the squares is: 9.9e-05lsqnonlin积分法functio

7、n KineticsEst1_int% 动力学参数辨识: 用积分法进行反应速率分析得到速率常数k和反应级数nclear allclcglobal CAmt = 0 20 40 60 120 180 300;CAm = 10 8 6 5 3 2 1'% 非线性拟合beta0 = 0.0053 1.39;tspan = 0 20 40 60 120 180 300;CA0 = 10;beta,resnorm,resid,exitflag,output,lambda,jacobian = . lsqnonlin(OptObjFunc,beta0,tspan,CA0,CAm)ci = nlpa

8、rci(beta,resid,jacobian) % 拟合效果图(实验与拟合的比较) t4plot CA4plot = ode45(KineticsEqs,tspan(1) tspan(end),CA0,beta);plot(t,CAm,'bo',t4plot,CA4plot,'k-')legend('Exp','Model')xlabel('时间, s')ylabel('浓度, mol/L') % 残差关于拟合值的残差图t CAc = ode45(KineticsEqs,tspan,CA0,bet

9、a);figureplot(CAc,resid,'*')xlabel('浓度拟合值（mol/L)')ylabel('残差R (mol/L)')refline(0,0) % 参数辨识结果fprintf('Estimated Parameters:n')fprintf('tk = %.4f ±%.4fn',beta(1),ci(1,2)-beta(1)fprintf('tn = %.2f ±%.2fn',beta(2),ci(2,2)-beta(2) % -function f =

10、OptObjFunc(beta,tspan,CA0,CAm)t CAc = ode45(KineticsEqs,tspan,CA0,beta);f = CAc - CAm; % -function dCAdt = KineticsEqs(t,CA,beta)dCAdt = -beta(1)*CAbeta(2); % k = beta(1), n = beta(2)输入>> KineticsEst1_int输出Local minimum possible.lsqnonlin stopped because the final change in the sum of squares

11、relative to its initial value is less than the default value of the function tolerance.<stopping criteria details>beta = 0.0047 1.4555resnorm = 0.0941resid = 0 -0.2408 0.1813 0.0303 -0.0169 -0.0435 0.0131exitflag = 3output = firstorderopt: 0.0041 iterations: 3 funcCount: 12 cgiterations: 0 alg

12、orithm: 'trust-region-reflective' message: 1x457 charlambda = lower: 2x1 double upper: 2x1 doublejacobian = (2,1) -394.7680 (3,1) -566.6235 (4,1) -629.8164 (5,1) -588.7875 (6,1) -478.0379 (7,1) -305.6976 (2,2) -4.0422 (3,2) -5.4804 (4,2) -5.7590 (5,2) -4.5555 (6,2) -3.1228 (7,2) -1.3857ci =

13、0.0031 0.0063 1.2700 1.6410Estimated Parameters:k = 0.0047 ± 0.0016n = 1.46 ± 0.192、lsqcurvefit微分法function Lsqcurvefit_Diff clear allclcglobal CAmt = 0 20 40 60 120 180 300;CAm = 10 8 6 5 3 2 1;knots = 3;K = 3; sp = spap2(knots,K,t,CAm);pp = fnder(sp); dCAdt = fnval(pp,t) rAm = dCAdt;ti =

14、linspace(t(1),t(end),200);CAi = fnval(sp,ti);plot(t,CAm,'ro',ti,CAi,'b-')xlabel('t')ylabel('C_A')legend('实验值','B样条拟合') beta0 = 0.0053 1.39;beta,resnorm,residual,exitflag,output,lambda,jacobian = . lsqcurvefit(OptObjFunc,beta0,rAm,CAm); ci = nlparci(bet

15、a,residual,jacobian); fprintf('Estimated Parameters:n')fprintf('tk = %.4f ± %.4fn',beta(1),ci(1,2)-beta(1)fprintf('tn = %.2f ± %.2fn',beta(2),ci(2,2)-beta(2)fprintf(' The sum of the squares is: %.1enn',sum(residual.2) figureplot(t,rAm,'ro',t,Rate(CAm

16、,beta),'b*')xlabel('t')ylabel('dC_Adt')legend('Experiment','Kinetic Model') % -function f = OptObjFunc(beta,rAm)global CAmrAc = Rate(CAm,beta);f = rAcrAm+CAm; % -function rA = Rate(CA,beta)rA = -beta(1)*CA.beta(2); 输入>> Lsqcurvefit _Diff输出dCAdt = -0.1258

17、 -0.0977 -0.0696 -0.0502 -0.0179 -0.0132 -0.0039Local minimum possible.lsqcurvefit stopped because the final change in the sum of squares relative to its initial value is less than the default value of the function tolerance.<stopping criteria details>Estimated Parameters:k = 0.0055 ± 0.0

18、025n = 1.37 ± 0.22 The sum of the squares is: 9.9e-05 lsqcurvefit积分法function Lsqcurvefit_int clear allclcglobal CAmt = 0 20 40 60 120 180 300;CAm = 10 8 6 5 3 2 1'beta0 = 0.0053 1.39;tspan = 0 20 40 60 120 180 300;CA0 = 10;beta,resnorm,resid,exitflag,output,lambda,jacobian = . lsqcurvefit(O

19、ptObjFunc,beta0,tspan,CAm,CA0)ci = nlparci(beta,resid,jacobian) t4plot CA4plot = ode45(KineticsEqs,tspan(1) tspan(end),CA0, ,beta);plot(t,CAm,'bo',t4plot,CA4plot,'k-')legend('Exp','Model')xlabel('时间t, s')ylabel('浓度C_A, mol/L') t CAc = ode45(KineticsEqs

20、,tspan,CA0,beta);figureplot(CAc,resid,'*')xlabel('浓度拟合值（mol/L）')ylabel('残差R (mol/L)')refline(0,0)fprintf('Estimated Parameters:n')fprintf('tk = %.4f ± %.4fn',beta(1),ci(1,2)-beta(1)fprintf('tn = %.2f ± %.2fn',beta(2),ci(2,2)-beta(2)% -functio

21、n f = OptObjFunc(beta,tspan,CA0) t CAc = ode45(KineticsEqs,tspan,CA0,beta);f = CAc ;% -function dCAdt = KineticsEqs(t,CA,beta)dCAdt = -beta(1)*CAbeta(2); 输入>> Lsqcurvefit _int输出beta = 0.0047 1.4555resnorm = 0.0941resid = 0 -0.2408 0.1813 0.0303 -0.0169 -0.0435 0.0131exitflag = 3output = firstorderopt: 0.0041 iterations: 3