Equazioni e sistemi nonlineari.

Alvise Sommariva

Università degli Studi di Padova Dipartimento di Matematica Pura e Applicata

21 novembre 2011

Equazioni non lineari: esempi.

▶ Risoluzione f (x) = 0 con x ∈ [a, b] ⊆ ℝ, f ∈ C([a, b]). ▶ Esempio 1: equazioni polinomiali p_{N}(x) = 0 con p_{N} polinomio

di grado N. Possibili problemi (nessuna soluzione in ℝ, soluzioni multiple).

▶ Esempio 2: f (x) = sin (x) − x. Soluzione unica poichè f decrescente (vedi derivata).

Equazioni non lineari: convergenza e ordine convergenza

▶ Metodo iterativo: genera successione {x_{n}}_{n∈ℕ} che si desidera

convergere a x^{∗} tale che f (x^{∗}) = 0. ▶ Ordine di convergenza del metodo iterativo: sia {x_{k}} una

successione convergente ad x^{∗} e sia e_{k} = x_{k} − x^{∗} l’errore al passo k. Se esiste un numero p > 0 e una costante C ≠ 0 tale che

lim_{k→∞} \frac{|e_{k+1}|}{|e_{k}|^{p}} = C

allora p è chiamato ordine di convergenza della successione e C è la costante asintotica di errore. Per p = 1 la convergenza di dice lineare, per p = 2 si dice quadratica.

Equazioni non lineari: esempio ordine di convergenza

Supponiamo

\frac{|e_{k+1}|}{|e_{k}|^{p}} = \frac{1}{10}, |e_{0}| = 1

o equivalentemente

|e_{k+1}| = \frac{1}{10} |e_{k}|^{p}, |e_{0}| = 1.

Il metodo ha ordine di convergenza p visto che

lim_{n→∞} \frac{|e_{k+1}|}{|e_{k}|^{p}} = \frac{1}{10},

▶ p = 1: e_{0} = 1, e_{1} = 1/10, e_{2} = 1/100, e_{3} = 1/1000, . . .. ▶ p = 2: e_{0} = 1, e_{1} = 1/10, e_{2} = 1/1000, e_{3} = 1/10^{7}, ….

Alvise Sommariva  Equazioni e sistemi nonlineari.  4/ 44

Equazioni non lineari: esempio ordine di convergenza

Figura: Grafico che illustra l’errore di un metodo con convergenza p = 1 (in rosso a rombi) e p = 2 (in nero a cerchietti), per C = 1/10 ed e_{0} = 1.

Alvise Sommariva  Equazioni e sistemi nonlineari.  5/ 44

Metodo bisezione: definizione

Sia f : [a, b] → ℝ una funzione continua e supponiamo f (a) · f (b) < 0. Il metodo di bisezione genera una successione di intervalli (a_{k},b_{k}) con

▶ f (a_{k}) · f (b_{k}) < 0,

▶ [a_{k}, b_{k}] ⊂ [a_{k−1}, b_{k−1}],

▶ |b_{k} − a_{k}| = \frac{1}{2}|b_{k−1} − a_{k−1}|.

Fissate due tolleranze ϵ_{1}, ϵ_{2} si arresta l’algoritmo quando

|b_{k} − a_{k}| ≤ ϵ_{1} oppure |f ((a_{k} + b_{k})/2)| ≤ ϵ_{2}.

Metodo bisezione: esempio

Operativamente dati a ≤ b con f (a) · f (b) < 0 calcola c = (a + b)/2. Se f (a) · f (c) > 0 sostituisce c ad a, viceversa sostituisce c a b, fermandosi se le condizioni d’arresto sono verificate.

Studiamo f (x) = 0 con f (x) = sin (x) − x. Osserviamo che se x < 0 allora f (x) > 0 altrimenti f (x) ≤ 0.

1. a = −0.30000000, b = 0.40000000 → c = 0.05000000.

2. a = −0.30000000, b = 0.05000000 → c = −0.12500000 3. a = −0.12500000, b = 0.05000000 → c = −0.03750000 4. a = −0.03750000, b = 0.05000000 → c = 0.00625000

5….

Metodo bisezione: alcuni fatti

▶ Non ha una velocità di convergenza nel senso della definizione sopra definita. Esempio: se applico bisezione per risolvere sin (x) − x = 0, in [a, b] = [−3, 2] la successione |e_{n+1}/e_{n}| alterna valori 1.5 e 1/6 e quindi non converge con ordine 1.

▶ Usa solo il segno della funzione. ▶ Se f (a) · f (b) < 0 e f ∈ C([a, b]) allora converge (sempre!!) a

un x^{∗} tale che f (x^{∗}) = 0. ▶ Il test del residuo può essere non adatto a funzioni piatte o

con picchi intorno al punto cui converge. ▶ Fissata una tolleranza ϵ, e due punti iniziali a, b tali che

f (a) · f (b) < 0 (con f ∈ C([a, b])) per avere un’errore assoluto |x_{n} − x^{∗}| sulla soluzione x^{∗} inferiore ad ϵ necessitano al piú

