You're probably already familiar with Euler's method for approximating the solution to a differential equation. We want to approximate a solution to \(y'(t) = f(t,y(t))\).

Let's take a closer look at Euler's method:

1. Recall that a taylor series for a continuously differentiable function \(g(x)\) centered at \(x=c\) is \begin{equation*} g(x) = g(c) + \frac{g'(c)}{1!}(x-c) + \frac{g''(c)}{2!}(x-c)^2 + \frac{g^{(3)}(c)}{3!}(x-c)^3 + \frac{g^{(4)}(c)}{4!}(x-c)^4 + \cdots \end{equation*}

2. Write a Taylor series for \(y(t)\) by replacing \(g\) with \(y\), \(x\) with \(t+h\) and \(c\) with \(t\) \begin{equation*} y(t+h) = \dots \end{equation*}

3. Since we know that \(y'(t) = f(t,y(t))\) we can rewrite the Taylor series as \begin{equation*} y(t+h) = \dots \end{equation*}

4. Now we can look at this as a difference equation that is ready made for numerical implementation with a loop: \begin{equation*} y_{n+1} = y_n + h f(t_n,y_n) + \mathcal{O}(h^2) \end{equation*}

5. The approximation error for Euler's method is actually not second order as the previous equation would suggest. Instead, if we rearrange to write Euler's formula as an approximation of the derivative we get \begin{equation*} \frac{y_{n+1}-y_n}{h} = f(t_n,y_n) + \mathcal{O}(h). \end{equation*} What does this formula tell you about the accuracy of Euler's method?

Write code to implement Euler's method for initial value problems. Your MATLAB function should accept as input: \(f(t,y)\), tmin, tmax, the number of grid points (the value of \(h = \Delta t\) should be calculated within your code), and an intial condition. The output should be vectors for \(t\) and \(y\).

function [t,y] = MyEuler1D(f,tmin,tmax,num_pts,IC)

Test your code on a first order differential equation where you know the answer and then test your code on the differential equation \begin{equation*} y' = -\frac{1}{3}y+\sin(t) \text{ where } y(0) = 1. \end{equation*}

Write code that implements a 2D version of Euler's method that will solve a system of two differential equations in two dependent variables. Test your code on the following problem by showing a time evolution plot (time on \(x\) and populations on \(y\)) as well as a phase plot (\(x\) on the \(x\) and \(y\) on the \(y\) with time understood implicitly):

The Lotka-Volterra Predator-Prey Model:

Let \(x(t)\) denote the number of rabbits (prey) and \(y(t)\) denote the number of foxes (prey) at time \(t\). The relationship between the species can be modeled by the classic 1920's Lotka-Volterra Model: \begin{equation*} \left\{ \begin{array}{ll} x' \amp = \alpha x - \beta xy \\ y' \amp = -\delta y + \gamma xy \end{array} \right. \end{equation*} wherer \(\alpha, \beta, \gamma,\) and \(\delta\) are positive constants. For this problems take \(\alpha \approx 1\), \(\beta \approx 0.05\), \(\gamma \approx 0.01\), and \(\delta \approx 1\). Be sure to explain the meaning of each of the parameters and each of the components of the model.

The Midpoint Method:

Now we begin the journey of creating better sovlers than Euler. The midpoint method is defined by first taking a half step with Euler's method to approximate a solution at time \(t_{n+1/2} \equiv (t_n + t_{n+1})/2\) and then taking a full step using the value of \(f\) at \(t_{n+1/2}\) and the approximate \(y_{n+1/2}\). \begin{align*} y_{n+1/2} \amp = y_n + \frac{h}{2} f(t_n, y_n)\\ y_{n+1} \amp = y_n + h f(t_{n+1/2}, y_{n+1/2}) \end{align*}

Note: Indexing by \(1/2\) in a computer is nonsense. Instead, we implement the midpoint method with: \begin{align*} y_{temp} \amp = y_n + \frac{h}{2} f(t_n, y_n)\\ y_{n+1} \amp = y_n + h f\left( \frac{t_n+t_{n+1}}{2}, y_{temp}\right) \end{align*}

Write MATLAB code to implement the midpoint method

function [t,y]=MyMidpointMethod(f,tmin,tmax,num_pts,IC)

Test your code against your Euler1D code on the same single variable ODE as before. You will likely see very little difference on a very small step size, but for a smaller number of points there will be a remarkable difference.

The Runge-Kutta Method:

Another method for approximating the solution to a first order initial value problem is to take several approximations and average them in a smart way. The Runge-Kutta 4 method is one (of many) such methods. In this method, each \(k_j\) is an approximation of the slope and we combine them in as a weighted average in the end. \begin{align*} k_1 \amp = f(t_n, y_n)\\ k_2 \amp = f(t_n + \frac{h}{2} , y_n + \frac{h}{2} k_1)\\ k_3 \amp = f(t_n + \frac{h}{2} , y_n + \frac{h}{2} k_2)\\ k_4 \amp = f(t_n + h , y_n + h k_3)\\ y_{n+1} \amp = y_n + \frac{h}{6} \left( k_1 + 2 k_2 + 2 k_3 + k_4 \right) \end{align*}

Write a MATLAB function that implements the Runge-Kutta 4 method in one dimension.

function [t,y]=MyRk4(f,tmin,tmax,num_pts,IC)

Test the problem on a known differential equation.

Modify your Runge-Kutta 4 code to work for two dependent variables. I'll get you started: We want to solve \begin{equation*} \left\{ \begin{array}{ll} x' \amp = f(t,x,y) \\ y' \amp = g(t,x,y) \end{array} \right. \end{equation*} and to do so we extend the Runge Kutta method as \begin{align*} k_1 \amp = f(t_n,x_n, y_n)\\ q_1 \amp = g(t_n,x_n, y_n)\\ k_2 \amp = f(t_n + \frac{h}{2} ,x_n+\frac{h}{2} k_1, y_n + \frac{h}{2} q_1)\\ q_2 \amp = g(t_n + \frac{h}{2} ,x_n+\frac{h}{2} k_1, y_n + \frac{h}{2} q_1)\\ k_3 \amp = \dots\\ q_3 \amp = \dots\\ k_4 \amp = \dots\\ q_4 \amp = \dots\\ x_{n+1} \amp = x_n + \frac{h}{6} \left( k_1 + 2 k_2 + 2 k_3 + k_4 \right)\\ y_{n+1} \amp = y_n + \frac{h}{6} \left( q_1 + 2 q_2 + 2 q_3 + q_4 \right) \end{align*}

Test your code on the predator prey model in problem 6.1.3.

Solving systems of ordinary differential equations would become challenging if we were to continue coding in the same way as in the previous problem. Write a MATLAB function that accepts any number of right-hand sides from a system of differential equations and then leverages the fact that MATLAB works very well with vectors to create Euler and Runge-Kutta solutions to these systems. Devise several systems to test your code (including 1D and 2D).

Consider problem 14 in Chapter 11 of the textbook (pages 296-298). Complete part (a) as instructed in the book. For part (b) you will have to create RK and Euler code that works for four dependent variables (or use your code from the previous problem). You do not have to do part (c) unless you want to teach yourself some of the build-in techniques for ODE solving that MATLAB has.

Consider problem 16 in Chapter 11 of the textbook (pages 298-299). Complete part (a) using both your MyEuler2D, and MyRK2D codes. Don't worry about part (b).

The major trouble with the RK (and midpoint) methods is that they take many more function evaluations than Euler's method. We can improve upon Euler's method in the following way:

We want to solve \(y' = f(t,y)\) so:

1. Approximate the derivative by looking forward in time(!) \begin{equation*} \frac{y_{n+1} - y_n}{h} \approx f(t_{n+1} , y_{n+1}) \end{equation*}

2. Rearrange to get the difference equation \begin{equation*} y_{n+1} = y_n + h f(t_{n+1},y_{n+1}). \end{equation*}

3. Notice that we will have \(t_{n+1}\) but we do not have \(y_{n+1}\). The major trouble is that \(y_{n+1}\) shows up on both sides of the equation. Can you think of a way to solve for it? … you have code that does this step!!!

4. This method is called the backward Euler method and is known as an implicit method since you need to solve a nonlinear equation at each step. The advantage (usually) is that you can take far fewer steps with reasonably little loss of accuracy.

Write MATLAB code to implement the backward Euler's method for a 1D initial value problem.

This problem is a spin off of the first problem in section 11.2.6 in the text as well as exercise #8 on page 295.

In this model there are two characters, Romeo and Juliet, whose affection is quantified on the scale from \(-5\) to \(5\), with \(-5\) being hysterical hatred, 0 being indifference, and 5 being ecstatic love.

The characters struggle with frustrated love due to the lack of reciprocity of their feelings. Mathematically,

• Romeo: “My feelings for Juliet decrease in proportion to her love for me.”

• Juliet: “My love for Romeo grows in proportion to his love for me.”

• Juliet's emotional swings lead to many sleepless nights, which consequently dampens her emotions.

This give rise to \begin{equation*} \left\{ \begin{array}{ll} \frac{dx}{dt} \amp = -\alpha y \\ \frac{dy}{dt} \amp = \beta x - \gamma y^2 \end{array} \right. \end{equation*} where \(x(t)\) is Romeo's love for Juliet and \(y(t)\) is Juliet's love for Romeo at time \(t\).

Your tasks:

1. First implement this 2D system with \(x(0) = 2\), \(y(0)=0\), \(\alpha=0.2\), \(\beta=0.8\), and \(\gamma=0.1\) for \(t \in [0,60]\). What is the fate of this pair's love under these assumptions?

2. Write MATLAB code that approximates the parameter \(\gamma\) that will result in Juliet having a feeling of indifference at \(t=30\). Your code should not need human supervision: you should be able to tell it that you're looking for indifference at \(t=30\) and turn it loose to find an approximation for \(\gamma\). Assume throughout this problem that \(\alpha=0.2\), \(\beta=0.8\), \(x(0)=2\), and \(y(0)=0\). Write a description for how your code works in your homework document and also submit your MATLAB file (along with any support files) demonstrating how it works. Hint: One way to do this problem is very similar to a bisection method for root finding.

   • Shoot two solutions with two different parameters.

   • Compare their results at the desired time against the result of shooting with the average value of the parameter.

   • Use the logic of the bisection method to make a new estimate of the parameter.

We wish to solve the boundary valued problem \(x'' + 4x = \sin(t)\) with initial condition \(x(0)=1\) and boundary condition \(x(1)=2\). Notice that you do not have the initial position and initial velocity as you normally would with a second order differential equation. Devise a method for finding a numerical solution to this problem.

Hint: First write the problem as a system of first order differential equations. Then think about how your bisection method code might help you.