Euler method

In mathematics and computational science, the Euler method (also called the forward Euler method) is a first-order numerical procedure for solving ordinary differential equations (ODEs) with a given initial value. It is the most basic explicit method for numerical integration of ordinary differential equations and is the simplest Runge–Kutta method. The Euler method is named after Leonhard Euler, who first proposed it in his book Institutionum calculi integralis (published 1768–1770).^[1]

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 method often serves as the basis to construct more complex methods, e.g., predictor–corrector method.

Given the initial value problem

$${\displaystyle y'=y,\quad y(0)=1,}$$

we would like to use the Euler method to approximate ${\displaystyle y(4)}$.^[5]

Using step size equal to 1 (h = 1)

The Euler method is

$${\displaystyle y_{n+1}=y_{n}+hf(t_{n},y_{n}).}$$

so first we must compute ${\displaystyle f(t_{0},y_{0})}$. In this simple differential equation, the function ${\displaystyle f}$ is defined by ${\displaystyle f(t,y)=y}$. We have

$${\displaystyle f(t_{0},y_{0})=f(0,1)=1.}$$

By doing the above step, we have found the slope of the line that is tangent to the solution curve at the point ${\displaystyle (0,1)}$. Recall that the slope is defined as the change in ${\displaystyle y}$ divided by the change in ${\displaystyle t}$, or ${\textstyle {\frac {\Delta y}{\Delta t}}}$.

The next step is to multiply the above value by the step size ${\displaystyle h}$, which we take equal to one here:

$${\displaystyle h\cdot f(y_{0})=1\cdot 1=1.}$$

Since the step size is the change in ${\displaystyle t}$, when we multiply the step size and the slope of the tangent, we get a change in ${\displaystyle y}$ value. This value is then added to the initial ${\displaystyle y}$ value to obtain the next value to be used for computations.

$${\displaystyle y_{0}+hf(y_{0})=y_{1}=1+1\cdot 1=2.}$$

The above steps should be repeated to find ${\displaystyle y_{2}}$, ${\displaystyle y_{3}}$ and ${\displaystyle y_{4}}$.

$${\displaystyle {\begin{aligned}y_{2}&=y_{1}+hf(y_{1})=2+1\cdot 2=4,\\y_{3}&=y_{2}+hf(y_{2})=4+1\cdot 4=8,\\y_{4}&=y_{3}+hf(y_{3})=8+1\cdot 8=16.\end{alineado}}}$$

Due to the repetitive nature of this algorithm, it can be helpful to organize computations in a chart form, as seen below, to avoid making errors.

The conclusion of this computation is that ${\displaystyle y_{4}=16}$. The exact solution of the differential equation is ${\displaystyle y(t)=e^{t}}$, so ${\displaystyle y(4)=e^{4}\approx 54.598}$. Although the approximation of the Euler method was not very precise in this specific case, particularly due to a large value step size ${\displaystyle h}$, its behaviour is qualitatively correct as the figure shows.

Using other step sizes

As suggested in the introduction, the Euler method is more accurate if the step size ${\displaystyle h}$ is smaller. The table below shows the result with different step sizes. The top row corresponds to the example in the previous section, and the second row is illustrated in the figure.

The error recorded in the last column of the table is the difference between the exact solution at ${\displaystyle t=4}$ and the Euler approximation. In the bottom of the table, the step size is half the step size in the previous row, and the error is also approximately half the error in the previous row. This suggests that the error is roughly proportional to the step size, at least for fairly small values of the step size. This is true in general, also for other equations; see the section Global truncation error for more details.

Other methods, such as the midpoint method also illustrated in the figures, behave more favourably: the global error of the midpoint method is roughly proportional to the square of the step size. For this reason, the Euler method is said to be a first-order method, while the midpoint method is second order.

We can extrapolate from the above table that the step size needed to get an answer that is correct to three decimal places is approximately 0.00001, meaning that we need 400,000 steps. This large number of steps entails a high computational cost. For this reason, higher-order methods are employed such as Runge–Kutta methods or linear multistep methods, especially if a high accuracy is desired.^[6]

For this third-order example, assume that the following information is given:

$${\displaystyle {\begin{aligned}&y'''+4ty''-t^{2}y'-(\cos {t})y=\sin {t}\\&t_{0}=0\\&y_{0}=y(t_{0})=2\\&y'_{0}=y'(t_{0})=-1\\&y''_{0}=y''(t_{0})=3\\&h=0.5\end{aligned}}}$$

A partir de esto podemos aislar y ''' para obtener la ecuación:

$${\displaystyle y'''=f{\left(t,y,y',y''\right)}=\sin {t}+\left(\cos {t}\right)y+t^{2}y'-4ty''}$$

Usando eso podemos obtener la solución para${\displaystyle {\vec {y}}_{1}}$: $${\displaystyle {\begin{aligned}{\vec {y}}_{1}&={\begin{pmatrix}y_{1}\\y_{1}'\\y_{1}''\end{pmatrix}}={\begin{pmatrix}y_{0}\\y'_{0}\\y''_{0}\end{pmatrix}}+h{\begin{pmatrix}y'_{0}\\y''_{0}\\f{\left(t_{0},y_{0},y'_{0},y''_{0}\right)}\end{pmatrix}}\\[1ex]&={\begin{pmatrix}2\\-1\\3\end{pmatrix}}+0.5\cdot {\begin{pmatrix}-1\\3\\\sin {0}+(\cos {0})\cdot 2+0^{2}\cdot (-1)-4\cdot 0\cdot 3\end{pmatrix}}\\[1ex]&={\begin{pmatrix}1.5\\0.5\\4\end{pmatrix}}\end{aligned}}}$$Y utilizando la solución para${\displaystyle {\vec {y}}_{1}}$, podemos obtener la solución para${\displaystyle {\vec {y}}_{2}}$:$${\displaystyle {\begin{aligned}{\vec {y}}_{2}&={\begin{pmatrix}y_{2}\\y_{2}'\\y_{2}''\end{pmatrix}}={\begin{pmatrix}y_{1}\\y'_{1}\\y''_{1}\end{pmatrix}}+h{\begin{pmatrix}y'_{1}\\y''_{1}\\f{\left(t_{1},y_{1},y'_{1},y''_{1}\right)}\end{pmatrix}}\\[1ex]&={\begin{pmatrix}1.5\\0.5\\4\end{pmatrix}}+0.5\cdot {\begin{pmatrix}0.5\\4\\\sin {0.5}+(\cos {0.5})\cdot 1.5+0.5^{2}\cdot 0.5-4\cdot 0.5\cdot 4\end{pmatrix}}\\[1ex]&={\begin{pmatrix}1.75\\2.5\\0.9604...\end{pmatrix}}\end{aligned}}}$$Podemos continuar este proceso utilizando la misma fórmula durante el tiempo que sea necesario para encontrar cualquiera.${\displaystyle {\vec {y}}_{i}}$deseado.

El método de Euler se puede derivar de varias maneras.

Esta línea de pensamiento puede continuarse para llegar a varios métodos lineales de múltiples pasos .

El error de truncamiento local del método de Euler es el error cometido en un solo paso. Es la diferencia entre la solución numérica después de un paso,${\displaystyle y_{1}}$y la solución exacta en ese momento${\displaystyle t_{1}=t_{0}+h}$La solución numérica viene dada por

$${\displaystyle y_{1}=y_{0}+hf(t_{0},y_{0}).}$$

Para obtener la solución exacta, utilizamos el desarrollo de Taylor mencionado en la sección Derivación anterior:

$${\displaystyle y(t_{0}+h)=y(t_{0})+hy'(t_{0})+{\tfrac {1}{2}}h^{2}y''(t_{0})+{\mathcal {O}}{\left(h^{3}\right)}.}$$

El error de truncamiento local (LTE) introducido por el método de Euler viene dado por la diferencia entre estas ecuaciones:

$${\displaystyle \mathrm {LTE} =y(t_{0}+h)-y_{1}={\tfrac {1}{2}}h^{2}y''(t_{0})+{\mathcal {O}}{\left(h^{3}\right)}.}$$

Este resultado es válido si${\displaystyle y}$tiene una tercera derivada acotada. ^{[ 10 ]}

Esto demuestra que para pequeños${\displaystyle h}$, el error de truncamiento local es aproximadamente proporcional a${\displaystyle h^{2}}$Esto hace que el método de Euler sea menos preciso que las técnicas de orden superior, como los métodos de Runge-Kutta y los métodos lineales de pasos múltiples , para los cuales el error de truncamiento local es proporcional a una potencia mayor del tamaño del paso.

Se puede obtener una formulación ligeramente diferente para el error de truncamiento local utilizando la forma de Lagrange para el término restante en el teorema de Taylor . Si${\displaystyle y}$tiene una segunda derivada continua, entonces existe una${\displaystyle \xi \in [t_{0},t_{0}+h]}$de tal manera que ^{[ 11 ]}

$${\displaystyle \mathrm {LTE} =y(t_{0}+h)-y_{1}={\tfrac {1}{2}}h^{2}y''(\xi ).}$$

En las expresiones anteriores para el error, la segunda derivada de la solución exacta desconocida${\displaystyle y}$puede ser reemplazado por una expresión que involucre el lado derecho de la ecuación diferencial. De hecho, se deduce de la ecuación${\displaystyle y'=f(t,y)}$que ^{[ 12 ]}

$${\displaystyle y''(t_{0})={\frac {\partial f}{\partial t}}\left(t_{0},y(t_{0})\right)+{\frac {\partial f}{\partial y}}\left(t_{0},y(t_{0})\right)\,f\left(t_{0},y(t_{0})\right).}$$

El error de truncamiento global es el error en un tiempo fijo.${\displaystyle t_{i}}$, después de la cantidad de pasos que el método necesite tomar para llegar a ese tiempo desde el tiempo inicial. El error de truncamiento global es el efecto acumulativo de los errores de truncamiento locales cometidos en cada paso. ^{[ 13 ]} El número de pasos se determina fácilmente como${\textstyle {\frac {t_{i}-t_{0}}{h}}}$, que es proporcional a${\textstyle {\frac {1}{h}}}$y el error cometido en cada paso es proporcional a${\displaystyle h^{2}}$(véase la sección anterior). Por lo tanto, cabe esperar que el error de truncamiento global sea proporcional a${\displaystyle h}$. ^{[ 14 ]}

Este razonamiento intuitivo se puede precisar. Si la solución${\displaystyle y}$tiene una segunda derivada acotada y${\displaystyle f}$es Lipschitz continua en su segundo argumento, entonces el error de truncamiento global (denotado como${\displaystyle |y(t_{i})-y_{i}|}$) está delimitado por

$${\displaystyle \left|y(t_{i})-y_{i}\right|\leq {\frac {hM}{2L}}\left(e^{L(t_{i}-t_{0})}-1\right)}$$

dónde${\displaystyle M}$es una cota superior para la segunda derivada de${\displaystyle y}$en el intervalo dado y${\displaystyle L}$es la constante de Lipschitz de${\displaystyle f}$. ^{[ 15 ]} O más simplemente, cuando${\displaystyle y'(t)=f(t,y)}$, el valor${\textstyle L={\text{max}}\left(|{\frac {d}{dy}}\left[f(t,y)\right]|\right)}$(de tal manera que${\displaystyle t}$se trata como una constante). En cambio,${\textstyle M=\max \left(\left|{\frac {d^{2}}{dt^{2}}}\left[y(t)\right]\right|\right)}$donde función${\displaystyle y(t)}$es la solución exacta que solo contiene la${\displaystyle t}$variable.

La forma precisa de esta cota tiene poca importancia práctica, ya que en la mayoría de los casos la cota sobreestima enormemente el error real cometido por el método de Euler. ^{[ 16 ]} Lo importante es que muestra que el error de truncamiento global es (aproximadamente) proporcional a${\displaystyle h}$Por esta razón, se dice que el método de Euler es de primer orden. ^{[ 17 ]}

The Euler method can also be numerically unstable, especially for stiff equations, meaning that the numerical solution grows very large for equations where the exact solution does not. This can be illustrated using the linear equation $${\displaystyle y'=-2.3y,\qquad y(0)=1.}$$ The exact solution is ${\displaystyle y(t)=e^{-2.3t}}$, which decays to zero as ${\displaystyle t\to \infty }$. However, if the Euler method is applied to this equation with step size ${\displaystyle h=1}$, then the numerical solution is qualitatively wrong: It oscillates and grows (see the figure). This is what it means to be unstable. If a smaller step size is used, for instance ${\displaystyle h=0.7}$, then the numerical solution does decay to zero.

If the Euler method is applied to the linear equation ${\displaystyle y'=ky}$, then the numerical solution is unstable if the product ${\displaystyle hk}$ is outside the region $${\displaystyle \left\{z\in \mathbf {C} \,{\big |}\,|z+1|\leq 1\right\},}$$ illustrated on the right. This region is called the (linear) stability region.^[18] In the example, ${\displaystyle k=-2.3}$, so if ${\displaystyle h=1}$ then ${\displaystyle hk=-2.3}$ which is outside the stability region, and thus the numerical solution is unstable.

This limitation — along with its slow convergence of error with ${\displaystyle h}$ — means that the Euler method is not often used, except as a simple example of numerical integration. Frequently models of physical systems contain terms representing fast-decaying elements (i.e. with large negative exponential arguments). Even when these are not of interest in the overall solution, the instability they can induce means that an exceptionally small timestep would be required if the Euler method is used.

In step ${\displaystyle n}$ of the Euler method, the rounding error is roughly of the magnitude ${\displaystyle \varepsilon y_{n}}$ where ${\displaystyle \varepsilon }$ is the machine epsilon. Assuming that the rounding errors are independent random variables, the expected total rounding error is proportional to ${\textstyle {\frac {\varepsilon }{\sqrt {h}}}}$.^[19] Thus, for extremely small values of the step size the truncation error will be small but the effect of rounding error may be big. Most of the effect of rounding error can be easily avoided if compensated summation is used in the formula for the Euler method.^[20]

A simple modification of the Euler method which eliminates the stability problems noted above is the backward Euler method: $${\displaystyle y_{n+1}=y_{n}+hf(t_{n+1},y_{n+1}).}$$ This differs from the (standard, or forward) Euler method in that the function ${\displaystyle f}$ is evaluated at the end point of the step, instead of the starting point. The backward Euler method is an implicit method, meaning that the formula for the backward Euler method has ${\displaystyle y_{n+1}}$ on both sides, so when applying the backward Euler method we have to solve an equation. This makes the implementation more costly.

Otras modificaciones del método de Euler que contribuyen a la estabilidad dan lugar al método de Euler exponencial o al método de Euler semiimplícito .

Los métodos más complejos pueden alcanzar un orden superior (y mayor precisión). Una posibilidad es utilizar más evaluaciones de la función. Esto se ilustra con el método del punto medio , que ya se menciona en este artículo: $${\displaystyle y_{n+1}=y_{n}+hf\left(t_{n}+{\tfrac {1}{2}}h,y_{n}+{\tfrac {1}{2}}hf(t_{n},y_{n})\right).}$$ Esto da lugar a la familia de métodos de Runge-Kutta .

La otra posibilidad es utilizar más valores pasados, como lo ilustra el método de dos pasos de Adams-Bashforth: $${\displaystyle y_{n+1}=y_{n}+{\tfrac {3}{2}}hf(t_{n},y_{n})-{\tfrac {1}{2}}hf(t_{n-1},y_{n-1}).}$$ Esto da lugar a la familia de métodos lineales multipaso . Existen otras modificaciones que utilizan técnicas de detección compresiva para minimizar el uso de memoria ^{[ 21 ].}

En la película Figuras ocultas , Katherine Johnson recurre al método de Euler para calcular la reentrada del astronauta John Glenn desde la órbita terrestre. ^{[ 22 ]}

• Método de Crank-Nicolson
• El descenso de gradiente también utiliza pasos finitos, en este caso para encontrar mínimos de funciones.
• Lista de métodos de Runge-Kutta
• Método lineal de pasos múltiples
• Integración numérica (para el cálculo de integrales definidas)
• Métodos numéricos para ecuaciones diferenciales ordinarias

• Contenido multimedia relacionado con el método de Euler en Wikimedia Commons.
• Implementaciones del método de Euler en diferentes lenguajes por Rosetta Code
• "Método de Euler" , Enciclopedia de Matemáticas , EMS Press , 2001 [1994]