# Matlab Assignment 2

### Due week of 4/23

This problem has two parts. In part one, we will extend what we know about Runge-Kutta methods to systems of equations. In the second, we will use pplane to examine some two-dimensional phase-planes, similar to our use of dfield to draw direction fields.

To do the second part of the assignment you need to download some additional programs, following the instructions here.

• If you chose to work with the MATLAB-based programs in assignment one, then you already have everything you need. The only difference in the instructions is that you type pplane2 at the MATLAB prompt (after the characters >> ).
• If you used the Java versions of the programs, available here, make sure you download and run the program named pplane.jar.

### Problem 1: Euler method for a second order equation

Runge-Kutta with two or more dependent variables is exactly the same as with one dependent variable. Suppose we have a system of two equations. Instead of calling the variables $y_1(t)$ and $y_2(t)$, let's call them $x(t)$ and $y(t)$. To solve a system of the form
\begin{aligned} \frac{dx}{dt}& = f(t,x,y)\\ \frac{dy}{dt}& = g(t,x,y), \end{aligned}
we define a sequence of times $t_k = k\cdot h$ and a sequence of approximations $(x_k,y_k)\approx(x(t_k),y(t_k))$, along with the obvious modification to the Euler updating rule:
\begin{aligned} x_{k+1}& = x_k + h\cdot f(t_k,x_k,y_k)\\ y_{k+1}& = y_k + h \cdot g(t_k,x_k,y_k). \end{aligned}

In the attached MATLAB-generated webpage, I solve the system
\begin{aligned} \frac{dx}{dt}& = t-\sin{x}+y\\ \frac{dt}{dt}& = t^2+x-\sin{y}, \end{aligned}
with initial condition $(x,y)=(3,2)$ from $t=0$ to $t=6$. Note that I plot the solution twice. The first graph shows both $x$ and $y$ as functions of $t$. The second graph shows the trajectory in the $(x,y)$ with initial condition $\theta=1$, $\theta'=0$.)-\$plane.

In Section 7.1 of Boyce and DiPrima, you learn how to convert a second-order equation
$x''(t)=f(t,x,x')$
into a system of equations second order equations. In this assignment you will be considering the forced damped pendulum equations:
$\frac{d^2\theta}{dt^2}+ \gamma \frac{d\theta}{dt} + \sin{\theta} = F_0 \cos{\omega t}.$
(We have chosen our units to make the problem as simple-looking as possible.)

1. (Pencil and paper) Rewrite this as a system of two first order equations, suitable for solving using the Euler method.
2. Modify the example program to solve the undamped, unforced system obtained by setting $\gamma=0$ and $F=0$ from $t=0$ to $t=20$ with initial condition $\theta=1$, $\theta'=0$. If you choose time step $h=0.1$, you will find a strange result: the amplitude of the pendulum's motion will grow slowly over time! (You know this shouldn't happen because the system conserves energy. That means the solutions should be periodic, and should look like closed curves trajectories are viewed in the $(x,y)$ plane. Instead, you should see a spiral.) Experiment with the value of $h$ to make the spiraling effect disappear, at least to the accuracy visible to the naked eye on your plots. Plot the trajectory and state how small a value of $h$ you used.
3. Using the same value of $h$, solve the system with the same initial conditions but $\gamma=0.1$, $F=0.1$, and $\omega=1$. Run the simulation up to $t=200$. Describe what you see.
4. Now increase the forcing to $F=1$. Describe what you see. This is an example of mathematical chaos.

### Problem 2: Using software to plot phase planes

Using the phase plane program described in the introduction, plot the phase plane for the Lotka-Volterra model
\begin{aligned} x' & = (a - b y)x \\ y' & = (cx -d)y \end{aligned}
Here $x(t)$ represents the population of a prey species, say mice, and $y(t)$ is the population of a predator species, for example owls. If there are no owls, then the population of of mice grows exponentially, and if there are no mice, then the owls die off. However if the initial populations are both positive, then interesting stuff happens. Plot the phase plane including a number of solutions of this system with $a=b=c=d=1$. You decide what is are good limits for your axes. Describe the solutions in a few words. Do you think that this model gives realistic behavior?