Euler integration method for solving differential equations

In mathematics there are several types of ordinary differential equations (ODE), like linear, separable, or exact differential equations, which are solved analytically, giving an exact solution. This means that there is a specific method to be applied in order to extract a general exact solution.

For example:
\[\frac{dy}{dx} = x \tag{1}\]

is a first order separable differential equation, which has the exact solution:

\[y(x) = \frac{x^2}{2} + C \tag{2}\]

where C is a constant.

In practice, most of the differential equation do not have a standard form and can not be solved with analytic methods, which means we can not find a general solution y(x). This is the case for most of the differential equations derived from physical models (mechanical, electrical, thermal, etc.).

In this case, we need to use numerical methods to be able to determine the solution of the differential equation. Bear in mind that with numerical methods:

• we get an approximation of the solution, not the exact solution
• the solution is calculated incrementally, step by step

One of the simplest integration method is the Euler integration method, named after the mathematician Leonhard Euler. The Euler method is a first-order method, which means that the local error (error per step) is proportional to the square of the step size, and the global error (error at a given time) is proportional to the step size.

The Euler integration method is also an explicit integration method, which means that the state of a system at a later time (next step) is calculated from the state of the system at the current time (current step).

\[y(t + \Delta t) = f(y(t)) \tag{3}\]

The Euler integration method is also called the polygonal integration method, because it approximates the solution of a differential equation with a series of connected lines (polygon).

Line equation

In order to have a better understanding of the Euler integration method, we need to recall the equation of a line:

\[y = m \cdot x + n \tag{4}\]

where:

m – is the slope of the line
n – is the offset
(x,y) – coordinates

If we apply a differentiation to the line equation (4), we get:

\[\frac{dy}{dx} = m \tag{5}\]

which means that the slope m of the line is equal with the differential of y(x).

Euler method

The Euler method gives an approximation for the solution of the differential equation:

\[\frac{dy}{dt} = f(t,y) \tag{6}\]

with the initial condition:

\[y(t_0) = y_0 \tag{7}\]

where t is continuous in the interval [a, b].

The Euler algorithm for differential equations integration is the following:

Step 1. Define the integration start parameters: N, a, b, h, t₀ and y₀.

N is the number of integration steps, it is defined by the user (e.g 10, 100, etc.). For a fixed integration interval, the higher the number of integration steps, the better the approximation of the exact solution. Very high steps implies high computing power. Usually there must be a compromise between the accuracy and the time taken to solve the integration.

a and b are the start and end of the integration interval. If, for example we want to approximate the solution of a differential equation between 0 and 1, then a = 0 and b = 1.

The size of the interval and the number of integration steps define the integration step size h. The smaller the step size, the better the approximation, the smaller the integration error. It is possible to directly define the step size, which will further determine the number of integration steps. The step size h is calculated as:

\[h = \frac{b-a}{N} \tag{8}\]

The initial conditions t₀, y₀ represent the solution (y₀) of the differential equation at a given time (t₀). Usually t₀ is equal with the start value of the integration interval a.

Step 2. Initialise the calculation loop index i = 1.

Step 3. (Loop) Calculate the function argument t_i and the function approximation w_i as:

\[ \begin{split}
t_{i} &= t_{0} + h \cdot i\\
w_{i} &= w_{i-1} + h \cdot f(t_{i-1}, w{i-1})
\end{split} \]

Note that the initial function approximation w₀ is equal with the initial solution y₀.

Step 4. If i < N, increment i = i + 1 and repeat Step 3.

Step 5. If i = N, the algorithm is complete and w_i will be the approximation of the solution y(t), for i = 1, 2, … N.

In each step of the iteration, the Euler approximation calculate the end point of a line. The starting point A₀ is known, it has the coordinates (t₀, y₀). The point A₁ is calculated based on the point A₀ and the slope f(t,y). The next points A_n are calculated based on the previous points A_n-1 and the slope. As we seen in the line equation, the slope is equal with the differential of y(t).

The is a direct link between the Euler approximation used in Step 3 and the line equation. In the picture below is depicted where every parameter of the line equation is found in the Euler approximation. This image summarises quite well how the Euler approximation (integration method) works.

Euler integration method example

Let’s apply the Euler integration and solve the following ordinary differential equation:

\[ \begin{split}
& \frac{dy}{dt} = \frac{1}{t^2}-\frac{y}{t}-y^2\\
& y(1) = -1\\
& 1 \leq t \leq 2
\end{split} \]

The Euler approximation must be performed in 10 and 30 steps. The exact solution of the equation is:

\[y = – \frac{1}{t}\]

We will use the exact solution to compare against the Euler approximation. For a better understanding, we are going to apply the method step-by-step (manual) and also using a Scilab and a C script.

Step-by-step (manual) method

First, we’ll define the integration start parameters: N, a, b, h, t₀ and y₀.

N = 10
a = 1
b = 2
h = 0.1
t₀ = 1
y₀ = -1

Second, we’ll write the expression of the slope f(t,w):

\[f(t,w) = \frac{1}{t^2}-\frac{w}{t}-w^2\]

