Math 0250 Computer Assignment 3

In this assignment, you'll write a short program to implement the ``Forward Euler Method'' to solve a first order differential equation. You should turn in answers to the questions enumerated at the end of this page. Be sure to give complete answers to the questions. Don't just turn in a bunch of numbers.

Recall that Euler's method solves an equation of the form

\begin{eqnarraystar}
x'(t) &=& f(t,x(t))\  x(t_0) &=& x_0\end{eqnarraystar}

by the recursive formulas

\begin{eqnarraystar}
t_{n+1} &=& t_n + h\  x_{n+1} &=& x_n + h f(t_n,x_n)\end{eqnarraystar}

starting with the known initial values t0 and x0. Here h is a fixed small positive number, and you hope that the table of values generated by the recipe above will be a good approximation to the exact solution if h is small enough. Here you'll test that hypothesis by trying it out on the equation

\begin{eqnarraystar}
x' &=& -x^2\ x(0)&=&1.\end{eqnarraystar}

Notice that this particular equation can be solved exactly:

x = 1/(t+1)

so you'll be able to compare the approximate solution coming from Euler's method with the exact solution.

You'll first want to write a short computer program to run Euler's method for you. Here's how to do it with Matlab or Octave. Fire up your favorite text editor (``Notepad'' will do if you're using Windows), and create the file em.m containing the following lines:

function [T,X]=em(h)

   t0=0;
   last_t=10;
   x0=1;

   T=[];
   X=[];
   x = x0;

   for t = t0:h:last_t
      X=[X,x];
      T=[T,t];
      x = x + h * (-x^2);
   end

Here's an explanation of what's going on in the program. The first line tells the system that you're defining a function called em that will receive a value h and return a pair of values T and X. Eventually, you'll call this function from the Matlab or Octave command line, giving it a value for h, the increment to use for Euler's Method, and it will return the lists of t and x values in T and X. It's set up this way so that you can easily pass in different values of h without having to edit the program itself.

The next three lines assign labels to the starting values of t and x. This is for convenience, in case you want to change them later.

Next, the lists to receive the t and x values are initialized to start out empty, and the first value of x is set to x0.

The for loop is where the real action is. The beginning line of the loop says to repeat the loop with values of t running from t0 (0 in this case) up to last_t (in this case 10) in increments of h. Each time through the loop, you append the current values of x and t to the lists X and T respectively, and apply the Euler recursion formula to calculate the next value of x.

The semicolons at the ends of lines suppress the output, so that you don't see a lot of junk on your screen as the calculations proceed.

Save the program, exit your editor, and then fire up Matlab or Octave. If you've started Matlab from the Windows ``Start'' menu, you may have to change your current directory so that Matlab can find your program. For example, if you saved your program on a floppy disk, you'd enter

cd a:
If you're using Matlab, enter
more on
to keep the output from scrolling of the top of the screen. At the prompt, enter
[T,X] = em(1);
to run Euler's method with h = 1, and store the results in T and X. You won't see the result, because screen output was suppressed by the semicolon. To see the values in the lists T and X, each arranged in a row, you can just enter those symbols without a trailing semicolon. A more useful way to display the results is to make a table by lining T and X up in adjacent columns. This will show the functional relationship between t and x. To make the table, enter
[T',X']
The two "primes" convert the rows to columns, and putting the columns inside square brackets, separated by commas, places them side by side and properly aligned.

You can also plot x as a function of t by entering

plot(T,X)
Neither your table nor your plot will be very interesting, nor will they give a good approximation to the exact solution. That's because the value you used for h was too big.

Now here are some things for you to try.

1.
Run Euler's method again, but this time with h = 0.1. Make a table showing the approximate and exact values of x when t = 0, 2, 4, 6, 8, 10. Also include a plot of the approximate solution, obtained by either transcribing by hand from the screen, or by printing a computer plot.
2.
The error in an approximation is the difference between the approximate value and the exact value. At t = 5 and 10, calculate the error in the Euler approximation with h = 1, 0.5, 0.1, and 0.01. You may want to enter format long to get more decimal places.


Frank Beatrous
2/2/1998