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}\]


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

Example of line equations

Image: Example of line equations

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, ht0 and y0.

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 t0, y0 represent the solution (y0) of the differential equation at a given time (t0). Usually t0 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 ti and the function approximation wi 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 w0 is equal with the initial solution y0.

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

Step 5. If i = N, the algorithm is complete and wi 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 A0 is known, it has the coordinates (t0, y0). The point A1 is calculated based on the point A0 and the slope f(t,y). The next points An are calculated based on the previous points An-1 and the slope. As we seen in the line equation, the slope is equal with the differential of y(t).

Graphical representation of Euler integration method

Image: Graphical representation of Euler integration method

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 approximation and line equation

Image: Euler approximation and line equation

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, ht0 and y0.

N = 10
a = 1
b = 2
h = 0.1
t0 = 1
y0 = -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:

StepTimeSlopeEuler approximationExact solutionAbsolute error
iti = t0 + h·if(ti-1, wi-1)wi = wi-1 + h·f(ti-1, wi-1)yi = -1/ti|yi – wi|

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.

#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]);
            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;

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);
            w(i) = w(i-1) + h * f(t(i-1),w(i-1));
z = [[t0 t']',[w0 w']']

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

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

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.

a = 1;
b = 2;
N = 10;
t0 = 1;
w0 = -1;
outEuler = eulerODE1(a,b,N,t0,w0);
time = outEuler(:,1);
legend('Euler approximation','Exact solution',2)

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

Euler integration method - example 1 (10 steps)

Image: Euler integration method – example 1 (10 steps)

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.

Euler integration method - example 1 (30 steps)

Image: Euler integration method – example 1 (30 steps)

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.

Leave a Reply

Ad Blocker Detected

Dear user, Our website provides free and high quality content by displaying ads to our visitors. Please support us by disabling your Ad blocker for our site. Thank you!