The iteration loop is going to be done in the table below:

Step	Time	Slope	Euler approximation	Exact solution	Absolute error
i	ti = t0 + h·i	f(ti-1, wi-1)	wi = wi-1 + h·f(ti-1, wi-1)	yi = -1/ti	|yi – wi|
1	1.1	1.0000000	-0.9000000	-0.9090909	0.0090909
2	1.2	0.8346281	-0.8165372	-0.8333333	0.0167961
3	1.3	0.7081591	-0.7457213	-0.7692308	0.0235095
4	1.4	0.6092475	-0.6847965	-0.7142857	0.0294892
5	1.5	0.5303982	-0.6317567	-0.6666667	0.0349100
6	1.6	0.4664990	-0.5851068	-0.6250000	0.0398932
7	1.7	0.4139668	-0.5437101	-0.5882353	0.0445252
8	1.8	0.3702295	-0.5066872	-0.5555556	0.0488684
9	1.9	0.3334030	-0.4733469	-0.5263158	0.0529689
10	2.0	0.3020810	-0.4431388	-0.5000000	0.0568612

As expected, there is a difference between the Euler approximation and the exact solution. The absolute error increases in each iteration step because at every step, the current approximation (i) is based on a previous approximation (i-1), which also has an error.

C script

The Euler method can be defined in any programming language. Below you can see the implementation in a C code.

#include<stdio.h>
#define N 10
float f (float, float);

int main(void)

{
    float a=1;
    float b=2;
    float t0=1;
    float w0=-1;
    float h=0;
    int i=0;
    float w[N]={0};
    float t[N]={0};
    
    h = (b-a)/N;
    printf("t: %f\t w: %f\n",t0,w0);
    
    for (i=0; i<N; i++)
    {
        t[i] = t0 + h * (i+1);
        if (i==0)
        {
            w[i] = w0 + h * f(t0,w0);
            printf("t: %f\t w: %f\n",t[i],w[i]);
        }
        else
        {
            w[i] = w[i-1] + h * f(t[i-1],w[i-1]);
            printf("t: %f\t w: %f\n",t[i],w[i]);
        }
    }
    return 0;
}

float f (float t, float y)
{
    float z;
    z = 1/(t*t) - y/t - y*y;
    return z;
}

Running the executable will output the following results:

t: 1.000000  w: -1.000000
t: 1.100000  w: -0.900000
t: 1.200000  w: -0.816537
t: 1.300000  w: -0.745721
t: 1.400000  w: -0.684796
t: 1.500000  w: -0.631757
t: 1.600000  w: -0.585107
t: 1.700000  w: -0.543710
t: 1.800000  w: -0.506687
t: 1.900000  w: -0.473347
t: 2.000000  w: -0.443139 

Scilab script

Using Scilab is a very easy and flexible way to experiment different integration step sizes and also give the possibility to plot the results. For an easier understanding we are going to define several Scilab function (*.sci), for:

• the slope function f(t,y)
• the Euler approximation w(t,y)
• the exact solution y(t)

The Scilab function for the slope function is going to be defined in a file f.sci, with the following content:

function [z]=f(t,y)
    z = 1/t^2 - y/t - y^2;
endfunction

The Scilab function for the Euler approximation is going to be defined in a file eulerODE1.sci, with the following content:

function [z]=eulerODE1(a,b,N,t0,w0)
    h = (b-a)/N;
    for i=1:N
        t(i) = t0 + h * i;
        if i==1
            w(i) = w0 + h * f(t0,w0);
        else
            w(i) = w(i-1) + h * f(t(i-1),w(i-1));
        end
    end
z = [[t0 t']',[w0 w']']
endfunction

The exact solution is going to be defined in a file y.sci, with the following content:

function [z]=y(t)
    z = -1./t;
endfunction

In a Scilab script, in our case named runEuler.sce, we are going to define the initial parameters of the Euler integration, call the *.sci functions and plot the results.

clear
exec('f.sci')
exec('eulerODE1.sci')
exec('y.sci')
clc
a = 1;
b = 2;
N = 10;
t0 = 1;
w0 = -1;
outEuler = eulerODE1(a,b,N,t0,w0);
time = outEuler(:,1);
plot(time,outEuler(:,2),'rs-')
plot(time,y(time),'bs-')
xlabel('t','FontSize',2)
ylabel('y(t)','FontSize',2)
legend('Euler approximation','Exact solution',2)
title('x-engineer.org','FontSize',2)
xgrid()

After running the script, the following plot is being generated:

It’s clearly visible that there is a significand difference between the exact solution of the differential equation and its Euler approximation. In order to improve the results, we are going to increase the number the integration steps N = 30 and run the Scilab script again.

As expected, the error between the exact solution and the Euler approximation is reduced. This is achieved by having 3 times more integration steps, which means more calculation power.

The way the Euler integration method is implemented in Scilab is very flexible. In order to try out other functions, all you need to do is change the f.sci and y.sci files with the corresponding functions, change the initial parameters and run the *.sce file.