近日点进动与升交点西移Matlab数值计算程序
MATLAB科学计算软件入门教程
MATLAB科学计算软件入门教程第一章：MATLAB基础知识MATLAB是一种专业的科学计算软件，具有强大的数学计算和数据分析能力。
在使用MATLAB进行科学计算前，我们需要先了解一些基本知识。
1.1 MATLAB界面打开MATLAB后，我们会看到一个主界面。
主界面中有命令窗口、当前文件夹窗口、工作空间窗口和编辑器窗口等基本功能区域。
1.2 MATLAB变量和数据类型MATLAB中的变量可以用来存储各种类型的数据，如数字、字符串、矩阵等。
常见的数据类型包括：double（双精度浮点数）、char（字符）、logical（逻辑值）等。
1.3 MATLAB基本操作在MATLAB中，可以使用基本的数学运算符进行加、减、乘、除等计算操作。
另外，还可以通过内置函数实现更复杂的数学运算。
例如，sin函数可以计算正弦值，sum函数可以计算矩阵元素的和等。
第二章：MATLAB矩阵和向量操作2.1 创建矩阵和向量在MATLAB中，可以使用方括号来创建矩阵和向量。
例如，使用[1，2；3，4]可以创建一个2x2的矩阵。
2.2 矩阵和向量的加减乘除运算MATLAB提供了丰富的矩阵和向量运算函数，可以进行加法、减法、乘法、除法等运算操作。
例如，可以使用矩阵相乘函数*来计算矩阵的乘法。
2.3 矩阵和向量的索引和切片在MATLAB中，可以使用索引和切片操作来获取矩阵和向量中的特定元素或子集。
例如，使用矩阵名加上行和列的索引可以获取矩阵中指定位置的元素。
第三章：MATLAB数据可视化3.1 绘制二维图形MATLAB提供了丰富的绘图函数，可以绘制二维曲线、散点图、柱状图、等高线图等。
例如，可以使用plot函数来绘制二维曲线。
3.2 绘制三维图形MATLAB还可以绘制三维图形，如三维曲线、三维散点图、三维曲面等。
例如，可以使用plot3函数来绘制三维曲线。
3.3 图像处理与显示MATLAB提供了图像处理和显示的函数，可以加载、编辑和保存图像。
Matlab中的数值计算方法简介
Matlab中的数值计算方法简介引言：Matlab是一种强大的数值计算软件，广泛应用于工程、科学、金融等领域。
它拥有丰富的数值计算方法库，可以帮助研究者和工程师解决各种数值计算问题。
本文将简要介绍几种常见的数值计算方法，并说明它们在Matlab中的实现和应用。
一、插值法插值法是一种通过已知数据点之间的插值，估计未知数据点的数值的方法。
常见的插值方法包括线性插值、拉格朗日插值和样条插值。
在Matlab中，我们可以使用interp1函数进行插值计算。
该函数可以根据给定的数据点，计算出在指定位置的插值结果。
我们可以通过设置插值的方法和插值节点的数目来调整插值的精度与计算效率。
二、数值积分数值积分是一种通过近似求解定积分的方法。
在Matlab中，我们可以使用quad和quadl函数进行数值积分。
这些函数可以自动选择合适的数值积分方法，并提供了较高的精度和计算效率。
我们只需提供被积函数和积分区间，即可获得近似的积分结果。
对于一些特殊形式的积分，如复杂函数或无穷积分，Matlab还提供了相应的函数供我们使用。
三、线性方程组求解线性方程组的求解是数值计算中的一个重要问题。
在实际应用中，我们经常会遇到大规模线性方程组的求解问题。
在Matlab中，我们可以使用矩阵运算功能和线性方程组求解函数来解决这类问题。
Matlab提供了一系列的求解函数，包括直接法和迭代法。
其中，直接法适用于小规模线性方程组，迭代法则适用于大规模线性方程组。
我们可以根据具体情况选择合适的方法和函数来求解线性方程组。
四、微分方程求解微分方程是许多科学和工程问题的数学模型，求解微分方程是数值计算中的常见任务。
在Matlab中，我们可以使用ode45函数来求解常微分方程的初值问题。
该函数采用龙格-库塔方法，对微分方程进行数值积分，并给出近似的解析结果。
对于偏微分方程和其他更复杂的微分方程问题，Matlab还提供了更多的求解函数和工具箱供我们使用。
五、最优化问题求解最优化问题是指在特定约束条件下，求解给定目标函数的最大值或最小值的问题。
实验一MATLAB基本操作及运算
实验一MATLAB基本操作及运算MATLAB是一种强大的数值计算和数据可视化工具，广泛应用于科学研究、工程设计、数据分析等领域。
本文将介绍MATLAB的基本操作和运算。
首先，我们需要了解MATLAB中的基本数据类型，包括数值型、字符型和逻辑型。
数值型可以是整数、实数、复数等；字符型用单引号或双引号包围字符；逻辑型用true和false表示。
MATLAB提供了各种数学运算函数，包括四则运算、三角函数、指数函数等。
例如，加法可以使用加号（+），减法可以使用减号（-），乘法可以使用乘号（*），除法可以使用除号（/）。
三角函数可以使用sin、cos、tan等函数，指数函数可以使用exp函数。
此外，还可以使用log 函数进行对数运算。
MATLAB还可以进行矩阵运算。
矩阵可以使用方括号（[]）表示，每一行用分号（;）分隔。
可以使用矩阵乘法运算符（*）进行矩阵相乘，使用点乘运算符（.）进行矩阵对应元素的运算。
矩阵还可以进行转置、逆运算等。
除了基本运算，MATLAB还提供了各种其他功能。
例如，可以使用plot函数进行数据可视化，使用subplot函数绘制多个图形。
可以使用for循环和while循环进行循环操作，使用if语句进行条件判断。
MATLAB还可以进行文件读写操作。
可以使用load函数从文件中加载数据，使用save函数将数据保存到文件中。
可以使用fopen函数打开文件，使用fclose函数关闭文件。
可以使用fprintf函数写入文本文件，使用fscanf函数读取文本文件。
还可以使用imread函数读取图像文件，使用imwrite函数保存图像文件。
MATLAB还具备向量化的能力。
向量化是指使用矩阵代替循环进行计算，能够提高代码的执行效率。
例如，可以使用点乘运算符（.）对矩阵的每个元素进行计算，而不是使用循环逐个计算。
使用向量化的方法，可以更加简洁地编写代码。
在MATLAB中还有很多强大的功能等待探索，例如符号计算、模拟仿真、深度学习等。
lambert轨道问题求解代码matlab
Lambert轨道问题是由瑞士数学家约翰·海因里希·兰伯特于1761年提出的一个数学问题，用于求解两个给定位置的天体之间的转移轨道。
这个问题在航天领域中具有重要意义，例如在飞行器的轨道设计、太空探测器的轨道调整等方面都有广泛的应用。
为了解决Lambert轨道问题，可以通过编写MATLAB代码来进行求解。
在编写MATLAB代码进行Lambert轨道问题求解时，主要需要考虑以下几个方面：1. 问题的数学建模：首先需要将Lambert轨道问题进行数学建模，将问题描述为数学表达式。
根据两个给定的位置点和转移时间，可以得到一个由开普勒方程定义的转移轨道。
通过求解这个开普勒方程，可以得到相应的轨道参数。
2. MATLAB代码实现：根据数学建模的结果，可以编写MATLAB代码来实现Lambert轨道问题的求解。
一般可以使用数值计算的方法，例如牛顿-拉夫逊法或者二分法等来求解开普勒方程，得到轨道参数。
3. 算法优化：在编写MATLAB代码的过程中，需要考虑算法的效率和准确性。
可以通过改进算法的迭代次数、初始值选择等方式来提高算法的求解速度和精度。
4. 代码测试与验证：编写MATLAB代码后，需要进行代码的测试与验证，验证代码的准确性和可靠性。
可以通过一些已知的实际案例进行验证，例如国际空间站的航天器轨道设计等。
编写MATLAB代码来求解Lambert轨道问题是一个复杂且具有挑战性的任务，需要对数学建模、数值计算以及编程技术具有较高的要求。
通过不断的学习和实践，可以不断提高自己的能力，进而解决实际的工程问题。
Lambert轨道问题是航天领域中一个重要的数学问题，由瑞士数学家约翰·海因里希·兰伯特于18世纪提出，用于求解两个给定位置的天体之间的转移轨道。
对于航天器的轨道设计、太空探测器的轨道调整等方面具有广泛的应用。
为了解决Lambert轨道问题，可以借助MATLAB等数学建模和编程工具来进行求解。
matlab点的轨迹 -回复
matlab点的轨迹-回复如何在Matlab中绘制点的轨迹在Matlab中，我们可以使用绘图函数来绘制点的轨迹。
点的轨迹是点在一段时间内运动形成的路径。
这对于分析和可视化运动的模式非常有用。
在本文中，我们将学习如何使用Matlab绘制点的轨迹。
首先，我们需要定义点的初始位置和速度。
假设我们的点的初始位置为(0, 0)，速度为(1, 1)。
我们可以使用Matlab的矩阵来表示点的位置和速度。
下面是一个例子：matlab定义点的初始位置和速度position = [0, 0];velocity = [1, 1];接下来，我们需要确定点的运动时间和步长。
时间用来模拟点的运动过程，步长用来控制时间的增量。
我们可以通过设置时间步长来调整点的轨迹的精度。
例如，如果我们将时间步长设置为0.1秒，那么点的位置将每0.1秒更新一次。
下面是一个例子：matlab定义时间和步长time = 0:0.1:10; 从0到10，步长为0.1然后，我们可以使用循环来模拟点的运动过程。
在每个时间步长，我们可以根据点的初始位置和速度来更新点的位置。
这样，我们就可以得到点在每个时间步长的位置。
下面是一个例子：matlab模拟点的运动过程for i = 1:length(time)更新点的位置position(i+1, :) = position(i, :) + velocity;end最后，我们可以使用Matlab的绘图函数来绘制点的轨迹。
我们可以使用plot函数将点的位置连接起来，从而形成轨迹。
下面是一个例子：matlab绘制点的轨迹plot(position(:, 1), position(:, 2));xlabel('X轴');ylabel('Y轴');title('点的轨迹');运行上述代码，我们就可以在Matlab中绘制出点的轨迹。
通过上面的步骤，我们可以使用Matlab绘制出点的轨迹。
学习使用MATLAB进行科学计算
学习使用MATLAB进行科学计算MATLAB是一种强大的科学计算软件，被广泛应用于科学研究和工程领域。
作为一名科学家或工程师，学习并熟练使用MATLAB可以帮助我们更高效地进行数据分析、模拟和可视化等工作。
在本文中，我将介绍一些MATLAB的基本概念和常用功能，帮助初学者快速入门。
首先，让我们从最基本的操作开始。
打开MATLAB后，你将看到一个命令窗口。
在这里，你可以输入各种命令，并立即获得结果。
试着输入一个简单的数学运算，如2+2，然后按下回车。
你将看到结果显示在命令窗口中。
除了进行简单的数学运算外，MATLAB还具备处理矩阵和向量的能力。
在MATLAB中，矩阵和向量可以用来表示和处理大量的数据。
你可以使用方括号来创建矩阵和向量，比如：A = [1, 2, 3; 4, 5, 6; 7, 8, 9];这是一个3x3的矩阵A，其中包含了一些数字。
你可以通过在命令窗口中输入矩阵的名称来查看其内容：AMATLAB将会显示矩阵A的内容。
你还可以通过索引来访问矩阵中的元素。
比如，通过输入A(1,2)，你将获得矩阵A中第一行第二列的元素。
在进行科学计算时，我们经常需要对数据进行统计分析。
MATLAB提供了许多有用的函数来执行这些任务。
例如，你可以使用mean函数来计算矩阵或向量的平均值，使用std函数来计算标准差，并使用hist函数来生成直方图。
试试看吧：data = [1, 2, 3, 4, 5];mean(data)std(data)hist(data)MATLAB将计算出数据的平均值、标准差，并生成直方图。
此外，MATLAB还具备强大的可视化功能，使我们能够更好地理解和呈现数据。
使用plot函数，你可以绘制函数曲线或数据点，使用imshow函数，你可以显示图像，使用surf函数，你可以创建3D曲面。
让我们看一些例子：x = [0:0.1:2*pi];y = sin(x);plot(x, y)在这个例子中，我们首先创建一个包含0到2π之间一系列值的向量x。
matlab数值求解常微分方程快速方法
MATLAB是一种用于科学计算和工程应用的高级编程语言和交互式环境。
它在数学建模、模拟和分析等方面有着广泛的应用。
在MATLAB 中，常微分方程的数值求解是一个常见的应用场景。
在实际工程问题中，通常需要对常微分方程进行数值求解来模拟系统的动态行为。
本文将介绍MATLAB中对常微分方程进行数值求解的快速方法。
1. 基本概念在MATLAB中，可以使用ode45函数来对常微分方程进行数值求解。
ode45是一种常用的Runge-Kutta法，它可以自适应地选取步长，并且具有较高的数值精度。
使用ode45函数可以方便地对各种类型的常微分方程进行求解，包括一阶、高阶、常系数和变系数的微分方程。
2. 函数调用要使用ode45函数进行常微分方程的数值求解，需要按照以下格式进行函数调用：[t, y] = ode45(odefun, tspan, y0)其中，odefun表示用于描述微分方程的函数，tspan表示求解的时间跨度，y0表示初值条件，t和y分别表示求解得到的时间序列和对应的解向量。
3. 示例演示为了更好地理解如何使用ode45函数进行常微分方程的数值求解，下面我们以一个具体的例子来进行演示。
考虑如下的一阶常微分方程：dy/dt = -2*y其中，y(0) = 1。
我们可以编写一个描述微分方程的函数odefun：function dydt = odefun(t, y)dydt = -2*y;按照上述的函数调用格式，使用ode45函数进行求解：tspan = [0 10];y0 = 1;[t, y] = ode45(odefun, tspan, y0);绘制出解曲线：plot(t, y);4. 高级用法除了基本的函数调用方式外，MATLAB中还提供了更多高级的方法来对常微分方程进行数值求解。
可以通过设定选项参数来控制数值求解的精度和稳定性，并且还可以对刚性微分方程进行求解。
5. 性能优化在实际工程应用中，常常需要对大规模的常微分方程进行数值求解。
MATLAB科学计算使用教程
MATLAB科学计算使用教程第一章：MATLAB入门MATLAB（Matrix Laboratory）是一种用于科学计算和技术计算的强大软件工具。
本章将介绍如何安装MATLAB，并进行初步的配置和设置。
同时还将介绍MATLAB的基本操作，如变量的定义和使用、基本数学运算、矩阵的创建和操作等。
第二章：数据处理与分析本章将介绍MATLAB在数据处理与分析方面的强大功能。
涵盖了数据的导入和导出、数据预处理、常用统计分析方法、数据可视化等内容。
具体包括：使用MATLAB读取和写入常见数据格式，例如Excel、CSV、TXT等；数据清洗和处理的常用方法，如缺失值处理、异常值检测等；常用统计分析方法的实现，如假设检验、方差分析等；数据可视化方法和技巧，如统计图表的绘制和优化。
第三章：信号处理与滤波本章将介绍MATLAB在信号处理和滤波方面的应用。
包括信号生成和操作、常用信号处理方法、数字滤波器设计等内容。
具体包括：使用MATLAB生成各类常用信号，如正弦信号、方波信号等；对信号进行时域和频域的分析；常用的信号处理方法，如时域滤波、频域滤波、小波变换等；数字滤波器的设计和实现。
第四章：图像处理与计算机视觉本章将介绍MATLAB在图像处理和计算机视觉方面的应用。
涵盖图像读取和显示、图像处理和增强、计算机视觉算法等内容。
具体包括：使用MATLAB读取和显示图像文件，如JPEG、PNG 等；图像的基本处理和增强，如灰度变换、滤波器应用、颜色空间转换等；图像分割和特征提取方法；计算机视觉算法的实现，如目标检测、图像识别等。
第五章：数学建模与优化本章将介绍MATLAB在数学建模与优化方面的应用。
包括数学建模的基本方法、优化问题和求解方法等。
具体包括：数学建模的基本步骤和实现思路，如问题分析、建立数学模型等；常见数学建模问题的解决方法，如线性规划、非线性规划等；优化问题的MATLAB求解方法，如线性规划求解器、遗传算法优化等。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
近日点进动与升交点西移Matlab数值计算程序这里只提供行星摄动轨道的数值计算程序。
分为方程部分m文件和画图部分m文件。
八大行星的轨道方程程序都可以根据模板(moban)程序调整得到。
保存名就是moban.m。
%后的绿色字体是注释部分不执行。
...是连接符合，连接上下行。
function ydot=moban(t,y)%模拟mb星相对太阳轨迹。
y(1)是Φ，y(2)是V，y(3)是θ,y(4)是R，y(5)是σ，y(6)是ψ。
%Gm(i),Sp(i),w1(i)中i=1水2金3地4火5木6土7天8海。
忽略小行星和柯依伯带天体，考虑卫星质量。
%木星卫星总Gm=2.66e13，土星卫星1.031e13,天王卫星6.1e11，海王卫星1.44e12。
%q5为σp,op为φp,Psi0为模拟行星黄经初值-其他行星黄经初值,Psi为Θp。
%初值举例：[op0,V0,pi/2,a,q5,0,op0]%梁彬彬编写于2011年4月5日GS=1.3271244e20;Gm=[2.203208e13,3.24858597e14,4.0350324e14,4.282838e13,1.267393645e17,3.795089544e16,5.7951591e15,6.83796713e15]; Sp=[5.790923e10,1.08209475e11,1.4959787e11,2.2794382e11,7.78340821e11,1.42666642e12,2.87065819e12,4.49839644e12]; q5=[0.12225995,0.059248274,2.67e-7,0.032283205,0.022766022,0.043388743,0.013485074,0.030890865];w1=[8.266758e-7,3.236395e-7,1.990987e-7,1.058577e-7,1.6784e-8,6.761258e-9,2.369788e-9,1.208208e-9];%行星公转角速度op0=[1.3518936,2.2968964,1.7966015,5.8652901,0.25706047,1.6161553,2.983715,0.78478315];op=op0+w1*t;w5=[6.856e-13,1.532e-12,1.334e-12,1.567e-12,-1.132e-12,1.597e-12,-2.345e-13,6.546e-14];%升交点西移角速度Psi0=-[0.84530995,1.3383157,0,0.86497713,1.7536005,1.9837835,1.291839,2.3000686];Psi=y(5)+Psi0+w5*t;ip=Sp.*(cos(op).*(cos(Psi).*cos(y(1))-sin(Psi).*sin(y(1)).*cos(y(5)))+...sin(op).*(cos(q5).*(sin(Psi).*cos(y(1))+cos(Psi).*cos(y(5)).*sin(y(1)))+sin(y(1)).*sin(y(5)).*sin(q5)));jp=-Sp.*(cos(op).*(cos(Psi).*sin(y(1))+sin(Psi).*cos(y(1)).*cos(y(5)))+...sin(op).*(cos(q5).*(sin(Psi).*sin(y(1))-cos(Psi).*cos(y(5)).*cos(y(1)))-cos(y(1)).*sin(y(5)).*sin(q5)));kp=Sp.*(sin(op).*(cos(y(5)).*sin(q5)-cos(Psi).*cos(q5).*sin(y(5)))+sin(Psi).*cos(op).*sin(y(5)));Dp2=Sp.^2-2*y(4).*ip+y(4)^2;R3=1./Dp2.^1.5-1./Sp.^3;f156=sum(Gm.*kp.*R3/y(2)/sin(y(3)));f2=sum(Gm.*(ip.*cos(y(3))+jp.*sin(y(3))).*R3-Gm*y(4)*cos(y(3))./Dp2.^1.5);f3=sum(Gm.*(ip.*sin(y(3))-jp.*cos(y(3))).*R3-Gm*y(4)*sin(y(3))./Dp2.^1.5)/y(2);ydot=[y(2)*sin(y(3))/y(4)-f156*sin(y(1))/tan(y(5));-GS/y(4)^2*cos(y(3))+f2;(GS/y(4)^2/y(2)-y(2)/y(4))*sin(y(3))-f3;y(2)*cos(y(3));f156*cos(y(1));f156*sin(y(1))/sin(y(5));y(2)*sin(y(3))/y(4)];各行星摄动方程形式是一样的，各种初值不一样。
第一块包括摄动体初值，太阳引力常数。
需把GS换成太阳与模拟行星的总引力常数，并把Gm、Sp、q5、w1、op0、w5各矩阵中模拟行星的那列数据去掉，再把Psi0模拟的行星那列移数据提出来就行了。
比如水星，改变矩阵中的第一列。
保存名为Mercury.m。
受摄相对轨道微分方程组(2.12)是六列表达式，由于(2.11)下面说明文字的原因ydot才成为七列的。
function ydot=Mercury(t,y)%模拟水星相对太阳轨迹。
y(1)是Φ，y(2)是V，y(3)是θ,y(4)是R，y(5)是σ，y(6)是ψ。
%Gm(i),Sp(i),w1(i)中i=1金2地3火4木5土6天7海。
忽略小行星和柯依伯带天体，考虑卫星质量。
%木星卫星总Gm=2.66e13，土星卫星1.031e13,天王卫星6.1e11，海王卫星1.44e12。
%q5为σp,op为φp,Psi0为水星黄经初值-行星黄经初值,Psi为Θp。
%初值举例：[1.3518936,3.88583e4,pi/2,6.98174e10,0.12225995,0,1.3518936]%梁彬彬编写于2011年4月5日GSM=1.3271246413e20;Gm=[3.24858597e14,4.0350324e14,4.282838e13,1.267393645e17,3.795089544e16,5.7951591e15,6.83796713e15];Sp=[1.08209475e11,1.4959787e11,2.2794382e11,7.78340821e11,1.42666642e12,2.87065819e12,4.49839644e12];q5=[0.059248274,2.67e-7,0.032283205,0.022766022,0.043388743,0.013485074,0.030890865];w1=[3.236395e-7,1.990987e-7,1.058577e-7,1.6784e-8,6.761258e-9,2.369788e-9,1.208208e-9];%行星公转角速度op0=[2.2968964,1.7966015,5.8652901,0.25706047,1.6161553,2.983715,0.78478315];op=op0+w1*t;w5=[1.532e-12,1.334e-12,1.567e-12,-1.132e-12,1.597e-12,-2.345e-13,6.546e-14];%升交点西移角速度Psi0=0.84530995-[1.3383157,0,0.86497713,1.7536005,1.9837835,1.291839,2.3000686];Psi=y(5)+Psi0+w5*t;ip=Sp.*(cos(op).*(cos(Psi).*cos(y(1))-sin(Psi).*sin(y(1)).*cos(y(5)))+...sin(op).*(cos(q5).*(sin(Psi).*cos(y(1))+cos(Psi).*cos(y(5)).*sin(y(1)))+sin(y(1)).*sin(y(5)).*sin(q5))); jp=-Sp.*(cos(op).*(cos(Psi).*sin(y(1))+sin(Psi).*cos(y(1)).*cos(y(5)))+...sin(op).*(cos(q5).*(sin(Psi).*sin(y(1))-cos(Psi).*cos(y(5)).*cos(y(1)))-cos(y(1)).*sin(y(5)).*sin(q5))); kp=Sp.*(sin(op).*(cos(y(5)).*sin(q5)-cos(Psi).*cos(q5).*sin(y(5)))+sin(Psi).*cos(op).*sin(y(5)));Dp2=Sp.^2-2*y(4).*ip+y(4)^2;R3=1./Dp2.^1.5-1./Sp.^3;f156=sum(Gm.*kp.*R3/y(2)/sin(y(3)));f2=sum(Gm.*(ip.*cos(y(3))+jp.*sin(y(3))).*R3-Gm*y(4)*cos(y(3))./Dp2.^1.5);f3=sum(Gm.*(ip.*sin(y(3))-jp.*cos(y(3))).*R3-Gm*y(4)*sin(y(3))./Dp2.^1.5)/y(2);ydot=[y(2)*sin(y(3))/y(4)-f156*sin(y(1))/tan(y(5));-GSM/y(4)^2*cos(y(3))+f2;(GSM/y(4)^2/y(2)-y(2)/y(4))*sin(y(3))-f3;y(2)*cos(y(3));f156*cos(y(1));f156*sin(y(1))/sin(y(5));y(2)*sin(y(3))/y(4)];轨道方程的数值计算采用高效的Adams-Bashforth-Moulton PECE法[ode113]计算。
而四阶五级Runge-Kutta-Felhberg法[ode45]效率较低，容易卡死。
计算出来要画图，为了找近日点精确一些，要加入插值程序。
有时还需要对比同一星体几种引力理论的不同处理结果，因此编写了以下程序。
保存名为P_orbit.m。
function P_orbit(fun,year,y0,P,varargin)%fun是轨道方程文件名；year是模拟儒略年份范围；y0是初值；P选择作图模式；varargin是对比轨道方程文件名。
%P=1时只作升交点西移和近日点进动图。