n = ceil(\frac{log(b − a) − log ϵ}{log 2})

iterazioni del metodo.

Alvise Sommariva  Equazioni e sistemi nonlineari.  8/ 44

Metodo Newton: interpretazione analitica

Il metodo di Newton richiede che ▶ f sia derivabile con continuità su un intervallo [a, b] di ℝ; ▶ f^{(1)}(x) ≠ 0 in [a, b].

Se x_{k} ∈ [a, b] è l’ultima approssimazione del metodo, allora dalla formula di Taylor centrata in x_{k} abbiamo per ξ_{k} ∈ I(x^{∗}, x_{k})

0 = f (x^{∗}) = f (x_{k}) + f′(x_{k})(x^{∗} − x_{k}) + f′′(ξ_{k})(x^{∗} − x_{k})^{2}/2

Tralasciando i termini di ordine superiore

0 ≈ f (x_{k}) + f′(x_{k})(x^{∗} − x_{k})

e quindi se f^{(1)}(x_{k}) ≠ 0, dopo facili conti

x^{∗} ≈ x_{k+1} := x_{k} − f (x_{k})/f^{(1)}(x_{k}).

Metodo Newton: interpretazione geometrica

Se ben definito, il metodo di Newton genera una successione {x_{k}} definita da

x_{k+1} := x_{k} − f (x_{k})/f^{(1)}(x_{k}).

Sia r la retta tangente al grafico della funzione f nel punto (x_{k},f (x_{k})). Allora x_{k+1} è l’intersezione della retta r con l’asse delle ascisse, cioè la retta di equazione y = 0.

Metodo Newton: interpretazione geometrica

Figura: Grafico che illustra geometricamente le iterazioni del metodo Newton per il calcolo dello zero di f (x) = x^{2} − 2.

Alvise Sommariva  Equazioni e sistemi nonlineari.  11/ 44

Metodo Newton: convergenza

Nel caso del metodo di Newton la convergenza non è in generale garantita. Esistono degli esempi in cui il metodo produce una successione non convergente. Ciò nonostante esistono molti teoremi che illustrano quando si è certi che il metodo converga.

▶ Alcuni sono detti di convergenza locale e dicono che se x_{0} è in un intorno I sufficientemente piccolo della soluzione x^{∗} allora il metodo converge ad x_{∗}. Usualmente non sono costruttivi e non permettono di definire chiaramente I .

▶ Altri sono detti di convergenza globale e dicono che se x_{0} appartiene a un ben definito intorno I di x^{∗} allora il metodo converge.

Metodo Newton: un teorema di convergenza locale

Uno zero x^{∗} si dice semplice se f (x^{∗}) = 0 e f^{(1)}(x^{∗}) ≠ 0.

Teorema. Sia x^{∗} ∈ (a, b) uno zero semplice di f : [a, b] → ℝ. Si supponga inoltre f ∈ C^{2}([a, b]). Allora per x_{0} ∈ [a, b] sufficientemente vicino a x^{∗} le iterazioni del metodo di Newton

x_{k+1} = x_{k} − f (x_{k})/f^{(1)}(x_{k}), k = 0, 1, 2, …

sono ben definite e convergono quadraticamente a x^{∗}.

Metodo Newton: un teorema di convergenza locale (dim.)

Traccia della dimostrazione.

Sia I_{1} intorno di x^{∗} in cui f^{(1)} ha segno costante (permanenza del segno).

1. Dalla formula di Taylor centrata in x^{∗} e valutata in x_{k} abbiamo per ξ_{k} ∈ I(x^{∗}, ξ_{k}), ricordando la succ. del metodo di Newton

0 = f (x^{∗}) = f (x_{k}) + f^{(1)}(x_{k})(x^{∗} − x_{k}) + f^{(2)}(ξ_{k})(x^{∗} − x_{k})^{2}/2
= f^{(1)}(x_{k})(x^{∗} − x_{k+1}) + f^{(2)}(ξ_{k})(x^{∗} − x_{k})^{2}/2  (1)

Posto e_{k} = |x^{∗} − x_{k}| ricaviamo per M = max_{x,y∈I_{1}} |f^{(2)}(ξ_{k})|/2 · |f^{(1)}(x_{k})|

e_{k+1} = \frac{|f^{(2)}(ξ_{k})|}{2 · |f^{(1)}(x_{k})|}e_{k}^{2} ≤ Me_{k}^{2}

Metodo Newton: un teorema di convergenza locale (dim.)

Sia I_{2} = [x^{∗} − δ, x^{∗} + δ] ⊆ I_{1} con Mδ < 1. Se x_{0} ∈ I_{2} allora e_{0} ≤ δ e

e_{1} ≤ Me_{0}^{2} ≤ Mδ^{2} ≤ δ

da cui per induzione ogni x_{k} ∈ I_{2}. Inoltre

