In this activity we will learn how to implement Euler's method in MATLAB to numerically solve an ordinary differential equation given an initial value.

Use the MATLAB Live Editor to edit and run this Live Script in your browser or your desktop.

- Read each section carefully.
- Run the code within a section before starting to read the next.
- To run the code from each section, click into the code section and then click on the Run Section button (from the toolstrip) or click on the blue stripe on left side of that section as shown below:

Remark: Run the code of each section from top to bottom, otherwise you may get an error.

- The end of a section is indicated with a thin line, like the next one -

Euler's method is one of the simplest numerical methods for approximating solutions of first-order initial-value problems

with initial condition .

The idea behind this method is to find approximate values for the solution at equally spaced numbers

where is the step size. The differential equation tells us that the slope at is

so the approximate value of the solution when is

as it is shown in the Figure below.

Similarly, we find that . In general, we have

.

Consider the Initial Value Problem (IVP)

where

.

Let's compute ten values of with a step size . That is, we need to compute from to .

Step 1) First, we need to define two arrays to store the values of every and . This can be done with the following code:

t = zeros(1, 10+1);

y = zeros(1, 10+1);

Observe that we wrote 10+1 instead of 11 just to emphasize that we are calculating the next ten values, since we already know and , the initial conditions.

The function zeros() creates an array of all zeros. This is just a placeholder to store all the values of and .

Step 2) Now we need to define the initial conditions and the step size :

t(1) = 0;

y(1) = 1;

dt = 0.2;

Remark 1: Note that we start with y(1) = 0, whereas in the theoretical work we write y(0) = 0. MATLAB does not allow you to put an index of zero. Keep this in mind in case otherwise you will get an error. For more details see: Array Indexing.

Step 3) To compute all the values and we can use a for loop. For example, if we want to compute 10 values we write:

for k = 1:10

t(k+1) = t(k) + dt;

y(k+1) = y(k) + dt * (1 - exp(-4 * t(k)) - 2 * y(k));

end

Remark 2: Here we are using the colon operator ":" which is one of the most useful operators in MATLAB. It is used to create vectors, subscript arrays, and specify for iterations. In particular, 1:10 creates a row vector, containing integers from 1 to 10 with an increment of one.

Step 4) Finally, we plot the result:

plot(t, y, 'ro-')

Here we are using the command

plot(x, y, LineSpec)

