Metodo di Eulero

Da Wikipedia, l'enciclopedia libera.
bussola Disambiguazione – Se stai cercando l'omonimo metodo di fattorizzazione, vedi Metodo di fattorizzazione di Eulero.
Illustrazione del metodo di Eulero. La curva sconosciuta è in blu, e la sua approssimazione polinomiale in rosso.

In matematica e informatica, il Metodo di Eulero è una procedura numerica del primo ordine per risolvere equazioni differenziali ordinarie (ODE) fornito un valore iniziale. Questo è il più basilare dei metodi espliciti per l'integrazione numerica di equazioni differenziali ordinarie ed è il più semplice metodo Runge-Kutta. Prende il nome da Leonhard Euler, il quale lo espose nel suo libro Institutionum calculi integralis (pubblicato nel 1768–70).[1]

Il metodo di Eulero è un metodo del primo ordine, il che significa che errori locali (errori per step) sono proporzionali alla radice dell'incremento, e l'errore globale (errore ad un tempo fornito) è proporzionale al valore dell'incremento. Il metodo di Eulero è spesso utilizzato come base di metodi più complessi.

Descrizione geometrica non formale[modifica | modifica sorgente]

Si consideri il problema del calcolo di una forma di una curva sconosciuta che inizia in un certo punto e soddisfa una data equazione differenziale. In questo caso, una equazione differenziale può essere immaginata come una formula tramite cui il coefficiente angolare della tangente alla curva può essere calcolata in ogni punto della curva, una volta che la posizione di quel punto è stata calcolata.

L'idea è che mentre la curva è inizialmente sconosciuta, il suo punto iniziale, che chiamiamo A_0, è noto (vedi l'immagine in alto a destra). Quindi, dall'equazione differenziale, il coefficiente angolare della curva nel punto A_0 può essere calcolato, e quindi, anche la retta tangente.

Si prenda un piccolo incremento sulla retta tangente fino ad un punto A_1. In questo intervallo il coefficiente angolare non cambia molto, quindi A_1 sarà vicino alla curva. Se assumiamo che A_1 sia ancora sulla curva, possiamo utilizzare lo stesso ragionamento applicato al punto A_0. Dopo diversi passi, una curva poligonale A_0A_1A_2A_3\dots viene generata. In generale, questa curva non diverge troppo dalla curva originale sconosciuta, e l'errore tra le due curve può essere diminuito se l'incremento è piccolo abbastanza e l'intervallo di calcolo è finito.[2]

Formulazione del metodo[modifica | modifica sorgente]

Supponiamo di voler approssimare la soluzione del problema del valore iniziale

y'(t) = f(t,y(t)), \qquad \qquad y(t_0)=y_0.

Si scelga un valore h per la dimensione di ogni incremento e si definisca t_n = t_0 + nh. Ora, uno step del metodo di Eulero da t_n a t_{n+1} = t_n + h è[3]

 y_{n+1} = y_n + hf(t_n,y_n).

Il valore di y_n è un'approssimazione della soluzione della ODE al tempo t_n: y_n \approx y(t_n). Il metodo di Eulero è esplicito, ovvero la soluzione y_{n+1} è una funzione esplicita di y_i per i \leq n.

