Matlab Assignment 2

NJIT Math 222 Spring 2018

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.

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 y1(t)y_1(t) and y2(t)y_2(t), let's call them x(t)x(t) and y(t)y(t). To solve a system of the form
dxdt=f(t,x,y)dydt=g(t,x,y),\begin{aligned} \frac{dx}{dt}& = f(t,x,y)\\ \frac{dy}{dt}& = g(t,x,y), \end{aligned}
we define a sequence of times tk=kht_k = k\cdot h and a sequence of approximations (xk,yk)(x(tk),y(tk))(x_k,y_k)\approx(x(t_k),y(t_k)), along with the obvious modification to the Euler updating rule:
xk+1=xk+hf(tk,xk,yk)yk+1=yk+hg(tk,xk,yk).\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
dxdt=tsinx+ydtdt=t2+xsiny,\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)(x,y)=(3,2) from t=0t=0 to t=6t=6. Note that I plot the solution twice. The first graph shows both xx and yy as functions of tt. The second graph shows the trajectory in the (x,y)(x,y) with initial condition θ=1\theta=1, θ=0\theta'=0.)-$plane.

Your Assignment

In Section 7.1 of Boyce and DiPrima, you learn how to convert a second-order equation
into a system of equations second order equations. In this assignment you will be considering the forced damped pendulum equations:
d2θdt2+γdθdt+sinθ=F0cosωt.\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 γ=0\gamma=0 and F=0F=0 from t=0t=0 to t=20t=20 with initial condition θ=1\theta=1, θ=0\theta'=0. If you choose time step h=0.1h=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)(x,y) plane. Instead, you should see a spiral.) Experiment with the value of hh 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 hh you used.
  3. Using the same value of hh, solve the system with the same initial conditions but γ=0.1\gamma=0.1, F=0.1F=0.1, and ω=1\omega=1. Run the simulation up to t=200t=200. Describe what you see.
  4. Now increase the forcing to F=1F=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
x=(aby)xy=(cxd)y\begin{aligned} x' & = (a - b y)x \\ y' & = (cx -d)y \end{aligned}
Here x(t)x(t) represents the population of a prey species, say mice, and y(t)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=1a=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?