\documentclass[minion]{homework}
\newcommand{\Reals}{\mathbb{R}}
\doclabel{Math F615: Homework 3}
\docdate{Due: February 6, 2019}
\usepackage{graphicx}

\newcommand{\bfx}{\mathbf{x}}
\newcommand{\bfv}{\mathbf{v}}

\begin{document}

\begin{problems}

\problem Recall the two step Adams-Bashforth method:
\[
u_{i+1} = u_i + \frac{h}{2}(3f_i-f_{i-1}).
\]
\begin{subproblems}
	\item Write down the stability polynomial $p(\rho)$ for this method applied to the problem $u'=\lambda u$.
	\item The equation $p(\rho)=0$ will involve the expresion $\lambda h$.  Solve for $\lambda h$ to write
	\[
       \lambda h = f(\rho)
	\]
	for some function $f$.
	\item Numerically determine values of $f(\rho)$ where $\rho$ lives on the unit circle of complex numbers.  These will 
	generate values of
	$\lambda h$ where the associated root of the characteristic
	polynomial has size one, and is therefore potentially on the 
	boundary of the stability region. 
	Generate a plot of the values of $f(\rho)$ as $\rho$
	varies around the unit circle to see the boundary of the 
	absolute stability region.
\end{subproblems}

\problem Recall that a linear multistep method (LMM) has the form
\[
\alpha_{k} u_{k+n}+\cdots +\alpha_1 u_{1+n} + \alpha_0 u_n = h(
\beta_{k} f_{k+n}+\cdots +\beta_1 f_{1+n} + \beta_0 f_n
)
\]
\begin{subproblems}
\item Show that the method is consistent if and only if
\begin{equation}
\begin{aligned}
\alpha_{k}+\cdots+\alpha_0 &= 0\\
k\alpha_{k}+\cdots+ 1\alpha_1+0\alpha_0 &= \beta_k+\cdots+\beta_0
\end{aligned}
\end{equation}
\item Use the previous result to show that every consistent one-step LMM
is zero stable.
\end{subproblems}

\problem Implement Euler's method and the Runge-Kutta RK4 method 
described in Table 1.3 for \textbf{vector} valued ODEs (i.e, systems).
Test your work against the IVP $\mathbf x' = A\mathbf x$ where
\[
A = \begin{pmatrix} 0 & -1 \\ 1 & 0 \end{pmatrix}
\]
with initial condition $\mathbf x_0 = (1,0)$.  Part of the exercise
is to compute the true solution of this system.  Hint: convert into
a second-order scalar ODE. \textbf{Verify that your code has the theoretical
order of convergence.}

\problem You are going to solve the 2-body problem of gravitation for the Earth-moon system.  Give two bodies with masses $m_1$ and $m_2$
at positions $x_1$ and $x_2$ in cartesian coordinates,
the force of body 2 on body 1 is
\[
F_{21} = G m_1 m_2 \frac{x_2-x_1}{|x_2-x_1|^3}
\]
where $G$ is the gravitational constant.  The force of body 1 on body 2
is the same, with the numbers 1 and 2 interchanged.

Newton's Laws then read
\begin{equation}
\begin{aligned}
\label{eq:2body}
m_1 \bfx_1'' &= F_{21}\\
m_2 \bfx_2'' &= F_{12}.
\end{aligned}
\end{equation}

We have the following physical constants:
\begin{align*}
m_{\text{Earth}} = 5.972 \times 10^{24} \mathrm{kg}\\
m_{\text{moon}} =  7.342 \times 10^{22} \mathrm{kg}\\
G = 6.67408 \times 10^{-11} \frac{ \mathrm{m}^3 }{ \mathrm{kg} \, \mathrm{s}^2}
\end{align*}

\begin{subproblems}
\item Convert system \eqref{eq:2body} into a first order system by
introducing the variables $\bfv_1=\bfx_1'$ and $\bfv_2=\bfx_2'$.
\item Writing $\bfx_i = [x_i,y_i]$ and $\bfv_i=[v_i,w_i]$, 
we will keep track of all of the scalar system variables
in vectors
\[
[x_1,y_1,x_2,y_2,v_1,w_1,v_2,w_2]^T.
\]
Using this convention, write a right-hand side 
function \texttt{z=earthmoon(t,u)} for the system you wrote 
down in part 1. 
\item Suppose at time $t=0$ the Earth is stationary and located
at the origin, the moon has position $x=3.565 \times 10^8$
and $y=0$ meters with velocity $1.09 \times 10^3 m/s$ in
the positive $y$ direction.

Use these initial conditions, and each of the solvers from problem 2, to generate an approximate solutions over a 40 day time interval (convert to seconds!) using daily (M=40) and then hourly (M=$40\times24$) time steps. That is, you are generating four different approximate
solutions.  Make basic plots of the computed trajectories.
 Describe in a few words what you see, and how these results relate to the local truncation errors of the schemes.

\item An important property of an isolated physical system is that
its energy is conserved.  For the 2-body problem, the energy
is 
\[
E(t) = \frac{m_1}{2}|\bfx_1'|^2 + \frac{m_2}{2}|\bfx_2'|^2 + U(\bfx_1,\bfx_2)
\]
where $U$ is the potential energy
\[
U(\bfx_1,\bfx_2) = -\frac{Gm_1m_2}{|\bfx_1-\bfx_2|}.
\]
If you have not seen this before, you should compute 
the derivative $E'(t)$ and show, using the differential equation, 
that it is zero (and hence $E(t)$ is constant).

Start by computing the value $E(0)$ exactly from the initial conditions. Then compute, for each of your four approximate
solutions, an array of scalar energy values E(t) computed at the 
$M+1$ sample times.  Plot the four energy curves from the four runs in one figure. Explain and comment. Describe an ``energy error norm'' which is small if the solution is of high accuracy, and report the values for the four runs.

\end{subproblems}

\problem Text: 1.10 a-d.  Note that your text uses different notation from what we've used in class.  My $h$ is your text's $k$. My $u_n$ is your text's $y_n$.
You may find it helpful to show that if $u'=f(t,u)$ then
\[
u''(t) = f_t(t,u) + f_u(t,u)f(t,u).
\]

\end{problems}
\end{document}