微分方程拟合

一种高质量拟合常微分方程参数的方法：低版本1stOpt和MATLAB联用作者: 月只蓝 (站内联系TA) 发布: 2014-07-010.引言就方程参数拟合问题而言，常见非线性方程（组）参数拟合和常微分方程（组）参数拟合两种。本帖主要讨论常微分方程参数的拟合问题。限于本人水平，存在纰漏甚至错误在所难免，欢迎大家批评指正，同时希望能抛砖引玉。1.影响拟合质量的因素Origin软件的Help文档总结了影响拟合质量的主要因素，总共5条，列出如下：（1）参数初值（2）模型（公式）的正确性或者适用性（3）迭代计算的中止的判断指标高低（4）数据点的充分性（5）数据点的噪声实际上，拟合算法的先进性同样影响拟合的质量，在此基础上，加上一条：（6）算法先进与否及其适用性本帖主要讨论第一个因素，即参数初值，对拟合质量的影响。在本体后续的讨论中可以看到，参数初值的选取，对最终拟合结果影响相当显著。2.本帖讨论的例子出处本帖讨论的例子出处见：/bbs/viewthread.php?tid=7595881&fpage=1直观起见，对上述例子中的问题表述如下：已知自变量t，因变量y，y和t满足一下函数式：dy/dt=0.001414/(1/(4.5727*104*k1*(1-0.02119*y)k2)+1/k3/45727)；（需要指出的是，该公式中参数可进一步合并，在本帖中并未作合并处理）已知 t=；y=；拟合出参数k1，k2，k3。3.结果与讨论运行MATLAB软件，调用lsqnonlin函数和ode45函数，参数初值按k1=k2=k3=1选取，程序代码如下：function KineticsEst1_Int_modified_by_Yuezhilanclear all;clcformat longtspan=-31; yexp=;k0=; %请注意这里，初值的选取y0=20.4122;lb=-*1e3;ub=; yy=;= . lsqnonlin(ObjFunc,k0,lb,ub,tspan,y0,yexp);fprintf(nn使用函数lsqnonlin()估计得到的参数值为:n)fprintf(t待拟合参数 k1 = %.6fn,k(1)fprintf(t待拟合参数 k2 = %.6fn,k(2)fprintf(t待拟合参数 k3 = %.6fn,k(3)fprintf( t残差平方和= %.6fnn,resnorm)ts=0:1:max(tspan);=ode45(KineticsEqs,ts,y0,k);= ode45(KineticsEqs,tspan,y0,k); y=XXsim(2:end);xexp=yexp;R2=1-sum(xexp-y).2)./sum(xexp-mean(y).2);fprintf(nt决定系数R-Square = %.6f,R2);figure(1)plot(ts,ys,b,tspan,yy,or),legend(计算值,实验值,Location,best);yr=y-yexp;figure(2)plot(tspan(2:end),yr,r*,),axis();figure(3)plot(yexp,y,ro,b-);%-function f = ObjFunc(k,tspan,y0,yexp) = ode45(KineticsEqs,tspan,y0,k) ;ysim = Xsim(2:end);size(ysim);size(yexp);f=ysim-yexp;%-function dydt = KineticsEqs(t,y,k)beta(1)=k(1);beta(2)=k(2);beta(3)=k(3);dydt = 0.001414/(1/(4.5727*104*beta(1)*(1-0.02119*y)beta(2)+1/beta(3)/45727);计算结果：使用函数lsqnonlin()估计得到的参数值为: 待拟合参数 k1 = 585.320105 待拟合参数 k2 = 14.003057 待拟合参数 k3 = 0.065162 残差平方和= 1.705667 决定系数R-Square = 0.945774Warning: Imaginary parts of complex X and/or Y arguments ignored In KineticsEst1_Int_modified_by_Yuezhilan at 30Warning: Imaginary parts of complex X and/or Y arguments ignored In KineticsEst1_Int_modified_by_Yuezhilan at 33Warning: Imaginary parts of complex X and/or Y arguments ignored In KineticsEst1_Int_modified_by_Yuezhilan at 35 由结果可见，决定系数R-Square = 0.945774，拟合质量一般，同时MATLAB警告结果含有复数。通常而言，在作拟合的时候，常赋予初值诸如（0 0 0）、（1 1 1）、（0 1 1）之类比较简单的数值以试探；更合理的方法，是事先对数据和公式进行观察，猜测初值的取值范围。不过，对于拟合公式较复杂的情况，通过观察实验数据很难判断初值的范围，比如在本例中，通过观察y随t的变化趋势可见（见附图1），基本上y随t递增，也就是说，dy/dt在t取值范围内应0，仅能大概判断k1、k2和k3大于0的可能性比较大。继续尝试其它初值的选取，结果列出如下（k1 k2 k3） 决定系数R-square（1 0 0） 0.958200（0 1 0） 0.000000（拟合失败）（0 0 1） 0.957533（1 1 0） 0.943656（1 0 1） 0.962539（0.5 1 1） 0.918617（1 0.5 1） 0.958019（1 1 0.5） 0.919118（0.5 0.5 0.5） 0.955740 由此可见，不同的初值取法，拟合结果差异显著，且在上述尝试的结果中，决定系数最高仅为0.962539。1stOpt软件处理拟合问题强大高效，使用者在操作层面上不需要赋予初值，即可获得极佳的拟合效果。低版本的1stOpt（比如1.5）软件并不支持常微分方程的拟合，但其对非线性方程的拟合已经足够强大。为了利用低版本1stOpt软件对代数方程强大的拟合能力，对本例问题初值的选取提出了一个解决方法，步骤如下：1.MATLAB计算出数据的数值导数dy/dt；2.记dy/dt=Y，原常微分方程化为Y=f（t，y，k1，k2，k3）的二元代数方程，1stOpt拟合出各个k值；3.将上述拟合出的k值作为初值，输入到上述MATLAB程序中。上述步骤1中求数值导数的MATLAB程序如下：function initial_value_determinationclear all;clct=;yexp=;knots=3;K=3;sp=spap2(knots,K,t,yexp);pp=fnder(sp);dydt=fnval(pp,t);计算结果如下：ans = 0 2.41696012708104220.412199999999999 1.000000000000000 1.99491516763736522.730319999999999 2.000000000000000 1.57287020819368924.556239999999999 3.000000000000000 1.15082524875001225.886790000000001 4.000000000000000 0.72878028930633526.690190000000001 5.000000000000000 0.30673532986265927.163340000000002 6.000000000000000 0.09421459315903727.458659999999998 7.000000000000000 0.09121807919546927.649190000000001 8.000000000000000 0.08822156523190127.576149999999998 9.000000000000000 0.08522505126833227.70635000000000010.000000000000000 0.07393137772835027.77302999999999911.000000000000000 0.06263770418836727.86195000000000012.000000000000000 0.05134403064838427.91911000000000013.000000000000000 0.04005035710840227.95721000000000014.000000000000000 0.02875668356841927.988969999999998上述步骤2中1stOpt代码和计算结果如下：Parameters k1,k2,k3;Variable t,dydt,y;Function dydt = 0.001414/(1/(4.5727*104*k1*(1-0.02119*y)k2)+1/k3/45727);Data; 0 2.41696012708104220.412199999999999 1.000000000000000 1.99491516763736522.730319999999999 2.000000000000000 1.57287020819368924.556239999999999 3.000000000000000 1.15082524875001225.886790000000001 4.000000000000000 0.72878028930633526.690190000000001 5.000000000000000 0.30673532986265927.163340000000002 6.000000000000000 0.09421459315903727.458659999999998 7.000000000000000 0.09121807919546927.649190000000001 8.000000000000000 0.08822156523190127.576149999999998 9.000000000000000 0.08522505126833227.70635000000000010.000000000000000 0.07393137772835027.77302999999999911.000000000000000 0.06263770418836727.86195000000000012.000000000000000 0.05134403064838427.91911000000000013.000000000000000 0.04005035710840227.95721000000000014.000000000000000 0.02875668356841927.988969999999998均方差(RMSE): 0.127523255347567残差平方和(SSE): 0.243932709816613相关系数(R): 0.98766224345764相关系数之平方(R2): 0.975476707151779决定系数(DC): 0.973324807678033卡方系数(Chi-Square): 0.341389826024371F统计(F-Statistic): 517.108255871652参数 最佳估算- -k1 58788769.7449094k2 26.9764629032194k3 0.0339