Il metodo di Eulero integra una equazione differenziale del primo ordine, ma ogni equazione differenziale di ordine N può essere rappresentata come una equazione differenziale del primo ordine: per modificare l'equazione

 y^{(N)}(t) = f(t, y(t), y'(t), \ldots, y^{(N-1)}(t)) ,

introduciamo variabili ausiliari z_1(t)=y(t), z_2(t)=y'(t),\ldots, z_N(t)=y^{(N-1)}(t) e otteniamo l'equazione equivalente

 \mathbf{z}'(t)
  = \begin{pmatrix} z_1'(t)\\ \vdots\\ z_{N-1}'(t)\\ z_N'(t) \end{pmatrix}
  = \begin{pmatrix} y'(t)\\ \vdots\\ y^{(N-1)}(t)\\ y^{(N)}(t) \end{pmatrix}
  = \begin{pmatrix} z_2(t)\\ \vdots\\ z_N(t)\\ f(t,z_1(t),\ldots,z_N(t)) \end{pmatrix}

Questo è un sistema del primo ordine nella variabile \mathbf{z}(t) e può essere approssimato col metodo di Eulero o con qualsiasi altro metodo per sistemi del primo ordine.[4]

Esempio[modifica | modifica sorgente]

Dato il valore iniziale del problema

y'=y, \quad y(0)=1,

vogliamo utilizzare il metodo di Eulero per approssimare y(4).[5]

Usando 1 come incremento (h = 1)[modifica | modifica sorgente]

Illustrazione della funzione y'=y, y(0)=1. In blu il risultato del metodo di Eulero; in verde il metodo del punto medio; in rosso, la funzione esatta, y=e^t. L'incremento è h = 1.0.

Il metodo di Eulero prevede:

 y_{n+1} = y_n + hf(t_n,y_n).  \qquad \qquad

Quindi innanzitutto dobbiamo calcolare f(t_0, y_0). In questa semplice equazione differenziale, la funzione f è definita da f(t,y) = y. Abbiamo

 f(t_0,y_0) = f(0,1) = 1. \qquad \qquad

Col passo precedente, abbiamo trovato il coefficiente angolare della retta tangente alla curva della soluzione al punto (0,1). Si ricordi che il coefficiente angolare è definito come l'incremento di y diviso per l'incremento di t, ovvero  \Delta y/ \Delta t.

Il passo successivo è moltiplicare il valore precedente per la dimensione dell'incremento h, che abbiamo definito essere uguale a 1:

 h \cdot f(y_0) = 1 \cdot 1 = 1. \qquad \qquad

Dato che la dimensione dell'incremento è la differenza di t, quando moltiplichiamo la dimensione dell'incremento e il coefficiente angolare della tangente, otteniamo la differenza del valore y. Questo valore è poi sommato al valore iniziale y per ottenere il prossimo valore da utilizzare nei calcoli.

 y_0 + hf(y_0) = y_1 = 1 + 1 \cdot 1 = 2. \qquad \qquad

Il passo precedente dev'essere ripetuto per trovare  y_2, y_3 e y_4.

 \begin{align}
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{align}

A causa della natura ripetitiva di questo algoritmo, può essere utile organizzare i calcoli sotto forma di grafico, come mostrato sotto, per evitare di commettere errori.

n y_n t_n f(t_n,y_n) h \Delta y y_{n+1}
0 1 0 1 1 1 2
1 2 1 2 1 2 4
2 4 2 4 1 4 8
3 8 3 8 1 8 16

La conclusione di questo calcolo è  y_4 = 16 . La soluzione esatta della equazione differenziale è  y(t) = e^t , quindi  y(4) = e^4 \approx 54.598 . Cioè, l'approssimazione del metodo di Eulero non è buona in questo caso. In ogni modo, come è mostrato dalla figura, il suo comportamento è qualitativamente corretto.

Usando altre dimensioni per l'incremento[modifica | modifica sorgente]

The same illustration for h = 0.25.

Come ricordato nell'introduzione, il metodo di Eulero è tanto più accurato minore è il passo h. La tabella sottostante mostra i diversi risultati sfruttando passi di diverse dimensioni. La prima riga corrisponde all'esempio della sezione precedente mentre la seconda riga è illustrata in figura.

passo risultato con Eulero errore
1 16 38.598
0.25 35.53 19.07
0.1 45.26 9.34
0.05 49.56 5.04
0.025 51.98 2.62
0.0125 53.26 1.34

L'errore riportato nell'ultima colonna della tabella indica la differenza tra la soluzione esatta per  t = 4 e l'approssimazione con il metodo di Eulero. Si può constatare che al dimezzarsi del passo corrisponde approssimativamente un dimezzamento dell'errore, questo suggerisce che l'errore è all'incirca proporzionale alla dimensione del passo, almeno per piccoli valori di  h . Questo è vero in generale ed anche per altre equazioni.

Altri metodi, come il metodo del punto medio, sono generalmente più precisi: sfruttando il punto medio l'errore è all'incirca proporzionale al quadrato della del passo. Per questo motivo il metodo di Eulero viene definito del primo ordine, mentre il metodo del punto medio del secondo ordine.

Dalla tabella è possibile affermare che il passo necessario per approssimare un risultato fino alla terza cifra decimale è circa 0.00001: il che significa che saranno necessari 400000 passi. Questo alto numero di passi implica un alto costo computazionale. Per questo motivo generalmente si preferiscono metodi numerici più precisi (cioè di ordine maggiore) come ad esempio il metodo di Runge-Kutta oppure i metodi lineari multistep, soprattutto se è necessaria grande accuratezza.[6]

Derivazione[modifica | modifica sorgente]

The Euler method can be derived in a number of ways. Firstly, there is the geometrical description mentioned above.

Un'altra possibilità è considerare la serie di Taylor della funzione y intorno t_0:

 y(t_0 + h) = y(t_0) + h y'(t_0) + \frac{1}{2}h^2 y''(t_0) + O(h^3).

The differential equation states that y'=f(t,y). If this is substituted in the Taylor expansion and the quadratic and higher-order terms are ignored, the Euler method arises.[7] The Taylor expansion is used below to analyze the error committed by the Euler method, and it can be extended to produce Runge–Kutta methods.

A closely related derivation is to substitute the forward finite difference formula for the derivative,

 y'(t_0) \approx \frac{y(t_0+h) - y(t_0)}{h}

in the differential equation y' = f(t,y). Again, this yields the Euler method.[8] A similar computation leads to the midpoint rule and the backward Euler method.

Finally, one can integrate the differential equation from t_0 to t_0 + h and apply the fundamental theorem of calculus to get:

 y(t_0+h) - y(t_0) = \int_{t_0}^{t_0+h} f(t,y(t)) \,\mathrm{d}t.

Now approximate the integral by the left-hand rectangle method (with only one rectangle):

 \int_{t_0}^{t_0+h} f(t,y(t)) \,\mathrm{d}t \approx h f(t_0, y(t_0)).

Combining both equations, one finds again the Euler method.[9] This line of thought can be continued to arrive at various linear multistep methods.

Local truncation error[modifica | modifica sorgente]

The local truncation error of the Euler method is error made in a single step. It is the difference between the numerical solution after one step, y_1, and the exact solution at time t_1 = t_0+h. The numerical solution is given by

 y_1 = y_0 + h f(t_0, y_0). \quad

For the exact solution, we use the Taylor expansion mentioned in the section Derivation above:

 y(t_0 + h) = y(t_0) + h y'(t_0) + \frac{1}{2}h^2 y''(t_0) + O(h^3).

The local truncation error (LTE) introduced by the Euler method is given by the difference between these equations:

 \mathrm{LTE} = y(t_0 + h) - y_1 = \frac{1}{2} h^2 y''(t_0) + O(h^3).

This result is valid if y has a bounded third derivative.[10]

This shows that for small h, the local truncation error is approximately proportional to h^2. This makes the Euler method less accurate (for small h) than other higher-order techniques such as Runge-Kutta methods and linear multistep methods, for which the local truncation error is proportial to a higher power of the step size.

A slightly different formulation for the local truncation error can be obtained by using the Lagrange form for the remainder term in Taylor's theorem. If y has a continuous second derivative, then there exists a \xi \in [t_0,t_0+h] such that

 \mathrm{LTE} = y(t_0 + h) - y_1 = \frac{1}{2} h^2 y''(\xi). [11]

In the above expressions for the error, the second derivative of the unknown exact solution y can be replaced by an expression involving the right-hand side of the differential equation. Indeed, it follows from the equation y'=f(t,y) that

y''(t_0) = {\partial f\over\partial t}(t_0, y(t_0)) + {\partial f\over\partial y}(t_0, y(t_0)) \, f(t_0, y(t_0)).[12]

Errore di troncamento globale[modifica | modifica sorgente]

The global truncation error is the error at a fixed time t, after however many steps the methods needs to take to reach that time from the initial time. The global truncation error is the cumulative effect of the local truncation errors committed in each step.[13] The number of steps is easily determined to be (t-t_0)/h, which is proportional to 1/h, and the error committed in each step is proportional to h^2 (see the previous section). Thus, it is to be expected that the global truncation error will be proportional to h.[14]

This intuitive reasoning can be made precise. If the solution y has a bounded second derivative and f is Lipschitz continuous in its second argument, then the global truncation error (GTE) is bounded by

 |\text{GTE}| \le \frac{hM}{2L}(e^{L(t-t_0)}-1) \qquad \qquad

where M is an upper bound on the second derivative of y on the given interval and L is the Lipschitz constant of f.[15]

The precise form of this bound of little practical importance, as in most cases the bound vastly overestimates the actual error committed by the Euler method.[16] What is important is that it shows that the global truncation error is (approximately) proportional to h. For this reason, the Euler method is said to be first order.[17]

Stabilità Numerica[modifica | modifica sorgente]

Solution of y' = -2.3y computed with the Euler method with step size h=1 (blue squares) and h=0.7 (red circles). The black curve shows the exact solution.

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

 y' = -2.3y, \qquad y(0) = 1.

The exact solution is y(t) = e^{-2.3t}, which decays to zero as t \to \infty. However, if the Euler method is applied to this equation with step size 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 h=0.7, then the numerical solution does decay to zero.

The pink disk shows the stability region for the Euler method.

If the Euler method is applied to the linear equation y' = k y, then the numerical solution is unstable if the product hk is outside the region

 \{ z \in \mathbf{C} \mid |z+1| \le 1 \},

illustrated on the right. This region is called the (linear) instability region.[18] In the example, k equals −2.3, so if h=1 then 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 h—means that the Euler method is not often used, except as a simple example of numerical integration.

Errori di arrotondamento[modifica | modifica sorgente]

The discussion up to now has ignored the consequences of rounding error. In step n of the Euler method, the rounding error is roughly of the magnitude εyn where ε is the machine epsilon. Assuming that the rounding errors are all of approximately the same size, the combined rounding error in N steps is roughly Nεy0 if all errors points in the same direction. Since the number of steps is inversely proportional to the step size h, the total rounding error is proportional to ε / h. In reality, however, it is extremely unlikely that all rounding errors point in the same direction. If instead it is assumed that the rounding errors are independent rounding variables, then the total rounding error is proportional to  \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]

