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

工程应用数学作业（lsqnonlin和lsqcurvefit） 材研1301 目的：学会使用lsqnonlin 和lsqcurvefit两种不同的最优化函数，并比较用微分法和积分法来估算微分动力学方程的参数。 过程：一般的动力学方程如下： ， 分别用微分法和积分法进行反应速率分析得到速率常数k和反应级数n。 1、lsqnonlin微分法 function KineticsEst1_Diff % 动力学参数辨识: 用微分法进行反应速率分析得到速率常数k和反应级数n clear all clc % 动力学数据 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处的导函数值 % 绘制浓度拟合曲线 ti = 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] = ... lsqnonlin(@OptObjFunc,beta0,[],[],[],rAm,CAm); ci = nlparci(beta,residual,jacobian); % 参数辨识结果 fprintf(Estimated Parameters:\n) fprintf(\tk = %.4f ± %.4f\n,beta(1),ci(1,2)-beta(1)) fprintf(\tn = %.2f ± %.2f\n,beta(2),ci(2,2)-beta(2)) fprintf( The sum of the squares is: %.1e\n\n,sum(residual.^2)) % 绘制反应速率拟合曲线 figure plot(t,rAm,ro,t,Rate(CAm,beta),b*) xlabel(t) ylabel(dC_Adt) legend(Experiment,Kinetic Model) % ------------------------------------------------------------------ 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*CA^n, 其中k=beta(1), n=beta(2) 输入 KineticsEst1_Diff 输出 dCAdt = -0.1258 -0.0977 -0.0696 -0.0502 -0.0179 -0.0132 -0.0039 Local minimum possible. lsqnonlin 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.0025 n = 1.37 ± 0.22 The sum of the squares is: 9.9e-05 lsqnonlin积分法 function KineticsEst1_int % 动力学参数辨识: 用积分法进行反应速率分析得到速率常数k和反应级数n clear all clc global CAm t = [0