Me_{k} ≤ M^{2}e_{k−1}^{2} = (Me_{k−1})^{2} ≤ … ≤ (Me_{0})^{2^{k}}

e visto che Me_{0} ≤ Mδ < 1 abbiamo e_{k} → 0, cioè il metodo di Newton converge quadraticamente. Infatti si vede facilmente che essendo ξ_{k} ∈ I(x^{∗}, x_{k}) e x_{k} → x^{∗} necessariamente f ∈ C^{(2)}([a, b])

lim_{k} \frac{e_{k+1}}{e_{k}^{2}} = lim_{k}|f^{(2)}(ξ_{k})|/(2 · |f^{(1)}(x_{k})|)
= lim_{k}|f^{(2)}(x^{∗})|/(2 · |f^{(1)}(x^{∗})|) ≠ 0

Metodo Newton: un teorema di convergenza

Teorema. Sia f ∈ C^{2}([a, b]), con [a, b] intervallo chiuso e limitato. Se

1. f (a) · f (b) < 0;

2. f^{(1)}(x) > 0, per ogni x ∈ [a, b];

3. f^{(2)}(x) > 0 allora le iterate x_{k} fornite dal metodo di Newton sono strettamente decrescenti e convergono all’unica soluzione x^{∗} in [a, b], per x_{0} = b.

Metodo Newton: un teorema di convergenza globale

Teorema. Sia f ∈ C^{2}([a, b]), con [a, b] intervallo chiuso e limitato. Se

1. f (a) · f (b) < 0;

2. f^{(1)}(x) ≠ 0, per ogni x ∈ [a, b];

3. f^{(2)}(x) ≥ 0 o f^{(2)}(x) ≤ 0, per ogni x ∈ [a, b];

4. |f (a)/f^{(1)}(a)| < b − a e |f (b)/f^{(1)}(b)| < b − a, allora il metodo di Newton converge all’unica soluzione x^{∗} in [a, b], per ogni x_{0} ∈ [a, b].

Metodo Newton: un teorema di convergenza globale (dim.)

Traccia della dimostrazione. Supponiamo che

1. f^{(2)}(x) ≥ 0 se x ∈ (a, b).

2. f^{(1)}(x) > 0 se x ∈ (a, b).

Allora:

1. La soluzione x^{∗} è unica poichè f cambia segno ed è crescente.

2. L’ultimo punto garantisce che x_{1} ∈ (a, b).

3. Da hp. 2 e f (x_{k}) + f^{(1)}(x_{k})(x^{∗} − x_{k}) = f^{(1)}(x_{k})(x^{∗} − x_{k+1}) :

0 = f (x^{∗}) = f (x_{k}) + f^{(1)}(x_{k})(x^{∗} − x_{k}) + f^{(2)}(ξ_{k})(x^{∗} − x_{k})^{2}/2
= f^{(1)}(x_{k})(x^{∗} − x_{k+1}) + f^{(2)}(ξ_{k})(x^{∗} − x_{k})^{2}/2

e quindi 0 ≥ f^{(1)}(x_{k})(x^{∗} − x_{k+1}) cioè x^{∗} ≤ x_{k+1}.

4. La formula di Newton x_{k+1} = x_{k} − f (x_{k})/f^{(1)}(x_{k}) mostra che

se x_{k} > x^{∗} allora x_{k+1} < x_{k} (cioè {x_{k}}_{k=1,...} decrescente).

Alvise Sommariva  Equazioni e sistemi nonlineari.  18/ 44

Metodo Newton: un teorema di convergenza globale (dim.)

Quindi, comunque sia scelto x_{0} ∈ [a, b], abbiamo che per k ≥ 1 si ha x_{k} ∈ [x^{∗}, b] e {x_{k}}_{k=1,...} decrescente e quindi converge a un certo L ∈ [x^{∗}, b]. Dalla formula del metodo di Newton, per continuità

L = lim_{k}x_{k+1} = lim_{k}(x_{k} − f (x_{k})/f^{(1)}(x_{k}))
= lim_{k}x_{k} − lim_{k}f (x_{k})/f^{(1)}(x_{k})
= lim_{k}x_{k} − lim_{k}f (x_{k})/ lim_{k}f^{(1)}(x_{k}) = L − f (L)/f^{(1)}(L)

da cui F (L)/f^{(1)}(L) = 0 ed essendo f^{(1)}(x) > 0 per x ∈ [a, b] abbiamo F (L) = 0. Ma x^{∗} è l’unico zero di f per cui L = x^{∗} e x_{k} → x^{∗}.

Metodo Newton: alcuni fatti

1. Il metodo di Newton non è sempre convergente.

2. Se converge, non è detto che l’ordine di convergenza sia p = 2 (conv. quadratica).

3. Se uno zero x^{∗} di f non è semplice allora la convergenza non è quadratica.

Metodi di punto fisso

Consideriamo i problemi di punto fisso

x = ϕ(x)