Modifications and extensions[modifica | modifica sorgente]

A simple modification of the Euler method which eliminates the stability problems noted in the previous section is the backward Euler method:

 y_{n+1} = y_n + h f(t_{n+1}, y_{n+1}).

This differs from the (standard, or forward) Euler method in that the function 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  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.

Other modifications of the Euler method that help with stability yield the exponential Euler method or the semi-implicit Euler method.

More complicated methods can achieve a higher order (and more accuracy). One possibility is to use more function evaluations. This is illustrated by the midpoint method which is already mentioned in this article:

 y_{n+1} = y_n + h f \Big( t_n + \tfrac12 h, y_n + \tfrac12 h f(t_n, y_n) \Big).

This leads to the family of Runge–Kutta methods.

The other possibility is to use more past values, as illustrated by the two-step Adams–Bashforth method:

 y_{n+1} = y_n + \tfrac32 h f(t_{n}, y_{n}) - \tfrac12 h f(t_{n-1}, y_{n-1}).

This leads to the family of linear multistep methods.

Voci correlate[modifica | modifica sorgente]

Note[modifica | modifica sorgente]

  1. ^ Butcher 2003, p. 45; Hairer, Nørsett & Wanner 1993, p. 35
  2. ^ Atkinson 1989, p. 342; Butcher 2003, p. 60
  3. ^ Butcher 2003, p. 45; Hairer, Nørsett & Wanner 1993, p. 36
  4. ^ Butcher 2003, p. 3; Hairer, Nørsett & Wanner 1993, p. 2
  5. ^ Vedi anche Atkinson 1989, p. 344
  6. ^ Hairer, Nørsett & Wanner 1993, p. 40
  7. ^ Atkinson 1989, p. 342; Hairer, Nørsett & Wanner 1993, p. 36
  8. ^ Atkinson 1989, p. 342
  9. ^ Atkinson 1989, p. 343
  10. ^ Butcher 2003, p. 60
  11. ^ Atkinson 1989, p. 342
  12. ^ Stoer & Bulirsch 2002, p. 474
  13. ^ Atkinson 1989, p. 344
  14. ^ Butcher 2003, p. 49
  15. ^ Atkinson 1989, p. 346; Lakoba 2012, equation (1.16)
  16. ^ Iserles 1996, p. 7
  17. ^ Butcher 2003, p. 63
  18. ^ Butcher 2003, p. 70; Iserles 1996, p. 57
  19. ^ Butcher 2003, pp. 74–75
  20. ^ Butcher 2003, pp. 75–78

Bibliografia[modifica | modifica sorgente]

Collegamenti esterni[modifica | modifica sorgente]

Template:Wikibooks