Skip to the content of the web site.

Laboratory 10

This laboratory will investigate two topics:

  1. Systems of IVPs and Higher-order IVPs (MATH 211), and
  2. Root Finding (MATH 215).

10.1 Systems of IVPs and Higher-order IVPs

We have seen initial-value problems where we have that the first-derivative of one function of one variable is described by a function of time and the value of the function:

tex:$$y^{(1)}(t) = f(t, y(t))$$.

How do we deal with a differential equation with a higher order, or a system of coupled differential equations.

So as to not alarm you, the beauty is that the code you have already written will be reused with only three small changes to three line of code.

10.1.1 Predator-Prey Models

Suppose you have rabbits: left to themselves, they will grow exponentially (assuming adequate food), but if you introduce foxes into the equation, the more foxes, the two will play one off another: an increase in rabbits will result in an increase in foxes, but at some point, the foxes will overpopulate and decimate the rabbits, at which point, the foxes die out, once again allowing the rabbits to begin increasing exponentially.

Suppose that Consider the following system of differential equations:

tex:$$\dot{r}(t) = r(t)(1 - 2f(t))$$
tex:$$\dot{f}(t) = -f(t)(1 - r(t))$$

Here tex:$$r(t)$$ represents the density of rabbit families per acre while tex:$$f(t)$$ represents the density of fox families per acre. If you begin with the density of rabbits of 0.5 families per acre and a density of foxes of 0.3 families per acre, the relatively low number of foxes allows the rabbit population to grow exponentially and the equally relative low number of rabbits causes the foxes to begin to die off. However, as the population of rabbits increases, the population of foxes rebounds until they begin to decimate the rabbits. This cycle is shown in Figure 1.


Figure 1. The population densities of foxes and rabbits over time.

To understand this, we note that at time tex:$$t = 0$$, tex:$$\dot{r}(0) = r(0)(1 - 2f(0)) = 0.5\cdot(1 - 2\cdot 0.3) = 0.4$$, and therefore the population density of rabbits is increasing. The fox population density, however, is decreasing: tex:$$\dot{f}(0) = -f(0)(1 - r(t)) = -0.3(1 - 0.5) = -0.15$$. Using Euler's method, we could therefore predict the population at time tex:$$t = 0.1$$: tex:$$r(0.1) \approx r(0) + h \dot{r}(0) = 0.5 + 0.1\cdot 0.4 = 0.54$$ and tex:$$f(0.1) \approx f(0) + h \dot{f}(0) = 0.3 + 0.1\cdot (-0.15) = 0.285$$. We could carry on and now predict tex:$$r(0.2)$$ and tex:$$f(0.2)$$.

The differential equation described above falls into a category of differential equations called predator-prey equations (see the wikipedia article on the Lotka-Volterra equation). There is no analytic solution for these equations and they must be solved numerically.

Another way of looking at the interaction of these two species is to plot tex:$$f(t)$$ versus tex:$$r(t)$$ over time. Such a plot is shown in Figure 2.


Figure 2. An animation of the population densities of foxes and rabbits over time.

Again, when both population densities are relatively small, the number of rabbits increases, but this increase is followed by an increase in foxes which in turn decreases the rabbit population density. Without a food supply, the foxes themselves begin to die off allowing the rabbits a chance to recover.

10.1.2 Lorenz Attractor

Another interesting differential equation are the Lorenz equations. You will have heard of the butterfly effect and this description is, as we will see, a direct response to the solution to this system of three differential equations:

tex:$$\dot{x}(t) = 10(y(t) - x(t))$$
tex:$$\dot{y}(t) = x(t)(28 - z(t)) - y(t)$$

tex:$$\dot{z}(t) = x(t)y(t) - \frac{8}{3}z(t)$$

with, for example, the initial values tex:$$x(0) = -5$$, tex:$$y(0) = 0$$, and tex:$$z(0) = 28$$.

The solutions to these coupled initial-value problems must be solved numerically. The solutions the these three are shown in Figure 3.


Figure 3. The solutions to the Lorenz equations as shown above with the given initial conditions (click for a larger image).

Again, like before, we can make a plot of the three points in space and watch how the three change over time. This is shown in Figure 4 which gives us the familiar butterfly.


Figure 4. A plot over time of the solutions tex:$$(x(t), y(t), z(t))$$.

This is reasonably straight-forward to implement: Modify dp45 so that it accepts a column vector as an initial condition. This requires three changes:

  • Change any reference to y_out( k ) to y_out( :, k ) and any reference to y_out( k + 1 ) to y_out( :, k + 1 ), and
  • Change abs( y1 - z1 ) to calculate norm( y1 - z1, Inf ).

Now, create the function

function [dy] = f10a( t, y )
    dy = [y(1)*(1 - 2*y(2)), -y(2)*(1 - y(1))]';
end

With a little bit of thought, you should notice that this is the predator-prey model described above. If you wanted, you could write it as:

function [dy] = f10a( t, y )
    r = y(1);    % rabbits
    f = y(2);    % foxes

    dy = [r*(1 - 2*f), -f*(1 - r)]';
end

Next, execute the following to give plots similar to Figure 1 and Figure 2:

