2
2
以万有引力的固定不动的施力质点
0
m
所在位置为坐标原点
O
, 建立直角
坐标系
Oxy
，质点的运动微分方程为
0
3
Gmm
mr r
r
= -
v v
&&
，分量方程为：
0 0
2 2 2 2
2 2 2 2
( ) ( )
Gm Gm
x y
x y
x y x y
x y x y
= - = -
+ +
+ +
&& &&
，
(1-4)
这两个方程都是二阶常微分方程，定义解矢量为
y
，令
(1) (2) (3) (4)y x y x y y y y= = = =
& &
， ， ，
(1-5)
可将方程组(1-4)化为：
( )
( )
0
3
2 2
2
0
3
2 2
2
(1)
(1) (2)
(2)
(1) (3)
(3)
(3) (4)
(4)
(1) (3)
Gm y
dy dy
y
dt dt
y y
Gm y
dy dy
y
dt dt
y y
×
= = -
+
×
= = -
+
，
，
(1-6)
(1)编写微分方程组函数文件 yxlcfun.m：
function ydot=yxlcfun(t,y,flag,p) %函数首行,p 为参量 Gm
0
ydot=[y(2); p*y(1)/sqrt(y(1).^2+y(3).^2).^3;
y(4); p*y(3)/sqrt(y(1).^2+y(3).^2).^3]; %建立微分方程组
(2)解微分方程的主程序 yxlc.m：
p=-1; %取 Gm
0
=1
y0=[-10 0.2 6 0.2;-25 0.5 5 0;-25 0.8 6 0]; %三组不同初始条件
plot(0,0, '*r') %画出 O 点
for i=1:3 %分别以不同初始条件解 3 次方程
[t,y]=ode45('yxlcfun',[0:0.1:300],y0(i,:),[ ],p);
hold on
axis([-25 25 -20 20]); %指定坐标范围
comet(y(:,1),y(:,3)) %绘出质点运动轨迹(x,y)
end %结束循环
解出的结果如图 1-7 所示：
剩余31页未读，继续阅读
