University of Trento, Department of Industrial Engineering
21.1 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. Come in tutto il capitolo poniamo e_k = x_k-\alpha.
Q-convergenza
Data una successione \{x_k\} convergente ad \alpha diremo che:
(A) Q-lineare. La successione converge Q-linearmente se esiste C<1 tale che, per ogni k sufficientemente grande,
\left\lvert e_{k+1}\right\rvert \leq C \left\lvert e_k\right\rvert, \qquad C<1 .
(B) Q-ordine p. La successione converge con Q-ordine p>1 se esiste C>0 tale che, per ogni k sufficientemente grande,
\left\lvert e_{k+1}\right\rvert \leq C \left\lvert e_k\right\rvert^p .
Qui non serve richiedere C<1: è però necessario che l’errore sia abbastanza piccolo, precisamente che C\left\lvert e_{k_0}\right\rvert^{p-1}<1 per qualche k_0, altrimenti la disuguaglianza non implica affatto la convergenza. È il motivo per cui tutti i teoremi di convergenza quadratica o cubica del capitolo sono locali: richiedono x_0 in un intorno sufficientemente piccolo di \alpha.
(A’) Q-superlineare. La successione converge Q-superlinearmente se, per ogni k sufficientemente grande,
cioè una successione con Q-ordine p>1 è Q-superlineare, e una successione Q-superlineare è definitivamente Q-lineare per qualunque C<1 fissato. Le implicazioni inverse sono false: ad esempio x_k = k^{-k} converge a zero Q-superlinearmente, perché x_{k+1}/x_k\to0, ma non ha Q-ordine p per nessun p>1. All’interno della classe (B), se p>p' allora il metodo di ordine p è asintoticamente più veloce del metodo di ordine p'.
Una lettura operativa: con Q-ordine p il numero di cifre corrette viene moltiplicato per p a ogni iterazione. Con p=2 (Newton) le cifre raddoppiano, con p=3 (Halley) triplicano; con convergenza lineare di fattore C si guadagnano invece -\log_{10}C cifre per iterazione, cioè circa 0.3 cifre per la bisezione.
R-convergenza
La Q-convergenza è un criterio piuttosto restrittivo, che non include le successioni la cui decrescita non è monotona. Prendiamo, ad esempio, la successione x_k = \dfrac{1}{2^k}, che converge Q-linearmente a 0. In questo caso, abbiamo:
e scegliendo C = \dfrac{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} \\[1mm]
0 & \text{se } k \text{ è pari}
\end{cases}
Questa successione, pur convergendo a zero almeno altrettanto in fretta della precedente, non è Q-convergente. Infatti, per k pari,
\dfrac{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.
L’idea della R-convergenza è allora di non chiedere nulla al rapporto fra errori consecutivi, ma di limitare l’errore con una successione maggiorante che converga a zero con la velocità voluta.
Data una successione \{x_k\} convergente ad \alpha diremo che:
(A) R-lineare. La successione converge R-linearmente se esiste una successione di numeri positivi \eta_k tale che
ovvero, equivalentemente, se \left\lvert e_k\right\rvert\leq M q^{p^{k}} per opportuni M>0 e 0<q<1. È esattamente la forma della tesi del Teorema 20.4.
(A’) R-superlineare. La successione converge R-superlinearmente se esiste una successione di numeri positivi \eta_k tale che
La successione y_k dell’esempio precedente, che non è Q-convergente, è R-linearmente convergente: basta prendere \eta_k = 2^{-k}. In generale la Q-convergenza di un dato ordine implica la R-convergenza dello stesso ordine, ma non viceversa: la R-convergenza descrive l’andamento medio dell’errore, la Q-convergenza quello di ogni singolo passo.
Ordine ed efficienza dei metodi del capitolo
Con queste definizioni, in base alle discussioni precedenti possiamo dire che:
il metodo dicotomico (Sezione 20.3) è R-lineare con fattore 1/2 (non è Q-lineare: l’errore del punto medio può occasionalmente crescere);
dove i metodi collegati da \simeq hanno lo stesso ordine e vanno confrontati sul fattore asintotico.
La Figura 21.1 riassume questa scala affiancandole le due quantità che l’ordinamento non mostra: il numero \nu di valutazioni richieste da una iterazione e l’indice di efficienza p^{1/\nu} discusso qui sotto.
Figura 21.1: I metodi del capitolo ordinati per velocità di convergenza per iterazione, con il costo di una iterazione e l’indice di efficienza
ImportanteL’ordine non basta: l’indice di efficienza
L’ordinamento precedente confronta l’ordine per iterazione e non tiene conto del costo di una iterazione. Il confronto corretto usa l’indice di efficienza
E = p^{1/\nu},
dove \nu è il numero di valutazioni di funzione (e derivate) richieste da una iterazione: è il fattore di cui si riduce l’errore per singola valutazione. Si ottiene
cioè, a parità di valutazioni, l’ordinamento si ribalta: i metodi che non usano derivate e riutilizzano i valori già calcolati (secanti, Müller) battono quelli di ordine più alto, e Halley è solo marginalmente migliore di Newton. Si vedano le tabelle di Sezione 20.13.
Per i metodi di ordine p=1 l’indice vale 1 e non discrimina: fra di essi il confronto va fatto sul fattore di riduzione asintotico (1/2 per la bisezione, dipendente da f per corde e regula falsi).
21.2 Verifica numerica dell’ordine
Per verificare la corretta implementazione di un metodo numerico è utile controllare che le iterazioni convergano con l’ordine teorico atteso: un errore di programmazione (un segno sbagliato, una derivata approssimata male) si manifesta quasi sempre come una perdita di ordine, anche quando il metodo continua a convergere. Supponiamo che un metodo converga con ordine p. Vicino alla soluzione possiamo aspettarci che
\left\lvert e_{k+1}\right\rvert \approx C \left\lvert e_k\right\rvert^p,
\qquad\text{e analogamente}\qquad
\left\lvert e_{k+2}\right\rvert \approx C \left\lvert e_{k+1}\right\rvert^p.
Facendo il rapporto tra queste due espressioni la costante C, che non è nota, si elimina:
e la stima di p si avvicina rapidamente al valore teorico 2. All’iterata 6, però, l’errore scende al livello della precisione di macchina, l’iterazione diventa stazionaria (x_7=x_6) e la stima dell’ordine perde ogni significato: il valore riportato nell’ultima riga non ha alcun contenuto.
Questo comportamento è esattamente quello previsto in Sezione 20.4. Qui \epsilon_f\approx\varepsilon_{\rm mach}\cdot2\approx4.4\cdot10^{-16} e \left\lvert f'(\sqrt 2)\right\rvert=2\sqrt2, quindi la 20.1 prevede un’incertezza di circa 1.6\cdot10^{-16}, cioè meno di una unità di ultima posizione: è il livello a cui l’iterazione si blocca. Si osservi inoltre che il punto fisso raggiunto, 1.4142135623730949, non è il numero macchina più vicino a \sqrt2, che è 1.4142135623730951: differisce da esso di una unità di ultima posizione. Infine, anche il valore di riferimento usato per calcolare e_k è il numero macchina sqrt(2), a sua volta arrotondato: errori inferiori a un ulp non sono misurabili in questo modo e vanno valutati in precisione estesa.
NotaMorale
Non si deve mai iterare “fino a che l’errore è zero”. Un metodo di ordine elevato raggiunge la precisione di macchina in pochissimi passi: ogni iterazione successiva è lavoro sprecato che, nei casi peggiori (si veda l’esempio della radice doppia perturbata più avanti), può addirittura peggiorare il risultato.
Stima dell’ordine senza conoscere la radice
La 21.1 ha un difetto che la rende inservibile proprio quando servirebbe: richiede \alpha. Su un problema di prova la radice è nota e la formula va benissimo, ma la verifica dell’ordine serve soprattutto quando si collauda un’implementazione nuova su un problema nuovo — cioè esattamente quando \alpha è l’incognita. Serve una stima che usi solo le iterate.
L’idea: sostituire l’errore con l’incremento
L’unica quantità osservabile è l’incremento fra due iterate successive,
Il fattore fra parentesi si controlla facilmente nei due regimi che ci interessano.
Se la convergenza ha ordine p>1, allora e_{k+1}/e_k \approx C\left\lvert e_k\right\rvert^{p-1}\to0, e quindi
\Delta_k = -e_k\bigl(1+O(\left\lvert e_k\right\rvert^{p-1})\bigr) :
l’incremento è l’errore, a meno di un fattore che tende a 1.
Se la convergenza è lineare con fattore C<1, allora e_{k+1}/e_k\to C e
\Delta_k = -e_k\,(1-C)\bigl(1+o(1)\bigr) :
l’incremento non è l’errore, ma gli è proporzionale, con costante 1-C.
In entrambi i casi \left\lvert\Delta_k\right\rvert = \gamma\left\lvert e_k\right\rvert\bigl(1+o(1)\bigr) con \gamma costante. È questo il punto decisivo: la 21.1 è un rapporto fra logaritmi di rapporti di errori, e in un rapporto la costante \gamma si semplifica. Si può dunque sostituire e con \Delta senza alterare il limite:
\boxed{\;
p \approx \dfrac{\log\bigl(\left\lvert\Delta_{k+1}\right\rvert/\left\lvert\Delta_{k+2}\right\rvert\bigr)}
{\log\bigl(\left\lvert\Delta_{k}\right\rvert/\left\lvert\Delta_{k+1}\right\rvert\bigr)}\;}
\tag{21.2}
ovvero, esplicitando gli incrementi in termini di iterate,
p \approx \dfrac{\log\bigl(\left\lvert x_{k+2}-x_{k+1}\right\rvert/\left\lvert x_{k+3}-x_{k+2}\right\rvert\bigr)}
{\log\bigl(\left\lvert x_{k+1}-x_{k}\right\rvert/\left\lvert x_{k+2}-x_{k+1}\right\rvert\bigr)} .
Rispetto alla 21.1 serve una iterata in più — quattro consecutive anziché tre — e nient’altro: nessuna informazione a priori, nessuna soluzione di riferimento. È la forma usata nei test automatici di validazione delle librerie numeriche, e la stessa idea si ritrova nello studio sperimentale dell’ordine dei metodi per equazioni differenziali, dove la soluzione esatta non è quasi mai disponibile.
Verifica
Mettiamo alla prova la formula su tre metodi applicati alla stessa equazione x^2-2=0, scelti in modo da coprire i tre regimi: Newton (p=2), secanti (p=\varphi\approx1.618) e il metodo delle corde con pendenza congelata m=f'(2)=4, che converge linearmente con fattore \left\lvert 1-f'(\alpha)/m\right\rvert=1-\sqrt2/2\approx0.293. Per ciascuno calcoliamo entrambe le stime, quella che usa \alpha e quella che non la usa.
f <-function(x) x^2-2; df <-function(x) 2*xalpha <-sqrt(2)itera <-function(passo, x0, n) { x <- x0; out <- xfor (k in1:n) { x <-passo(x); out <-c(out, x) } out}secanti <-function(x0, x1, n) { out <-c(x0, x1)for (k in1:n) { x2 <- x1 -f(x1)*(x1-x0)/(f(x1)-f(x0)) out <-c(out, x2); x0 <- x1; x1 <- x2 } out}# stima che usa la radice esatta: eq. @eq-stima-ordinep_con <-function(x, alpha) { e <-abs(x - alpha); n <-length(e); p <-rep(NA_real_, n)for (k in3:n)if (all(e[(k-2):k] >0)) p[k] <-log(e[k-1]/e[k]) /log(e[k-2]/e[k-1]) p}# stima che usa i soli incrementi: eq. @eq-stima-ordine-deltap_senza <-function(x) { d <-abs(diff(x)); n <-length(x); p <-rep(NA_real_, n)if (n >=4)for (k in3:(n-1))if (all(d[(k-2):k] >0)) p[k] <-log(d[k-1]/d[k]) /log(d[k-2]/d[k-1]) p}confronta <-function(x, didascalia, righe =seq_along(x)) { knitr::kable(data.frame(k =seq_along(x),`$x_k$`=sprintf("%.15f", x),`$p$ con $\\alpha$`=sprintf("%.4f", p_con(x, alpha)),`$p$ dai soli $\\Delta$`=sprintf("%.4f", p_senza(x)),check.names =FALSE)[righe, ],row.names =FALSE, align ="r", caption = didascalia)}confronta(itera(function(x) x -f(x)/df(x), 2, 6), "Newton: ordine atteso $p=2$")
Newton: ordine atteso p=2
k
x_k
p con \alpha
p dai soli \Delta
1
2.000000000000000
NA
NA
2
1.500000000000000
NA
NA
3
1.416666666666667
1.8503
1.9681
4
1.414215686274510
1.9839
1.9995
5
1.414213562374690
1.9998
2.0000
6
1.414213562373095
NA
0.6296
7
1.414213562373095
NA
NA
confronta(secanti(1, 2, 6), "Secanti: ordine atteso $p=\\varphi\\approx1.618$")
Secanti: ordine atteso p=\varphi\approx1.618
k
x_k
p con \alpha
p dai soli \Delta
1
1.000000000000000
NA
NA
2
2.000000000000000
NA
NA
3
1.333333333333333
-5.7130
5.6789
4
1.400000000000000
0.8782
0.6585
5
1.414634146341463
2.0246
2.3375
6
1.414211438474870
1.5023
1.4935
7
1.414213562057320
1.6666
1.6650
8
1.414213562373095
1.6075
NA
confronta(itera(function(x) x -f(x)/4, 2, 8),"Corde con $m=4$: ordine atteso $p=1$", righe =1:9)
Corde con m=4: ordine atteso p=1
k
x_k
p con \alpha
p dai soli \Delta
1
2.000000000000000
NA
NA
2
1.500000000000000
NA
NA
3
1.437500000000000
0.6788
0.6375
4
1.420898437500000
0.9571
0.9458
5
1.416160345077515
0.9885
0.9852
6
1.414782814334998
0.9967
0.9958
7
1.414380211400584
0.9990
0.9988
8
1.414262365800194
0.9997
0.9996
9
1.414227855970503
0.9999
NA
Le tabelle dicono tre cose.
Le due stime coincidono non appena il metodo entra nel regime asintotico: per Newton entrambe danno 2.00, per le secanti entrambe si assestano attorno a 1.61, per le corde entrambe danno 1.00. La colonna che richiede \alpha non aggiunge nulla.
La formula funziona anche per p=1. È la conseguenza della cancellazione della costante \gamma=1-C vista sopra, e non è ovvia: la stima resta valida per i metodi lineari, dove misura correttamente p=1 (il fattoreC, se interessa, si legge direttamente dal rapporto \left\lvert\Delta_{k+1}\right\rvert/\left\lvert\Delta_k\right\rvert).
Nel transitorio iniziale le stime non significano nulla. Per le secanti le prime righe danno valori negativi o maggiori di 5: le iterate sono ancora lontane dalla radice e la relazione \left\lvert e_{k+1}\right\rvert\approx C\left\lvert e_k\right\rvert^p, che è asintotica, non vale.
AvvisoCome leggere una stima dell’ordine
Non fidarsi di un singolo numero. Si cerca un plateau: tre o quattro stime consecutive che si stabilizzano sullo stesso valore. Un valore isolato può essere un artefatto del transitorio o dell’arrotondamento.
Attenzione all’ultima riga. Sul plateau di arrotondamento la stima basata sull’errore smette semplicemente di essere definita, perché e_k diventa esattamente zero; quella basata sugli incrementi, invece, continua a produrre un numero. Nella tabella di Newton l’ultima riga ne è l’esempio: 0.63, un valore perfettamente plausibile e del tutto privo di significato. È il caso più insidioso, perché non dà alcun segnale d’allarme.
Finestra utile. Le stime attendibili sono quelle calcolate quando \left\lvert\Delta_k\right\rvert è già piccolo rispetto a \left\lvert x_k\right\rvert ma ancora molto maggiore di \varepsilon_{\rm mach}\left\lvert x_k\right\rvert: fra il transitorio e il plateau di Sezione 20.4. Con un metodo di ordine alto questa finestra può essere larga due o tre iterazioni soltanto, ed è il motivo per cui conviene fare la verifica con una tolleranza lasca, non spinta al limite della macchina.
Variante: stimare prima la radice
Un’alternativa è usare le iterate per stimare \alpha, e poi applicare la 21.1 come se la radice fosse nota. Lo strumento naturale è l’estrapolazione \Delta^2 di Aitken:
che elimina il termine dominante dell’errore e converge ad \alpha più in fretta della successione di partenza.
aitken <-function(x) { n <-length(x); a <-rep(NA_real_, n)for (k in1:(n-2)) { d1 <- x[k+1]-x[k]; d2 <- x[k+2]-x[k+1]; den <- d2 - d1if (is.finite(den) && den !=0) a[k] <- x[k] - d1^2/den } a}xn <-itera(function(x) x -f(x)/df(x), 2, 6)knitr::kable(data.frame(k =seq_along(xn),`$\\hat\\alpha_k$ (Aitken)`=sprintf("%.15f", aitken(xn)),`$\\lvert\\hat\\alpha_k-\\alpha\\rvert$`=sprintf("%.1e", abs(aitken(xn)-alpha)),check.names =FALSE),row.names =FALSE, align ="r",caption ="Stima della radice per estrapolazione di Aitken")
Stima della radice per estrapolazione di Aitken
k
\hat\alpha_k (Aitken)
\lvert\hat\alpha_k-\alpha\rvert
1
1.400000000000000
1.4e-02
2
1.414141414141414
7.2e-05
3
1.414213560532626
1.8e-09
4
1.414213562373095
0.0e+00
5
1.414213562373095
2.2e-16
6
NA
NA
7
NA
NA
Già alla terza riga \hat\alpha ha una decina di cifre corrette, e usandola al posto di \alpha nella 21.1 si riottengono le stesse stime della colonna «p con \alpha». Il prezzo è una maggiore sensibilità all’arrotondamento: il denominatore \Delta_{k+1}-\Delta_k è una differenza di differenze, e si annulla per cancellazione non appena le iterate si avvicinano alla precisione di macchina. Per la sola verifica dell’ordine la 21.2 è più robusta e va preferita; l’estrapolazione di Aitken resta utile quando serve davvero un valore accurato di \alpha a partire da una successione che converge lentamente.
21.3 Radici multiple e rallentamento della convergenza
Quando si risolvono equazioni non lineari f(x) = 0 con il metodo di Newton, è noto che il metodo ha convergenza quadratica se la radice è semplice. Tuttavia, in presenza di radici multiple, la convergenza diventa lineare, causando un significativo rallentamento. Vediamo come gestire queste situazioni e recuperare un ordine superiore.
Rallentamento della convergenza
Supponiamo che la radice r abbia molteplicità m, cioè che
f(x) = (x - r)^m g(x), \qquad g(r) \neq 0 .
Il fattore di riduzione dell’errore si calcola esplicitamente. Derivando,
Il metodo di Newton su una radice di molteplicità m ha dunque convergenza lineare con fattore di riduzione esattamente 1 - 1/m: per m=2 l’errore si dimezza a ogni passo, per m=3 si riduce di un terzo, e per m grande il metodo diventa inutilizzabile.
Si noti che il rallentamento si somma al problema di condizionamento discusso in Sezione 20.4. Su una radice di molteplicità m l’accuratezza raggiungibile è comunque limitata a O(\epsilon_f^{1/m}), quindi le iterazioni in più non servono a nulla.
Esempio in R
Usiamo un unico driver che restituisce la storia delle iterate, con la salvaguardia necessaria quando la derivata si annulla (cosa che su una radice multipla accade esattamente nella radice).
newton_storia <-function(f, df, x0, n, m =1) { x <- x0 out <-numeric(0)for (k inseq_len(n)) { d <-df(x)if (!is.finite(d) || d ==0) break# derivata nulla: non si può proseguire x <- x - m *f(x) / d out <-c(out, x) } out}
Consideriamo f(x) = (x-2)^2, che ha una radice doppia in x=2. In questo caso il passo di Newton si calcola a mano,
Il rapporto \left\lvert e_{k+1}\right\rvert/\left\lvert e_k\right\rvert vale costantemente 1/2: servono circa 3.3 iterazioni per cifra decimale, come per la bisezione. Su una radice semplice le cifre corrette sarebbero invece raddoppiate a ogni passo.
Una radice doppia perturbata
Una seconda difficoltà, distinta dalla precedente, si presenta quando la radice doppia viene perturbata. La funzione
non ha radici reali: il suo minimo vale \epsilon>0. Il metodo di Newton si avvicina comunque a x=2 finché il passo è dominato dal termine (x-2)/2, ma appena (x-2)^2 scende sotto \epsilon il passo diventa \approx-\epsilon/\bigl(2(x-2)\bigr), cioè un rapporto fra due quantità piccole, e le iterate si allontanano di nuovo dal minimo, oscillando in modo irregolare.
epsilon <-1e-8f_eps <-function(x) (x -2)^2+ epsilondf_eps <-function(x) 2* (x -2)st3 <-newton_storia(f_eps, df_eps, x0 =4, n =30)tail(data.frame(k =seq_along(st3), x = st3, f =f_eps(st3)), 12)
Le iterate vagano in un intorno di ampiezza \approx\sqrt{\epsilon}=10^{-4} senza mai stabilizzarsi. Questo illustra due cose: che su radici multiple o quasi multiple un criterio d’arresto basato su \left\lvert f(x)\right\rvert si attiva su un intero intervallo di ampiezza O(\sqrt{\epsilon}) (coerentemente con la 20.2), e che occorre sempre un limite al numero di iterazioni.
Correzione del metodo di Newton per radici multiple
Se la molteplicità m è nota, basta moltiplicare il passo per m:
Lo schema equivale ad applicare il metodo di Newton a f^{1/m}, che ha in r una radice semplice, e recupera quindi la convergenza quadratica. Per f(x)=(x-2)^2 la correzione individua la radice in un solo passo, dopodiché la derivata si annulla e il driver si arresta:
newton_storia(f, df, x0 =4, n =10, m =2)
[1] 2
Se la molteplicità non è nota, si può invece applicare il metodo di Newton alla funzione
u(x) = \dfrac{f(x)}{f'(x)},
che ha in r una radice semplice qualunque sia m. Lo schema risultante,
converge quadraticamente senza conoscere m, al prezzo del calcolo di f'' e di una maggiore sensibilità agli errori di arrotondamento (il denominatore è una differenza di quantità che tendono entrambe a zero). Si confronti la 21.5 con il passo di Halley dato dalla 20.22. Le due formule differiscono solo per il fattore 1/2 davanti a f f''.
Una terza possibilità, quando si sospetta una radice multipla ma non se ne conosce l’ordine, è stimare m dal rapporto osservato fra errori consecutivi: se \left\lvert e_{k+1}\right\rvert/\left\lvert e_k\right\rvert\to C, allora m\approx1/(1-C).
Riepilogo su ordine e radici multiple
Q-convergenza descrive il comportamento di ogni singolo passo, R-convergenza l’andamento medio dell’errore: la prima implica la seconda, non viceversa. La bisezione è R-lineare, le secanti sono R-convergenti di ordine \varphi, Newton è Q-quadratico, Halley Q-cubico.
L’ordine teorico va verificato numericamente con la 21.1 (o la sua variante senza \alpha): la perdita di ordine è il sintomo più affidabile di un’implementazione sbagliata.
La stima dell’ordine ha senso solo finché l’errore resta sopra il livello di arrotondamento: al di sotto, i numeri ottenuti sono privi di significato.
Su una radice di molteplicità m Newton degrada a convergenza lineare con fattore 1-1/m; si recupera l’ordine due o con la correzione 21.4, se m è noto, o applicando Newton a f/f'21.5 altrimenti.
Su radici multiple l’accuratezza raggiungibile è comunque limitata a O(\epsilon_f^{1/m}): nessuna correzione dello schema può superare questo limite, che dipende solo dal problema.
21.4 Esercizi
Bisezione. Quante iterazioni servono per localizzare con errore minore di 10^{-12} una radice contenuta in [1,2]? E in [0,100]? Verificare sperimentalmente la formula k_{\min}=\lceil\log_2((b-a)/\tau)\rceil.
Robustezza della bisezione. Si implementi la bisezione con il test f(c)f(a)<0 anziché con il confronto dei segni e la si applichi a f(x)=10^{-200}(x-1) su [0,2] in doppia precisione. Che cosa succede, e perché?
Ordine di Newton quando f''(\alpha)=0. Si applichi Newton a f(x)=x^3-8 partendo da x_0=3 e poi a f(x)=x^3+x partendo da x_0=1, stimando p con la 21.1 in entrambi i casi. Perché nel primo si ottiene p=2 e nel secondo p=3? Si giustifichi con la costante asintotica f''(\alpha)/(2f'(\alpha)).
Condizione di Fourier. Verificare che f(x)=e^x-2 soddisfa le ipotesi del Teorema 20.3 su [0,1] e che Newton partendo da x_0=1 genera una successione decrescente. Che cosa accade partendo da x_0=0?
Fallimento di Newton. Determinare numericamente, con una bisezione sulla funzione \phi(x_0)=x_1(x_0)+x_0, il valore critico x^\ast per cui Newton applicato a \arctan x entra in un ciclo di periodo due. Confrontare con il valore 1.3917452002707\ldots riportato in Sezione 20.6.4.
Secanti contro Newton. Risolvere \cos x = x con entrambi i metodi, contando le valutazioni di funzione necessarie per scendere sotto 10^{-14}. Il risultato conferma l’indice di efficienza p^{1/\nu}?
Regula falsi e varianti. Applicare la regula falsi a f(x)=x^{10}-1 su [0,1.3] e osservare la stagnazione dell’estremo destro. Implementare poi la variante Illinois (si dimezza f(a_k) ogni volta che a_k non viene aggiornato) e confrontare il numero di iterazioni.
Radici multiple. Per f(x)=(x-1)^3(x-2) si stimi numericamente la molteplicità della radice x=1 a partire dal rapporto osservato fra errori consecutivi, usando m\approx1/(1-C), e si verifichi che 21.4 con il valore stimato recupera la convergenza quadratica.
Accuratezza raggiungibile. Valutare p(x)=(x-1)^5 nella forma espansa in 1\pm10^{-3} e confrontare l’ampiezza della zona di rumore con la previsione O(\epsilon_f^{1/5}) della 20.2.
Parabola tangente e Halley. Implementare entrambi i metodi per f(x)=x^3-2x-5 e verificare numericamente che l’ordine è 3 per entrambi. Quale dei due è più robusto quando f'' è grande, e perché?