hold off
[ts, pts] = dp45( @f10a, [0, 6], [0.5 0.3]', 0.1, 1e-8 );
plot( pts(2,:), pts(1,:) )     % rabbits on the vertical axis, foxes on the horizontal
plot( ts, pts(1,:) )           % plot population density of rabbits over time
hold on
plot( ts, pts(2,:) )           % plot population density of foxes over time

Similarly, write a function [dv] = f10c( t, v ) which, in a similar manner, implements the three coupled Lorenz equations. In this case, v(1), v(2) and v(3) will represent tex:$$x(t)$$, tex:$$y(t)$$ and tex:$$z(t)$$, respectively. When you run the following code, you should get an output similar to Figure 5.

[ts, pts] = dp45( @f10c, [0, 100], [-5 0 28]', 0.1, 1e-5 );
size( pts )
  ans =
             3       34145
plot3( pts(1,:), pts(2,:), pts(3,:) )


Figure 5. The output in Matlab of the given commands.

This is your first use of the plot3 command in Matlab for plotting points in three dimensions. You will notice a few entries in the toolbar, highlighted in Figure 6.


Figure 6. Five tools for 3-dimensional plots: Zoom In, Zoom Out, Pan, Rotate 3D, Data Cursor, and Brush/Select Data.

These five tools do, in order:

  • Zoom In: self explanatory,
  • Zoom Out: self explanatory,
  • Pan: move the image in the visual view of the plot,
  • Rotate 3D: rotate the points in three dimensions,
  • Data Cursor: determine the x-, y- and z-values of a point in the image, and
  • Change the colour of a portion of a plot.

Assignment: use the pan tool to move the image around, then rotate the image to change the angle, and finally select the Zoom In tool and click twice on one region of the plot to zoom in on the plot. Save this image and print it. Each student must do their own plot.

10.2 Root Finding

Given a real-valued function of a real variable tex:$$g$$, suppose we want to find a root of the function, that is, a point tex:$$r$$ such that tex:$$g(r) = 0$$.

We will look at three techniques for approximating

10.2.3 Newton's Method

You may have already seen Newton's method.

tex:$$g(r) = g(x_k) + g'(x_k)(r - x_k) + \frac{1}{2}g''(\xi)(r - x_k)^2$$.

First, if tex:$$r$$ is a root, then tex:$$g(r) = 0$$:

tex:$$0 = g(x_k) + g'(x_k)(r - x_k) + \frac{1}{2}g''(\xi)(r - x_k)^2$$.

Next, assuming that the derivative at the point tex:$$x_k$$ is not zero, that is, tex:$$g'(x_k) \ne 0$$, we can divide the equation by that derivative:

tex:$$0 = \frac{g(x_k)}{g'(x_k)} + r - x_k + \frac{1}{2}\frac{g''(\xi)}{g'(x_k)}(r - x_k)^2$$.

The next step is to isolate the error term for tex:$$x_k$$ on one side and group the other:

tex:$$r - \left( x_k - \frac{g(x_k)}{g'(x_k)} \right) = -\frac{1}{2}\frac{g''(\xi)}{g'(x_k)}(r - x_k)^2$$.

The expression in the parentheses is the formula for Newton's method, and thus, this equals tex:$$x_{k + 1}$$:

tex:$$r - x_{k + 1} = -\frac{1}{2}\frac{g''(\xi)}{g'(x_k)}(r - x_k)^2$$.

That is, the error of the tex:$$(k + 1)^{\rm st}$$ approximation is a constant multiplied by the square of the error of the tex:$$k^{\rm th}$$ approximation. Thus, the claim is, if the error is tex:$$0.001$$ at one step, the error will be on the order of tex:$$0.000\,001$$ with the next step, and on the order of tex:$$0.000\, 000\, 001$$ with the next.

Now, assuming that we are close to the root and the second derivative is continuous, then it is a reasonably good approximation that

tex:$$-\frac{1}{2}\frac{g''(\xi)}{g'(x_k)} \approx -\frac{1}{2}\frac{g''(r)}{g'(r)}$$

and we will use this to demonstrate how well Newton's method works.

Consider the function tex:$$g:x \mapsto x \cos(x) + 2$$ and let us start with tex:$$x_0 = \pi$$.

   tex:$$k$$   tex:$$x_k$$Absolute ErrorEstimated Absolute Error
tex:$$\left|\frac{1}{2} \frac{g''(r)}{g'(r)}(r - x_{k - 1})^2\right|$$
03.141592653589793tex:$$6.4284 \times 10^{-1}$$
12.000000000000000tex:$$4.9876 \times 10^{-1}$$tex:$$7.2016 \times 10^{-1}$$
22.522524071586095tex:$$2.3768 \times 10^{-2}$$tex:$$4.3351 \times 10^{-2}$$
32.498648774013250tex:$$1.0698 \times 10^{-4}$$tex:$$9.8451 \times 10^{-5}$$
42.498755760658375tex:$$1.9940 \times 10^{-9}$$tex:$$1.9948 \times 10^{-9}$$
52.498755762652415tex:$$4.4409 \times 10^{-16}$$tex:$$6.9294 \times 10^{-19}$$
62.498755762652415