Numerical methods matlab 3e asu

MATLAB sessions: Laboratory 3

MAT 275 Laboratory 3 Numerical Solutions by Euler and Improved Euler Methods

(scalar equations)

In this session we look at basic numerical methods to help us understand the fundamentals of numerical approximations. Our objective is as follows.

1. Implement Euler’s method as well as an improved version to numerically solve an IVP.

2. Compare the accuracy and efficiency of the methods with methods readily available in MATLAB.

3. Apply the methods to specific problems and investigate potential pitfalls of the methods.

Instructions: For your lab write-up follow the instructions of LAB 1.

Euler’s Method

To derive Euler’s method start from y(t0) = y0 and consider a Taylor expansion at t1 = t0 + h:

y(t1) = y(t0) + y ′(t0)(t1 − t0) + . . .

= y0 + hf(t0, y(t0)) + . . .

= y0 + hf(t0, y0) + . . .

For small enough h we get an approximation y1 for y(t1) by suppressing the . . ., namely

y1 = y0 + hf(t0, y0) (L3.1)

The iteration (L3.1) is repeated to obtain y2 ≃ y(t2), . . . such that

yn+1 = yn + hf(tn, yn) tn+1 = tn + h

Geometrically, the approximation made is equivalent to replacing the solution curve by the tangent line at (t0, y0). From the figure we have

f(t0, y0) = f(t0, y(t0)) = y ′(t0) = tan θ =

y1 − y0 h

,

from which (L3.1) follows.

.

........... ........... ........... ..........s

s y0

y1

y(t1)

t0 t1

θ

h

� � � � � �

As an example consider the IVP

y′ = 2y = f(t, y) with y(0) = 3.

Note that here f does not explicitly depend on t (the ODE is called autonomous), but does implicitly through y = y(t). To apply Euler’s method we start with the initial condition and select a step size h. Since we are constructing arrays t and y without dimensionalizing them first it is best to clear these names in case they have been used already in the same MATLAB work session.

>> clear t y % no comma between t and y! type help clear for more info

>> y(1)=3; t(1)=0; h=0.1;

c⃝2011 Stefania Tracogna, SoMSS, ASU

MATLAB sessions: Laboratory 3

Since f is simple enough we may use the inline syntax:

>> f=inline(’2*y’,’t’,’y’)

f =

Inline function:

f(t,y) = 2*y

Note that the initialization y(1)=3 should not be interpreted as “the value of y at 1 is 3”, but rather “the first value in the array y is 3”. In other words the 1 in y(1) is an index, not a time value! Unfortunately, MATLAB indices in arrays must be positive (a legacy from Fortran...). The successive approximations at increasing values of t are then obtained as follows:

>> y(2)=y(1)+h*f(t(1),y(1)), t(2)=t(1)+h,

y =

3.0000 3.6000

t =

0 0.1000

>> y(3)=y(2)+h*f(t(2),y(2)), t(3)=t(2)+h,

y =

3.0000 3.6000 4.3200

t =

0 0.1000 0.2000

The arrays y and t are now 1× 3 row arrays. The dimension increases as new values of y and t are computed. It is obviously better to implement these steps in a for loop.

y(1)=3; t(1)=0; h = 0.1;

for n = 1:5

y(n+1)= y(n)+h*f(t(n),y(n));

t(n+1) = t(n)+h;

end

Note that the output in each command has been suppressed (with a ;). The list of computed y values vs t values can be output at the end only by typing:

>> [t(:),y(:)] % same as [t’,y’] here

ans =

0 3.0000

0.1000 3.6000

0.2000 4.3200

0.3000 5.1840

0.4000 6.2208

0.5000 7.4650

The next step is to write a function file that takes in input the function defining the ODE, the time span [t0, tfinal], the initial condition and the number of steps used. Note that, if we are given the number

of steps used, N , then the stepsize h can be easily computed using h = tfinal − t0

N .

The following function implements these ideas.

euler.m

function [t,y] = euler(f,tspan,y0,N)

% Solves the IVP y’ = f(t,y), y(t0) = y0 in the time interval tspan = [t0,tf]

% using Euler’s method with N time steps.

% Input:

c⃝2011 Stefania Tracogna, SoMSS, ASU

MATLAB sessions: Laboratory 3

% f = name of inline function or function M-file that evaluates the ODE

% (if not an inline function, use: euler(@f,tspan,y0,N))

% For a system, the f must be given as column vector.

% tspan = [t0, tf] where t0 = initial time value and tf = final time value

% y0 = initial value of the dependent variable. If solving a system,

% initial conditions must be given as a vector.

% N = number of steps used.

% Output:

% t = vector of time values where the solution was computed

% y = vector of computed solution values.

m = length(y0);

t0 = tspan(1);

tf = tspan(2);

h = (tf-t0)/N; % evaluate the time step size

t = linspace(t0,tf,N+1); % create the vector of t values

y = zeros(m,N+1); % allocate memory for the output y

y(:,1) = y0’; % set initial condition

for n=1:N

y(:,n+1) = y(:,n) + h*f(t(n),y(:,n)); % implement Euler’s method

end

t = t’; y = y’; % change t and y from row to column vectors

end

Remark: You should notice that the code above is slightly different from the first one we wrote (in particular, note the use of “:” when creating the output y). Although the two codes are equivalent in the scalar case, only the second one will work also for systems of Differential Equations.

We can implement the function with, say, N = 50 by typing the following commands:

>> [t,y] = euler(f,[0,.5],3,50); % use @f if defined in separate function

>> [t,y]

ans =

0 3.0000

0.0100 3.0600

0.0200 3.1212

0.0300 3.1836

0.0400 3.2473

0.0500 3.3122

: :

0.4500 7.3136

0.4600 7.4598

0.4700 7.6090

0.4800 7.7612

0.4900 7.9164

0.5000 8.0748

An even longer output is obtained for larger values of N . To compare the two approximations with N = 5 and N = 50 we plot both approximations on the same figure, together with the exact solution y(t) = 3e2t. Note that we can no longer call the output simply [t,y] because every time we call the function we will lose the output from the previous computations. Thus we will call [t5,y5] the output from Euler with N = 5 and [t50,y50] the output with N = 50. The exact solution of the IVP is y = 3e2t. We will store the exact solution in the vector y.

>> [t5,y5] = euler(f,[0,.5],3,5); % use @f if defined in separate function

c⃝2011 Stefania Tracogna, SoMSS, ASU

MATLAB sessions: Laboratory 3

>> [t50,y50] = euler(f,[0,.5],3,50);

>> t = linspace(0,.5,100); y = 3*exp(2*t); % evaluate the exact solution

>> plot(t5,y5,’ro-’,t50,y50,’bx-’,t,y,’k-’); axis tight;

>> legend(’Euler N = 5’,’Euler N = 50’,’Exact’,2);

0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 3

3.5

4

4.5

5

5.5

6

6.5

7

7.5

8 Euler N = 5 Euler N = 50 Exact

Figure L3a: Euler’s method applied to y′ = 2y, y(0) = 3 with step h = 0.1 and h = 0.01, compared to the exact solution.

IMPORTANT REMARK When using 5 steps of size h = 0.1 the approximation of the exact value y(0.5) = 3e ≃ 8.1548 is stored in y5(6). On the other hand when using 50 intervals of size h = 0.01, it is stored in y50(51). To avoid confusion we can reference these values by y5(end) and y50(end), respectively. The exact value of the function at t = 0.5 can be retrieved using y(end).

We can find the error (exact - approximation) at t = 0.5 by entering

>> e5 = y(end) - y5(end) % error with N = 5

e5 =

0.6899

>> e50 = y(end) - y50(end) % error with N = 50

e50 =

0.0801

and we can calculate the ratio of the errors to see how the approximation improves when the number of steps increases by a factor 10

>> ratio = e5/e50

ratio =

8.6148

EXERCISES

1. (a) If you haven’t already done so, enter the following commands:

>> f=inline(’2*y’,’t’,’y’);

>> t=linspace(0,.5,100); y=3*exp(2*t); % defines the exact solution of the ODE

>> [t50,y50]=euler(f,[0,.5],3, 50); % solves the ODE using Euler with 50 steps

Determine the Euler’s approximation for N = 500 and N = 5000 and fill in the following table with the values of the approximations, errors and ratios of consecutive errors at t = 0.5.

c⃝2011 Stefania Tracogna, SoMSS, ASU

MATLAB sessions: Laboratory 3

Some of the values have already been entered based on the computations we did above. Include the table in your report, as well as the MATLAB commands used to find the entries.

N approximation error ratio

5 7.4650 0.6899

50 0.0801 8.6148

500

5000

(b) Examine the last column. How does the ratio of consecutive errors relate to the number of steps used? Your answer to this question should confirm the fact that Euler’s method is of “ order h ”, that is, every time the stepsize is decreased by a factor k, the error is also reduced (approximately) by the same factor.

