1、工程應(yīng)用數(shù)學(xué)作業(yè)（lsqnonlin和lsqcurvefit）材研1301目的：學(xué)會(huì)使用lsqnonlin 和lsqcurvefit兩種不同的最優(yōu)化函數(shù)，并比較用微分法和積分法來(lái)估算微分動(dòng)力學(xué)方程的參數(shù)。過(guò)程：一般的動(dòng)力學(xué)方程如下：，分別用微分法和積分法進(jìn)行反應(yīng)速率分析得到速率常數(shù)k和反應(yīng)級(jí)數(shù)n。1、lsqnonlin微分法function KineticsEst1_Diff % 動(dòng)力學(xué)參數(shù)辨識(shí): 用微分法進(jìn)行反應(yīng)速率分析得到速率常數(shù)k和反應(yīng)級(jí)數(shù)nclear allclc % 動(dòng)力學(xué)數(shù)據(jù)t = 0 20 40 60 120 180 300;CAm = 10 8 6 5 3 2 1; % 用最小

2、二乘樣條擬合法計(jì)算微分dCA/dt-使用不經(jīng)過(guò)實(shí)驗(yàn)點(diǎn)的B樣條插值函數(shù)knots = 3;K = 3; % 三次B樣條sp = spap2(knots,K,t,CAm);pp = fnder(sp); % 計(jì)算B樣條函數(shù)的導(dǎo)函數(shù)dCAdt = fnval(pp,t) % 計(jì)算t處的導(dǎo)函數(shù)值 % 繪制濃度擬合曲線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、;實(shí)驗(yàn)值','B樣條擬合') % 非線性擬合beta0 = 0.0053 1.39;beta,resnorm,residual,exitflag,output,lambda,jacobian = . lsqnonlin(OptObjFunc,beta0,rAm,CAm); ci = nlparci(beta,residual,jacobian); % 參數(shù)辨識(shí)結(jié)果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) % 繪制反應(yīng)速率擬合曲線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% 動(dòng)力學(xué)參數(shù)辨識(shí): 用積分法進(jìn)行反應(yīng)速率分析得到速率常數(shù)k和反應(yīng)級(jí)數(shù)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) % 擬合效果圖(實(shí)驗(yàn)與擬合的比較) t4plot CA4plot = ode45(KineticsEqs,tspan(1) tspan(end),CA0,beta);plot(t,CAm,'bo',t4plot,CA4plot,'k-')legend('Exp','Model')xlabel('時(shí)間, s')ylabel('濃度, mol/L') % 殘差關(guān)于擬合值的殘差圖t CAc = ode45(KineticsEqs,tspan,CA0,bet

9、a);figureplot(CAc,resid,'*')xlabel('濃度擬合值（mol/L)')ylabel('殘差R (mol/L)')refline(0,0) % 參數(shù)辨識(shí)結(jié)果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('實(shí)驗(yàn)值','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('時(shí)間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