3.1 Il metodo di Bisezione

---

Il metodo di bisezione è un metodo numerico che consiste nel dimezzare di volta in volta l'Intervallo di separazione e usare il punto medio di tale intervallo come approssimazione della soluzione $\xi$.
Consente di costruire una successione {$x_k$} di approssimazioni di $\xi$, ottenute assumendo $x_k$ come punto medio di un intervallo $[a_{k-1},b_{k-1}]$ con $k\in {\mathbb{N}}^{+}$.

☑️ Ipotesi

Algoritmo

Individuato unIntervallo di separazione $[a_0,b_0]$ ne calcolo il punto medio che chiamerò $x_1$:
$x_1=\frac{b_0+a_0}{2}$
Ora controllo in quale tra gli intervalli $[a_0,x_1]$ e $[x_1,b_0]$ appartiene la soluzione $\xi$, verificando che il prodotto $f(a_{k-1})\cdot f(x_k)$ oppure il prodotto $f(b_{k-1})\cdot f(x_k)$ sia negativo.

Una volta individuato l'intervallo, ripeto il procedimento usando quello come nuovo intervallo di separazione. Supponiamo ad esempio che la soluzione cada nell'intervallo $[a_0,x_1]$, questo sarà proprio il nuovo intervallo di separazione $[a_1,b_1]=[a_0,x_1]$.

Chiamo $x_2$ il punto medio tra $a_1$ e $b_1$:
$x_2=\frac{a_1+b_1}{2}$
Continuando a ripetere questo processo, ammesso che il metodo converga, $x_k$ sarà man mano sempre più vicino alla soluzione $\xi$.

• Input: $a=a_0,b=b_0$ per $k=1,2,...$
  • $x_k=\frac{b_{k-1}-a_{k-1}}{2}$
• se:
  • $f(a_{k-1})\cdot f(x_k)<0\Rightarrow [a_k,b_k]=[a_{k-1},x_k]$
  • $f(b_{k-1})\cdot f(x_k)<0\Rightarrow [a_k,b_k]=[x_k,b_{k-1}]$

Errori

Errore di troncamento

Nell'applicazione del metodo di bisezione, compiamo un Errore di troncamento dovuto quindi al numero di iterazioni che compiamo

Convergenza

Calcolo l'errore $e_k$ tenendo presente che esso dovrà per forza essere minore della metà dell'ampiezza dell'intervallo.
A sua volta l'ampiezza dell'intervallo $[a_{k-1},b_{k-1}]$ sarà la metà dell'intervallo della iterazione predente $[a_{k-2},b_{k-2}]$ e così via

$
\begin{align}
|e_{k}| = |\xi - x_{k}| \le & \frac{b_{k-1} - a_{k-1}}{2} \

&\le \frac{1}{2} \frac{b_{k-2}-a_{k-2}}{2} \

\vdots \

&\le \frac{1}{2^{k}} (b_{0}- a_{0})

\end{align}
$$
\Rightarrow |e_{k}| = |\xi - x_{k}| \le \frac{1}{2^{k}} (b_{0} - a_{0})
$\text{<a class="internal-link" target="" data-note-icon="" href="/universita/triennale/2-anno/2-semestre/analisi-numerica/appunti/01-introduzione-ai-metodi-numerici/\#convergenza">Converge</a>?}$
\lim_{k \to \infty} |e_{k}| = \lim_{k \to \infty} \left|\frac{1}{2^{k}}(b_{0}-a_{0}) \right| = 0
$

Ordine di convergenza

Stima iterazioni necessarie

Per ottenere l'approssimazione della soluzione con una tolleranza scelta $\epsilon$ sono necessarie almeno un certo numero $k$ di iterazioni. Questo numero può essere stimato a partire dalla relazione:

$|e_k|=\frac{1}{2^k}(b-a)\le \epsilon$

Risolvendo per $k$ otteniamo infatti:

$
\begin

2^{k} &\ge \frac{b-a}{\varepsilon} \

log_{2}\left( 2^{k} \right) &\ge log_{2}\left( \frac{b-a}{\varepsilon} \right) \

k &\ge log_{2}\left( \frac{b-a}{\varepsilon} \right)

\end{align}
$
Tenendo presente che $k$ deve per definizione essere un numero intero positivo.

Criterio di arresto

Criterio di arresto a posteriori

Posso interrompere l'algoritmo quando l'errore al passo k è minore di una tolleranza scelta $\epsilon$.
$|e_k|\le \epsilon$

Implementazione in Matlab

function [soluzioneApprox,erroreStimato, iterazioniNecessarie] = MetodoBisezioneNonLineari(estremoInferiore, estremoSuperiore, funzione, tolleranza, numeroMaxIterazioni)

