Queste note sono basate sugli appunti fatti con Gianmarco Manzini negli anni 1995-2005

Zeri di funzione

Appunti di calcolo numerico

Autore/Autrice
Affiliazione

Enrico Bertolazzi

University of Trento, Department of Industrial Engineering

Introduzione

Un problema comune in molte applicazioni consiste nel trovare una soluzione x tale che:

f(x) = 0,

dove f: \Bbb{R} \to \Bbb{R} è una funzione continua. Un numero reale \alpha che soddisfa f(\alpha) = 0 è chiamato radice della funzione. Se f(x) è una funzione lineare, ossia della forma:

f(x) = a x + b,

dove a \neq 0 e b sono costanti assegnate, la radice si calcola facilmente come \alpha = -b/a.

Se invece f(x) ha una forma non lineare, ad esempio:

f(x) = a x^{n} - b,

e se a e b hanno lo stesso segno o se n è dispari, si può determinare una radice reale con \alpha = \sqrt[n]{b/a}. Tuttavia, nei problemi pratici, la funzione f(x) è spesso non lineare e non esiste una formula esplicita per trovare facilmente una radice. In questi casi, è necessario utilizzare metodi numerici per approssimare le radici.

I metodi numerici per approssimare le soluzioni di problemi del tipo f(x) = 0 si suddividono principalmente in due categorie:

  • Metodi iterativi a un punto
  • Metodi iterativi a due punti

Metodi iterativi a due punti

Questi metodi, come il metodo della bisezione (o dicotomia) e il metodo delle corde (o falsa posizione), determinano a ogni passo due estremi di un intervallo che contiene una radice. La lunghezza di questo intervallo si riduce progressivamente, consentendo di approssimare la radice con la precisione desiderata.

Metodi iterativi a un punto

I metodi a un punto generano una successione di valori che converge verso la radice. Esempi di questi metodi includono il metodo delle secanti, il metodo di Newton1-Raphson2 e il metodo delle iterazioni a punto fisso.

In alcuni casi particolari, come nel caso di polinomi, esistono tecniche speciali che non rientrano nelle due categorie principali sopra descritte.

Metodo di bisezione (o dicotomico)

Il metodo di bisezione è una semplice applicazione del seguente teorema, che enunciamo senza dimostrazione.

Teorema 1 (Teorema: Esistenza di uno zero per cambio di segno) Sia f:[a,b] \to \Bbb{R} una funzione continua su [a,b] e tale che f(a)f(b) < 0. Allora esiste almeno un punto \alpha \in (a,b) tale che f(\alpha) = 0.

Se f(x) è continua e a_0 e b_0 sono due punti tali che f(a_0)f(b_0) < 0, significa che f(a_0) e f(b_0) hanno segno opposto. Il teorema Teorema 1 garantisce l’esistenza di almeno una radice \alpha nell’intervallo (a_0, b_0).

Ora consideriamo il punto medio dell’intervallo:

c_0 = \frac{a_0 + b_0}{2}.

Si possono verificare solo tre possibili situazioni:

  • Caso 1: f(c_0) = 0; in tal caso, c_0 è una radice esatta.

  • Caso 2: f(c_0)f(a_0) < 0; in questo caso, esiste almeno una radice nell’intervallo (a_0, c_0).

  • Caso 3: f(c_0)f(b_0) < 0; in questo caso, esiste almeno una radice nell’intervallo (c_0, b_0).

Escludendo il primo caso in cui la radice è già stata trovata, è sempre possibile ridurre l’intervallo di ricerca dimezzandolo. Ripetendo questo procedimento, il metodo dicotomico consente, data una tolleranza \epsilon, di trovare un intervallo di lunghezza inferiore a \epsilon che contiene una radice.

L’algoritmo può essere descritto come segue:

Algoritmo: Metodo di Bisezione

Input: a, b, tolleranza \epsilon
Output: Approssimazione della radice

  • Inizializza a_0 \gets a, b_0 \gets b, i \gets 0
  • While (b_i - a_i) \geq \epsilon do:
    • Calcola il punto medio: c \gets \dfrac{a_i + b_i}{2}
    • If f(c) = 0, restituisci c (la radice esatta è trovata)
    • If f(c)f(a_i) < 0 Then
      • Allora esiste almeno una radice in (a_i, c):
      • Aggiorna: a_{i+1} \gets a_i, b_{i+1} \gets c
    • Else
      • esiste almeno una radice in (c, b_i):
      • Aggiorna: a_{i+1} \gets c, b_{i+1} \gets b_i
    • End If
    • Incrementa i \gets i + 1
  • End while

Alla fine, c contiene l’approssimazione della radice con precisione \epsilon.
Restituisci c.

È possibile calcolare il numero minimo di passi necessari per approssimare una radice con l’accuratezza desiderata. Consideriamo che, ad ogni iterazione, la lunghezza dell’intervallo si dimezza. Dunque, per l’intervallo all’iterazione i abbiamo:

\begin{aligned} b_i - a_i &= \frac{1}{2}(b_{i-1} - a_{i-1}), \\ &\vdots \\ &= 2^{-i}(b_0 - a_0), \\ &= 2^{-i}(b - a), \end{aligned}

Imponendo che b_i - a_i < \epsilon, si ottiene:

2^{-i}(b - a) < \epsilon,\qquad \Rightarrow \qquad \frac{b - a}{\epsilon} < 2^i.

Applicando il logaritmo base 2 ad entrambi i membri:

i > \log_2\left(\frac{b - a}{\epsilon}\right) = \log_2 (b - a) - \log_2 (\epsilon).

La parte intera del risultato fornisce il numero minimo di iterazioni necessarie per garantire che l’errore sia minore o uguale a \epsilon.

Metodo delle false posizioni (Regula Falsi)

Sia data una funzione f:[a,b] \to \Bbb{R} continua nell’intervallo [a,b] e tale che f(a)f(b) < 0. In queste condizioni, il teorema degli zeri garantisce che esiste almeno uno zero di f(x) nell’intervallo [a,b].

Per stimare la posizione della radice, costruiamo la retta che interpola i punti:

\begin{pmatrix} a \\ f(a) \end{pmatrix}, \qquad \begin{pmatrix} b \\ f(b) \end{pmatrix},

ossia l’equazione della retta diventa:

y = f(a) + (x-a) \dfrac{f(b)-f(a)}{b-a}, \tag{1}

Figura 1: Costruzione della regula falsi

Lo zero della retta interpolante, ottenuto risolvendo l’equazione Equazione 1, fornisce la seguente stima della radice \hat{x} della funzione f(x):

0 = f(a) + (\hat{x}-a) \dfrac{f(b)-f(a)}{b-a},\qquad \Rightarrow \qquad \hat{x} = \dfrac{af(b)-bf(a)}{f(b)-f(a)}.

A questo punto, possiamo applicare un procedimento simile a quello del metodo di bisezione, dividendo l’intervallo in [a,\hat{x}] o [\hat{x},b], a seconda del segno di f(\hat{x}). Tuttavia, la regula falsi non garantisce necessariamente che gli intervalli costruiti diventino sempre piu piccoli come mostrato in Figura 2.

Figura 2: Andamento della regula falsi

Infatti, il metodo genera una successione di due punti x_{i} e a_i, dove x_i è un’approssimazione della radice e a_i è tale che f(x_{i})f(a_i) < 0 dove solo x_i converge alla radice.

Il metodo delle false posizioni può essere riassunto nel seguente algoritmo.

Algoritmo: Regula Falsi

Input: a, b

Output: Approssimazione della radice

  • Inizializza: a_0 \gets a, x_0 \gets b, i \gets 0
  • While \lvert f(x_{i}) \rvert \geq \epsilon do:
    • Calcola la nuova approssimazione della radice: \hat{x} \gets \dfrac{a_i f(x_{i}) - x_{i} f(a_i)}{f(x_{i}) - f(a_i)}
    • If f(x_{i}) f(\hat{x}) < 0 then:
      • Aggiorna x_{i+1} \gets \hat{x}, a_{i+1} \gets a_i
    • Else:
      • Aggiorna x_{i+1} \gets x_{i}, a_{i+1} \gets \hat{x}
    • End If
    • Incrementa: i \gets i+1
  • End While
  • Restituisci x_{i} come approssimazione della radice.

Convergenza della regula falsi

Per discutere la convergenza della regula falsi, assumiamo che la funzione f(x) abbia una derivata seconda regolare su tutto l’intervallo di definizione e consideriamo le seguenti ipotesi:

  1. x_0 < a_0;
  2. f(x_0) < 0 e f(a_0) > 0;
  3. f''(x) \geq 0 per ogni x \in [x_0, a_0].

