Page 315 - Applied Numerical Methods Using MATLAB

304 ORDINARY DIFFERENTIAL EQUATIONS

function dx = df_sat(t,x)
global G Me Re
dx = zeros(size(x));
r = sqrt(sum(x(1:2).^2));
if r <= Re, return; end % when colliding against the earth surface
GMr3 = G*Me/r^3;
dx(1) = x(3); dx(2) = x(4); dx(3) = -GMr3*x(1); dx(4) = -GMr3*x(2);
%nm6p05.m to solve a nonlinear d.e. on the orbit of a satellite
clear, clf
global G Me Re
G = 6.67e-11; Me = 5.97e24; Re = 64e5;
f = ’df_sat’; ;
t0 = 0; T = 24*60*60; tf = T; N = 2000;
R = 4.223e7;
v20s = [3071 3500 2000];
for iter = 1:length(v20s)
x10 = R; x20 = 0; v10 = 0; v20 = v20s(iter);
x0 = [x10 x20 v10 v20]; tol = 1e-6;
[tR,xR] = ode_RK4(f,[t0 tf],x0,N);
[t45,x45] = ode45(????????????);
[t23s,x23s] = ode23s(f,[t0 tf],x0);
plot(xR(:,1),xR(:,2),’b’, x45(:,1),x45(:,2),’k.’, ????????????)
[t45,x45] = ode45(f,[t0 tf],x0,odeset(’RelTol’,tol));
[t23s,x23s] = ode23s(?????????????????????????????????);
plot(xR(:,1),xR(:,2),’b’, x45(:,1),x45(:,2),’k.’, ????????????)
end

7
(i) (x 10 ,x 20 ) = (4.223 × 10 , 0)[m] and (x 30 ,x 40 ) = (v 10 ,v 20 ) =
(0, 3071)[m/s].
7
(ii) (x 10 ,x 20 ) = (4.223 × 10 , 0)[m] and (x 30 ,x 40 ) = (v 10 ,v 20 ) =
(0, 3500)[m/s].
7
(iii) (x 10 ,x 20 ) = (4.223 × 10 , 0)[m] and (x 30 ,x 40 ) = (v 10 ,v 20 ) =
(0, 2000)[m/s].
Run the program and check if the plotting results are as depicted in
Fig. P6.5.
(b) In Fig. P6.5, we see that the “ode23s()” solution path differs from
the others for case (ii) and the “ode45()”and “ode23s()” paths differ
from the “ode_RK4()” path for case (iii). But, we do not know which
one is more accurate. In order to find which one is the closest to the
true solution, apply the two routines “ode45()”and “ode23s()” with
smaller relative error tolerance of tol = 1e-6 to find the paths for the
three cases. Which one do you think is the closest to the true solution
among the paths obtained in (a)?
(cf) The purpose of this problem is not to compare the several MATLAB routines,
but to warn the users of the danger of abusing them. With smaller number
of steps (N) (i.e., larger step size), the routine “ode_RK4()” will also deviate
much from the true solution. The MATLAB built-in routines have too many
good features to be mentioned here. Note that setting the parameters such as