小行星轨道模型(论文)

-1-小行星轨道模型摘要围绕小行星运动观测点来推测小行星运动轨迹，通过开普勒第一定律，小行星运动轨迹为椭圆。将实际问题转化为数学问题，利用二次曲线，积分法求椭圆方程和椭圆周长。一、问题重述据报道：2013年2月16日有一颗直径大约50米的小行星与地球擦肩而过，小行星撞击地球危险可能再度引起公众的关注。要确定一颗小行星绕太阳运行的轨道，需要在轨道平面内建立以太阳为原点的空间直角坐标系，然后在不同时刻对小行星进行观测，以确定其轨道。已知在5个不同时刻对某颗小行星进行了5次观测，表B给出了相应的观测数据。表B：某小行星的5次观测数据（单位：天文单位）12345X坐标5.7646.2866.7597.1687.480Y坐标0.6481.2021.8322.5263.360注：一个天文单位等于地球到太阳的平均距离，即11101.4959787米。确定这颗小行星的轨道，如椭圆的半长轴、半短轴、半焦距、近日点、远日点，以及椭圆轨道的周长等。二、问题分析根据开普勒第一定律知，小行星的轨迹是以太阳为焦点椭圆，求解小行星轨迹，将小行星轨迹问题转化为数学模型问题，现建立椭圆方程（221234522210axaxyayaxay）利用5个点确定二次曲线的一般方程，并求出椭圆的重要参数。因此将资料锁定在线性代数方程组理论和椭圆的有关概念上。即用5个点的坐标数据分别带入椭圆的一般方程的线性代数方程组𝑎1𝑥2+2𝑎2𝑥𝑦+𝑎3𝑦2+2𝑎4𝑥+2𝑎5𝑦+1=0，该方程组的系数为A，右端顶为b，则：{𝑎1𝑥12+2𝑎2𝑥1𝑦1+𝑎3𝑦12+2𝑎4𝑥1+2𝑎5𝑦1+1=0𝑎1𝑥22+2𝑎2𝑥2𝑦2+𝑎3𝑦22+2𝑎4𝑥2+2𝑎5𝑦2+1=0𝑎1𝑥32+2𝑎2𝑥3𝑦3+𝑎3𝑦32+2𝑎4𝑥3+2𝑎5𝑦3+1=0𝑎1𝑥42+2𝑎2𝑥4𝑦4+𝑎3𝑦42+2𝑎4𝑥4+2𝑎5𝑦4+1=0𝑎1𝑥52+2𝑎2𝑥5𝑦5+𝑎3𝑦52+2𝑎4𝑥5+2𝑎5𝑦5+1=0；将这一包含5个未知数的线性方程组，写成矩阵的形式-2-22111111122222222222323333322424444422525555512221222122212221222axxyyxyaxxyyxyaxxyyxyaxxyyxyaxxyyxy。求解这一线性方程组，即可得到曲线方程的系数。三、模型的建立与求解根据上述建立的模型，对数据进行处理，利用MATLAB对数据进行处理编程得：图表1求解结果：图表2则小行星运行轨迹的一般方程为：220.05080.07020.03810.45310.264310xxyyxy[1]轨迹图如图所示：-3-图表3方案一：利用椭圆一般方程与椭圆标准方程之间的关系求解椭圆参数椭圆的一般方程:A𝑥2+𝐵𝑥𝑦+𝐶𝑦2+𝐷𝑥+𝐸𝑦+1=0椭圆几何中心:𝑋𝑐=𝐵𝐸−2𝐶𝐷4𝐴𝐶−𝐵2𝑌𝑐=𝐵𝐷−2𝐴𝐸4𝐴𝐶−𝐵2长短半轴分别为:𝑎2=2(𝐴𝑋𝑐2+𝐶𝑌𝑐2+𝐵𝑋𝑐𝑌𝑐−1)𝐴+𝐶−√(𝐴−𝐶)2+𝐵2𝑏2=2(𝐴𝑋𝑐2+𝐶𝑌𝑐2+𝐵𝑋𝑐𝑌𝑐−1)𝐴+𝐶+√(𝐴−𝐶)2+𝐵2离心率：𝑒=𝑐𝑎=√𝑎2−𝑏2𝑎故可编辑程序求解（只需做简单的四则运算）：-4-求解结果为：a=2.4631,b=0.8154,c=2.3242,e=0.9436于是椭圆长半轴a=2.4631,短半轴b=4.2147,半焦距c=2.3242。小行星近日点距离h1=a-c=0.1389，远日点距离h2=a+c=4.7873。最后计算面积、周长：根据椭圆面积公式abS则有S=6.3095由于椭圆周长计算较为复杂，且对于小星星来说要求精度较高，故采用220.128()4[4]()2.8ababCababab椭圆周长（a+b）-4此公式误差为十万分之三此时π取值为3.1416则利用MATLAB编程求解则求解得10.9611C椭圆周长最终结果：-5-椭圆标准方程为：22221;2.46310.8154xy长半轴：a=2.4631；短半轴：b=0.8154；半焦距：c=2.3242；离心率：e=0.9436；面积：S=6.3095；周长：L=10.9611；近日点：h1=0.1389；远日点：h2=4.7873。方案二：将二次曲线转化为椭圆的标准方程形式：22221xyab利用二次曲线理论求解，可得到椭圆旋转和平移两种变换后的方程如下：2212||0||DXYC123122352345,1aaaaaCDaaaaaaa12,为C的特征值。故，椭圆长半轴：21DaC；椭圆短半轴：22DbC；椭圆半焦距：22cab;离心率cea；面积abS;利用MATLAB求解：图表4特征值求解图表5求解结果则所有特征值120.0088,0.0801；7.0286e-04C同理有：0.05080.03510.0381D=0.03510.03810.1321D=0.00100.22650.13211.0000则。-6-故椭圆长半轴1DaC=12.7097,短半轴2DbC=4.2127。c=11.9912。则小行星近日点距离h1=a-c=0.7185，远日点距离h2=a+c=24.7009。则偏心率cea=0.9435，面积S=168.2076，椭圆周长的计算：将托圆周长计算转化为数学模型，已知椭圆长半轴a,偏心率e，求其周长。椭圆的方程为：xasinybcos，其微分为dxacosdybsind，弧元为2222222222cossin(1sin)sindsdxdyabdabd22222221sin1sinabcadadaa周长为/222041sinCaed其中e是偏心率。对于一个偏心率，可直接计算上述积分。第二类完全椭圆积分定义为/2220()1sinEkkxdx周长可表示为4()CaEe取长半轴为a为周长的单位，周长可表示为/2220*41sinCCeda可利用数值积分指令quadl直接对上式积分，也可利用符号积分指令int对上式积分。将参数e2向量和积分变量φ化为矩阵，则能利用梯形求和指令trapz求周长。用MATLAB的函数ellipke可直接计算第二类完全椭圆积分之值，不过，第二类完全椭圆积分定义为/220()1sinEmmxdx因此椭圆周长为24()CaEe将周长的椭圆积分公式按二项式展开得-7-π/2210(21)!!4[1(sin)]d2!(21)nnnnCaenn其中(2n–1)!!=1·3…·(2n–1)。利用积分π/220(21)!!πsind2!2nnnxxn可得用无穷级数表示的椭圆周长211(21)!!2π{1[]}212!nnnnCaenn图表6故椭圆的周长为L=56.5684做出标准椭圆图像：图表7四、模型评价-8-由以上两种模型方案解得不同结果，故可根据实际天文观察对模型进行取舍则取用选取方案二进行求解，此模型更加贴近实际解决问题，兼顾求解的复杂程度，在误差允许范围内可以将两者互为补充，互相印证。可拓展性强，由于本模型充分考虑了二次曲线利用积分法求解椭圆周长，结果更为精确，与公式法相比更具有适应性，可迁移到其它相关模型求解缺点：由于对天文知识的缺乏，对于天文学上如何计算分析求解了解甚少，以上模型仅适用于数学分析，考虑到天体其它星体对小行星轨迹的影响，在此没有找出相关资料作为解答。五参考文献[1]张志勇.精通MATLAB6.5[M].北京：北京航空航天大学出版社，2003[2]周祖逵.椭圆周长近似公式[J]．数学通报，1995[3]姜启源.数学建模第四版[M].北京：高等教育出版社，2011-9-