Con queste condizioni, è possibile dimostrare il seguente teorema.

Teorema 2 (Teorema: Convergenza della regula falsi) Sotto le ipotesi a–c, la regula falsi converge, e inoltre:

x_k \leq x_{k+1} < a_k = a_0, \qquad \text{per } k = 1, 2, \ldots.

Dimostrazione. Osserviamo innanzitutto che l’ipotesi f''(x) \geq 0 in [x_0, a_0] implica che la funzione f(x) è convessa su [x_0, a_0], ovvero:

f(t x + (1-t) y) \leq t f(x) + (1-t) f(y), \qquad \forall x, y \in [x_0, a_0].

Inoltre, poiché f(x_0) < 0 e f(a_0) > 0, esiste almeno una radice nell’intervallo [x_0, a_0]. Definiamo:

\eta = \sup \{ x \in [x_0, a_0] \,|\, f(x) < 0 \}.

Dalla continuità di f(x), segue che f(\eta) = 0. Ora, calcoliamo x_1:

x_1 = \frac{a_0 f(x_0) - x_0 f(a_0)}{f(x_0) - f(a_0)} = x_0 + (a_0 - x_0) \frac{f(x_0)}{f(x_0) - f(a_0)}.

Poiché f(x_0) < 0, f(a_0) > 0 e x_0 < a_0, segue che x_1 > x_0. Dalla convessità di f(x), possiamo affermare che:

\begin{aligned} f(x_1) &= f\left(x_0 \frac{f(a_0)}{f(a_0) - f(x_0)} + a_0 \left(1 - \frac{f(a_0)}{f(a_0) - f(x_0)} \right)\right), \\ &\leq \frac{f(a_0)}{f(a_0) - f(x_0)} f(x_0) + \left(1 - \frac{f(a_0)}{f(a_0) - f(x_0)}\right) f(a_0) = 0, \end{aligned}

e quindi x_1 \leq \eta. Procedendo in modo analogo, otteniamo:

x_0 \leq x_1 \leq x_2 \leq \ldots \leq x_k \leq \eta < a_0 = a_1 = \cdots = a_k.

Di conseguenza, x_k, per k = 0, 1, \ldots, è una successione monotona crescente e, dunque, convergente. Indichiamo il limite di questa successione con \mu. Si ha quindi:

\mu = \lim_{k \to \infty} x_k \leq \eta.

Inoltre, deve valere:

\mu = \frac{a_0 f(\mu) - \mu f(a_0)}{f(\mu) - f(a_0)},

da cui segue che:

(\mu - a_0) f(\mu) = 0.

Poiché \mu - a_0 \leq \eta - a_0 < 0, ne segue che f(\mu) = 0, dimostrando così la convergenza del metodo.

Metodo di Newton-Raphson

Consideriamo una funzione continua e differenziabile f : [a, b] \to \Bbb{R}. L’obiettivo è trovare una radice della funzione, ovvero un punto \alpha tale che f(\alpha) = 0. Se x_0 non è già uno zero di f(x), possiamo stimare la posizione della radice approssimando la funzione f(x) attraverso la retta tangente al grafico di f(x) nel punto (x_0, f(x_0)). L’equazione di questa retta tangente è:

y = f(x_0) + (x - x_0) f'(x_0). \tag{2}

Questa equazione rappresenta una linea retta con pendenza f'(x_0), che passa per il punto (x_0, f(x_0)). Lo zero di questa retta, cioè il punto in cui y = 0, fornisce una prima approssimazione della radice di f(x).

Figura 3

Costruzione grafica del metodo di Newton-Raphson

Per trovare il valore di x_1, lo zero della retta tangente, imponiamo y = 0 nell’equazione della retta:

0 = f(x_0) + (x_1 - x_0) f'(x_0),

da cui si ricava:

