Eulreexample

Lecture 30 Euler Methods
Numerical Solution of an IVP

Suppose we wish to numerically solve the initial value problem y = f (t, y), on an interval of time [a, b]. By a numerical solution, we must mean an approximation of the solution at a nite number of points, i.e. (t0 , y0 ), (t1 , y1 ), (t2 , y2 ), . . . , (tn , yn ), where t0 = a and tn = b. The rst of these points is exactly the initial value. If we take n steps as above, and the steps are evenly spaced, then the time change in each step is h= ba , n (30.2) y(a) = y0 , (30.1)
and the times ti are given simply by ti = a + ih. This leaves the most important part of nding a numerical solution: determining y1 , y2 , . . . , yn in a way that is as consistent as possible with (30.1). To do this, rst write the dierential equation in the indexed notation yi f (ti , yi ), (30.3)
and will then replace the derivative y by a dierence. There are many ways we might carry this out and in the next section we study the simplest.
The Euler Method

The most straight forward approach is to replace yi in (30.3) by its forward dierence approximation. This gives yi+1 yi = f (ti , yi ). h Rearranging this gives us a way to obtain yi+1 from yi known as Eulers method: yi+1 = yi + hf (ti , yi ). (30.4)
With this formula, we can start from (t0 , y0 ) and compute all the subsequent approximations (ti , yi ). This is very easy to implement, as you can see from the following program (which can be downloaded from the class web site):
109
110
LECTURE 30. EULER METHODS
function [T , Y] = myeuler(f,tspan,y0,n) % function [T , Y] = myeuler(f,tspan,y0,n) % Solves dy/dt = f(t,y) with initial condition y(a) = y0 % on the interval [a,b] using n steps of Euler s method. % Inputs: f -- a function f(t,y) that returns a column vector of the same % length as y % tspan -- a vector [a,b] with the start and end times % y0 -- a column vector of the initial values, y(a) = y0 % n -- number of steps to use % Outputs: T -- a n+1 column vector containing the times % Y -- a (n+1) by d matrix where d is the length of y % Y(t,i) gives the ith component of y at time t % parse starting and ending points a = tspan(1); b = tspan(2); h = (b-a)/n; % step size t = a; y = y0; % t and y are the current variables T = a; Y = y0; % T and Y will record all steps for i = 1:n y = y + h*f(t,y); t = a + i*h; T = [T; t]; Y = [Y; y]; end To use this program we need a function, such as the vector function for the pendulum: > f = inline([y(2);-.1*y(2)-sin(y(1)) + sin(t)],t,y) Then type: > [T Y] = myeuler(f,[0 20],[1;-1.5],5); Here [0 20] is the time span you want to consider, [1;-1.5] is the initial value of the vector y and 5 is the number of steps. The output T contains times and Y contains values of the vector as the times. Try: > size(T) > size(Y) Since the rst coordinate of the vector is the angle, we only plot its values: > theta = Y(:,1); > plot(T,theta) In this plot it is clear that n = 5 is not adequate to represent the function. Type: > hold on then redo the above with 5 replaced by 10. Next try 20, 40, 80, and 200. As you can see the graph becomes increasingly better as n increases. We can compare these calculations with Matlabs built-in function with the commands: > [T Y]= ode45(f,[0 20],[1;-1.5]); > theta = Y(:,1); > plot(T,theta,r)
111
The problem with the Euler method

You can think of the Euler method as nding a linear approximate solution to the initial value problem on each time interval. An obvious shortcoming of the method is that it makes the approximation based on information at the beginning of the time interval only. This problem is illustrated well by the following IVP: x + x = 0, x(0) = 1, x(0) = 0 . (30.5)
You can easily check that the exact solution of this IVP is x(t) = cos(t). If we make the standard change of variables y1 = x, then we get y1 = y2 , y2 = y1 . y2 = x,
Then the solution should be y1 (t) = cos(t), y2 (t) = sin(t). If we then plot the solution in the (y1 , y2 ) plane, we should get exactly a unit circle. We can solve this IVP with Eulers method: > g = inline([y(2);-y(1)],t,y) > [T Y] = myeuler(g,[0 4*pi],[1;0],20) > y1 = Y(:,1); > y2 = Y(:,2); > plot(y1,y2) As you can see the approximate solution goes far from the true solution. Even if you increase the number of steps, the Euler solution will eventually drift outward away from the circle because it does not take into account the curvature of the solution.
The Modied Euler Method

An idea which is similar to the idea behind the trapezoid method would be to consider f at both the beginning and end of the time step and take the average of the two. Doing this produces the Modied (or Improved) Euler method represented by the following equations: k1 = hf (ti , yi ) k2 = hf (ti + h, yi + k1 ) 1 yi+1 = yi + (k1 + k2 ) . 2 (30.6)
Here k1 captures the information at the beginning of the time step (same as Euler), while k2 is the information at the end of the time step. The following program implements the Modied method. It may be downloaded from the class web site.
112
LECTURE 30. EULER METHODS
function [T , Y] = mymodeuler(f,tspan,y0,n) % Solves dy/dt = f(t,y) with initial condition y(a) = y0 % on the interval [a,b] using n steps of the modified Eulers method. % Inputs: f -- a function f(t,y) that returns a column vector of the same % length as y % tspan -- a vector [a,b] with the start and end times % y0 -- a column vector of the initial values, y(a) = y0 % n -- number of steps to use % Outputs: T -- a n+1 column vector contianing the times % Y -- a (n+1) by d matrix where d is the length of y % Y(t,i) gives the ith component of y at time t % parse starting and ending points a = tspan(1); b = tspan(2); h = (b-a)/n; % step size t = a; y = y0; % t and y are the current variables T = a; Y = y0; % T and Y will record all steps for i = 1:n k1 = h*f(t,y); k2 = h*f(t+h,y+k1); y = y + .5*(k1+k2); t = a + i*h; T = [T; t]; Y = [Y; y]; end We can test this program on the IVP above: > [T Ym] = mymodeuler(g,[0 4*pi],[1;0],20) > ym1 = Ym(:,1); > ym2 = Ym(:,2); > plot(ym1,ym2)
Exercises
30.1 Download the les myeuler.m and mymodeuler.m from the class web site. 1. Type the following commands: > hold on > f = inline(sin(t)*cos(x),t,x) > [T Y]=myeuler(f,[0,12],.1,10); > plot(T,Y) Position the plot window so that it can always be seen and type: > [T Y]=myeuler(f,[0,12],.1,20); > plot(T,Y) Continue to increase the last number in the above until the graph stops changing (as far as you can see). Record this number and print the nal graph. Type hold off and kill the plot window. 2. Follow the same procedure using mymodeuler.m . 3. Describe what you observed. In particular compare how fast the two methods converge as n is increased (h is decreased).

Eulreexample

Uploaded by

Copyright:

Available Formats

Eulreexample

Uploaded by

Document Information

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Eulreexample

Uploaded by

Copyright:

Available Formats

Lecture 30 Euler Methods

Numerical Solution of an IVP

The Euler Method

LECTURE 30. EULER METHODS

The problem with the Euler method

The Modied Euler Method

LECTURE 30. EULER METHODS

You might also like