# MATHEMATICS

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