微分方程拟合

1、一种高质量拟合常微分方程参数的方法：低版本1stOpt和MATLAB联用作者: 月只蓝 (站内联系TA)    发布: 2014-07-010.引言就方程参数拟合问题而言，常见非线性方程（组）参数拟合和常微分方程（组）参数拟合两种。本帖主要讨论常微分方程参数的拟合问题。限于本人水平，存在纰漏甚至错误在所难免，欢迎大家批评指正，同时希望能抛砖引玉。1.影响拟合质量的因素Origin软件的Help文档总结了影响拟合质量的主要因素，总共5条，列出如下：（1）参数初值（2）模型（公式）的正确性或者适用性（3）迭代计算的中止的判断指标高低（4）数据点的充分性（5）数据点的噪声实际上，

2、拟合算法的先进性同样影响拟合的质量，在此基础上，加上一条：（6）算法先进与否及其适用性本帖主要讨论第一个因素，即参数初值，对拟合质量的影响。在本体后续的讨论中可以看到，参数初值的选取，对最终拟合结果影响相当显著。2.本帖讨论的例子出处本帖讨论的例子出处见：直观起见，对上述例子中的问题表述如下：已知自变量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.结果与讨论运行MAT

3、LAB软件，调用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使用函数lsqnonl

4、in()估计得到的参数值为: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

5、-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-

6、');%-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-

7、0.02119*y)beta(2)+1/beta(3)/45727);计算结果：使用函数lsqnonlin()估计得到的参数值为:        待拟合参数 k1 = 585.320105        待拟合参数 k2 = 14.003057        待拟合参数 k3 = 0.065162        残差平方和= 1.705667        决定系数R-Square = 0.945

8、774Warning: 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 ignore

9、d > 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，仅

10、能大概判断k1、k2和k3大于0的可能性比较大。继续尝试其它初值的选取，结果列出如下（k1 k2 k3）                                              决定系数R-square（1 0 0）       

11、                                           0.958200（0 1 0）                      

12、                             0.000000（拟合失败）（0 0 1）                                    

13、;               0.957533（1 1 0）                                                  

14、0;0.943656（1 0 1）                                                   0.962539（0.5 1 1）          

15、                                      0.918617（1 0.5 1）                           

16、60;                    0.958019（1 1 0.5）                                             

17、   0.919118（0.5 0.5 0.5）                                          0.955740                 

18、0;                                            由此可见，不同的初值取法，拟合结果差异显著，且在上述尝试的结果中，决定系数最高仅为0.962539。1stOpt软件处理拟合问题强大高效，使用者在操作层面上不需要赋予初值，即可获得极佳的拟合效果。低版本的1s

19、tOpt（比如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

20、=3;K=3;sp=spap2(knots,K,t,yexp);pp=fnder(sp);dydt=fnval(pp,t);计算结果如下：ans =                   0   2.416960127081042  20.412199999999999   1.000000000000000   1.994915167637365  22.73031999999

21、9999   2.000000000000000   1.572870208193689  24.556239999999999   3.000000000000000   1.150825248750012  25.886790000000001   4.000000000000000   0.728780289306335  26.690190000000001   5.00000000000

22、0000   0.306735329862659  27.163340000000002   6.000000000000000   0.094214593159037  27.458659999999998   7.000000000000000   0.091218079195469  27.649190000000001   8.000000000000000   0.08822156523

23、1901  27.576149999999998   9.000000000000000   0.085225051268332  27.706350000000000  10.000000000000000   0.073931377728350  27.773029999999999  11.000000000000000   0.062637704188367  27.86195000000

24、0000  12.000000000000000   0.051344030648384  27.919110000000000  13.000000000000000   0.040050357108402  27.957210000000000  14.000000000000000   0.028756683568419  27.988969999999998上述步骤2中1stOpt代码和计算结果如下：Para

25、meters 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.416960127081042  20.412199999999999   1.0000000000

26、00000   1.994915167637365  22.730319999999999   2.000000000000000   1.572870208193689  24.556239999999999   3.000000000000000   1.150825248750012  25.886790000000001   4.000000000000000   0.7287802893

27、06335  26.690190000000001   5.000000000000000   0.306735329862659  27.163340000000002   6.000000000000000   0.094214593159037  27.458659999999998   7.000000000000000   0.091218079195469  27.6491900000

28、00001   8.000000000000000   0.088221565231901  27.576149999999998   9.000000000000000   0.085225051268332  27.706350000000000  10.000000000000000   0.073931377728350  27.773029999999999  11.000000000000000   0.062637704188367  27.861950000000000  12.000000000000000   0.051344030648384  27.919110000000000  13.000000000000000   0.040050357108402  27.957210000000000  14.000000000000000   0.0287