Question

1. (Coding problem) Consider the satellite equations GMx (x2 + y2)3/2' GMY (x2 + y2)3/2" with the initial condition x(0) = 0, x'(0) = v, y(0) = 20000m + R, y'(0) = 0, where G = 6.67408 x 10-11mºkg-15-2, M = 5.974 x 1024 kg, R= 6378.1 x 10°m. The initial condition is such that at time t = 0, the satellite is 20km above our head with a velocity v in the direction perpendicular to vertical. If the velocity v is too small, the satellite's orbit will crash onto the Earth (when (x2 + y2) falls below the Earth’s radius R). This problem is to numerically determine the minimal velocity v so that the satellite does not crash onto the Earth. (1) First, choose a large value of v = 10,000m/s, and use the Runge-Kutta-4 method to compute the solution to time t = 10 hours. Verity in this case that the distance of the satellite from the center of the Earth V (x2 + y2) is always larger than R. Choose your time step h carefully so that the solution at t = 10 hours has relative accuracy less than 10-6. Then plot your solution in the (x,y) plane. (2) Now, gradually decrease your velocity v, and for each v, check its full orbit to decide if it crashes onto the Earth or not. From this, determine the minimum velocity Vmin so that the satellite does not crash onto the Earth.

Answer 1

(1) The system can be first converted into a system of four ordinary differential equations and then can be solved using ode45 (Runke-Kutta 4th order) of MATLAB by the following code:

G=6.67408*10^(-11);
M=5.974*10^(24);
R=6378.1*10^3;
v=10000;
f = @(t,x) [x(2);(-G*M*x(1))/((x(1))^2+(x(3))^2)^(3/2);x(4);(-G*M*x(1))/((x(3))^2+(x(3))^2)^(3/2)]; % The system %
[t,xa] = ode45(f,[0 10],[0 v 20000+R 0]);
distance=sqrt((xa(85,1))^2+(xa(85,3))^2) % distance of the satelite from the centre of earth %

plot(xa(:,1),xa(:,3))
xlabel('x(t)'), ylabel('y(t)')

You can see the distance we get is 6398880 which is larger than R, i.e., 6378100. The plot (x,y) is given below :

6.3981 X 106 6.3981 6.3981 6.3981 6.3981 y(t) 6.3981 6.3981 6.3981 6.3981 6.3981 0 1 2 3 4 6 7 8 9 10 5 X(t) 104

(b) To determine the v_{min}, one can change v in the above MATLAB code and check the distance for each of this v. Then see for which it is less than R.