2006年第三届研究生数学建模竞赛B题优秀论文.pdf
- 文档编号:14648634
- 上传时间:2023-06-25
- 格式:PDF
- 页数:57
- 大小:1.99MB
2006年第三届研究生数学建模竞赛B题优秀论文.pdf
《2006年第三届研究生数学建模竞赛B题优秀论文.pdf》由会员分享,可在线阅读,更多相关《2006年第三届研究生数学建模竞赛B题优秀论文.pdf(57页珍藏版)》请在冰点文库上搜索。
目录目录一、问题的提出.1二、问题的分析.2三、符号约定.3四、模型的假设.4五、问题二体轨道力学和运动方程的模型.45.1二体问题和运动方程:
.45.2进一步简化的运动方程:
.55.3轨道的几何方程的进一步探究:
.5六、问题-1解析解模型.76.1解析解模型的分析:
.76.2解析解模型的求解:
.8七、问题-1求最小组数的最小二乘法模型.97.1、离散格式下的求最小m值的算法.97.2、利用解析解求最小m值的算法.127.3模型的求解与检验.13八、问题-3高精度差商的最小二乘法模型.168.1解析解模型的求解和讨论.168.2考虑了高精度差商的最小二乘法模型的分析和求解.18九、问题-3变分伴随同化模型.279.1问题的分析.279.2问题的求解:
.329.3对模型误差敏感度的检验:
.37十、问题-4中央差商的最小二乘法模型.3910.1分析和模型修改.3910.2结果讨论.4010.3对模型误差敏感度的检验:
.42十一、问题-4变分伴随同化模型.4511.1问题的分析.4511.2问题的求解:
.4511.3对模型误差敏感度的检验:
.48十二、模型的推广.4912.1背景知识:
.4912.2考虑捕获量这个外界因素的影响:
.5012.3考虑更为广泛的捕食系统:
.51十三、模型的优缺点分析.54十四、关于本问题的几点讨论.54十五、模型的扩展.55一、问题的提出一、问题的提出包括“神舟六号”载人航天宇宙飞船、人造地球卫星等航天器围绕地球在轨运行的过程中,要受到很多力的作用,其中主要的是地球万有引力和航天器发动机作用力。
一:
一:
考虑航天器在仅受到地球万有引力、航天器自身发动机作用力的作用下作平面运动,将地球和航天器视为质点,试建立航天器运动的数学模型(只要列出模型,不要求解)。
显然这样的数学模型在精度上是远远不能满足实际需要的,在其他要求精确制导等有关高科技的实际问题中,我们都面临着类似的问题:
我们必须建立高精度的数学模型,必须高精度地估计模型中的大批参数,因为只有这样的数学模型才能解决实际问题,而不会出现差之毫厘,结果却失之千里的情况。
这时所建立数学模型的精度就成了数学模型的生命线。
例如上述问题中的航天器还要受到地球质量分布不均匀所引起的摄动力,大气阻力,日、月及其它星球的摄动引力的影响,以及航天器发动机为调整航天器自身姿态运作时作用力的影响。
这样不但数学模型十分复杂,而且在这些数学模型中还要涉及到许多重要的参数,如地球的引力场模型就有许多待定参数。
不仅如此,在对航天器进行测量时,还涉及到观测站的地理位置以及设备的系统误差等参数。
为此人们要设法利用长期积累的丰富的观测资料,高精度确定这些重要的参数。
由于航天器的问题太复杂,下面本题仅考虑较简单的确定高精度参数问题。
假设有一个生态系统,其中含有两种生物,即:
A生物和B生物,其中A生物是捕食者,B生物是被捕食者。
假设t时刻捕食者A的数目为xt,被捕食者B数目为,它们之间满足以下变化规律:
yt1234xtxtytytytxt初始条件为:
0506xtyt其中为模型的待定参数。
1kk6通过对此生态系统的观测,可以得到相关的观测数据。
观测数据的格式依次为:
观测时刻、A生物数目、B生物数目jt)(jtx)(jty二:
二:
请利用有关数据,解决以下问题:
1)在观测数据无误差的情况下,若已知215,求其它5个参数?
有关数据见数据文件:
DATA1.TXT1,3,4,5,6kk2)在观测数据无误差的情况下,若2也未知,问至少需要多少组观测数据,才能确定参数1kk6?
有关数据见数据文件:
DATA1.TXT3)在观测资料有误差(时间变量不含有误差)的情况下,请分别利用观测数据DATA2.TXT和DATA3.TXT,确定参数1kk6在某种意义下的最优解,并与仿真结果比较,进而改进你们的数学模型。
4)假设连观测资料的时间变量也含有误差,试利用数据DATA4.TXT,建立数学模型,确定参数在某种意义下的最优解。
1kk6二、问题的分析二、问题的分析问题问题中建立二体轨道力学和运动方程的模型来解决此问题。
分别分析了二体问题和运动方程,进一步简化的运动方程,以及讨论了轨道的几何方程和轨道可能的形状。
问题问题实际包含了四个小问:
(1)由于2已知,我们主要是用求解解析式的方法来直接求解出)6,5,4,3,1(kk
(2)采用最小二乘法来求解。
假设有m组数据,分解析解和离散差分格式这两种情况来讨论。
离散情况下:
首先将原微分方程组按某种差分格式离散化,然后将组观测数据离散后的差分格式,构造以参数为未知变量的二元一次线性方程组,根据最小二乘法的原理,结合矩阵的相关理论,讨论解存在、唯一、稳定的条件,从而确定求解参数需要观测资料的组数m。
解析解的利用:
从原始方程组求出包含未知参数的解析解,同样将m组数据代入解析解,构造出超定的二元一次方程组,然后用同
(1)的方法讨论,最终确定最小的。
(3)先用的是高精度差商的最小二乘模型,求出最优解,用仿真解与观测资料进行比较,效果非常好。
由于此问要考虑观测资料的误差,所以就必须对模型的误差敏感度进行讨论,我们发现高精度差商的最小二乘模型的敏感度也是比较好的。
然后我们又用了变分伴随同化模型,在精确度和敏感度上面做了进一步的改进。
(4)与(3)所用方法基本一致,也是先用中央差商的最小二乘模型求最优解,再用变分伴随同化模型作进一步的改进。
虽然此时不仅要考虑观测资料的误差还要考虑时间误差,但是模型的稳定性还是很好。
miyxii,.,1,mm三、符号约定三、符号约定)(tx-时刻捕食者A的数目t)(ty-被捕食者B数目)6,.,2,1(ii-种群模型中的待定参数jt-观测时间序列m-观测资料的数据组数gF-地球和航天器的万有引力G-万有引力常数M-地球质量m-航天器的质量r-地球与航天器之间的距离K-积分常数A-系数矩阵b-线性方程组右端的向量J-目标函数iilk,-中间迭代值h-等距的时间步长randn-0到1之间的随机数N-用来表示数据精度的正数0-参数的初值kd-搜索方向k-搜索步长w-权重系数QP,-伴随模式中的伴随变量obsx,-观测资料数据obsymeanx,-,meanyAB种群的平均值四、模型的假设四、模型的假设1考虑航天器在仅受到地球万有引力、航天器自身发动机作用力的作用下作平面运动;2将地球和航天器视为质点;3假设题目给定的数据资料完全正确合理,具有实际意义;五、问题二体轨道力学和运动方程的模型五、问题二体轨道力学和运动方程的模型主要用此模型来解决题目中问题。
5.1二体问题和运动方程:
二体问题和运动方程:
无论是航天器还是天然天体,它们在运动中的任何给定时刻,均会受到多个周围天体的万有引力作用,甚至还有其他的力,例如阻力、推力和太阳辐射压力等的作用。
本文中已经进行了简化,只考虑地球和航天器两个质体,只考虑地球万有引力和航天器发动机作用力。
分析的第一步应选择一个适于物体运动的坐标系。
这件事做起来并不容易,因为所选的任何坐标系其惯性特征都存在某种程度的不确性。
为不失一般性,假定存在某个合适的惯性坐标系ZYXO,分析所受得力:
2FFFg(其中rrGMmFg3)所以:
Fmvdtd)(把对时间的导数展开,得到:
Fdtdmvdtdvm(1-1)航天器可能不断排出某些质量以产生推力,所以dtdmv不为零,各项除以,可得航天器的一般运动方程为:
m,0000rrrrmmrmFrtt(1-2)r和r分别为航天器相对于惯性坐标系ZYXO的速度矢量和加速度矢量。
和分别是物体的质量和质量随时间变化。
mm5.2进一步简化的运动方程:
进一步简化的运动方程:
首先,作两个简化假设:
(1)物体为球对称的,这样就可以把物体看作质量集中在其中心。
(2)除了沿两物体中心连线作用的引力外,没有其他外力和内力作用。
由5.1可知此时方程简化为:
rrmMGr3)(1-3)又因为航天器的质量m将比天体的质量M小得多,因此有:
GMmMG)(所以运动方程又转化成:
03rrr(其中GM)(1-4)5.3轨道的几何方程的进一步探究:
轨道的几何方程的进一步探究:
根据5.2得到的二体的简化方程(1-4),该方程在形式上十分简单,但已经完全确定了相应轨道的形状和大小。
将方程式(1-4)两边同时与h叉乘,有rhrhrrhr33(1-5)考虑到h守恒和矢量运算规则:
)()()(cbacabcba及rrrr,所以hrhrhrhrdtd)()()(233rrdtdrrrvrrvrrrhr于是,可以将式(5)改写为)()(rrdtdhrdtd两边积分得Brrhr这里B是积分常矢量。
用r点乘该式就得到标量方程Brrrrhrr因为2,aaacbacba总成立,所以得到vrBrhcos2式中,v为常矢量B和矢径r之间的夹角。
解出r得轨道的几何方程为vBhrcos)/(1/2(1-6)令,/2hp/Be,则上式即为veprcos1(1-7)显然,轨道的几何方程是一个圆锥曲线的极坐标方程,中心引力体质心即为极坐标的原点,位于一焦点上,极角v为r与圆锥曲线上离焦点最近的一点与焦点连线间的夹角,常数p户称为“半正焦弦”,常数称为“偏心率”,它确定了方程式(7)表示的圆锥曲线的类型。
至此,可以把航天器的轨道运动总结如下:
(1)圆锥曲线族(圆、椭圆、抛物线、双曲线)为二体问题中的航天器惟一可能的运动轨道。
(2)中心引力体中心必定为圆锥曲线轨道的一个焦点。
(3)当航天器沿着圆锥曲线轨道运动时,其比机械能(单位质量的动能和势能之和)保持不变。
然而,动能和势能这两种形式的能量却可以相互转换。
这意味着,当航天器高度增高(r增加)时,其速度必定变慢,当r减少时,速度加快。
正是以这种方式,使保持不变。
(4)航天器绕中心引力体运动,当r和v沿轨道变化时,比角动量h保持不变。
(5)轨道运动总是处在一个固定于惯性空间的平面内。
六、问题六、问题-1解析解模型解析解模型问题要求当观测数据无误差时,并且已知2,求其他五个参数。
6.1解析解模型的分析:
解析解模型的分析:
由题目所给的微分方程组模型:
)()()()(4321txtydtdytytxdtdx将以上方程组的自变量t消去,得到的y与x之间的函数关系为:
xyxxyydxdy2143采用分离变量方法得到其通解为:
Kxxyy4321lnln(2-1)其中K为任意常数,如果给定了初始条件,则),(00yxK为定值。
即:
04030201lnlnxxyyK(2-2)由上面两式可得:
0)()ln()()ln(04030201xxxxyyyy(2-3)其中已知2.02,(2-3)式可以化为:
)(2.0)()ln()ln(0040301yyxxxxyy任意时刻的观测数据都应该满足上式。
假设有)(,iiyx1m组观测数据(包括初始时刻),利用上式就可以得到如下形式的方程组:
1133mmBA其中为待求参数,,T),(43113AB为系数矩阵。
通过求解上述方程组就可以得到要求的参数6.2解析解模型的求解:
解析解模型的求解:
利用data1.txt中提供的6组观测数据,可以得到式中的系数矩阵,求解得到:
1535,BA-21,2.02,213,14又因为所以:
)()(0605tytx105,606。
这样利用初值,由(2-2)式可得到:
82.13K,这样就可以求出关于x和y的关系式:
082.13ln122.0ln2xxyy(2-4)我们就可以求出到1.0t5.0t时间仿真的),(YX,跟真实的进行比较,如下图所示:
),(YX图1解析解模型仿真结果从图上看,可以观测值和仿真值吻合很好,所以我们这个)6,5,4,3,2,1(kK参数求得还是很好的。
七、问题七、问题-1求最小组数的最小二乘法模型求最小组数的最小二乘法模型分两种情形进行讨论。
7.1、离散格式下的求最小、离散格式下的求最小m值的算法值的算法原微分方程组为:
1234xtxtytytytxt0506xtyt将方程组左边的导数项离散化,采用向后差分格式:
).2(.)()0()()0()()()()()()()()()()(161511431121mityytxxtttytytytxtytttxtxtytxtxiiiiiiiiiiiiii将i看作未知变量,可以写出此方程组的矩阵形式:
bA,其中10000001000000)()()()(0000)()(000000)()()()(0000)()(2222222mmmmmtytxtytxtytytytxtytxtxtxA6)22(nA由分块矩阵组成,左上的两个分块阵都有n行数据(,即将第一个数据作为初始值),b同样可以写成如下格式:
1mn1)22(1112121212)()()()()()()()()()(nmmmmmmmmtytxtttytytttytytttxtxtttxtxb而,T654321因此问题转化为要求bA存在唯一稳定的非零解时至少需要多少组数据,即最小为多少时可以求出m。
一般情况下,bA所表示的是一个超定的线性方程组(方程个数大于未知数的个数)。
在最小二乘的意义下,可以求出该线性方程组的一组最优解,设为p,使其满足误差min!
2bAJPStep1:
证明“证明“min!
2bAJP”等价于“”等价于“bAAATpT”证明:
证明:
记为欧氏空间的内积,构造泛函:
bAbAJPP,21可以推出0,bAAAbAAJTpTppp利用p的任意性,得出bAAATpT,得证。
Step2:
证明证明bAAATpT可解可解证明:
证明:
设的相应齐次方程为bAAATpT0yAAT,只需证明0,ybAT,即垂直于bATy。
又因为0,yAAyAAybATppTT其中,为6R中内积,于是bAAATpT,得证。
Step3:
证明:
证明bAAATpT的解不唯一的解不唯一证明:
证明:
如果min!
2bAJP存在一个解,而*x0A有非零解,则有:
xmin!
)(2*bxtxAJ仍然是成立的。
不妨记该问题的最小二乘解集为,从该解集中取一个代表元作为问题的解,xbAxxxmin(广义逆)。
由Step1的证明,方程组有非零解的条件是0AAT,此时解为:
。
bAbAAATTp11)(若对m从2开始往上搜索,求出每种情况下的AAT和,并计算出每个值对应的模拟误差,由此综合分析比较,使得行列式m0AAT同时模拟误差在容许的范围之内的m即为所求。
7.2、利用解析解求最小、利用解析解求最小m值的算法值的算法首先来求原种群模型微分方程组的解析解,将原微分方程组改写为;)()(4321xydtdyyxdtdx两式相比,化简为dyyydxxx2143两边同时积分,移项得Kxxyy4321lnln结合初始条件)1(),1(65yx与初始值有关其中,,K为积分常数Kconstxxyy)1()1(ln)1
(2)1(ln431将m组数据分别代入该解析解,构造出包含未知参数的二元一次线性方程组,其矢量形式记为:
11bA,此时系数矩阵和等式右边的向量变为如下:
1000010000)()(00)(ln)(ln00)()(00)(ln)(ln2222mmmmtxtxtxtxtytytytyA)()(211tytxKKb方程组有非零解的条件是011AAT,此时K可以赋值为1,它并不影响行列式的求解。
此时解为:
11111)(bbAAATTp11A。
同上,对从2开始往m上搜索,求出每种情况下的AAT和,并计算出每个值对应的模拟误差,由此综合分析比较,使得行列式m0AAT同时模拟误差在容许的范围之内的即为所求。
m7.3模型的求解与检验模型的求解与检验根据上面的算法,利用data1.txt的对上述两种情况进行检验。
7.3.1、离散情况下求最小、离散情况下求最小m值值表2-1是m分别取2,3,等时对应的系数行列式AAT,二元一次线性方程组的最小二乘解i和对应误差J的情况。
由表分析可知,m=2时,AAT的数量级是10的-31次方,近似的可以认为行列式的值为零,因此=2对应的线性方程组无解。
当3后,行列式均是大于零的,即此种情况下线性方程组的是存在的,但不是每个存在的解都是稳定的,为找到综合意义上的最优解,还要比较通过各解模拟的模拟误差,模拟误差越小,说明其对应的解越稳定,此时的m就是要找的最小组数。
mm而表2-1中m=3与=4、5比较知(m=6时由于差分格式的误差和data1.txt数据的原因波动较大,可暂不比较),m=3时误差最小,说明提供三组数据时线性方程组即存在唯一,而且能在一定程度上体现出最小二乘法的稳定性。
所以,在离散的情况下至少需要3组数据就可以确定参数m。
表1离散情况下,不同m值时的AAT、和误差AAT123456模拟误差6m6.122562e+150.9119200.0252151.67910.00938110603.752589e+25m9.520669e+81.347800-0.9068116.239-4.21310010601.824217e+104m8.733238e+81.621800-0.8285316.545-4.32040010604.308100e+53m5.188582e+60.400190-0.74352359.37-104.710010601.968029e+5数据组值数参2m-3.343748e-310.2979900.000000.0000-74.1270010604.387772e+37.3.2、利用解析解求最小、利用解析解求最小m值值利用微分方程组的解析解和初始条件构造矢量方程11bA,此时该方程等式右边的向量中包含有一个积分常数1bK,虽然是常数,但其值是未知的。
对此式稍作变化:
解析解两边同时除以K,相当于将解分量14同时缩小了K倍,则有)1(),1(1lnln654321yxxyyx对应得向量变为:
1)22
(1)1()1(11nyxb此时亦可以求出不同m值时对应得AAT、和误差J,如表2。
(模拟时为减少误差,采用改进的欧拉方法-预报校正法)表2连续情况下,不同m值时对应得AAT、和误差JAAT123456模拟误差3m1.965255e+62.1707-0.21707-13.0251.085410601.802354e+65m2.108208e+42.1706-0.217-13.0251.085410607.678974e+24m2.504550e-10-3.23982.4833-2.5602-0.1634810605.735110e+63m5.415568e-279.9162-2.9495-12.359-0.3943310601.804139e+142m3.503917e-47-3.402960-1.1774e-015010604.229816e+11组据数值数参由上表,行列式AAT在=2,3,4时均很小,可以看作是零,则此三种情况下原二元线性方程组无解或者解不惟一,从=5开始mmAAT是大于零的数,最小的观测数据组数必须大于或等于5。
比较=5,6两种情况的模拟误差知,前者的误差要远小于后者,也即=5时,解的稳定性比较好。
由此可以得出结论:
连续的情况下至少需要5组观测数据才能完全确定参数mm。
表2-2的数据分析还在一定程度上检验了最小二乘法的稳定性。
当观测数据组数很小时,如=2,3,参数解变化剧烈,行列式近似为零,方程组是欠定的。
误差也很大;当观测数据增多到一定时,如=5,6,参数解相差很小,基本趋于稳定了。
当然,数据越多,参数解就会越精确。
mm7.3.3、两种情况的比较、两种情况的比较综合上述两者情况可以得出,至少需要3组数据才能确定参数i。
图2和图3是连续情况下用data1.txt的前5组数据确定参数,然后模拟并与观测数据的比较情况。
模拟时采用改进的欧拉方法预报校正方法。
由于数据较少,模拟具有一定的误差,但能够反映出总体的趋势。
总的来说离散情况下的模拟效果要比连续情况下稍差。
图2用data1.txt前5组观测数据确定的参数进行模拟对比(解析解情况)图3用data1.txt前5组观测数据确定的参数进行模拟对比(解析解情况)八、问题八、问题-3高精度差商的最小二乘法模型高精度差商的最小二乘法模型在用高精度差商的最小二乘法模型求解的之前,我们还是用问题一的解析解方法试着求解了一下,结果求不出来,但是在理论上还是有指导作用的,在这里也先把它写出来,起到抛砖引玉的作用.8.1解析解模型的求解和讨论解析解模型的求解和讨论8.1.1解析解模型的求解:
解析解模型的求解:
对题目所给的微分方程组模型:
)()()()(4321txtydtdytytxdtdx在6.1节求出其通解为:
Kxxyy4321lnln(12)其中K为任意常数,如果给定了初始条件,则),(00yxK为定值。
即:
04030201lnlnxxyyK由于K为任意常数,所以可以将方程(12)两边同除K,可得:
1lnln4321xKxKyKyK(13)将K1,K2,K3,K4重新看成新的待求未知数,设为,。
然后用data2中的数据带入,1Z2Z3Z4Zylny,xlnx求出系数阵,这样就是解AIAZ的矩阵了。
可以求出:
0.141111Z,-0.00707342Z,-0.716643Z,0.0716914Z代入(13)式,这样就可以求出关于x和y的关系式:
1071691.0ln71664.00070734.0ln14111.0xxyy这就是求出的解析解模型。
8.1.2解析解模型的讨论:
解析解模型的讨论:
根据前面求出的解析解模型,我们根据实际给出的X值,就可以求出相应的Y值,将这时的求出的和实际的进行比较,画图如下:
),(YX),(YX图4求出的和实际的比较图),(YX),(YX求解值和实际值吻合的很好,可见我们这个解析解的模型还是比较准确的。
但是由于K的值并不知道,所以此题并不能求出)6,5,4,3,2,1(kK,但是其在理论上有很好的指导意义!
8.2考虑了高精度差商的最小二乘法模型的分析和求解考虑了高精度差商的最小二乘法模型的分析和求解问题三在观测资料有误差(时间变量不含有误差)的情况下,同样可以通过最小二乘方法来解决。
在此主要将不同于问题二的部分加以讨论。
8.2.1、分析与模型变化、分析与模型变化此时观测数据比data1.txt增多,可以将问题二的模型基础上,将一阶微商设置成比一阶差商精度更高的差商格式,过程如下:
由泰勒级数:
32)(!
3)()(!
2)()()()(xxfxxfxxfxfxxf可以写出如下的变形式:
32)(!
3)()(!
2)()()()(xxfxxfxxfxfxxf32)2(!
3)()2(!
2)
(2)()()2(xxfxxfxxfxfxxf32)2(!
3)()2(!
2)
(2)()()2(xxfxxfxxfxfxxf从上面四个展式出发,可以构造出具有4阶精度的一阶导数的差分表达式:
)(12)2()(8)(8)2(xfxxxfxxfxxfxxf所以得到xxxfxxfxxfxxfxf12)2()(8)(8)2()(利用上面的结果,则此时离散形式下的矢量方程等式右端的向量可以改写成:
1)22(2121122112)()(12)()(8)(8)(12)()(8)(8)(niiiiiiiitytxttytytytyttxtxtxtxb;4;2,4,3mnmi其中,(说明:
该方法要求有m4组数据,虽然精度很高,但只适合于观测数据较多的情况,比如data2.txt和data3.txt。
)原微分方程组对应的龙格库塔方法下的差分方程方程组为:
),()
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 2006 第三 研究生 数学 建模 竞赛 优秀论文