in cui supponiamo che ϕ : Ω → Ω sia una funzione continua in Ω sia un intervallo compatto di ℝ. Notiamo subito che ogni problema di tipo f (x) = 0 si può ovviamente riscrivere in questa forma, posto ϕ(x) = f (x) + x. Nel metodo di punto fisso, si definisce la successione

x_{k+1} = ϕ(x_{k})

Teorema di punto fisso

Teorema. Si supponga che ϕ : Ω → Ω e ϕ′ siano funzione continue in Ω = [a, b] (intv. compatto). Inoltre sia λ = max_{x∈Ω} |ϕ′(x)| < 1. Allora:

1. Esiste uno ed un solo α tale che ϕ(α) = α in Ω;

2. Per ogni punto iniziale x_{0} ∈ Ω le iterate x_{n+1} = ϕ(x_{n}) convergono ad α;

3. Si ha

|α − x_{n}| ≤ \frac{λ^{n}}{1 − λ}|x_{0} − x_{1}|, n ≥ 0;

4. Infine

lim_{n→∞} \frac{α − x_{n+1}}{α − x_{n}} = ϕ′(α)

Metodi di punto fisso

Consideriamo un generico problema di punto fisso x = ϕ(x), e supponiamo

1. ϕ di classe C^{(2)},

2. ϕ(α) = α,

3. ϕ′(α) = 0, ϕ′′(α) ≠ 0,

4. x_{n} → α. Allora per ξ_{n} ∈ I(α, x_{n}),

x_{n+1} = ϕ(x_{n}) = ϕ(α) + ϕ′(α)(x_{n} − α) + ϕ′′(ξ_{n})(x_{n} − α)^{2}/2
= α + ϕ′′(ξ_{n})(x_{n} − α)^{2}/2

e quindi facilmente portando α a primo membro, per e_{n} = |x_{n} − α|

e_{n+1} = |ϕ′′(ξ_{n})|e_{n}^{2}/2 ⇒ lim_{n→∞}e_{n+1}/e_{n}^{2} = lim_{k}|ϕ′′(ξ_{n})|/2 = |ϕ′′(α)|/2

poichè ξ_{n} ∈ I(α, x_{n}) e x_{n} → α, e ϕ′′ ∈ C^{2}([a, b]) per hp.1. Quindi la convergenza è quadratica.

Alvise Sommariva  Equazioni e sistemi nonlineari.  23/ 44

Metodi di punto fisso e Newton

Il metodo di Newton è un particolare metodo di punto fisso. Infatti posto

ϕ(x) = x − f (x)/f′(x)

abbiamo che x_{k+1} = ϕ(x_{k}) = x_{k} − f (x_{k})/f′(x_{k}). Osserviamo che se f (α) = 0 e f′(α) = 0 allora

ϕ(α) = α − f (α)/f′(α) = α

cioè lo zero della f è pure sol. del problema di pto fisso e viceversa. Quindi se f ∈ C^{2}(Ω), f′(x) ≠ 0, ϕ : Ω → Ω otteniamo

ϕ′(x) = \frac{f (x)f′′(x)}{[f′(x)]^{2}}

e se f′(α) ≠ 0 (radice semplice) abbiamo ϕ′(α) = 0 da cui deduciamo nuovamente che la convergenza è quadratica.

Cenni alla risoluzione numerica di sistemi nonlineari

Problema. Data una funzione F : Ω ⊆ ℝ^{n} → ℝ^{n} si tratta di trovare gli x^{∗} ∈ Ω tali che F (x^{∗}) = 0. Esempio:

