如何进行微分方程组数值解的拟合？

maplev12

已知关于t,x(t),p(t)的微分方程组eq1,eq2，以及一组数据t1,x1,p1，求参数A,B,p0,a，拟合已知的这组数据，使数值解的残差最小。方程组及数据如下：

1. eq1 := diff(x(t), t) = 2*A*(sinh(p(t)+q*x(t))-x(t)*cosh(p(t)+q*x(t))):
   
2. eq2 := diff(p(t), t) = -2*B*(p0*sinh(a*x(t))+p(t)*cosh(a*x(t))):
   
3. t1:=[1981,1982,1983,1984,1985,1986,1987,1988,1989,1990,1991,1992,1993,1994,1995,1996,1997,1998,1999,2000,2001,2002,2003,2004,2005,2006,2007,2008,2009]:
   
4. x1:=[0.196462844,0.23113261,0.239836251,0.254280917,0.282074055,0.303702732,0.314438263,0.316017811,0.259551387,0.241967162,0.256846413,0.300113175,0.369964846,0.353602032,0.329298272,0.321924659,0.315817924,0.33655726,0.332913476,0.331782975,0.339368402,0.361496949,0.409111121,0.440819131,0.480019876,0.508510525,0.516623829,0.550329299,0.659601287]:
   
5. p1:=[0.052431095,0.090568707,0.108529793,0.151760212,0.134661926,0.088463767,0.115832638,0.112805837,0.040632032,0.038390249,0.091789367,0.142407064,0.139643147,0.130806818,0.109249803,0.100085234,0.092970338,0.078333469,0.076198357,0.084312795,0.083003176,0.09082068,0.100253788,0.100850399,0.113100352,0.126765342,0.14162395,0.096346682,0.091130745]:

复制代码

maplelab

还有一个q.

maplev12

对……q~那得如何去做啊~~也解不出解析解啊~~

maplelab

q也做参数吗?
像A,B,p0,a一样吗?

maplev12

对呀……A,B,p0,a,q这么五个参数~~

shamohu

用1stOpt到可以很方便求解这类微分方程拟合问题，只是这道题结果不好，不知是模型公式的问题还是数据的问题。

maplev12

是数据不太好…话说难道Maple开发者如此懒惰，这样的问题都不写好函数包…DETools也成摆设了吗～
如果 麻婆 实在不好解的话～楼上能帮我把1stOpt的解决代码贴出来吗…我对这个神器没接触过～

shamohu

1stOpt求解代码（注意1stOpt不区分大小写，因此a变成了a1）：

1. 
   
2. Variable t,x,p;
   
3. ODEFunction x' = 2*A*(sinh(p+q*x)-x*cosh(p+q*x));
   
4.             p' = -2*B*(p0*sinh(a1*x)+p*cosh(a1*x));
   
5. Data;
   
6. t=[1981,1982,1983,1984,1985,1986,1987,1988,1989,1990,1991,1992,1993,1994,1995,1996,1997,1998,1999,2000,2001,2002,2003,2004,2005,2006,2007,2008,2009];
   
7. x=[0.196462844,0.23113261,0.239836251,0.254280917,0.282074055,0.303702732,0.314438263,0.316017811,0.259551387,0.241967162,0.256846413,0.300113175,0.369964846,0.353602032,0.329298272,0.321924659,0.315817924,0.33655726,0.332913476,0.331782975,0.339368402,0.361496949,0.409111121,0.440819131,0.480019876,0.508510525,0.516623829,0.550329299,0.659601287];
   
8. p=[0.052431095,0.090568707,0.108529793,0.151760212,0.134661926,0.088463767,0.115832638,0.112805837,0.040632032,0.038390249,0.091789367,0.142407064,0.139643147,0.130806818,0.109249803,0.100085234,0.092970338,0.078333469,0.076198357,0.084312795,0.083003176,0.09082068,0.100253788,0.100850399,0.113100352,0.126765342,0.14162395,0.096346682,0.091130745];
   

复制代码

结果：
均方差(RMSE) 0.0337939616362601
残差平方和(SSE) 0.0662378468982349
相关系数(R) 0.991413214890611
决定系数(DC) 0.982388463456317

a -0.198549475971414
q 0.598621875594106
b 0.0231141231043744
p0 0.124394097698298
a1 -7.53318015151156

p的拟合结果很差。

maplev12

太感谢了…我再重新整理数据试试～效果别差劲的太过分就行～