Question

What's the f(t,y) in this question?

Euler's method requires f(t,y). What's the f(t,y) function in this assignment? q(t) and p(t) are vectors. How am I supposed to apply them in Euler's method?

One of science's great achievements was the discovery of Kepler's laws for planetary motion; in particular, that the planets follow closed elliptical orbits around the sun. In this assignment you will compare two different numerical methods for solving Kepler's problem for the motion of a simple solar system consisting of two planets (the two-body problem). For a system of two planets, we may assume one is fixed at the origin with the motion of the other planet being in a 2D plane. Let q(t) = 91 ), 22(t) /' p(t) = ( pi(t)) p2(t) ? be the position and momentum vectors of the moving planet, respectively. Kepler's laws give the following ordinary differential equation for q(t) and p(t): q'(t) = p(t), p'(t) = - (1) (91(t)2 +42(t)2)3/29(t). a) Write a code that implements Euler's method for this problem for time 0 <t<T = 200 and stepsize h = 0.0005. Use the initial conditions 1+e 91(0) = 1 – €, 92(0) = 0, P1(0) = 0, P2(0) = Vite e=0.6.

Answer 1

The Matlab scripts eulFwdSys.m, eulFwdSysSolver.m, dydtSys.m are posted at the end along with the generated plot. Please run the script eulFwdSysSolver.m to generate the desired results and plot.
Q' (t) = Pct) so .. (1) a;' = Ⓡp aź = 3 .. (2) and, P' (t) - 132 (t). or, b, c - a, B'= - = ...(4) ...4) So, the eqns (1), (2), (3) and (4) Forms the system Of first order differential equations az - Þe, 9210) = 1-e az' = pz , 9(0) = 0 ' = - 92 , B, (0) = Vite 19;2+2;2775' ? Vice where e=0.6

Thus, the Euler method is applied on the system of first order equations, wher ty (Ps, Pz, 91, 92, t ) = Þg tę (ba, pa, 93, 92, t) = þz tz (Ps, bz, 91,9,, t) = -as (9?+9,2,312 and, fa (Ps, bz,96,92, +) = - 92 (9: +972)302 therefore, FCB, ģit) to given as and, the dependent variable vector is I al 42 | = g in Matlab) where 9. = y(1), 92 = 4(2), by = 4(3), Đ, = Y(4), The Matlab scripts along with the generated plot is posted below.

eulFwdSys.m

function [t, y] = eulFwd (dydt, tSpan, initCond, timeStep)
% eulode: Euler ODE solver
% [t,y] = eulode (dydt, tSpan, initCond, timeStep):
% uses Euler's method to integrate an ODE

% input:
% dydt = name of the function M−file that evaluates the ODE
% tspan = [ti, tf] where ti and tf = initial and
% final values of the independent variable t
% initCond = initial value of dependent variables
% timeStep = step size of the independent variable

% output:
% t = vector of independent variable
% yp = vector of solution for dependent variable

tInitial = tSpan(1);
tFinal   = tSpan(2);
t = (tInitial: timeStep :tFinal)';

y(1, :) = initCond;
yp(1, :) = y(1, :);

for i = 1: length(t)-1   
    y(i + 1, :) = y(i, :) + dydt(t(i), y(i, :)) * timeStep;
    yp(i + 1, :) = y(i + 1, :);
end
end

eulFwdSysSolver.m

clear, close all; clc

tSpan = [0 200];
timeStep = 5e-4;
ec = 0.6;
y0 = [ 1-ec, 0, 0, sqrt((1+ec)/(1-ec))];

[t, y] = eulFwd(@dydtSys, tSpan, y0, timeStep);

figure
plot(t, y(:, 1), t, y(:, 2), t, y(:, 3), t, y(:, 4))
xlabel( '\bf t' )
ylabel( '\bf y' )
legend('q_1', 'q_2', 'p_1', 'p_2')
grid on
grid minor

dydtSys.m

function rhs = dydtSys(t, y)
% governing ODE dydt = rhs
% dydt = dydt(t, y)

% input:
% x = independent variable
% y = dependent variable

% output:
% dydt = rhs of the ODE

dydt1 = y(3);
dydt2 = y(4);
dydt3 = -y(1)/((y(1))^2+(y(2))^2)^(3/2);
dydt4 = -y(2)/((y(1))^2+(y(2))^2)^(3/2);

rhs = [dydt1, dydt2, dydt3, dydt4];

end
0.5 KM -0.5 – 20 40 60 100 140 160 180 200