\{ \begin{matrix}x − \frac{1}{2} cos (y) = 0 \\ y − \frac{1}{2} sin (x) = 0 \end{matrix}  (2)

Soluzione (unica!):

x = 0.48640515466592, y = 0.23372550195872.

Generalizzazione punto fisso e Teorema di Banach

Teorema. Sia (X , d) uno spazio metrico completo ed M un sottoinsieme non vuoto e chiuso di X . Se ϕ : M → M è L contrattiva cioè d(ϕ(x), ϕ(y)) ≤ Ld(x, y) con L < 1 allora

1. l’equazione x = ϕ(x) ha una ed una sola soluzione x^{∗};

2. la successione x_{k+1} = ϕ(x_{k}) (detta di Banach o di punto fisso) converge alla soluzione x^{∗} per ogni scelta del punto iniziale x_{0};

3. per k = 0, 1, 2, … si ha

d(x_{k}, x^{∗}) ≤ L^{k}(1 − L)^{−1}d(x_{0}, x_{1}), a priori

d(x_{k+1}, x^{∗}) ≤ L(1 − L)^{−1}d(x_{k+1}, x_{k}), a posteriori

4. per k = 0, 1, 2, … si ha

d(x_{k+1}, x^{∗}) ≤ Ld(x_{k}, x^{∗})

Punto fisso per sistemi e Teorema di Banach

Teorema. Sia ℝ^{d} (dotato della norma euclidea ∥x∥ = \sqrt{∑_{k=1}^{d} x_{k}^{d}}) ed M un sottoinsieme non vuoto e chiuso di ℝ^{d}. Se ϕ : M → M è L contrattiva cioè ∥ϕ(x) − ϕ(y∥ ≤ L∥x − y∥ con L < 1 allora

1. l’equazione x = ϕ(x) ha una ed una sola soluzione x^{∗};

2. la successione x_{k+1} = ϕ(x_{k}) (detta di Banach o di punto fisso) converge alla soluzione x^{∗} per ogni scelta del punto iniziale x_{0};

3. per k = 0, 1, 2, … si ha

∥x_{k} − x^{∗}∥ ≤ L^{k}(1 − L)^{−1}∥x_{0} − x_{1}∥, a priori

∥x_{k+1} − x^{∗}∥ ≤ L(1 − L)^{−1}∥x_{k+1} − x_{k}∥, a posteriori

4. per k = 0, 1, 2, … si ha

∥x_{k+1} − x^{∗}∥ ≤ L∥x_{k} − x^{∗}∥

Esempio punto fisso per sistemi

Esempio: il sistema (2) è eqv. a risolvere il problema di punto fisso

\{ \begin{matrix}x = \frac{1}{2} cos (y) \\ y = \frac{1}{2} sin (x) \end{matrix}  (3)

che si scrive come v = ϕ(v) con v = (x, y) e

ϕ(v) = (\frac{1}{2} cos (y), \frac{1}{2} sin (x)) .

Si dimostra che ϕ è L contrazione con L = 1/2. Per applicare il teorema di punto fisso per sistemi basta porre

M = [−1, 1] × [−1, 1] ⊂ ℝ^{2} (quadrato unitario)

e osservare che

ϕ(M) ⊆ [−1/2, 1/2] × [−1/2, 1/2] ⊂ M

in quanto −1 ≤ cos (x) ≤ 1, −1 ≤ sin (x) ≤ 1.

Alvise Sommariva  Equazioni e sistemi nonlineari.  28/ 44

L contrattiva (per esperti)

Mostriamo che ϕ : ℝ^{2} → ℝ^{2} è contrattiva. Dalla Taylor in ℝ^{d}

ϕ(v) = ϕ(v_{0}) + (ϕ′(ξ)) · (v − v_{0}), ξ ∈ S

con S il segmento che congiunge v con v_{0}, ϕ′ la matrice Jacobiana. Nel nostro esempio quindi,

(ϕ′)_{j,k}(v) = ( \begin{matrix}0 & (−1/2) sin (v_{1}) \\ (+1/2) cos (v_{2}) & 0 \end{matrix} )

Essendo ∥A∥_{∞} = max_{j} ∑_{i} |a_{i,j}| si ha

∥(ϕ′)(v)∥_{∞} = max_{j} ∑_{i}| ((ϕ′)(v))_{i,j} |
= max (|(−1/2) sin (v_{1})|, |(+1/2) cos (v_{2})|) ≤ \frac{1}{2}

Quindi

∥ϕ(v) − ϕ(v_{0})∥_{∞} ≤ ∥(ϕ′(ξ))∥_{∞}∥(v − v_{0})∥_{∞}, ξ ∈ S da cui ϕ : ℝ^{2} → ℝ^{2} è una L-contrazione con L = 1/2.

Alvise Sommariva  Equazioni e sistemi nonlineari.  29/ 44

Esempio punto fisso per sistemi

Successione punto fisso (convergente alla soluzione di (3)): x_{k+1} = \frac{1}{2} cos (y_{k}), y_{k+1} = \frac{1}{2} sin (x_{k})

\begin{matrix}x_{k} & y_{k} \\ 0.00000000000000 & 0.00000000000000 \\ 0.50000000000000 & 0.00000000000000 \\ 0.50000000000000 & 0.23971276930210 \\ 0.48570310513698 & 0.23971276930210 \\ 0.48570310513698 & 0.23341513183103 \\ 0.48644107261515 & 0.23341513183103 \\ … & … \\ 0.48640515466592 & 0.23372550195872 \end{matrix}

Alvise Sommariva  Equazioni e sistemi nonlineari.  30/ 44

Metodo Newton per sistemi

Supponiamo di dover risolvere f (v) = 0 con f : Ω ⊆ ℝ^{d} → ℝ^{d}. La successione di Newton viene descritta come

\{ \begin{matrix}f′(v_{k})) · h_{k+1} = −f (v_{k}), \\ v_{k+1} := v_{k} + h_{k+1} \end{matrix}  (4)

dove f′ è la matrice Jacobiana di f cioè se v = (v_{1}, …, v_{d}) allora

(f′(v_{k})))_{i,j} = \frac{∂f_{i}(v)}{∂v_{j}}.

Si noti che ad ogni iterazione si deve risolvere il sistema lineare f′(v_{k})) · h_{k+1} = −f (v_{k}).

Alvise Sommariva  Equazioni e sistemi nonlineari.  31/ 44

Esempio Newton per sistemi

Consideriamo il precedente esempio, come f (x) = 0. Allora posto v = (x, y) si ha f (x, y) = (f_{1}(x, y), f_{2}(x, y)) con

f_{1}(x, y) = x − \frac{1}{2} cos (y)
f_{2}(x, y) = y − \frac{1}{2} sin (x).

Inoltre

f′(x, y) = ( \begin{matrix}1 & \frac{1}{2} sin (y) \\ −\frac{1}{2} cos (x) & 1 \end{matrix} )  (5)
Alvise Sommariva  Equazioni e sistemi nonlineari.  32/ 44

Esempio Newton per sistemi

Successione Newton:

\{ \begin{matrix}f′(v_{k})) · h_{k+1} = −f (v_{k}), \\ v_{k+1} := v_{k} + h_{k+1} \end{matrix}  (6)

\begin{matrix}x_{k} & y_{k} \\ 0.000000000000000 & 0.000000000000000 \\ 0.500000000000000 & 0.250000000000000 \\ 0.486463513355199 & 0.233773076987732 \\ 0.486405155145713 & 0.233725502568820 \\ 0.486405154665921 & 0.233725501958721 \end{matrix}

Alvise Sommariva  Equazioni e sistemi nonlineari.  33/ 44

Bisezione in Matlab/Octave: demobisezione2011.m

Salviamo in demobisezione2011.m il seguente file

% f=i n l i n e ( ’ ( x−1) . ˆ 3 ’ ) ; a=0; b =1.5;
f=inline ( ’ x .ˆ 2 −2 ’ ) ; a=1; b=2;
tolres =10ˆ (−15) ; tolintv =10ˆ (−15) ; maxit =10000;

[ aa , bb]=bisezione ( a , b , tolintv , tolres , maxit , f ) ;
format long e ; [ aa bb ]

▶ inline: definisce funzione da valutare; ▶ bisezione2011: funzione bisezione (da fare); ▶ a,b: intervallo iniziale bisezione (input); ▶ tolres,tolintv,maxit: tolleranze bisezione (input); ▶ f: passaggio di funzione come variabile (input); ▶ aa,bb: vettori intervalli analizzati da bisez. (output); Salviamo in demobisezione2011.m il seguente file

% f=i n l i n e ( ’ ( x−1) . ˆ 3 ’ ) ; a=0; b =1.5;
f=inline ( ’ x .ˆ 2 −2 ’ ) ; a=1; b=2;
tolres =10ˆ (−15) ; tolintv =10ˆ (−15) ; maxit =10000;

[ aa , bb]=bisezione ( a , b , tolintv , tolres , maxit , f ) ;
format long e ; [ aa bb ]

▶ inline: definisce funzione da valutare; ▶ bisezione2011: funzione bisezione (da fare); ▶ a,b: intervallo iniziale bisezione (input); ▶ tolres,tolintv,maxit: tolleranze bisezione (input); ▶ f: passaggio di funzione come variabile (input); ▶ aa,bb: vettori intervalli analizzati da bisez. (output); Salviamo in demobisezione2011.m il seguente file

% f=i n l i n e ( ’ ( x−1) . ˆ 3 ’ ) ; a=0; b =1.5;
f=inline ( ’ x .ˆ 2 −2 ’ ) ; a=1; b=2;
tolres =10ˆ (−15) ; tolintv =10ˆ (−15) ; maxit =10000;

[ aa , bb]=bisezione ( a , b , tolintv , tolres , maxit , f ) ;
format long e ; [ aa bb ]

▶ inline: definisce funzione da valutare; ▶ bisezione2011: funzione bisezione (da fare); ▶ a,b: intervallo iniziale bisezione (input); ▶ tolres,tolintv,maxit: tolleranze bisezione (input); ▶ f: passaggio di funzione come variabile (input); ▶ aa,bb: vettori intervalli analizzati da bisez. (output); Salviamo in demobisezione2011.m il seguente file

% f=i n l i n e ( ’ ( x−1) . ˆ 3 ’ ) ; a=0; b =1.5;
f=inline ( ’ x .ˆ 2 −2 ’ ) ; a=1; b=2;
tolres =10ˆ (−15) ; tolintv =10ˆ (−15) ; maxit =10000;

[ aa , bb]=bisezione ( a , b , tolintv , tolres , maxit , f ) ;
format long e ; [ aa bb ]

▶ inline: definisce funzione da valutare; ▶ bisezione2011: funzione bisezione (da fare); ▶ a,b: intervallo iniziale bisezione (input); ▶ tolres,tolintv,maxit: tolleranze bisezione (input); ▶ f: passaggio di funzione come variabile (input); ▶ aa,bb: vettori intervalli analizzati da bisez. (output); Salviamo in demobisezione2011.m il seguente file

% f=i n l i n e ( ’ ( x−1) . ˆ 3 ’ ) ; a=0; b =1.5;
f=inline ( ’ x .ˆ 2 −2 ’ ) ; a=1; b=2;
tolres =10ˆ (−15) ; tolintv =10ˆ (−15) ; maxit =10000;

[ aa , bb]=bisezione ( a , b , tolintv , tolres , maxit , f ) ;
format long e ; [ aa bb ]

▶ inline: definisce funzione da valutare; ▶ bisezione2011: funzione bisezione (da fare); ▶ a,b: intervallo iniziale bisezione (input); ▶ tolres,tolintv,maxit: tolleranze bisezione (input); ▶ f: passaggio di funzione come variabile (input); ▶ aa,bb: vettori intervalli analizzati da bisez. (output); Salviamo in demobisezione2011.m il seguente file

% f=i n l i n e ( ’ ( x−1) . ˆ 3 ’ ) ; a=0; b =1.5;
f=inline ( ’ x .ˆ 2 −2 ’ ) ; a=1; b=2;
tolres =10ˆ (−15) ; tolintv =10ˆ (−15) ; maxit =10000;

[ aa , bb]=bisezione ( a , b , tolintv , tolres , maxit , f ) ;
format long e ; [ aa bb ]

▶ inline: definisce funzione da valutare; ▶ bisezione2011: funzione bisezione (da fare); ▶ a,b: intervallo iniziale bisezione (input); ▶ tolres,tolintv,maxit: tolleranze bisezione (input); ▶ f: passaggio di funzione come variabile (input); ▶ aa,bb: vettori intervalli analizzati da bisez. (output);

Bisezione in Matlab/Octave: criterio arresto

Il criterio del residuo |f (c)| < tol non è spesso accettabile:

▶ funzione piatta: si pensi a risolvere l’equazione f (x) = 0 con f (x) = 10^{−50} · x; il residuo |f (1)| = 10^{−50} ma 1 è molto distante da x^{∗} = 0.

▶ funzione ripida: si pensi a risolvere l’equazione f (x) = 10^{50} · x = 0; il residuo |f (10^{−20})| = 10^{30} seppure 10^{−20} sia molto vicino a x^{∗} = 0.

Bisezione in Matlab/Octave: criterio arresto

Siano a < b e c = (a + b)/2. Diciamo residuo pesato |f (c) · w| con

w := (\frac{f (b) − f (a)}{b − a})^{−1} .

Si vede subito che w^{−1} è un rapporto incrementale.

▶ funzione piatta: si pensi a risolvere l’equazione f (x) = 0 con f (x) = 10^{−50} · x; w è grande (w = 10^{50}) e quindi |f (1)| = 10^{−50} ma |f (1)w| = 1;

▶ funzione ripida: si pensi a risolvere l’equazione f (x) = 0 con f (x) = 10^{50} · x; w è piccolo (w = 10^{−50}) e quindi |f (10^{−20})| = 10^{30} ma |f (10^{−20})w| = 10^{30} · 10^{−50} = 10^{−20}

Per f più generali, il test del residuo pesato |f (c)w| < tol prova ad adattare situazioni in cui il test del residuo |f (c)| < tol non sia affidabile.

Bisezione in Matlab/Octave: bisezione.m

f u n c t i o n [ aa , bb , ko]=bisezione ( a , b , tinv , tres , maxit , f )

i f b < a s=b ; b=a ; a=s ; end % A g g i u s t a e r r o r i u t e n t e .
ko=0; fa=f e v a l ( f , a ) ; fb=f e v a l ( f , b ) ; aa=[a ] ; bb=[b ] ;
i f fa∗fb > 0 ko=1; r e t u r n ; end
i f fa = 0 aa=[aa ; a ] ; bb=[bb ; a ] ; r e t u r n ; end % a s o l .
i f fb = 0 aa=[aa ; b ] ; bb=[bb ; b ] ; r e t u r n ; end % b s o l .

f o r index =1: maxit

c=(a+b ) / 2 ; fc=f e v a l ( f , c ) ;
w=(b−a ) /( fb−fa ) ; wres=abs ( fc∗w ) ; % Residuo p e s a t o .
semilung=(b−a ) / 2 ; % Semilung . i n t e r v a l l o .

     i f ( wres<tres ) | ( semilung<tinv ) | ( fc==0)

aa=[aa ; c ] ; bb=[bb ; c ] ; r e t u r n ;% OK e x i t .

end

i f s i g n ( fc ) = s i g n ( fa ) % Aggiornamento .

         aa=[aa ; c ] ; bb=[bb ; b ] ; a=c ; fa=fc ;

e l s e

         aa=[aa ; a ] ; bb=[bb ; c ] ; b=c ; fb=fc ;

end end

k0=1;
Alvise Sommariva  Equazioni e sistemi nonlineari.  37/ 44

Bisezione in Matlab/Octave: esempio demobisezione

Usiamo demobisezione per il calcolo della soluzione positiva di x^{2} − 2 = 0, cioè \sqrt{2} = 1.41421356237309 …

>> demobisezione2011

ans =
1.000000000000000 e+00  2.000000000000000 e+00
1.000000000000000 e+00  1.500000000000000 e+00
1.250000000000000 e+00  1.500000000000000 e+00
1.375000000000000 e+00  1.500000000000000 e+00

…

1.414213562373092 e+00  1.414213562373099 e+00
1.414213562373096 e+00  1.414213562373096 e+00
>> l e n g t h ( aa ) % ITERAZIONI COMPIUTE.
ans = 49
>> s q r t ( 2 )
ans = 1.414213562373095 e+00

>>

Newton in Matlab/Octave: demonewton2011.m

Il file demonewton2011.m usa Newton per calcolo \sqrt{2}.

% f=i n l i n e ( ’ ( x−1) . ˆ 3 ’ ) ; x0 =1.5;
f=inline ( ’ x .ˆ 2 −2 ’ ) ; df=inline ( ’ 2∗x ’ ) ; x0=1;
tolstep =10ˆ (−15) ; maxit =1000;

[ x , ko]=newton2011 ( x0 , tolstep , maxit , f , df ) ;

format long e ; x

▶ Il comando inline permette di definire le funzioni

f (x) = x^{2} − 2 e f^{(1)}(x) = 2x, in forma vettoriale (si noti il .). ▶ Il valore da cui parte il metodo di Newton è x0 (vicino a x^{∗}). ▶ tolstep definisce la tolleranza del criterio di arresto dello step

|x_{k+1} − x_{k}| ≤ tolstep

uscendo al massimo dopo maxit iterazioni. ▶ newton2011 è la routine che esegue il metodo di Newton. In

output vettore iterazioni x, flag ko (0: ok, 1: ko).

Alvise Sommariva  Equazioni e sistemi nonlineari.  39/ 44

Newton in Matlab/Octave: newton2011.m

f u n c t i o n [ x , ko]=newton2011 ( x0 , tol , kmax , f , df )
x_old=x0 ; ko=0; x=[x0 ] ;

step=realmax ; k=0;
w h i l e ( abs ( step ) > tol ) & ( k < kmax )
k=k+1; fx=f e v a l ( f , x_old ) ;
i f fx = 0 r e t u r n ; end

     dfx=f e v a l ( df , x_old ) ;

i f dfx = 0 ko=1; r e t u r n ; end

     step=−fx/dfx ; x_new=x_old+step ;
     x=[x ; x_new ] ; x_old=x_new ;

end

ko=1;

Newton in Matlab/Octave: esempio radice quadrata.

>> % NEWTON. CALCOLO RADICE QUADRATA.

>> demonewton2011
x =

1.000000000000000 e+00
1.500000000000000 e+00
1.416666666666667 e+00
1.414215686274510 e+00
1.414213562374690 e+00
1.414213562373095 e+00
1.414213562373095 e+00
>> d i f f ( x )
ans =
5.000000000000000 e−01
−8.333333333333326e−02
−2.450980392156854e−03
−2.123899820016817e−06
−1.594724352571575e−12
−2.220446049250313e−16

>>

Applicando il metodo di bisezione nell’intervallo [1, 2] sarebbero servite 49 iterazioni.

Alvise Sommariva  Equazioni e sistemi nonlineari.  41/ 44

Newton in Matlab/Octave: esempio (x − 1)^{3}.

Risolviamo con il metodo di Newton l’equazione (x − 1)^{3} = 0 la cui soluzione è 1 con molteplicità 3 (zero multiplo!). Ci si aspetta convergenza solo lineare. Quale punto iniziale poniamo x_{0} = 1.5. Il metodo offre la soluzione dopo 83 iterazioni.

>> demonewton2011
x =

1.500000000000000 e+00
1.333333333333333 e+00
1.222222222222222 e+00
1.148148148148148 e+00
1.098765432098765 e+00

…

1.000000000000003 e+00
1.000000000000002 e+00
>> l e n g t h ( x )
ans =

83 >>

Bisezione in Matlab/Octave: esempio (x − 1)^{3}.

Sullo stesso esempio studiamo bisezione, con a = 0 e b = 1.5. Il metodo offre la soluzione dopo 46 iterazioni.

>> demobisezione2011

ans =
0.000000000000000 e+00  1.500000000000000 e+00
7.500000000000000 e−01  1.500000000000000 e+00
7.500000000000000 e−01  1.125000000000000 e+00
9.375000000000000 e−01  1.125000000000000 e+00
9.375000000000000 e−01  1.031250000000000 e+00
9.843750000000000 e−01  1.031250000000000 e+00
9.843750000000000 e−01  1.007812500000000 e+00

…

9.999999999999432 e−01  1.000000000000028 e+00
9.999999999999858 e−01  1.000000000000028 e+00
>> l e n g t h ( a )
ans =

46 >>

Alvise Sommariva  Equazioni e sistemi nonlineari.  43/ 44

Riferimenti bibliografici

▶ K. Atkinson, W. Han, Elementary Numerical Analysis, Wiley,

(2004). ▶ A. Quarteroni, F. Saleri, Introduzione al Calcolo Scientifico,

Springer, (2002).