抛体运动的Runge—Kutta解

达瑞 杨桂芝

文章主题:抛体运动 空气阻力 二体运动

文章内容：第4卷第2期1994年2月哲里术畜牧学院学报21.4.2.1994.一抛体运动的--解'姿缸.(畜牧学院基础部)邹牧岐(白城市工业交通学校)蕉茎(蒙医学院基础部)林长圣(民族师范学院物理系)03摘要研究抛体运动的基本方法有略去空气阻力的怍甩.亦不考虑地茸自转的影响.在惯性参考系中分析抛体运动的抛物线理论,以及将物体与地球视为二体运动的慵圆理论.本文既考虑空气阻力的作甩.又考虑到地球自转的影响,给出了抛体运动的--敷值解.得到了抛体运动过程中的速度,位置及孰遗.揭示了空气阻力及地茸自转对抛体运动影响的规律,数值结果与解析比较表明.方法不但精度高而且筒便?亦可应甩其数值分析更为复杂的抛体问题.关键词嗟旦,——方法苗^力,二侔适萄抛体运动是物理学中的一个典型问题,但实际问题相当复杂,第一.空气阻力作用的影响.在求解抛体运动问题时利用建立的主要适用于亚音速情况的()=形式.而马也夫斯基和萨布斯基则建立形式为()=的定律.其中—().并且与的值都是根据速度的间隔选择的.考虑到空气阻力作用不连续性.亚西切,德泽尔斯基等通常将标准函数()和()之间的经验关系以表格的形式给出.其中()一()/;第二.地球自转对抛体的影响.地球在空间完成复杂的运动.其进动周期为2600年,章动周期为18.6年.由于地球进动与章动的角速度很小.故一般不予考虑,但地球自转的角速度.=7.29211×10～/一般不能忽略.第三,地球表面形状的影响,我们所称地球外形是与非扰动状态的海平面相重合以引力势面为界的外形(这里的非扰动是指没有涨潮,落潮以及其它扰动),在这个意义下,地球的外形为一椭球体,其半长轴=6378988.半短轴为6356912.椭圆体轴在数值上相差42852.对于这样一个椭球俸,地球的引力势应由的引力定律确定:其中=(景:=』字'.__..._....'.....''..'..................'.●''一√?'):一..而其展开式可以由包含多项的形式表示为:(,1')..(.)于是地球引力势可表示为:=_((+)?收穑日期994一一∞10哲里木畜牧学院学报第4卷式中为地球半径,…:为无量纲系数.当然,如果地球外形为一球形,且质量均匀分布,刚=/,而且在距地心距离为处的引力加速度,与地球表面处引力加建度.之比为/.=(/);第四.地球自转对重力的影响,我们知道视重—.一,其中.为地球对一●物体的引力矢量.为物体的惯性离心力矢量,其大小为=0.引入克列罗函敷3—0.00528001后,地球表面上任意一纬度处的重力加速度与赤道上的重力加速度,之间满足—(14-日!);第五,地球的自转,使得研究相对地表面作措体运动的问题变成丁非惯性系中的力学问题,因此就必须考虑割力的作用,练上所述,考虑诸多田素的影殉而建立起来的抛体运动微分方程相当复杂,无法求其分析解.阻对具体问题可进行简化,然后寻求其解,最简单的方法为抛物线理论,该方法的基本思想是,略去空气阻力等渚因素的影啊.只考虑重力的作用,并在惯性参考系中研究抛体运动.其次是椭圆理论,椭圆理论在考虑影响困素上与抛物线理论相同,只不过是将抛物体与地球之间的运动视为二体问题.还有假定阻力满足()一+形式的阻力解析函数法,引入函数()傲变换,从而可以分量变换求解的虚速度法以及予先假定运动特性的解析函数为已知的运动特性解析法等等:述方法对某些特殊问题可以成功地进行分析求解,但还远离实际问题本身.因此研究和发展了求解抛体运动阿题的数值方法,诸如矗方法,方法,方法等等.本文则给出了数值求解在空气阻力和地球自转影响共同作用下的抛体运动的—方法,研究了抛体运动规律,并与解析解进行比较.●抛体的运动微分方程1.1空气阻力作用抛俸的运动微分方程将地球视为一惯性参考系,不考虑地球自转及其形状对抛体运动的影响,在北纬(南纬&;)处建立一坐标系,其轴沿径线向南,轴沿纬线向东,轴自地心指向天顶,设空气阻力满足:()=于是有斋—一一-鍪=一毒—一一亲一其中=/,为抛体质量1.2力作用下抛体的运动微分方程将地球视为以转动的非惯性参考系,不考虑空气阻力,地球形状等因素的影响,在北纬处建立坐标系,则其抛体的运动微分方程为嚣一2鲁:—一2(五+.囊)(2)靠=.…一1.3空气强力与力共同作用下抛体的运动微分方程建立与情形1完全相同的坐标系,则抛体的运动微分方程为第2期达瑞等:抛体运动的—解1≥一一-警2一_一…～2(5五—.五)尝;一…耋一—…2--方法(3)设三元联立一阶微舒方程组形式如:蒙=(,,,)(,,,)塞:(,..)则可有:;.+_】一+叫=.+一】=一其中,,.分别为每步迭代时的增量+自.时的+...起,逐步计算.且计算公式相同,在数值计算过程中取且则(+22+2+.)言+22+2十?)(2:+2+.)1;(,,,,^):一(.+/',/2.—/2.+/2.+/2)一(;+2+.+±/2,+:/2+./2)【.;(.+,.+,+,.五十)1(;,,,%):;(,+/,+/2.+/2,+/2),一(.+△/2,.+2/,+:/,+2/)【.(,++;土,+,,.+3);(.,,,,%)±=(.+/2,.+/2.+/2.+/2)=(.+△/2,+/2,+2/2,.+2/2)△【'一(,+△,.+,,+)若令(,)=且:1则我们可以很容易地证明--方法的精度.设:=;+++=(;,1)=—(+,+/2)一(1+/2)3一(+/',/2,+/2)=(1+△/2+'/4)112哲里木畜牧学院学报第4卷=(.+△.+,.△)=(--+2<>14一)所以+等++等+等+++等等等故其误差数量级为0(△).具有四级精确度3数值求群辽摧的程序化对三元二阶的微分方程组(),(2),(3)与六元一阶的微分方程等效,一般地有:—=(.1.'..)=面=(,.._,.)亲一=.,,,..)则—计算公式为一1=,一△'+———=+.一1一.+()…一().+()(,+(囊).=(囊).+若将.....分别记为.:....则有'一吉(一