Matlab实验报告开普勒方程

Matlab实验报告西安交通大学车辆01班教师：阮小娥10015005郝杰2011.06.22一、实验问题：开普勒方程近似解与方程求根天文学中，有一类著名的方程—开普勒方程,x=qsinx+a(0q1,a为常数)，开普勒方程是用来确定行星在其运动轨道上的位置的。如何求解方程并使其解达到一定的精度要求呢？开普勒方程近似解与方程求根。二、问题的分析、（1）理论知识：开普勒方程是一个超越方程，很难得出严格的分析解，但是，已经证明这个方程存在唯一解。如果已知某一作椭圆运动的天体的轨道要素，利用二体问题的关系式可以得到任意给定时刻t时的平近点角M，而后采用图解法、数值法或近似迭代法求解开普勒方程得出偏近点角E，再利用二体问题的其他积分而得到t时刻天体在轨道上的坐标和速度。对于抛物线轨道和双曲线轨道也有相应的开普勒方程。首先为了求解开普勒方程，通常将求函数f（x）的零点或方程f（x）=0的根的问题。其次，在平面坐标系中绘制函数f（x）的曲线，大致了解f（x）的零点情况和位置，如果有多个零点，则要分区逐个求解。二体问题运动方程的一个积分。它反映天体在其轨道上的位置与时间t的函数关系。对于椭圆轨道，开普勒方程可以表示为E－esinE＝M，式中E为偏近点角，M为平近点角，都是从椭圆轨道的近日点开始起算，沿逆时针方向为正，E和M都是确定天体在椭圆轨道上的运动和位置的基本量。如果定义天体在轨道上运动的平均角速度为n，天体过近日点的时刻为τ，则对任一给定时刻t，天体从近日点出发所走过的角度就是平近点角M＝n（t－τ）。这样，开普勒方程给出了天体在轨道上运动的位置与时间t的关系。三、实验目的：通过对开普勒方程求根问题的解决，熟悉非线性方程的几种数值方法，并进行试验。数学建模及程序设计及结果展示：Number1.用二分法求方程x=0.5sinx+1的近似根（误差小于0.00001）程序设计：f=inline('x-0.5*sin(x)-1');a=1;b=2;dlt=1.0e-5;k=1;whileabs(b-a)dltc=(a+b)/2;iff(c)==0break;elseiff(c)*f(b)0a=c;elseb=c;endfprintf('k=%d,x=%.5f\n',k,c);k=k+1;end实验结果：Number2：用“切线法”求方程x=0.5sinx+1的近似根（误差小于0.00001）.程序设计：f=inline('x-0.5*sin(x)-1');df=inline('1-0.5*cos(x)');d2f=inline('0.5*sin(x)');a=1;b=2;dlt=1.0e-5;iff(a)*d2f(a)0x0=a;elsex0=b;endm=min(abs(df(a)),abs(df(b)));k=0;whileabs(f(x0))m*dltk=k+1;x1=x0-f(x0)/df(x0);x0=x1;fprintf('k=%dx=%.5f\n',k,x0);end实验结果：Number3：求方程组的近似解.程序设计：functionf=group1(x)f=[sin(x(1))+x(2)+x(3)^2*exp(x(1))-4;x(1)+x(2)*x(3);x(1)*x(2)*x(3)+2];►[x,fval]=fsolve('group1',[1,1,1])x=1.4142-1.37011.0322fval=1.0e-012*0.11550.0007-0.0071实验结果：Number4求方程组的近似解.程序设计：functionf=group2(x)f=[9*x(2)^2-12*x(1)-54*x(2)+61;x(1)*x(2)-2*x(1)+1];[x,fval]=fsolve('group2',[0,0])x=1.09021.0828fval=1.0e-011*0.1044-0.0324[x,fval]=fsolve('group2',[-2,2])x=-1.56822.6377fval=1.0e-007*0.32180.0260实验结果：实验进一步拓展与收获：（1）：二分法一般地，对于函数f(x),如果存在实数c,当x=c时，若f(c)=0,那么把x=c叫做函数f(x)的零点。解方程即要求f(x)的所有零点。假定f(x)在区间（x，y）上连续先找到a、b属于区间（x，y），使f(a)，f(b)异号，说明在区间(a,b)内一定有零点，然后求f[(a+b)/2],现在假设f(a)0,f(b)0,ab①如果f[(a+b)/2]=0，该点就是零点，如果f[(a+b)/2]0,则在区间（(a+b)/2，b)内有零点，(a+b)/2=a，从①开始继续使用中点函数值判断。如果f[(a+b)/2]0，则在区间(a,(a+b)/2)内有零点，(a+b)/2=b，从①开始继续使用中点函数值判断。这样就可以不断接近零点。通过每次把f(x)的零点所在小区间收缩一半的方法，使区间的两个端点逐步迫近函数的零点，以求得零点的近似值，这种方法叫做二分法。从以上可以看出，每次运算后，区间长度减少一半，是线形收敛。另外，二分法不能计算复根和重根。（2）：一般迭代法非线性方程一般很难求得数值解，而且在实际应用中也没有必要求得精确的数值解，往往只要求得满足一定精度要求的近似解。常用的非线性方程求解方法主要有两种：搜索法（直接法）和迭代法。直接迭代法就最简单的迭代法。直接迭代法的解思路为使用某个固定公式反复校正根的近似值，使其逐步精确，直到得出满足精度要求的结果位置。非线性方程的一般形式为f(x)=0，用迭代法求解时，首先需要将一般形式的方程改写为以下形式x=g(x)（1）上式左右两端都含有未知的x，我们用一个估计值x0带入方程右端求出x，将求出的x写为x1，方程变为x1=g(x0)。再将x1带入右端，又可求得x2，这样可得xk+1=g(xk)（2）上式即迭代格式，由上式反复计算可得到一个数列x0，x1，x2，…，xk，…。如果此数列有极限，这个极限就是方程x=g(x)的根。所得数列的极限存在时，称迭代格式收敛，反之，则称迭代格式发散。此时无法通过迭代求解。如果两次迭代计算的偏差小于规定的允许误差ε，即满足收敛判据∣xk+1-xk∣ε（3）则终止计算。如果g(x)有连续的一阶导数g’(x)，若满足∣g’(x*)∣1，则对任意初值x0均收敛。g’(x)取值不同，迭代过程的收敛情况不同。直接迭代法的具体步骤：①给定初值x0、计算精度；②用迭代格式xk+1=g(xk)进行迭代计算；③判断迭代结果是否满足收敛判据，如果满足，终止计算并输出结果，否则返回步骤②。（3）：切线法：切线法又称为牛顿法,是一种一般情况下具有二阶收敛速度的非线性方程的数值解法.具体方法如下:设x*是方程f(x)=0的根,又x0为x*附近的一个值,将f(x)在x0附近做泰勒展开:f(x)=f(x0)+(x-x0)f'(x0)+1/2(x-x0)²f''(ξ)其中ξ在x和x0之间令x=x*,则:0=f(x*)=f(x0)+(x*-x0)f'(x0)+1/2(x*-x0)²f''(ξ)去掉x*-x0的二次项得到:f(x0)+x*f'(x0)-x0f'(x0)≈0即x*≈x0-f(x0)/f'(x0)令x1=x0-f(x0)/f'(x0)并由此构成一个递推式x[k+1]=x[k]-f(x[k])/f'(x[k])([]表示下标)可以证明,当f(x)∈C[a,b]且满足以下条件时,由以上递推式产生的序列最后收敛到f(x)=0在[a,b]上的唯一根(1)f(a)f(b)0(2)f'(x)≠0,(x∈[a,b])(3)f''(x)在[a,b]上恒正或恒负(4)初值x0∈[a,b]应满足f(x0)f''(x0)0计算实例:1.求解f(x)=x-cosx=0的实根由零点定理知f(x)=0在(0,π/2)内有实根f'(x)=1+sinx,由迭代公式有:x[n+1]=x[n]-(x[n]-cosx[n])/(1+sinx[n])取x0=π/4得到:x1=0.73936133x2=0.739085178x3=0.739085133x4=0.739085133所以x=0.739085133.....2.任意数开n次方为了说明的方便,在此就常见的开3次方作较详细的说明,对于其他的可以类比计算设x=³√A则x³=A所以x³-A=0采用递推公式x[n+1]=x[n]-(x[n]³-A)/(3x[n]²)([]表示下标)即可求出³√A的任意精度近似值.初值x[0]一般取与³√A接近的整数.举例求³√28,取x[0]=3,迭代结果如下:x[1]=3.037037037037037x[2]=3.036589037977101x[3]=3.036588971875664x[4]=3.036588971875663x[5]=3.036588971875663从上面可以看出,只要迭代4次即可求出15位精度的近似值五：实验心得与体会：Matlab实验工具真的是一种非常实用的工具，一般我们用人脑难以及计算的数据，他都可以帮我们解决掉。但是Matlab难在使用该工具之前，你要将问题的解决方法分析清楚，首先要确定问题的类型，然后根据已知的理论进行数学建模，建好摸之后，还要将数学模型套用Matlab语言将其编写出来，这样这一类问题便都将迎刃而解。所以，只要分析好一个问题，便是解决了一类的很多问题。实验当中用到了一般迭代法和二分法，学习了这两种问题的解决方法后，有利于我们以后解决其他类似的复杂问题。