x_1 = x_0 - \dfrac{f(x_0)}{f'(x_0)}.

Questa espressione fornisce una prima approssimazione x_1 della radice di f(x). Ripetendo il procedimento utilizzando x_1, otteniamo un nuovo punto x_2, che dovrebbe essere un’approssimazione ancora più accurata della radice.

Questo suggerisce un processo iterativo, dove partendo da un’ipotesi iniziale x_0, possiamo generare una successione di punti x_1, x_2, x_3, \ldots, ognuno dei quali rappresenta una stima migliore della radice di f(x). La formula iterativa è la seguente:

\boxed{ x_{k+1} = x_k - \dfrac{f(x_k)}{f'(x_k)}, \qquad k = 0, 1, 2, \ldots } \tag{3}

Questa procedura è nota come metodo di Newton-Raphson. Ogni iterazione si basa sul calcolo della tangente alla curva di f(x) nel punto corrente x_k, e il punto in cui questa tangente incontra l’asse x diventa il nuovo punto x_{k+1}. Se la funzione è sufficientemente regolare e l’approssimazione iniziale x_0 è abbastanza vicina alla radice, il metodo converge rapidamente verso lo zero di f(x).

Vale il seguente teorema:

Teorema 3 (Convergenza del metodo di Newton per radici semplici) Sia f \in C^2([a,b]) e \alpha una radice semplice di f(x), cioè tale che f'(\alpha) \neq 0. Allora, se la stima iniziale x_0 si trova sufficientemente vicino a \alpha, ossia x_0 \in (\alpha - \epsilon, \alpha + \epsilon) per un \epsilon abbastanza piccolo, la successione generata dal metodo di Newton-Raphson:

x_{k+1} = x_k - \frac{f(x_k)}{f'(x_k)}

soddisferà le seguenti condizioni:

  • x_k \in (\alpha - \epsilon, \alpha + \epsilon) per k = 1, 2, 3, \ldots.
  • La convergenza sarà quadratica, ossia: \left\lvert x_{k+1} - \alpha\right\rvert \leq C \left\lvert x_k - \alpha\right\rvert^2 per k = 1, 2, 3, \ldots, dove C è una costante positiva.

Dimostrazione. Consideriamo lo sviluppo di f(x) attorno alla radice \alpha tramite la serie di Taylor:

f(\alpha) = f(x_k) + f'(x_k)(\alpha - x_k) + \frac{f''(\eta_k)}{2}(x_k - \alpha)^2,

dove \eta_k è un punto compreso tra x_k e \alpha. Utilizzando l’iterazione del metodo di Newton-Raphson: Assumendo f'(x_k)\neq 0 e sapendo che f(\alpha)=0 possiamo scivere

0 = \dfrac{f(x_k)}{f'(x_k)} + \alpha - x_k + \frac{f''(\eta_k)}{2f'(x_k)}(x_k - \alpha)^2, \tag{4}

utilizzando l’espressione per il passo di Newton:

x_{k+1} = x_k - \frac{f(x_k)}{f'(x_k)},

nella Equazione 4 possiamo eliminare \frac{f(x_k)}{f'(x_k)} ottenendo:

\alpha - x_{k+1} = \frac{f''(\eta_k)}{2f'(x_k)}(x_k - \alpha)^2.

Prendendo il valore assoluto, otteniamo:

\left\lvert x_{k+1} - \alpha\right\rvert \leq \frac{M}{2 \left\lvert f'(x_k)\right\rvert} \left\lvert x_k - \alpha\right\rvert^2,

dove M = \sup_{\eta \in [a,b]} \left\lvert f''(\eta)\right\rvert. Poiché \alpha è una radice semplice, esiste un \epsilon tale che 2 \left\lvert f'(x_k)\right\rvert > \left\lvert f'(\alpha)\right\rvert per ogni x_k \in (\alpha - \epsilon, \alpha + \epsilon). Questo implica che:

\left\lvert x_{k+1} - \alpha\right\rvert \leq \frac{M}{\left\lvert f'(\alpha)\right\rvert} \left\lvert x_k - \alpha\right\rvert^2.

Riducendo eventualmente \epsilon, possiamo garantire che:

\frac{M}{\left\lvert f'(\alpha)\right\rvert} < 1,

quindi la successione \{x_k\} rimane in (\alpha - \epsilon, \alpha + \epsilon) e converge a \alpha con velocità quadratica. Definendo C = \frac{M}{\left\lvert f'(\alpha)\right\rvert}, otteniamo:

\left\lvert x_{k+1} - \alpha\right\rvert \leq C \left\lvert x_k - \alpha\right\rvert^2,

che conclude la dimostrazione.

Metodo delle Secanti

Anziché controllare ripetutamente il segno del prodotto f(x_{i}) f(a_i), come nel metodo della Regula Falsi, possiamo utilizzare i valori delle ultime due iterate per calcolare una nuova approssimazione della radice. Questo porta al metodo delle secanti, che si esprime nella seguente forma:

\boxed{ x_{i+1} = \dfrac{x_{i-1} f(x_{i}) - x_{i} f(x_{i-1})}{f(x_{i}) - f(x_{i-1})}, \qquad i = 1, 2, \ldots }

Nel metodo delle secanti, quando si raggiunge la convergenza, si ha che f(x_{i}) \approx f(x_{i-1}), il che può causare problemi di cancellazione numerica. Al contrario, nel metodo della Regula Falsi questo non si verifica, poiché è garantito che f(x_{i}) f(a_i) < 0.

Inoltre, se x_{i} non è sufficientemente vicino alla radice, non è garantito che il metodo converga. Pertanto, è utile esaminare la convergenza locale di questo schema. Sia \alpha una radice di f(x). Sottraendo \alpha da entrambi i membri della formula sopra, otteniamo:

\begin{aligned} x_{i+1} - \alpha &= \dfrac{x_{i-1} f(x_{i}) - x_{i} f(x_{i-1})}{f(x_{i}) - f(x_{i-1})} - \alpha \\ &= \dfrac{x_{i-1} f(x_{i}) - x_{i} f(x_{i-1}) - \alpha (f(x_{i}) - f(x_{i-1}))}{f(x_{i}) - f(x_{i-1})} \\ &= \dfrac{(x_{i-1} - \alpha) f(x_{i}) - (x_{i} - \alpha) f(x_{i-1})}{f(x_{i}) - f(x_{i-1})} \\ &= (x_{i} - \alpha)(x_{i-1} - \alpha) \dfrac{\dfrac{f(x_{i})}{x_{i} - \alpha} - \dfrac{f(x_{i-1})}{x_{i-1} - \alpha}}{f(x_{i}) - f(x_{i-1})}. \end{aligned}

Inoltre, poiché f(\alpha) = 0, possiamo scrivere:

\begin{aligned} \dfrac{\dfrac{f(x_{i})}{x_{i} - \alpha} - \dfrac{f(x_{i-1})}{x_{i-1} - \alpha}}{f(x_{i}) - f(x_{i-1})} &= \dfrac{\dfrac{f(x_{i}) - f(\alpha)}{x_{i} - \alpha} - \dfrac{f(x_{i-1}) - f(\alpha)}{x_{i-1} - \alpha}}{f(x_{i}) - f(x_{i-1})} \\ &= \dfrac{\dfrac{f(x_{i}) - f(\alpha)}{x_{i} - \alpha} - \dfrac{f(x_{i-1}) - f(\alpha)}{x_{i-1} - \alpha}}{x_{i} - x_{i-1}} \left( \dfrac{f(x_{i}) - f(x_{i-1})}{x_{i} - x_{i-1}} \right)^{-1}, \end{aligned} \tag{5}

dove possiamo applicare il teorema di Lagrange3, che afferma che:

\dfrac{f(x_{i}) - f(x_{i-1})}{x_{i} - x_{i-1}} = f'(\eta_i), \qquad \eta \in I[x_{i-1}, x_{i}], \tag{6}

dove I[a, b, c, \ldots] denota il più piccolo intervallo chiuso contenente a, b, c, \ldots.

Lemma 1 Se f \in C^{2}, allora vale la seguente uguaglianza:

\dfrac{\dfrac{f(\alpha + h) - f(\alpha)}{h} - \dfrac{f(\alpha - k) - f(\alpha)}{k}}{h + k} = \dfrac{1}{2} f''(\zeta),

dove \zeta \in [\alpha - k, \alpha + h].

Dimostrazione. Consideriamo la funzione

G(t) := \dfrac{\dfrac{f(\alpha + th) - f(\alpha)}{h} - \dfrac{f(\alpha - tk) - f(\alpha)}{k}}{h + k}.

La funzione

H(t) := G(t) - G(1) t^{2}

si annulla per t = 0 e t = 1. Pertanto, secondo il teorema di Rolle4, esiste un punto \eta \in (0,1) tale che H'(\eta) = 0.

Calcoliamo H'(t):

H'(t) = G'(t) - 2G(1)t,

dove

G'(t) = \dfrac{f'(\alpha + th) - f'(\alpha - tk)}{h + k}.

Valutando G'(\eta), otteniamo:

\begin{aligned} 0 &= G'(\eta) - 2G(1)\eta, \\ G(1) &= \dfrac{1}{2\eta} G'(\eta), \\ &= \dfrac{f'(\alpha + \eta h) - f'(\alpha - \eta k)}{2\eta(h + k)} \\ &= \dfrac{1}{2} f''(\zeta). \end{aligned}

Concludiamo quindi che la relazione è valida.

Ponendo h=x_{i}-\alpha e k=\alpha-x_{i-1} nel lemma Lemma 1 otteniamo

\dfrac{\dfrac{f(x_{i})-f(\alpha)}{x_{i}-\alpha}-\dfrac{f(x_{i-1})-f(\alpha)}{x_{i-1}-\alpha}} {x_{i}-x_{i-1}} = \dfrac{1}{2} f''(\zeta), \tag{7}

sostituendo Equazione 5 ed Equazione 6 in Equazione 7 otteniamo

x_{i+1}-\alpha = (x_{i}-\alpha)(x_{i-1}-\alpha) \dfrac{f''(\zeta_{i})}{2f'(\eta_{i})},\qquad \zeta_{i},\eta_{i}\in I[x_{i-1},x_{i},\alpha] \tag{8}

questa equazione è utile per dimostrare il seguente teorema:

Teorema 4 (Convergenza del Metodo delle Secanti) Se \alpha è uno zero semplice di f \in C^{2}([a,b]) (cioè f'(\alpha) \neq 0), allora esiste un intervallo J = (x - \epsilon, x + \epsilon) tale che, per ogni scelta di x_0, x_1 \in J, il metodo delle secanti converge. Inoltre, si ha:

\left\lvert x_{i} - \alpha\right\rvert \leq C q^{p^{i}}, \quad 0 < q < 1, \quad p = \dfrac{1 + \sqrt{5}}{2} \approx 1.618034.

Dimostrazione. Poiché \alpha è uno zero semplice e f'(x) è continua, esiste un \epsilon > 0 tale che:

\left\lvert f'(x)\right\rvert > \dfrac{1}{2} \left\lvert f'(\alpha)\right\rvert, \quad \forall x \in (\alpha - \epsilon, \alpha + \epsilon).

Inoltre, poiché f''(x) è continua, possiamo trovare una costante K per cui:

\left\lvert f''(x)\right\rvert \leq K, \quad \forall x \in (\alpha - \epsilon, \alpha + \epsilon).

Definiamo J = (x - \epsilon, x + \epsilon) e poniamo M = \dfrac{K}{\left\lvert f'(\alpha)\right\rvert}. Allora, per ogni \eta \in J e \zeta \in J, possiamo affermare:

\left\lvert\dfrac{1}{2} \dfrac{f''(\zeta)}{f'(\eta)}\right\rvert \leq \left\lvert\dfrac{f''(\zeta)}{f'(\alpha)}\right\rvert \leq M.

Supponendo che \epsilon M < 1 (altrimenti basta ridurre \epsilon), otteniamo che x_{i} \in J per i = 1, 2, \ldots. Infatti, dalla relazione in Equazione 8, abbiamo:

\left\lvert x_{i+1} - \alpha\right\rvert \leq \left\lvert x_{i} - \alpha\right\rvert \left\lvert x_{i-1} - \alpha\right\rvert M \leq \epsilon^{2} M < \epsilon, \quad \implies \quad x_{i+1} \in J.

Definiamo e_{i} = M \left\lvert x_{i} - \alpha\right\rvert; dalla relazione precedente otteniamo:

\dfrac{e_{i+1}}{M} \leq \dfrac{e_{i}}{M} \dfrac{e_{i-1}}{M} M, \quad \implies \quad e_{i+1} \leq e_{i} e_{i-1}.

Poniamo E_{0} \geq e_{0} e E_{1} \geq e_{1}, con E_{i+1} = E_{i} E_{i-1}. Allora, per i = 0, 1, 2, \ldots, abbiamo che e_{i} \leq E_{i}. Infatti, se ciò vale fino a k, allora:

e_{k+1} \leq e_{k} e_{k-1} \leq E_{k} E_{k-1} = E_{k+1}.

Ponendo E_{i} = \exp(c z^{i}), otteniamo:

\exp(c z^{i+1}) = \exp(c z^{i}) \exp(c z^{i-1}) = \exp(c (z^{i} + z^{i-1})).

Prendendo il logaritmo e dividendo per c z^{i-1}, arriviamo all’equazione:

z^{2} = z + 1, \quad \implies \quad z_{1,2} = \dfrac{1 \pm \sqrt{5}}{2}.

Scegliendo p = \dfrac{1 + \sqrt{5}}{2}, possiamo esprimere la soluzione particolare della relazione di ricorrenza E_{k+1} = E_{k} E_{k-1} come:

E_{0} = \exp(c), \quad E_{1} = \exp(cp), \quad E_{k} = \exp(cp^{k}).

Assumendo che la costante c soddisfi le seguenti condizioni:

\begin{aligned} e_{0} &\leq E_{0} = \exp(c), & \quad c \geq \log(e_{0}), \\ e_{1} &\leq E_{1} = \exp(cp), & \quad c \geq \dfrac{\log(e_{1})}{p}, \end{aligned}

definiamo q = \exp(c). Da ciò segue che:

e_{i} \leq \exp(cp^{i}), \quad \implies \quad M \left\lvert x_{i} - \alpha\right\rvert \leq q^{p^{i}}.

Infine, osserviamo che:

e_{i} = \left\lvert x_{i} - \alpha\right\rvert M < \epsilon M < 1,

perciò possiamo scegliere c < 0 e q = \exp(c) < 1.

Metodo delle Parabole Tangenti

Il metodo delle parabole tangenti, noto anche come metodo di Müller, è un algoritmo numerico utilizzato per trovare approssimazioni agli zeri di una funzione non lineare. A differenza del metodo delle secanti, che utilizza una retta, il metodo di Müller usa una parabola come interpolante tra tre punti distinti per migliorare l’approssimazione della radice.

Descrizione del Metodo

Il metodo parte da tre punti distinti x_{k-2}, x_{k-1}, x_k, e costruisce una parabola che interpola la funzione f(x) in questi punti. La parabola ha la seguente forma generale:

P(x) = a(x - x_k)^2 + b(x - x_k) + c

I coefficienti a, b e c sono determinati imponendo che la parabola passi attraverso i punti (x_{k-2}, f(x_{k-2})), (x_{k-1}, f(x_{k-1})), e (x_k, f(x_k)).

Costruzione del metodo di Müller

Consideriamo l’espansione di Taylor attorno alla radice \alpha a partire da un punto x_0:

0 = f(\alpha) = f(x_0) + h f'(x_0) + \frac{h^2}{2} f''(x_0) + \mathcal{O}(h^3).

Trascurando il termine \mathcal{O}(h^3), otteniamo la seguente condizione su h:

0 = f(x_0) + h f'(x_0) + \frac{h^2}{2} f''(x_0), \tag{9}

risolvendo per h, otteniamo:

h_{1,2} = \frac{-f'(x_0) \pm \sqrt{f'(x_0)^2 - 2f(x_0) f''(x_0)}}{f''(x_0)}. \tag{10}

Osserviamo che:

\begin{aligned} \sqrt{f'(x_0)^2 - 2f(x_0) f''(x_0)} &= \sqrt{f'(x_0)^2 \left(1 - 2 \frac{f(x_0) f''(x_0)}{f'(x_0)^2}\right)} \\ &= \left\lvert f'(x_0)\right\rvert \sqrt{1 - 2 \frac{f(x_0) f''(x_0)}{f'(x_0)^2}}. \end{aligned} \tag{11}

Poiché f'(x_0) = \mathrm{sign}(f'(x_0)) \left\lvert f'(x_0)\right\rvert, possiamo riscrivere l’equazione precedente come:

\sqrt{f'(x_0)^2 - 2f(x_0) f''(x_0)} = \mathrm{sign}(f'(x_0)) f'(x_0) \sqrt{1 - 2 \frac{f(x_0) f''(x_0)}{f'(x_0)^2}}. \tag{12}

Combinando l’equazione Equazione 10 con Equazione 12, otteniamo:

h_{1,2} = -\frac{f'(x_0)}{f''(x_0)} \left(1 \pm \mathrm{sign}(f'(x_0)) \sqrt{1 - 2 \frac{f(x_0) f''(x_0)}{f'(x_0)^2}}\right).

Scelgo la soluzione che produce il valore assoluto di h più piccolo, il che corrisponde a \pm \mathrm{sign}(f'(x_0)) = -1. Quindi:

h = -\frac{f'(x_0)}{f''(x_0)} \left(1 - \sqrt{1 - 2 \frac{f(x_0) f''(x_0)}{f'(x_0)^2}}\right). \tag{13}

Quando f(x_0) \approx 0, abbiamo \sqrt{1 - 2 \frac{f(x_0) f''(x_0)}{f'(x_0)^2}} \approx 1. Per evitare significativi errori di floating point, possiamo riscrivere l’espressione usando l’identità:

1 - \sqrt{A} = \frac{1 - \sqrt{A}}{1 + \sqrt{A}} (1 + \sqrt{A}) = \frac{1 - A}{1 + \sqrt{A}}.

Applicandola otteniamo:

1 - \sqrt{1 - 2 \frac{f(x_0) f''(x_0)}{f'(x_0)^2}} = \frac{2 \frac{f(x_0) f''(x_0)}{f'(x_0)^2}}{1 + \sqrt{1 - 2 \frac{f(x_0) f''(x_0)}{f'(x_0)^2}}}.

Sostituendo questo risultato nell’equazione (h), otteniamo:

h = -\frac{1}{f'(x_0)} \frac{2 f(x_0)}{1 + \sqrt{1 - 2 \frac{f(x_0) f''(x_0)}{f'(x_0)^2}}}.

Il che porta al metodo numerico:

\boxed{ x_{k+1} = x_k - \frac{f(x_k)}{f'(x_k)} \left( \dfrac{2}{1 + \sqrt{1 - 2 \dfrac{f(x_k) f''(x_k)}{f'(x_k)^2}}} \right) }. \tag{14}

  • Ordine di Convergenza: Il metodo di Müller ha un ordine di convergenza circa 1.84. Questo è superiore al metodo delle secanti (1.62), ma inferiore rispetto al metodo di Newton (2).

  • Instabilità numerica: Se i punti scelti sono troppo vicini tra loro, il metodo può essere soggetto a instabilità numeriche e errori di arrotondamento.

Esempio Numerico

Un semplice esempio di applicazione del metodo di Müller potrebbe essere la ricerca di uno zero della funzione f(x). Partendo da tre valori iniziali vicini alla radice stimata, il metodo approssima progressivamente il valore della radice con una parabola.

# Esempio in R
f <- function(x) x^3 - 2*x - 5

# Inizializzazione dei punti
x0 <- 3; x1 <- 4; x2 <- 2

# Iterazione del metodo di Muller
for (i in 1:5) {
  h1 <- x1 - x0
  h2 <- x2 - x1
  delta1 <- (f(x1) - f(x0)) / h1
  delta2 <- (f(x2) - f(x1)) / h2
  a <- (delta2 - delta1) / (h2 + h1)
  b <- a*h2 + delta2
  c <- f(x2)
  radice <- sqrt(b^2 - 4*a*c)

  # Scelta del segno per garantire stabilità
  if (abs(b + radice) > abs(b - radice)) {
    denominatore <- b + radice
  } else {
    denominatore <- b - radice
  }

  # Nuovo punto
  dx <- -2*c / denominatore
  x3 <- x2 + dx

  # Aggiornamento dei punti
  x0 <- x1
  x1 <- x2
  x2 <- x3

  print(x3) # Stampa dell'approssimazione della radice
}
[1] 2.111111
[1] 2.094815
[1] 2.094552
[1] 2.094551
[1] 2.094551

Metodo di Halley

Il metodo presentato in Equazione 14 presenta una limitazione: se il termine sotto radice è negativo, non possiamo applicarlo, poiché otterremmo numeri complessi. In tal caso, la parabola tangente non interseca l’asse x. Per risolvere questo problema, utilizziamo lo sviluppo di Taylor della funzione g(z) = \sqrt{1 + z}:

g(z) = g(0) + z g'(0) + \mathcal{O}(z^2), \quad g'(z) = \frac{1}{2\sqrt{1+z}}.

da cui g'(0)=1/2 e

\sqrt{1+z} = 1 + \frac{z}{2}+ \mathcal{O}(z^2),

Da questo sviluppo, otteniamo l’approssimazione:

\sqrt{1 - 2\frac{f(x_0)f''(x_0)}{f'(x_0)^2}} \approx 1-\frac{f(x_0)f''(x_0)}{f'(x_0)^2}.

Sostituendo questa approssimazione nel passo di Müller in Equazione 14 e semplificando, giungiamo al seguente metodo numerico:

\boxed{ x_{k+1} = x_k - \frac{f(x_k)f'(x_k)}{f'(x_k)^2 - \frac{1}{2}f(x_k)f''(x_k)} } \tag{15}

Questo è il Metodo di Halley.

Vale il seguente teorema:

Teorema 5 (convergenza del Metodo di Halley) Sia f \in C^{3}([a,b]) e \alpha una radice semplice di f(x), cioè f'(\alpha) \neq 0. Se x_0 \in (\alpha - \epsilon, \alpha + \epsilon) con \epsilon sufficientemente piccolo, allora la successione generata dallo schema Equazione 15 soddisferà le seguenti proprietà:

  1. x_k \in (\alpha - \epsilon, \alpha + \epsilon) per k = 1, 2, 3, \ldots

  2. \left\lvert x_{k+1} - \alpha\right\rvert \leq C \left\lvert x_k - \alpha\right\rvert^{2} per k = 1, 2, 3, \ldots

Dimostrazione. Sviluppando con la serie di Taylor

\begin{aligned} 0 = f(\alpha) &= f(x_k) + f'(\zeta_k)(\alpha -x_k), \\ 0 = f(\alpha) &= f(x_k) + f'(x_k)(\alpha -x_k) + f''(\omega_k)\dfrac{(x_k-\alpha)^{2}}{2}, \\ 0 = f(\alpha) &= f(x_k) + f'(x_k)(\alpha -x_k) + f''(x_k)\dfrac{(x_k-\alpha)^{2}}{2} + f'''(\eta_{k})\dfrac{(x_k-\alpha)^{3}}{6}, \end{aligned}

definendo \epsilon_k = x_k-\alpha

\begin{aligned} {\color{cyan}f(x_k)} &= f'(\zeta_k)\epsilon_k, \\ {\color{red}f(x_k)} &= f'(x_k)\epsilon_k - f''(\omega_k)\dfrac{\epsilon_k^{2}}{2}, \\ {\color{blue}f(x_k)} &= f'(x_k)\epsilon_k - f''(x_k)\dfrac{\epsilon_k^{2}}{2} + f'''(\eta_{k})\dfrac{\epsilon_k^{3}}{6}, \end{aligned}

e utilizzando Equazione 15:

\begin{aligned} \epsilon_{k+1} = x_{k+1}-\alpha &= x_k-\alpha -\dfrac{f(x_k)f'(x_k)} {f'(x_k)^2-\frac{1}{2}{\color{cyan}f(x_k)}f''(x_k)} \\ &= \dfrac{ (x_k-\alpha)\left(f'(x_k)^2-\frac{1}{2}{\color{red}f(x_k)}f''(x_k)\right) -{\color{blue}f(x_k)}f'(x_k) } {f'(x_k)^2-\frac{1}{2}{\color{cyan}f(x_k)}f''(x_k)} \\ &= \dfrac{\epsilon_k\left(f'(x_k)^2-\frac{1}{2}{\color{red}f(x_k)}f''(x_k)\right)-f'(x_k) \left({\color{blue}f'(x_k)\epsilon_k - f''(x_k)\dfrac{\epsilon_k^{2}}{2} + f'''(\eta_{k})\dfrac{\epsilon_k^{3}}{6}}\right)} {f'(x_k)^2-\frac{1}{2}{\color{cyan}f(x_k)}f''(x_k)} \\ &= \dfrac{ -\dfrac{\epsilon_k}{2}{\color{red}f(x_k)}f''(x_k) +f'(x_k) f''(x_k)\dfrac{\epsilon_k^{2}}{2} -f'(x_k)f'''(\eta_{k})\dfrac{\epsilon_k^{3}}{6} } {f'(x_k)^2-\frac{1}{2}{\color{cyan}f(x_k)}f''(x_k)} \\ &= \dfrac{ -\dfrac{\epsilon_k}{2} \left({\color{red} f'(x_k)\epsilon_k - f''(\omega_k)\dfrac{\epsilon_k^{2}}{2} }\right) f''(x_k) +f'(x_k) f''(x_k)\dfrac{\epsilon_k^{2}}{2} -f'(x_k)f'''(\eta_{k})\dfrac{\epsilon_k^{3}}{6} } {f'(x_k)^2-\frac{1}{2}{\color{cyan}f(x_k)}f''(x_k)} \\ &= \dfrac{ \frac{1}{4} f''(x_k)f''(\omega_k) -\frac{1}{6}f'(x_k)f'''(\eta_{k}) } {f'(x_k)^2-\frac{1}{2}f'(\zeta_k)f''(x_k)\epsilon_k}\epsilon_k^{3} \end{aligned}

da cui

\epsilon_{k+1} = \dfrac{3f''(x_k)f''(\omega_k)-2f'(x_k)f'''(\eta_{k})} {2f'(x_k)^2-f'(\zeta_k)f''(x_k)\epsilon_k} \frac{\epsilon_k^3}{6}.

Prendendo in valore assoluto entrambi i membri e ponendo

\begin{aligned} M_0 &= \sup_{\eta\in[a,b]}\left\lvert f(\eta)\right\rvert \\ M_1 &= \sup_{\eta\in[a,b]}\left\lvert f'(\eta)\right\rvert \\ M_2 &= \sup_{\eta\in[a,b]}\left\lvert f''(\eta)\right\rvert \\ M_3 &= \sup_{\eta\in[a,b]}\left\lvert f'''(\eta)\right\rvert \\ M &= \frac{1}{6}\left[3\sup_{\eta\in[a,b]}\left\lvert f''(\eta)\right\rvert^2 +2\left(\sup_{\eta\in[a,b]}\left\lvert f(\eta)\right\rvert\right) \left(\sup_{\eta\in[a,b]}\left\lvert f'''(\eta)\right\rvert\right)\right] \\ &= \frac{1}{6}\left[3M_2^2+2M_0M_3\right] \end{aligned}

otteniamo (se \epsilon_k è sufficientemente piccolo)

\left\lvert\epsilon_{k+1}\right\rvert \leq \left\lvert\epsilon_k\right\rvert^{3}\dfrac{M}{2\left\lvert f'(x_k)\right\rvert-M_1M_2\left\lvert\epsilon_k\right\rvert}.

Poiché \alpha è supposto zero semplice si ha che \left\lvert f'(\alpha)\right\rvert>0. Inoltre, per la continuità di f', esisterà \epsilon tale che 2\left\lvert f'(x)\right\rvert>(3/2)\left\lvert f'(\alpha)\right\rvert per ogni x\in(\alpha-\epsilon,\alpha+\epsilon) e M_1M_2\left\lvert\epsilon_k\right\rvert<\left\lvert f'(\alpha)\right\rvert/2. Se x_k\in(\alpha-\epsilon,\alpha+\epsilon) allora avremo

\left\lvert\epsilon_{k+1}\right\rvert \leq \left\lvert\epsilon_k\right\rvert^{3}\dfrac{M}{2\left\lvert f'(x_k)\right\rvert-M_1M_2\left\lvert\epsilon_k\right\rvert} \leq \left\lvert\epsilon_k\right\rvert^{3}\dfrac{M}{\left\lvert f'(\alpha)\right\rvert}. \tag{16}

Eventualmente riducendo \epsilon possiamo supporre \epsilon\dfrac{M}{\left\lvert f'(\alpha)\right\rvert}<1, cosicché dalla Equazione 16 segue il punto 1,

\left\lvert x_{k+1}-\alpha\right\rvert \leq \epsilon, \quad\implies\quad x_{k+1} \in (\alpha-\epsilon,\alpha+\epsilon).

Inoltre ponendo C=M/\left\lvert f'(\alpha)\right\rvert si trova verificato il punto 2.

Velocità di Convergenza

Consideriamo una successione \{x_0, x_1, x_2, \ldots\} che converge a \alpha, ossia:

\alpha = \lim_{k \to \infty} x_k.

Un aspetto fondamentale è determinare la velocità con cui la successione converge a \alpha. Questa analisi permette di confrontare l’efficacia di diversi metodi numerici nel risolvere lo stesso problema.

Q-convergenza

Da una successione \{x_0,x_1,x_2,\ldots\} convergente ad \alpha diremo che

(A)\quad Q-lineare:

La succesione converge Q-linearmente se per ogni k:

\left\lvert x_{k+1}-\alpha\right\rvert \leq C \left\lvert x_k-\alpha\right\rvert, \qquad C<1

(B)\quad Q-ordine p:

La succesione converge con Q-ordine p se per ogni k:

\left\lvert x_{k+1}-\alpha\right\rvert \leq C \left\lvert x_k-\alpha\right\rvert^p

(non serve limitare C<1)

(A')\quad Q-superlineare:

La succesione converge Q-superlinearmente se per ogni k:

\left\lvert x_{k+1}-\alpha\right\rvert \leq C_k \left\lvert x_k-\alpha\right\rvert \qquad\textrm{con}\qquad \lim_{k\to\infty} C_k = 0.

Naturamente confrontando le velocità di convergenza vale

(A)\leq(A')\leq (B)

e nella classe (B) se p>p' allora il meotodo di ordine-p e più veloce del metodo di ordine p'.

R-Convergenza

La Q-convergenza è un criterio piuttosto restrittivo che non include le successioni con convergenza non monotona. Prendiamo, ad esempio, la successione x_k = \frac{1}{2^k}, che è Q-linearmente convergente a 0. In questo caso, abbiamo:

\left\lvert x_{k+1} - 0\right\rvert = \frac{1}{2^{k+1}} \leq C \left\lvert x_k - 0\right\rvert = \frac{C}{2^k},

scegliendo C = \frac{1}{2}, la disuguaglianza risulta vera per ogni k.

Tuttavia, consideriamo la successione:

y_k = \begin{cases} \dfrac{1}{2^k} & \text{se } k \text{ è dispari} \\ 0 & \text{se } k \text{ è pari} \end{cases}

Questa successione, pur essendo convergente, non è Q-convergente. Infatti, per k pari, abbiamo:

\frac{1}{2^{k+1}} = \left\lvert y_{k+1} - 0\right\rvert \not\leq C \left\lvert y_k - 0\right\rvert = C \cdot 0 = 0, \quad \text{per ogni } C > 0.

Da una successione \{x_0,x_1,x_2,\ldots\} convergente ad \alpha diremo che

(A)\quad R-lineare:

La succesione converge R-linearmente se esiste successione y_k take che \left\lvert x_k-\alpha\right\rvert\leq M\left\lvert y_k-\alpha\right\rvert e per ogni k:

\left\lvert y_{k+1}-\alpha\right\rvert \leq C \left\lvert y_k-\alpha\right\rvert, \qquad C<1

(B)\quad R-ordine p:

La succesione converge con R-ordine p se esiste successione y_k take che \left\lvert x_k-\alpha\right\rvert\leq M\left\lvert y_k-\alpha\right\rvert e per ogni k:

\left\lvert y_{k+1}-\alpha\right\rvert \leq C \left\lvert y_k-\alpha\right\rvert^p

(non serve limitare C<1)

(A')\quad R-superlineare:

La succesione converge Q-superlinearmente se esiste successione y_k take che \left\lvert x_k-\alpha\right\rvert\leq M\left\lvert y_k-\alpha\right\rvert e per ogni k:

\left\lvert y_{k+1}-\alpha\right\rvert \leq C_k \left\lvert y_k-\alpha\right\rvert \qquad\textrm{con}\qquad \lim_{k\to\infty} C_k = 0.

Con queste definizioni in base alle diuscussioni precedenti possiamo dire che:

  • Il metodo dicotomico è R-lineare
  • Il metodo di Newton è Q-quadratico, cioè Q-convergente di ordine p=2
  • Il metodo di Halley è Q-cubico, cioè Q-convergente di ordine p=3
  • Il metodo delel secanti è R-convergente di ordine p=(1+\sqrt{5})/2=1.618\ldots

e quindo ordinando i metodi in base alla velocità di convergenza

[\textrm{dicotomico}]\leq [\textrm{secanti}]\leq [\textrm{Newton}]\leq [\textrm{Halley}]

Verifica Numerica dell’Ordine

Per verificare la corretta implementazione di un metodo numerico, è utile controllare che le iterazioni convergano alla soluzione con l’ordine teorico atteso. Supponiamo che un metodo numerico converga con ordine p. Vicino alla soluzione, possiamo aspettarci che:

\left\lvert x_{k+1} - \alpha\right\rvert \approx C \left\lvert x_k - \alpha\right\rvert^p.

In modo analogo, vale anche:

\left\lvert x_{k+2} - \alpha\right\rvert \approx C \left\lvert x_{k+1} - \alpha\right\rvert^p.

Facendo il rapporto tra queste due espressioni otteniamo:

\dfrac{\left\lvert x_{k+1} - \alpha\right\rvert}{\left\lvert x_{k+2} - \alpha\right\rvert} \approx \dfrac{\left\lvert x_k - \alpha\right\rvert^p}{\left\lvert x_{k+1} - \alpha\right\rvert^p}.

Prendendo il logaritmo di entrambi i membri, otteniamo:

\log \left( \dfrac{\left\lvert x_{k+1} - \alpha\right\rvert}{\left\lvert x_{k+2} - \alpha\right\rvert} \right) \approx p \log \left( \dfrac{\left\lvert x_k - \alpha\right\rvert}{\left\lvert x_{k+1} - \alpha\right\rvert} \right),

da cui è possibile ricavare una stima numerica per l’ordine p del metodo.

\boxed{ p \approx\dfrac{\log\dfrac{\left\lvert x_{k+1}-\alpha\right\rvert}{\left\lvert x_{k+2}-\alpha\right\rvert} }{\log\dfrac{\left\lvert x_{k}-\alpha\right\rvert}{\left\lvert x_{k+1}-\alpha\right\rvert}}, }

Esempio con metodo Newton calcolo \sqrt{2} (usando f(x)=x^2-2) e stima ordine:

x_{k+1} = x_k - \dfrac{x_k^2-2}{2x_k}=\dfrac{x_k}{2}+\dfrac{1}{x_k}

le iterate sono state calcolate in doppia precisione.

Tabella 1: Iterate
iter x_k f(x_k) e_k=\left\lvert x_k-\sqrt{2}\right\rvert p
1 2 2 0.585786438 —
2 \dfrac{3}{2} \dfrac{1}{4} 0.0857864376269 —
3 \dfrac{17}{12} \dfrac{1}{144} 0.00245310429357 1.85
4 \dfrac{577}{408} \dfrac{1}{166464} 2.1239\cdot 10^{-6} 1.98
5 1.4142135623746899106 4.51\cdot 10^{-12} 1.594861\cdot 10^{-12} 1.999
6 1.4142135623730950488 -4.77\cdot 10^{-21} -1.68872421\cdot 10^{-21} 1.465

si può notare che quando l’errore e molto piccolo gli errori di arrotondamento floating point diventano dominanti mancando cifre per calcolare la differenza tra a soluzione esatta e approssimata. Quindi la stima dell’errore non è più buona (vedi iterata 6).

Radici multiple e rallentamento della convergenza

Quando ci si trova a risolvere equazioni non lineari del tipo f(x) = 0 utilizzando il metodo di Newton, è noto che il metodo presenta convergenza quadratica se la radice è semplice. Tuttavia, in presenza di radici multiple, la convergenza diventa lineare, causando un significativo rallentamento. Esploriamo come gestire queste situazioni e recuperare una convergenza di ordine superiore.

Rallentamento della Convergenza con Radici Multiple

Quando una radice r ha molteplicità m, cioè quando f(x) = (x - r)^m g(x) con g(r) \neq 0, la convergenza del metodo di Newton rallenta. La ragione è che, vicino a una radice multipla, f(x) e le sue derivate si annullano in modo più rapido. In particolare, il metodo di Newton per una radice multipla mostra convergenza lineare con un fattore di riduzione proporzionale a 1 - 1/m.

Esempio in R

Supponiamo di avere una funzione con radice doppia f(x) = (x - 2)^2 + \epsilon, dove \epsilon è un piccolo termine aggiunto. Vediamo come si comporta il metodo di Newton.

newton_method <- function(f, df, initial_guess, tol = 1e-6, max_iter = 100) {
  x <- initial_guess
  for (i in 1:max_iter) {
    x_new <- x - f(x) / df(x)  # Calcola il nuovo punto
    cat(sprintf("Iterazione %d: x = %g\n", i, x_new))
    if ( abs(x_new - x) < tol || abs(f(x_new)) < tol ) {
      return(x_new)  # Restituisce la soluzione se la differenza è inferiore alla tolleranza
    }
    x <- x_new
  }
  warning("La radice non è stata trovata dopo il numero massimo di iterazioni.")
  return(NULL)  # Se non si trova la radice
}
# Definizione della funzione e derivata
epsilon <- 1e-8
f <- function(x) (x - 2)^2 + epsilon
df <- function(x) 2 * (x - 2)

# Applichiamo il metodo di Newton
result <- newton_method(f, df, initial_guess = 4, tol=1e-4)
Iterazione 1: x = 3
Iterazione 2: x = 2.5
Iterazione 3: x = 2.25
Iterazione 4: x = 2.125
Iterazione 5: x = 2.0625
Iterazione 6: x = 2.03125
Iterazione 7: x = 2.01562
Iterazione 8: x = 2.00781
print(result)
[1] 2.007812

da notare che ho usato una tolleranza molto alta. Se uso una tolleranza più spinta quando arrivo vicoino alla radice ho dei rapporti \approx 0 / \approx 0 e mi riallontano dalla radice:

result <- newton_method(f, df, initial_guess = 4, tol=1e-8, max_iter=40)
Iterazione 1: x = 3
Iterazione 2: x = 2.5
Iterazione 3: x = 2.25
Iterazione 4: x = 2.125
Iterazione 5: x = 2.0625
Iterazione 6: x = 2.03125
Iterazione 7: x = 2.01562
Iterazione 8: x = 2.00781
Iterazione 9: x = 2.00391
Iterazione 10: x = 2.00195
Iterazione 11: x = 2.00097
Iterazione 12: x = 2.00048
Iterazione 13: x = 2.00023
Iterazione 14: x = 2.00009
Iterazione 15: x = 1.99999
Iterazione 16: x = 2.00074
Iterazione 17: x = 2.00036
Iterazione 18: x = 2.00017
Iterazione 19: x = 2.00005
Iterazione 20: x = 1.99993
Iterazione 21: x = 2.00004
Iterazione 22: x = 1.9999
Iterazione 23: x = 2
Iterazione 24: x = 2.0018
Iterazione 25: x = 2.0009
Iterazione 26: x = 2.00044
Iterazione 27: x = 2.00021
Iterazione 28: x = 2.00008
Iterazione 29: x = 1.99998
Iterazione 30: x = 2.00023
Iterazione 31: x = 2.00009
Iterazione 32: x = 1.99999
Iterazione 33: x = 2.00057
Iterazione 34: x = 2.00028
Iterazione 35: x = 2.00012
Iterazione 36: x = 2.00002
Iterazione 37: x = 1.99975
Iterazione 38: x = 1.99989
Iterazione 39: x = 1.99999
Iterazione 40: x = 2.00079
Warning in newton_method(f, df, initial_guess = 4, tol = 1e-08, max_iter = 40):
La radice non è stata trovata dopo il numero massimo di iterazioni.
print(result)
NULL

Correzione del Metodo di Newton per Radici Multiple

Se conosciamo la molteplicità della radice m, possiamo modificare il metodo di Newton per recuperare la convergenza quadratica.

newton_multiple <- function(f, df, x0, m, tol = 1e-4, max_iter = 100) {
  x <- x0
  iter <- 0
  repeat {
    x_new <- x - m * f(x) / df(x)
    iter <- iter + 1
    cat(sprintf("Iterazione %d: x = %g\n", iter, x_new))
    if ( abs(x_new - x) < tol || abs(f(x_new)) < tol || iter >= max_iter) break
    x <- x_new
  }
  return(list(root = x, iterations = iter))
}

# Applichiamo il metodo corretto per la radice multipla
result_multiple <- newton_multiple(f, df, x0 = 4, m = 2)
Iterazione 1: x = 2
print(result_multiple)
$root
[1] 4

$iterations
[1] 1

Iterazioni di punto fisso

Vogliamo ora dare una formulazione più generale alle iterazioni del tipo Newton-Raphson.

Definizione 1 (Punto fisso) Data una funzione g:\Bbb{R}\mapsto \Bbb{R} diremo che \alpha è un punto fisso se vale:

\alpha = g(\alpha).

Cercare gli zeri di una funzione f(x) è equivalente a cercare i punti fissi della funzione g(x) = x-f(x).

Un modo molto naturale per cercare i punti fissi è quello di usare delle iterazioni. Ad esempio dato un punto iniziale x_0 costruendo la successione

\boxed{ x_{k+1} =g(x_k),\qquad k=1,2,\ldots } \tag{17}

possiamo chiederci quando la successione \{x_{i}\}_{i=0}^{\infty} è convergente e se converge a \alpha. Ad esempio lo schema di Newton~Equazione 3 può essere messo nella forma Equazione 17 ponendo

g(x) = x-\dfrac{f(x)}{f'(x)},

fondamentale per lo studio di schemi iterativi della forma Equazione 17 è il seguente teorema:

Teorema 6 (Convergenza punto fisso) Sia g(x) funzione continua e \alpha un suo punto fisso, cioè \alpha=g(\alpha). Sia g(x) soddisfacente la seguente condizione

\left\lvert g(x)-g(x')\right\rvert \leq L \left\lvert x-x'\right\rvert, \tag{18}

per ogni x, x' in un intervallo chiuso I\equiv[\alpha-\rho,\alpha+\rho], dove la costante L soddisfa:

0 \leq L < 1,

allora

  • (i)\quad Posto x_0\in I allora tutte le iterate x_k definite da Equazione 17 stanno in I, cioè \alpha-\rho \leq x_k \leq \alpha+\rho, \qquad k=1,2,\ldots

  • (ii)\quad (esistenza) le iterate convergono ad \alpha, \lim_{k\mapsto\infty} x_k = \alpha.

  • (iii)\quad (unicità) \alpha è l’unico punto fisso di g(x) in [x_0-\rho,x_0+\rho].

Dimostrazione.

  • (i)\quad basta osservare che \begin{aligned} \left\lvert\alpha-x_k\right\rvert & = \left\lvert g(\alpha)-g(x_{k-1})\right\rvert, \\ & \leq L \left\lvert\alpha-x_{k-1}\right\rvert = L \left\lvert g(\alpha)-g(x_{k-2})\right\rvert, \\ & \leq L^{2} \left\lvert\alpha-x_{k-2}\right\rvert = L^{2} \left\lvert g(\alpha)-g(x_{k-3})\right\rvert,\\ & \vdots \\ & \leq L^{k} \left\lvert\alpha-x_0\right\rvert < \rho, \end{aligned} \tag{19} e quindi x_k è nell’intervallo [\alpha-\rho,\alpha+\rho].

  • (ii)\quad sempre dalla Equazione 19 e dal fatto che 0\leq L<1 segue che \lim_{k\mapsto\infty}\left\lvert\alpha-x_k\right\rvert=0.

  • (iii)\quad siano \alpha e \beta due punti fissi allora dalla Equazione 18 abbiamo \left\lvert\alpha-\beta\right\rvert = \left\lvert g(\alpha)-g(\beta)\right\rvert \leq L \left\lvert\alpha-\beta\right\rvert < \left\lvert\alpha-\beta\right\rvert. Questa contraddizione implica che \alpha=\beta.

Osservazione

Se la funzione g(x) è differenziabile possiamo scrivere

\begin{aligned} g(x)-g(x') &= g'(\eta)(x-x'),\qquad \eta \in (x,x') \\ \left\lvert g(x)-g(x')\right\rvert &\leq \left\lvert g'(\eta)\right\rvert\left\lvert x-x'\right\rvert, \end{aligned}

e quindi se abbiamo

\left\lvert g'(\zeta)\right\rvert\leq L < 1, \qquad \zeta\in [\alpha-\rho,\alpha+\rho],

dove \alpha è un punto fisso di g(x) allora per il teorema precedente le iterazioni Equazione 17 convergono ad \alpha.

Esempio

Considerando il metodo di Newton-Raphson come uno schema a punto fisso abbiamo

g(x) = x-\dfrac{f(x)}{f'(x)}, \qquad g'(x) = \dfrac{f(x)f''(x)}{(f'(x))^{2}},

Se \alpha è una radice semplice di f(x) cioè f(\alpha)\neq 0 allora

g'(\alpha) = \dfrac{f(\alpha)f''(\alpha)}{(f'(\alpha))^{2}} = 0,

se f(x) è una funzione C^{2} allora g'(x) è continua nell’intorno di \alpha e scegliendo un intorno sufficientemente piccolo di \alpha avremo

\left\lvert g'(x)\right\rvert \leq L < 1, \qquad x\in [\alpha-\rho,\alpha+\rho]

e quindi per il teorema precedente lo schema di Newton-Raphson è convergente perlomeno purché si parti da una approssimazione sufficientemente vicina alla radice e la radice sia semplice.

Osservazione

Se \alpha è un punto fisso di g(x) e inoltre g'(\alpha)=0 è chiaro che nella iterazione Equazione 17 più ci avviciniamo alla radice \alpha più piccola è la costante L e quindi più rapidamente ci avviciniamo al punto fisso. Supponiamo ora che g''(x) esista e sia continua, allora sviluppando con Taylor attorno al punto fisso \alpha la funzione g(x) otteniamo:

g(x) = g(\alpha) + 0 + \dfrac{(x-\alpha)^{2}}{2}g''(\eta),

e da questa

\left\lvert x_{k+1}-\alpha\right\rvert = \left\lvert g(x_k)-g(\alpha)\right\rvert = \left\lvert\dfrac{(x_k-\alpha)^{2}}{2}g''(\eta_k)\right\rvert\leq \left\lvert\dfrac{g''(\eta_k)}{2}\right\rvert\left\lvert x_k-\alpha\right\rvert^{2}. \tag{20}

Da questa equazione segue che l’errore alla iterata successiva è proporzionale al quadrato dell’errore precedente. In tal caso diremo che il metodo è del secondo ordine. Quindi in particolare il metodo di Newton-Raphson è del secondo ordine nel caso di radici semplici.

Sia ora I\equiv [\alpha-\rho,\alpha+\rho] un intervallo in cui lo schema è convergente e

\left\lvert\dfrac{g''(x)}{2}\right\rvert \leq M, \qquad x\in I

allora dalla Equazione 20

\begin{aligned} \left\lvert x_k-\alpha\right\rvert & \leq M\left\lvert x_{k-1}-\alpha\right\rvert^{2}, \\ & \leq M\left( M\left\lvert x_{k-2}-\alpha\right\rvert^{2}\right)^{2}= M\cdot M^{2} \left\lvert x_{k-2}-\alpha\right\rvert^{4}, \\ & \leq M\cdot M^{2} \left( M \left\lvert x_{k-3}-\alpha\right\rvert^{2}\right)^{4} = M\cdot M^{2}\cdot M^{4} \left\lvert x_{k-3}-\alpha\right\rvert^{8}, \\ & \vdots \\ & \leq M^{1+2+4+8+\cdots+2^{k-1}} \left\lvert x_0-\alpha\right\rvert^{2^{k}}, \end{aligned}

osserviamo che 1+2+4+\cdots+2^{k-1}=2^{k}-1 e quindi

\left\lvert x_k-\alpha\right\rvert \leq \left(M\left\lvert x_0-\alpha\right\rvert\right)^{2^{k-1}}\left\lvert x_0-\alpha\right\rvert.

Scegliamo x_0 sufficientemente vicino ad \alpha in modo che

M\left\lvert x_0-\alpha\right\rvert < 1,

allora possiamo calcolare il numero di iterazioni necessarie per ridurre l’errore iniziale \left\lvert x_0-\alpha\right\rvert di ad esempio 10^{-m} ponendo

\left(M\left\lvert x_0-\alpha\right\rvert\right)^{2^{k-1}} \approx 10^{-m}, \tag{21}

e prendendo il logaritmo da entrambe le parti

k \approx \dfrac{1}{\log_{10}{2}}\log\left( \dfrac{m}{-\log_{10}(M\left\lvert x_0-\alpha\right\rvert)} \right),

osserviamo che sempre dalla Equazione 21 otteniamo

m \approx 2^{k} \log_{10}\left(1/(M\left\lvert x_0-\alpha\right\rvert)\right), \tag{22}

m è proporzionale al numero di cifre esatte nella approssimazione di \alpha. quindi la Equazione 22 significa che ad ogni iterazione un metodo al secondo ordine circa raddoppia le cifre esatte.

Osservazione

Supponiamo che sia \alpha un punto fisso di g(x) e inoltre g\in C^{p} con

g'(\alpha) = g''(\alpha)= \cdots = g^{(p-1)}(\alpha) = 0,

per il teorema di Taylor

g(x) = g(\alpha) + \dfrac{(x-\alpha)^{p}}{p!}g^{(p)}(\eta),

anche qui avremo

\left\lvert x_{k+1} -\alpha\right\rvert = \left\lvert g(x_k)-g(\alpha)\right\rvert \leq \dfrac{\left\lvert g^{(p)}(\eta)\right\rvert}{p!} \left\lvert x-\alpha\right\rvert^{p}.

Da questa equazionesegue che l’errore alla iterata successiva è proporzionale alla p-esima potenza dell’errore precedente. In tal caso diremo che il metodo è di ordine p.

Esempio

Nel metodo di Newton-Raphson se \alpha è una radice multipla di f(x) cioè

f(x) = (x-\alpha)^{n} h(x), \qquad n>1

dove h(\alpha)\neq 0 segue che

\begin{aligned} f'(x) &= n(x-\alpha)^{n-1} h(x)+(x-\alpha)^{n} h'(x) \\ f''(x) &= n(n-1)(x-\alpha)^{n-2} h(x)+2n(x-\alpha)^{n-1} h'(x)+(x-\alpha)^{n} h''(x), \end{aligned}

e di conseguenza

\begin{aligned} g'(x) &= \dfrac{f(x)f''(x)}{(f'(x))^{2}} = \dfrac{h(x)\left( n(n-1)h(x)+2n(x-\alpha) h'(x) +(x-\alpha)^{2} h''(x) \right)}{\left(nh(x)+(x-\alpha) h'(x)\right)^{2} }, \\ g'(\alpha) &= \dfrac{n(n-1)h(\alpha)^{2}}{n^{2}h(\alpha)^{2}}= 1-\dfrac{1}{n}, \end{aligned}

quindi

\left\lvert g'(\alpha) \right\rvert = 1-\dfrac{1}{n} < 1

e lo schema iterativo converge localmente anche se ora la convergenza è solo al primo ordine.

Elenco metodi numerici per zeri di funzione

Newton-Raphson

Ordine 2

\boxed{ x_{k+1} = x_k - \dfrac{f(x_k)}{f'(x_k)},\qquad k=1,2,\ldots }

Newton per radici multiple. Basta applicare Newton alla funzione F(x) = f(x)/f'(x) ottenendo:

\boxed{ x_{k+1} = x_k - \dfrac{F(x_k)}{F'(x_k)}=x_k -\dfrac{f(x_k)f'(x_k)}{f'(x_k)^2-f(x_k)f''(x_k)} ,\qquad k=1,2,\ldots }

  • R. L.Burden & J. D. Faires
    Numerical Analysis (9th ed.)
    Brooks Cole (2010)

Secanti

Ordine \frac{1+\sqrt{5}}{2}=1.618\ldots

\boxed{ x_{k+1} = \dfrac{f(x_k)x_{k-1}-f(x_{k-1})x_k} {f(x_k)-f(x_{k-1})}, \qquad k=1,2,\ldots }

Secanti (variante Brent)

\boxed{ x_{k+1} = x_k - \dfrac{f(x_k)}{D_k},\quad D_k=\begin{cases} \frac{f(x_k)-f(x_{k-1})}{x_k-x_{k-1}} & \left\lvert x_k-x_{k-1}\right\rvert > h \\ \frac{f(x_k)-f(x_k-h)}{h} & \left\lvert x_k-x_{k-1}\right\rvert \leq h \end{cases} \qquad k=1,2,\ldots }

Parabole tangenti o metodo di Muller

Ordine 3

\boxed{ x_{k+1} = x_k -2\dfrac{f(x_k}{f'(x_k)}\bigg/ \left( 1+\sqrt{1-2\dfrac{f(x_k)f''(x_k)}{f'(x_k)^2}} \right) \qquad k=1,2,\ldots }

  • D. E. Müller
    A method for solving algebraic equations using an automatic computer
    Mathematical Tables and Other Aids to Computation, Vol. 10, No. 56, 208-215, (1956)

Metodo di Halley

Ordine 3

\boxed{ x_{k+1} = x_k -\dfrac{f(x_k)f'(x_k)}{f'(x_k)^2-\frac{1}{2}f(x_k)f''(x_k)} \qquad k=1,2,\ldots }

  • J. F. Traub
    • Iterative Methods for the Solution of Equations*
      Prentice-Hall, 1964.

Metodo di Steffensen

Ordine 2

\boxed{ x_{k+1} = x_k -\dfrac{f(x_k)^2}{f(x_k+f(x_k))-f(x_k)} \qquad k=1,2,\ldots }

  • J. F. Steffensen
    Remarks on Iteration
    Skandinavisk Aktuarietidskrift, 1933, 16(1), 64-72.

Metodo delle corde tangenti

\boxed{ z = x_k -\dfrac{f(x_k)}{f'(x_k)},\qquad x_{k+1} = z -\dfrac{f(z)}{f'(x_k)}, \qquad k=1,2,\ldots }

Metodo di Ostrowski

Ordine 4

\boxed{ z = x_k -\dfrac{f(x_k)}{f'(x_k)},\qquad x_{k+1} = z -\dfrac{f(z)(z-x_k)}{2f(z)-f(x_k)}, \qquad k=1,2,\ldots }

  • A. M. Ostrowski
    Solution of equations and systems of equations
    Academic Press, New York, 1960