% Questa funzione serve ad applicare il metodo di bisezione a una funzione
% assegnata per ricavarne la soluzione approssimata.
% Input:
% estermoInferiore, estremoSuperiore: Estremi dell'intervallo di separazione
% funzione: la funzione, definita come anonymous function, a cui
% applicare il metodo
% tolleranza: La tolleranza con la quale si vuole trovare la soluzione
% numeroMaxIterazioni: Il numero max di iterazioni che si vogliono
% provare prima di interrompere il calcolo. (Se non specificato è
% inizializzato a 100

% ------------

% Outpu:
% soluzioneApprox: la soluzione approssimata
% erroreStimato: la media degli errori calcolati come b-a e come f(xn)
% iterazioniNecessarie: Il numero di iterazioni necessarie alla stima
% dell'errore
% Cerco la radice di una funzione f nell'intervallo [a,b] con il metodo di
% bisezione
% Suppongo di aver già individuato gli estremi degli intervalli di
% separazione
% Sfrutteremo un criterio di arresto basato sulla tolleranza:
% f = Funzione della quale vogliamo cercare la radice
% a = Estremo inferiore in cui è stata isolata la radice
% b = Estremo superiore in cui è stata isolata la radice
% toll = Tolleranza (\tollilon)
% numeroMaxIterazioni = Massimo numero di iterazioni

f = funzione;

% So che nell'intervallo definito sotto sicuremante ci sarà una soluzione
a = estremoInferiore;
b = estremoSuperiore;

% Definisco la tolleranza
toll = tolleranza;
maxIter = 100;
maxIter = numeroMaxIterazioni;

% Vogliamo calcolare:
% xn: Approx della radice
% err1 = |x(n) - x(n-1)|
% err2 = f(xn)
% iter = Numero di iterazioni

format long

% Inizializzo i parametri necessari
iter = 0; err1 = b-a; err2 = toll+1; x0 = a;

% Verifico che la radice esista
if f(a)*f(b) > 0
	error('Non ci sono soluzioni nellintervallo [a,b]') % Con questo comando blocco l'esecuzione dello script
elseif f(a)*f(b) == 0
	if f(a) == 0
		xn = a;
	elseif f(b) == 0
		xn = b;
	end
	return % Return perché non voglio andare avanti se queste condizioni del secondo if sono rispettate
end

% Calcolo successione
while (err1 > toll) && iter < maxIter
	xn = (a+b)/2;
	if f(a)*f(xn) < 0
		b = xn;
	elseif f(xn)*f(b) < 0
		a = xn;
end

iter = iter + 1;
err1 = abs(xn-x0);
err2 = abs(f(xn));
x0 = xn;
end

soluzioneApprox = xn;
erroreStimato = abs((err1 + err2)/2);
iterazioniNecessarie = iter;

end

3.2 Il metodo delle Tangenti

---

Il metodo delle tangenti, detto anche metodo di Newton sfrutta il Polinomio di Taylor per approssimare il valore delle funzioni. In questo caso è necessario scegliere un'approssimazione iniziale $x_0$.

Si tratta di fatto di un metodo di linearizzazione. Infatti, come vedremo, interromperemo il Polinomio di Taylor al primo ordine, confondendo quindi il valore della funzione con il valore della tangente.

☑️ Ipotesi

Per l'applicabilità del metodo di newton, la funzione deve rispettare le seguenti ipotesi:

Algoritmo

Scelgo un'approssimazione iniziale $x_0$.

Approssimo la funzione $f(x)$ con un polinomio di Taylor al secondo ordine:
$f(x)=f(x_0)+f^{\prime }(x_0)(x-x_0)+\frac{1}{2}{f}^{″}(\tau )(x-x_0{)}^2\text{ con }\tau \in [x_0,x]$
Se cerco la radice di $f(x)$ impongo la funzione approssimata uguale a 0:
$f(x)=f(x_0)+f^{\prime }(x_0)(x-x_0)\stackrel{!}{=}0$
da cui
$x_1=x_0-\frac{f(x_0)}{f^{\prime }(x_0)}$
motivo per cui dobbiamo imporre anche $f^{\prime }(x_k)\ne 0$ tra le #☑️ Ipotesi.

Dopodiché approssimo nuovamente la funzione con Taylor usando stavolta $x_1$ al posto di $x_0$:
$f(x)=f(x_1)+f^{\prime }(x_1)(x-x_1)\stackrel{!}{=}0$
da cui
$x_2=x_1-\frac{f(x_1)}{f^{\prime }(x_1)}$
Continuo ad iterare questo procedimento. Avrò così, in definitiva:
$x_k=x_{k-1}-\frac{f(x_{k-1})}{f^{\prime }(x_{k-1})}$
Il metodo convergerà per alcuni valori di $x_0\in J\subset I$ (Ossia scegliendo $x_0$ abbastanza vicino a $\xi$).

Convergenza

---

Estremo di Fourier

Se $f^{″}(x)\ne 0\mathrm{\forall }x\in I$ si può dimostrare che si può scegliere uno dei due estremi di $I$ come approssimazione iniziale $x_0$.

Supponiamo di avere una funzione monotona, senza massimi o minimi locali, a concavità fissa.
L'estremo di Fourier di un intervallo è quello dei due tale che il prodotto $f(a)\cdot f^{″}(a)>0$ dove $a$ è l'estremo.

Ordine di convergenza

Per calcolare l'ordine di convergenza devo calcolare il limite:
$\underset{k\to \mathrm{\infty }}{lim}\frac{|e_{k+1}|}{|e_k{|}^p}=\text{C}>0p\ge 1$
Calcolo $|e_{k+1}|$:

$|e_{k+1}|=\xi -x_{k+1}=(\xi -\frac{f(\xi )}{f^{\prime }(\xi )})-x_{k+1}$
essendo $\frac{f(\xi )}{f^{\prime }(\xi )}=0$ perché $f(\xi )=0$.

$
\begin{align}
|e_{k+1}| = \xi - x_{k+1} &= \left( \xi - \frac{f(\xi)}{f'(\xi)} \right) - x_{k+1} = \

&= \left( \xi - \frac{f(\xi)}{f'(\xi)} \right) - \left( x_{k} - \frac{f(x_{k})}{f'(x_{k})} \right) = \

&= (\xi - x_{k}) - \left( \frac{f(\xi)}{f'(\xi)} - \frac{f(x_{k})}{f'(x_{k})} \right) = \

\end{align}
$
Sviluppando $f(x_k)$ e $f^{\prime }(x_k)$ in serie di Taylor e ricordando $\xi -x_k=e_k$ ottengo:

$
\begin

|e_{k+1}| &= \left|
e_{k} - \left( \frac{f(\xi)}{f'(\xi)} - \frac{f(\xi) + f'(\xi)(-e_{k}) + \frac{1}{2} f''(\xi)e_{k}^{2} }{f'(\xi)} \right)
\right| \

&= \frac{1}{2} \left| \frac{f''(\xi)}{f'(\xi)} \right| |e_{k}^{2}|
\end{align}
$

Si ha quindi che

$\frac{|e_{k+1}|}{|e_k{|}^p}\stackrel{p=2}{\approx }\frac{1}{2}\left|\frac{f^{″}(\xi )}{f^{\prime }(\xi )}\right|$
quindi

$\underset{k\to \mathrm{\infty }}{lim}\frac{|e_{k+1}|}{|e_k{|}^2}=\frac{1}{2}\left|\frac{f^{″}(\xi )}{f^{\prime }(\xi )}\right|$
Allora l'ordine di convergenza del metodo delle tangenti è $p=2$.

Implementazione in Matlab

function xk = tangenti(x0,f, df, N_max)

% Input:
%   x0: approx iniziale
%   f, df: funzione e derivata
%   N_max: numero di iterazioni
% 
% Output:
%   xk: approx dello zero

for i = 1:N_max
    xk = x0 - f(x0) / df(x0);
    errk = abs(xk - x0);
    x0 = xk;
    fprintf('k = %3i x_k = %14.12f err=%9.3e f=%18.12f f''=%18.12f \n', i, xk, errk, f(x0), df(x0));
end

end

Metodo di Newton in ${\mathbb{R}}^n$

Al link 03.4 Metodo di Newton in n dimensioni

Se in $\mathbb{R}$ possiamo riscrivere una funzione sviluppandola in un Polinomio di Taylor, lo stesso possiamo fare con ${\mathbb{R}}^n$.

Sia $F(X):{\mathbb{R}}^n\to {\mathbb{R}}^n$, posso scrivere
$F(X)=F(X^{(k)})+J_F(X^{(k)})(X-X^{(k)})+...$

dove
$
\boldsymbol J_F(\boldsymbol X) =
\begin

\frac{\partial f_1}{\partial x_1} & \frac{\partial f_1}{\partial x_2} & \dots &\frac{\partial f_1}{\partial x_n} \

\frac{\partial f_2}{\partial x_1} & \frac{\partial f_2}{\partial x_2} & \dots &\frac{\partial f_2}{\partial x_n} \

\vdots & \vdots & \ddots & \vdots \

\frac{\partial f_n}{\partial x_1} & \frac{\partial f_n}{\partial x_2} & \dots &\frac{\partial f_n}{\partial x_n} \

\end{bmatrix}
$

3.3 Metodo delle Secanti

---

Il metodo delle Secanti consiste nell'approssimare localmente la funzione con la retta secante passante per i due estremi dell'intervallo $I$.

☑️ Ipotesi

Algoritmo

Scelgo 2 approssimazioni iniziali: $x_0,x_1$
Poi, per $k=2,3,...$ calcolo la soluzione approssimata come radice della retta secante passante per $x_0$ e $x_1$:
$x_k=x_{k-1}-f(x_{k-1})\frac{x_{k-1}-x_{k-2}}{f(x_{k-1})-f(x_{k-2})}$

Convergenza

Ordine di convergenza

L'ordine di convergenza del metodo delle secanti è
$p=\frac{1+\sqrt{5}}{2}\approx 1.62$
Essendo maggiore di 1 ma minore di 2 si dice che ha ordine di convergenza sopralineare.

Implementazione in Matlab

function [xk,errk, iter] = secanti_tol(x0,x1, f, N_max, tol)
% Input: (x0, x1, f, N_max, tol)
% x0, x1: Approx iniziali
%   f: Funzione di cui trovare gli zeri
%   N_max: Numero di iterazioni
%   tol: Tolleranza richiesta
%
% Output:
%   xk: approx dello zero
%   errk: Errore alla k-esima iterata
%   iter: Numero di iterazioni effettuate

errk  = tol + 1;
iter = 0;

while (iter < N_max) && (errk > tol)
    xk = x1 - f(x1)* (x1-x0) / (f(x1) - f(x0));
    errk = abs(xk-x1);
    x0 = x1;
    x1 = xk;
    iter = iter + 1;
end

fprintf('xk = %f', xk)
end

Metodo di Newton in ${\mathbb{R}}^n$

Se in $\mathbb{R}$ possiamo riscrivere una funzione sviluppandola in un Polinomio di Taylor, lo stesso possiamo fare con ${\mathbb{R}}^n$.

Sia $F(X):{\mathbb{R}}^n\to {\mathbb{R}}^n$, posso scrivere
$F(X)=F(X^{(k)})+J_F(X^{(k)})(X-X^{(k)})+...$

dove al posto della derivata, uso la matrice Jacobiana:
${\mathit{J}}_F(\mathit{X})=\left[\begin{array}{cccc}\frac{\partial f_1}{\partial x_1}& \frac{\partial f_1}{\partial x_2}& \dots & \frac{\partial f_1}{\partial x_n}\\ \frac{\partial f_2}{\partial x_1}& \frac{\partial f_2}{\partial x_2}& \dots & \frac{\partial f_2}{\partial x_n}\\ ⋮& ⋮& \ddots & ⋮\\ \frac{\partial f_n}{\partial x_1}& \frac{\partial f_n}{\partial x_2}& \dots & \frac{\partial f_n}{\partial x_n}\end{array}\right]$

Come metodo iterativo diventa l'equazione. La nuova approssimazione si trova risolvendo questo sistema

$\mathit{F}\left({\mathit{X}}^{(k)}\right)+{\mathit{J}}_F({\mathit{X}}^{(k)})\underset{\stackrel{\mathrm{△}}{=}{\mathit{Y}}^{(k+1)}}{\underbrace{({\mathit{X}}^{(k+1)}-{\mathit{X}}^{(k)})}}=0$
Ho ottenuto un problema lineare della forma:

${\mathit{J}}_F({\mathit{X}}^{(k)}){\mathit{Y}}^{(k+1)}=-\mathit{F}({\mathit{X}}^{(k)})$

Posso quindi scrivere la soluzione approssimata ${\mathit{X}}^{(k+1)}$ come:
${\mathit{X}}^{(k+1)}={\mathit{X}}^{(k)}+{\mathit{Y}}^{(k+1)}$
Se ${\mathit{J}}_F({\mathit{X}}^{(k)})$ è invertibile, allora
${\mathit{X}}^{(k+1)}=\underset{={\mathbf{\Phi }}_N\left({\mathit{X}}^{(k)}\right)}{\underbrace{{\mathit{X}}^{(k)}-{\left(J_F\left({\mathit{X}}^{(k)}\right)\right)}^{-1}\mathit{F}\left(X^{(k)}\right)}}$
Mi fermo quando
$\mathrm{dist}({\mathit{X}}^{(k+1)},{\mathit{X}}^{(k)})\le \epsilon$