Math 753/853 HW7 ODEs

Problem 1

Implement three Julia functions for numerical solution of 1st-order IVPs

\begin{align*} y' &= f(t,y) \\ y(t_0) & = y_0 \end{align*}

for forward Euler, Runge-Kutta 2nd order (midpoint method), and Runge-Kutta 4th order. The functions should take inputs f, y0, t0, Δt, N and returns $N+1$-dimensional vectors t and y, where t is the vector of discrete times $t_n = t_0 + n \Delta t$ and y is the vector of approximate values of $y$ at those times, $y(t_0), y(t_1), y(t_2), \ldots$.

For forwardeuler and rungekutta2, feel free to copy code from the lecture notebook.

In [ ]:
function forwardeuler(f, y0, t0, Δt, N)

function rungekutta2(f, y0, t0, Δt, N) # a.k.a. midpoint method

function rungekutta4(f, y0, t0, Δt, N) 

Problem 2

Compare the three numerical solution methods on the initial value problem

\begin{align*} y' &= t y + t^3 \\ y(0) &= 1 \end{align*}

to the true solution $y(t) = [2 + y(0)] \, e^{t^2/2} - t^2 - 2$ over the interval $0 \leq t \leq 2$. Your comparison should consist of two plots. The first plot should show $y(t)$ versus $t$ for the true solution and the three numerical solutions. The second plot should show the global errors of the three methods at the last time step versus the time step $\Delta t$.

For guidance on the plots and the structure of the computations to produce them, look at the lecture notebook for 2016-12-02.

In [ ]:

Problem 3

Revise your three numerical IVP solvers to work when $y$ and $f$ are vector-valued. If $d$ is the dimension of $y$ and $f$, the y output should now be an $d \times (N+1)$ matrix, whose columns are the $d$-dimensional vectors $y(t_0), y(t_1), y(t_2), \ldots$.

In [ ]:

Problem 4

The equations of motion for a small satellite orbiting a much more massive central body are

\begin{align*} x'' &= -\frac{Gmx}{(x^2+y^2)^{3/2}} \\ y'' &= -\frac{Gmy}{(x^2+y^2)^{3/2}} \end{align*}

In these equations, $(x,y)$ is the position of the satellite in the plane of the orbit, $G$ is the universal gravitational constant, and $m$ is the mass of the central body at the origin.

(a) Convert this system of two second-order ODEs to a system of four first-order ODEs in the four variables $x, y, v_x, v_y$, where $v_x = x'$ and $v_y = y'$.

(b) Write a Julia function f(t,x) that maps the vector $x = [x, y, v_x, v_y]$ into $x' = [x', y', v_x', v_y']$ according to the system of four first-order ODEs you derived in (a). Since this is a math class and not a physics class, you can set the physical constant $G$ and $m$ to 1.

(c) Compute numerical solutions to the system of ODEs using the initial condition

\begin{align*} x(0) &= 1 \\ y(0) &= 0 \\ v_x(0) &= 0 \\ v_y(0) &= 0.6 \end{align*}

using Forward Euler, 2nd-order Runge-Kutta, and 4th-order Runge-Kutta, using $\Delta t = 0.04$ over the interval $0 \leq t \leq 10$

(d) Make a plot of the $x(t), y(t)$ orbit of the satellite as computed using the three time-stepping algorithms. Plot Forward Euler with a blue line, 2nd-order Runge-Kutta in red, and 4th-order Runge-Kutta in green. Make sure the axes are equispaced in $x$ and $y$ by using axis("equal"), and crop the plot to $-1 \leq x \leq 1.5$ and $-1.5 \leq y \leq 1$ using xlim(-1, 1.5); ylim(-1.5, 1).

(e) Comment on your results. Describe the differences between the three approximate orbits. Which of them makes most sense from a physical perspective, and why?

In [ ]: