MATHEMATICS

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.