(c) Recall the geometrical interpretation of Euler’s method based on the tangent line. Using this geometrical interpretation, can you explain why the Euler approximations underestimate the solution in this particular example?

2. Consider the IVP y′ = −2y, y(0) = 3

for 0 ≤ t ≤ 10. The exact solution of the ODE is y = 3e−2t. The goal of this exercise is to visualize how Euler’s method is related to the slope field of the differ- ential equation. In order to do this we will plot the direction field together with the approximations and the exact solution.

(a) To plot the slope field we will use the MATLAB commands meshgrid and quiver. Enter the following commands:

t = 0:.45:10; y = -30:6:42 ; % define grid of values in t and y direction

[T,Y]=meshgrid(t,y); % creates 2d matrices of points in the ty-plane

dT = ones(size(T)); % dt=1 for all points

dY = -2*Y; % dy = -2*y; this is the ODE

quiver(T,Y,dT,dY) % draw arrows (t,y)->(t+dt, t+dy)

axis tight % adjust look

hold on

After running these commands you should get the graph of the slope field.

(b) Use linspace to generate a vector of 200 t-values between 0 and 10. Evaluate the solution y at these t-values and plot it in black together with the direction field (use ’linewidth’,2).

(c) Enter the function defining the ODE as inline. Use euler.m with N = 8 to determine the approximation to the solution. Plot the approxi- mated points in red (use ’ro-’,’linewidth’,2 as line-style in your plot), together with the exact solution and the direction field.

You should get Figure L3b. Based on the slope field and the geometrical meaning of Euler’s method explain why the approximations are so inaccurate for this particular value of the stepsize.

(d) Open a new figure by typing figure. Plot again the direction field but in a different window: t = 0:.4:10; y = -1:0.4:3; Repeat part (b) and repeat part (c) but this time with N = 16. You should get Figure L3c. Comment on the result from a geometrical point of view.

c⃝2011 Stefania Tracogna, SoMSS, ASU

MATLAB sessions: Laboratory 3

0 2 4 6 8 10 −30

−20

−10

0

10

20

30

40

Euler’s approximation with N=8

Figure L3b: Euler’s method applied to y′ = −2y, y(0) = 3 N = 8, compared to the exact solution.

0 2 4 6 8 10 −1

−0.5

0

0.5

1

1.5

2

2.5

3 Euler’s approximation with N=16

Figure L3c: Euler’s method applied to y′ = −2y, y(0) = 3 N = 16, compared to the exact solution.

Note: Brief explanations of the commands quiver and meshgrid are included in Appendix A. In Appendix B we describe the Graphical User Interface dfield8 for plotting slope fields.

Improved Euler’s Method

The improved Euler’s method is obtained by using an average of f values, i.e., a slope average:

y(t1) ≃ y(t0) + h

2 (f(t0, y(t0)) + f(t1, y(t1))) (L3.2)

Then substitute y(t0) = y0 and define y1 ≃ y(t1):

y1 = y0 + h

2 (f(t0, y0) + f(t1, y1)) . (L3.3)

The equation (L3.3) defines the trapezoidal method. Unfortunately, this formula defines y1 only implic- itly, i.e., y1 appears on both sides of the equality so that an equation must be solved to obtain y1. To avoid this problem, and since we already have made an approximation to get (L3.3), we replace y1 on the right-hand side by the approximation one would obtain by simply applying Euler’s method from

c⃝2011 Stefania Tracogna, SoMSS, ASU

MATLAB sessions: Laboratory 3

(t0, y0). The resulting quantity

y1 = y0 + h

2

( f(t0, y0) + f(t1, y0 + hf(t0, y0)︸ ︷︷ ︸

Euler y1 from (L3.1)

)

) (L3.4)

with t1 = t0 + h is the improved Euler approximation. This approximation can be thought of as a correction to the Euler approximation. The iteration (L3.4) is then repeated to obtain y2 ≃ y(t2), . . ., i.e.,

f1 = f(tn, yn) f2 = f(tn + h, yn + hf1)

yn+1 = yn + h 2 (f1 + f2)

tn+1 = tn + h

Minor modifications are made to the function euler.m to implement the improved Euler method.

EXERCISES

3. Modify the M-file euler.m to implement the algorithm for Improved Euler. Call the new file impeuler.m (include the file in your report). Test your code for the ODE y′ = 2y, y(0) = 3. First enter the function f(t, y) = 2y as inline function and then enter the following:

