Assignment 3 MATH4Q3 4 February 2000


This assignment is due in the locker in the basement of the Burke Science Building by 15:00 on Friday 18 February 2000.

1.
Milne's fourth-order predictor-corrector method for the initial value problem y'=f(x,y), y(0)=y0 is given by the formulas
y(P)n+1 = $\displaystyle y_{n-3} + \frac{4h}{3} [2f_n - f_{n-1} + 2 f_{n-2}] + \frac{14}{45} h^5 y^{(5)}(\xi)$  
y(C)n+1 = $\displaystyle y_{n-1} + \frac{h}{3}[f(x_{n+1},y^{(P)}_{n+1}) + 4 f_n + f_{n-1}] - \frac{1}{90} h^5 y^{(5)}(\xi)$  

where the terms proportional to h5 are the local truncation errors.
(a)
By applying the corrector of Milne's method to the model problem $y'=\lambda y$, show that it leads to a differential equation of the form

 
a yn+1 + b yn + c yn-1 = 0 (1)

where a,b,c are constants depending on $\lambda$ and h.
(b)
Assume that the solution of (1) has the form $y_n=\theta^n$ where $\theta\ne 0$ is a constant, and determine the two roots of (1), $\theta_1$, $\theta_2$. The stability region for Milne's method is determined by $\vert\theta_1\vert\le 1$and $\vert\theta_2\vert\le 1$. What happens to the stability of the method when $h\lambda \rightarrow 0$?
(c)
Assuming that the two fifth-order derivatives in the error terms are the same, find an estimate of the local truncation error en+1 for Milne's fourth-order method in terms of y(P)n+1 and y(C)n+1. How could this estimate be used in practice?
2.
(a)
Write Matlab functions to integrate the initial value problem

\begin{displaymath}y' = f(x,y), \mbox{\hspace{1em}}y(a)=y_0, \end{displaymath}

on an interval [a,b] using:
  • Euler's method
  • Modified Euler
  • Improved Euler
  • RK4
It is suggested that you implement, for example, Improved Euler as
[x,y]=eulerimp('f',a,y0,b,stepsize), where (x,y)=(xn,yn) is the computed solution.
(b)
Use the four Matlab function you wrote in (a) with stepsize h=1/40 and 1/80 to solve the following ODE:

 \begin{displaymath}y'=-y^3/2, \mbox{\hspace{1em}}y(0)=1 \mbox{\hspace{1em}}\mbox{with} \mbox{\hspace{1em}}y_{exact}(x) = \frac{1}{\sqrt{x+1}}. \end{displaymath} (2)

Calculate the absolute errors (i.e. |yexact(x) - ycomputed(x)|) at x=1 to see if they are roughly reduced by a half for Euler's method, by a quarter for Modified Euler and Improved Euler, and by 1/16 for RK4, as the stepsize h is halved from 1/40 to 1/80. Discuss the performance of the methods by comparing the computation time used by each method to calculate y(1) with stepsize h=1/40 and 1/80 (use the Matlab functions tic and toc to do the timing).
3.
(a)
Let g(z) be a user-defined function with possibly complex-valued coefficients and variables. Then the set of all such z in the complex plane that satisfy |g(z)|=1 is in general a curve or consists of several curves. Write a Matlab script to plot the curves(s) (in the complex plane) determined by the equation |g(z)|=1, where $g(z)=\sum_{k=1}^{10} z^{k}/k!$.
(b)
Recall that one of the RK3 methods is as follows:
yn+1 = $\displaystyle y_n + \frac{1}{6}(k_1+4k_2+ k_3), \mbox{\hspace{1em}}\mbox{where}$  
k1 = hf(xn,yn)  
k2 = $\displaystyle hf(x_n+\frac{1}{2}h, y_n + \frac{1}{2} k_1),$  
k3 = hf(xn+h, yn-k1+2k2).  

Show that the stability region of the above RK3 method is determined by $\vert E(\lambda h)\vert \le 1$ where

\begin{displaymath}E(\lambda h) = 1 + \lambda h + \frac{(\lambda h)^2}{2!} + \frac{(\lambda h)^3}{3!}. \end{displaymath}

(c)
In general, the stability region of a p-th order Runge-Kutta method is determined by $\vert E(\lambda h)\vert \le 1$ with $E(\lambda h) = \sum_{k=0}^{p} \frac{(\lambda h)^k}{k!}$. Write Matlab scripts (similar to that used in (a)) to plot the stability regions for p=1,2,3,4,5 by letting $\lambda h = z$. Are the stability regions larger for larger p?
4.
(a)
Modify your Runge-Kutta program of question 2(a) to solve a second-order ordinary differential equation.
(b)
Duffing's equation

\begin{displaymath}\frac{d^2 y}{dt^2} + y - \epsilon y^3 = 0, \mbox{\hspace{1em}}y(0)=1, \mbox{\hspace{1em}} \frac{dy}{dt}(0) = 0, \end{displaymath}

describes an oscillator with a slightly nonlinear restoring force. The size of the nonlinear part of the restoring force is given by the parameter $\epsilon$. Use your program from (a) with stepsize h=0.1 to solve Duffing's equation with $\epsilon=0.05$ from t=0 to $t=6\pi$. Plot the results compared with the linear oscillator ( $\epsilon=0$). Can the small term $\epsilon y^3$ be neglected? What is the effect of the nonlinearity?

About this document ...

This document was generated using the LaTeX2HTML translator Version 98.1p1 release (March 2nd, 1998)

Copyright © 1993, 1994, 1995, 1996, 1997, Nikos Drakos, Computer Based Learning Unit, University of Leeds.

The command line arguments were:
latex2html -split 0 -no_navigation 4q3-3-00.

The translation was initiated by Nicholas Kevlahan on 2000-02-18

Nicholas Kevlahan
2000-02-18