以天体问题为例求解常微分方程

《以天体问题为例求解常微分方程》由会员上传分享，免费在线阅读，更多相关内容在工程资料-天天文库。

1、以天体问题为例求解常微分方程1、Eulermethod用以对初值问题求解??=?(?,?),y(??)=??,在节点处用差商近似代替导数，则有：?????+?−?(??)=y’(??)=,可导出Euler法的数值解法，???+?????+hf(??，??),k=0,1,2…利用初值不断进行迭代，就可以求出微分方程的数值解。高阶微分方程：一个n阶方程(n)(n1)yf(t,y',y",,y)(n1)yy,yy',,yy设12n，可将上式化为一阶方程组y1'y2y'y23y'yn1ny'f(t,y,y,,y)n12n，

2、再利用Euler法求解。单体问题：以力心为原点，建立直角坐标系，根据牛顿运动定律和万有引力定律，GMm可得质点的运动微分方程为F=ma=mr=-,两个方向上的分力为r3rMmMmFx=−x,Fy=−y,其中r=x2+y2,两个分运动的微分方程分别为r3r3mx=Fx,my=Fy,将二阶微分方程组转化为一阶微分方程组：xvxGMxvx1.5x^2y^2yvyGMyvy1.5x^2y^2利用Euler法，设置初值，可得到轨道图x1=x+h*vx;y1=y+h*vy;r1=sqrt(x1^2+y1^2);vx1=v

3、x+h*(-3*x/(r^3));vy1=vy+h*(-3*y/(r^3));h=0.01;h=0.001h=0.0005优点：思路简单，适合于快速编程缺点：只有一阶精度（局部截断误差），当步数增多时，误差会因积累而越来越大那么该如何改进？？？1、向后欧拉法对于刚性微分方程（吸引解被附近快速改变的解所包围），用Euler法，误差很大，改用向后Euler法，误差较小，并且采用较大的步长，提高效率yxk+1=yxk+hf(xk+1，yxk+1)y(xo)=yo用某种办法使得斜率在右边点上是正确的。但是这种方法并不能直接求得yxk+1的迭代公式，对于单体，xk+1=x

4、k+hvxk+1GMxk+1vxk+1=vxk+h(−(x2+Y2)^1.5）k+1k+1根据不动点的迭代和Newton法可求出xk+12、梯形法hyxk+1=yxk+(f(xk，yk)+f(xk+1，yxk+hf(xk,yxk))2y(xo)=yo斜率用左端点的斜率与利用欧拉法给出的右端点的斜率的平均值来表示，提高了精度，二阶对于单体问题，写出迭代方程：GMxkvxk+1=vxk+h−x2+Y21.5kkhxk+1=xk+2(vxk+vxk+1)hGMxk+1GMxkvxk+1=vxk+2（−(x2+Y2)^1.5−(x2+Y2)^1.5）k+1k+1kkh=

5、0.01对比欧拉法，可发现梯形法精度更高一些2、Taylor方法h2hkyx=yx+hf(x，y)+f′(x,y)+…+fk−1(x,y)k+1kkkkkkk2k!k阶Taylor方法的截断误差为o（hk+1），是K阶的h2当K=2时，yx=yx+hf(x，y)+f′(x,y)，此即为我们熟悉k+1kkkkk2的蛙跳算法h2xit+1=xit+hvit+Fi(t)2mhvit+1=vit+[Fit+1+Fi(t)]2位移采用泰勒，速度采用梯形法优点：精度较高缺点：需手动去求倒数蛙跳法求单体问题：h=0.01当初速度过大，将脱离中心天体，飞离出去V=53、Rung

6、e-Kutta法中点方法：hhyxk+1=yxk+hf(xk+，yxk+f(xk,yk))22y(xo)=yo中点方法、梯形方法都是二阶Runge-Kutta族的一员当h=0.02时，h=0.01时二体问题：一质量为m的质点绕质量为M的质点运动，且M也会运动，则其运动微分方程分别为：GMm（−）GMm（−）r2r1r2r1F=ma1=mr1=,F=ma2=mr2=-,且r2rr2rGMm（x2−x1）GMm（y2−y1）对m:两个分运动方程为：mx1=my1=r2rr2rGMm（x2−x1）对M：两个分运动方程为：Mx2=-，r2rGMm（y2−y1）Mx1=-

7、r2r以下图形均是以蛙跳算法画出的下图为以中心天体为圆心建立的动坐标系下的环绕天体的运动轨道：三体：四体：