where LineSpec is a string (i.e., some text surrounded by apostrophes (')) that specifies the line style, marker symbol and colour of the graph. Some other values you can use are the following:

In this case, 'ro-' creates the plot using the color Red, with marker Circle and line style Solid. For more details see: plot().

Now we are ready to compute and and plot them. Run this section to see the result of the following code:

t = zeros(1, 10+1);

y = zeros(1, 10+1);

t(1) = 0;

y(1) = 1;

dt = 0.2;

for k = 1:10

t(k+1) = t(k) + dt;

y(k+1) = y(k) + dt * (1 - exp(-4 * t(k)) - 2 * y(k));

end

plot(t, y,'ro-')

t(end) % or t(11)

y(end) % or y(11)

Remark 3: Observe that we added in the last two lines t(end) and y(end). This gives you the last values from each array. More details here: Last array index. You can easily access any other value by replacing "end" by any other integer from 1 to 10.

An anonymous function is a function that is not written in an external file (usually named M-file), rather is associated with a function handle. Anonymous functions work as regular functions, accepting inputs and returning outputs. As an example to create an anonymous function with the handle sqr which squares a number. For example, run this section to see the output of the following code:

sqr = @(x) x.^2;

sqr(5)

(1 - exp(-4 * t(k)) - 2 * y(k))

to define the ODE. This can be easily replaced by

f(t(k), y(k))

using the following anonymous function:

f = @(t,y) 1 - exp(-4 * t) - 2 * y;

clear % This command just removes variables to free memory

% Define anonymous function

f = @(t,y) 1 - exp(-4 * t) - 2 * y;

t = zeros(1, 10+1);

y = zeros(1, 10+1);

t(1) = 0;

y(1) = 1;

dt = 0.2;

for k = 1:10

t(k+1) = t(k) + dt;

y(k+1) = y(k) + dt * f(t(k), y(k)); % Here we use the anonymous function

end

plot(t, y, 'bo-')

Modify the anonymous function, try for example and re-run this section to see the result.

We can also implement Euler's method by using pre-defined funcitons. In order to do that we must create our own MATLAB function

In MATLAB we must declare a new function before using it in our code. In the case of Live Scripts, we can add it at the very end of our document, see Appendix for more details.

Using the pre-defined function named my_euler, we can easily compute as many values as we want with any given step size. In this case, my_euler function requires five inputs:

- f = anonymous function,
- t_0 = initial value of t,
- y_0 = initial value of y,
- dt = step size,
- t_f = final value of t.

Let's compute from to using the same ODE with a step size . Run this section to see the result of the following code:

clear % This command just removes variables to free memory

% Define initial conditions and step size

t_0 = 0;

y_0 = 1;

dt = 0.05;

% Define the last value 't' we want to compute

t_f = 2;

% Define anonymous function

f = @(t,y) 1 - exp(-4 * t) - 2 * y;

% Apply Eulers' method with pre-define function

my_euler(f, t_0, y_0, dt, t_f);

Play around the previous code. Modify some parameters and re-run this section to see the changes.

Let's practice what we just learned.

We know that the general solution of

.

is given by

.

Using the initial condition of the IVP, , we can easily determine the value and hence the particular solution for this IVP. (You should check this is true in your notebook)

Compare the general solution with Euler's method for different values of .

If you run this section, you will plot the general solution for over the domain . Compare the general solution with Euler's method for the following values of the step size:

You can use any of the methods described above. What do you notice? What do you wonder? Write your code in the indicated space below:

t = linspace(0, 2.5, 1000); % 1000 points between 0 and 2.5

y = 1/2 + exp(-4 * t)/2; % Define general solution

plot(t, y, '-k', 'LineWidth', 2) % plot general solution

hold on

%%% Your code must go after this line %%%

%%% Your code must go before this line %%%

hold off

Remark 3: To help us plot the general solution we are using the command linspace(x1, x2, n) which generates n points between x1 and x2 (including the end points). The spacing between them is

(x2-x1)/(n-1).

This command is similar to the colon operator ":", but gives direct control over the number of points and always includes the endpoints.

Now consider the IVP

, .

In your notebook, calculate the exact solution. Then use MATLAB to plot the general solution. Then use Euler's method to plot the approximate solutions in the following cases:

- step size and last value of .
- step size and last value of .

What is the approximate value of for ? And for ? Compare these values with those obtained using the exact solution.

Write your code here:

Hint: You can create anonymous functions for and for the general solution.

function [ tval, yval ] = my_euler(f, t_0, y_0, dt, t_f)

% 'my_euler' function inputs:

% f = anonymous function

% t_0 = initial value of t

% y_0 = initial value of y

% dt = step size

% t_f = final value of t

% 'my_euler' outputs:

% plots (t, y) points

% Clear the 't' and 'y' variables

clear t y

% Calculate how many steps it takes to get from t_0 to t_f

nsteps = ceil((t_f -t_0)/dt);

% To lear more about the ceil function check:

% https://au.mathworks.com/help/matlab/ref/ceil.html

% initialise the storage vectors.

t = zeros(1, nsteps + 1);

y = zeros(1, nsteps + 1);

t(1) = t_0;

y(1) = y_0;

for k = 1:nsteps

% calculate the next t-value

t(k+1) = t(k) + dt;

% calculate the corresponding approximation for y

y(k+1) = y(k) + dt * f(t(k), y(k));

end

% Plot the solution with blue points

% and connected by red segments

plot(t, y,'bo', t, y,'r-', 'LineWidth', 2)

xlabel('t')

ylabel('y')

% define the function outputs:

tval = t(end); % or t(nsteps+1)

yval = y(end); % or y(nsteps+1)

end