>> [t5,y5] = impeuler(f,[0,.5],3,5); % use @f if defined in separate function

>> [t5,y5]

ans =

0 3.0000

0.1000 3.6600

0.2000 4.4652

0.3000 5.4475

0.4000 6.6460

0.5000 8.1081

Compare your output with the one above. Note that the improved Euler approximation with 5 steps is already more accurate than the Euler approximation with 50 steps! (hence the “improved”)

4. Consider the IVP y′ = 2y, y(0) = 3 0 ≤ t ≤ 0.5

(a) Determine the Improved Euler’s approximation for N = 50, N = 500 and N = 5000. Fill in the following table with the values of the approximations, errors and ratios of consecutive errors at t = 0.5. One of the values has already been entered based on the computations we did above. Recall that the exact solution to the ODE is y = 3e2t. Include the table in your report, as well as the MATLAB commands used to find the entries.

N approximation error ratio

5 8.1081

50

500

5000

c⃝2011 Stefania Tracogna, SoMSS, ASU

MATLAB sessions: Laboratory 3

(b) Examine the last column. How does the ratio of the errors relate to the number of steps used? Your answer to this question should confirm the fact that Improved Euler’s method is of “ order h2 ”, that is, every time the stepsize is decreased by a factor k, the error is reduced (approximately) by a factor of k2.

5. Repeat Problem 2 for Improved Euler. Compare the results with the ones obtained with Euler’s method.

c⃝2011 Stefania Tracogna, SoMSS, ASU

MATLAB sessions: Laboratory 3

APPENDIX A: The commands meshgrid and quiver

The function meshgrid: This command is especially important because it is used also for plotting 3D surfaces. Suppose we want to plot a 3D function z = x2 − y2 over the domain 0 ≤ x ≤ 4 and 0 ≤ y ≤ 4. To do so, we first take several points on the domain, say 25 points. We can create two matrices X and Y , each of size 5×5, and write the xy- coordinates of each point in these matrices. We can then evaluate z at these xy values and plot the resulting surface. Creating the two matrices X and Y is much easier with the meshgrid command:

x = 0:4; % create a vector x = [0, 1, 2, 3, 4]

y = -4:2:4; % create a vector y = [-4, -2, 0, 2, 4]

[X, Y] = meshgrid(x,y); % create a grid of 25 points and store

% - their coordinates in X and Y

Try entering plot(X,Y,’o’,’linewidth’,2) to see a picture of the 25 points. Then, just for fun try entering

Z = X.^2 - Y.^2;

surf(X,Y,Z);

to get a picture of the surface representing z = x2 − y2. (The graph will be quite rough because of the limited number of points, but you should nonetheless recognize a “saddle”). Thus, in our code for the direction field, the meshgrid command is used to generate the points in the xy plane at which we want to draw the arrows representing the slope or tangent line to the solution at the point. Type help meshgrid for more information.

The function quiver: The command quiver(X,Y,U,V) draws vectors (arrows) specified by U and V at the points specified by X and Y . In our examples the vectors drawn are (dT, dY) where dT =1 and dY = -2*y, thus these vectors are in the direction of the slope at the given point. The arrows are automatically scaled so as not to overlap. The matrices X and Y are built using the meshgrid command as explained above. Type help quiver for more explanation.

APPENDIX B: Plotting direction fields with dfield8

The Graphical User Interface dfield8 from J. Polking (available at http://math.rice.edu/∼dfield) can also be used to plot the direction field and selected solution curves. After downloading the file, enter dfield8 in the command window. A new window will pop up, one that looks like Figure L3d

Figure L3d: The dfield8 setup GUI

c⃝2011 Stefania Tracogna, SoMSS, ASU

MATLAB sessions: Laboratory 3

ODEs with up to two parameters can be entered in the Setup window. For instance, we can enter the ODE y′ = −2y, identify the independent variable as t and choose the display window’s dimension (0 ≤ t ≤ 10 and −4 ≤ y ≤ 4).

Click on Proceed to see the slope field. A sample display window is shown in Figure L3e.

Figure L3e: The dfield8 direction field plot

Solutions curves can be drawn by clicking on the display window. Specific initial conditions can be entered by clicking on Options ->Keyboard Input. A new window will arise prompting the user to input initial values for the independent and dependent variables as an initial condition for the IVP.

More information about various options are available at http://math.rice.edu/∼dfield.

c⃝2011 Stefania Tracogna, SoMSS, ASU