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

以天体问题为例求解常微分方程 1、Euler method 用以对初值问题求解 = (, ),y( )= ,在节点处用差商近似代替导数，则有： −( ) y’( )= + ,可导出Euler 法的数值解法， = + +hf( ， ),k=0,1,2…利用初值不断进行迭代，就可以求出微分 方程的数值解。 n 高阶微分方程：一个 阶方程 y (n) f (t,y ,y , ,y (n1) ) 设y 1 y ,y 2 y , ,y n y (n1) ，可将上式化为一阶方程组 y y  1 2  y 2 y 3    y n1 y n   y f (t,y , y , ,y )  n 1 2 n ，再利用Euler 法求解。 单体问题： 以力心为原点，建立直角坐标系，根据牛顿运动定律和万有引力定律， 可得质点的运动微分方程为F=ma=mr =-GMm ,两个方向上的分力为 3 r r Mm Mm Fx=− x, Fy=− y,其中r= x 2 + y 2, 两个分运动的微分方程分别为 3 3 r r mx =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=vx+h*(-3*x/(r^3)); vy1=vy+h*(-3*y/(r^3)); h=0.01; h=0.001 h=0.0005 优点：思路简单，适合于快速编程 缺点：只有一阶精度（局部截断误差），当步数增多时，误差会因积 累而越来越大 那么该如何改进？？？ 1、向后欧拉法 对于刚性微分方程（吸引解被附近快速改变的解所包围），用Euler 法，误差很大，改用向后Euler 法，误差较小，并且采用较大的步长， 提高效率 y xk +1 = y xk +hf(xk +1 ，y xk +1 ) y(xo )=yo 用某种办法使得斜率在右边点上是正确的。 但是这种方法并不能直接求得y xk +1 的迭代公式， 对于单体，xk +1 = xk + hvxk +1