Francesco Mugelli

Il Calcolo Numerico: un paio di applicazioni

Le formule di quadratura: come si calcolano gli integrali
Soluzione di equazioni non lineari
Il metodo di bisezione
Il metodo delle corde
Bibliografia


1. Le formule di quadratura: come si calcolano gli integrali.

Vogliamo risolvere il problema del calcolo approssimato di

supponendo di conoscere la funzione f(x) e che questa sia continua. Supponiamo poi di essere in grado di calcolare i valori di f(x) con l'accuratezza desiderata. Il primo metodo che viene in mente è suddividere l'intervallo di integrazione in un certo numero di intervallini,

a=x0< x1< ... < xn-1< xn=b,

ed approssimare il valore dell'area del sottografico della funzione con quella dell'unione dei rettangoli aventi per base l'intervallo [xi-1, xi] e altezza f(xi).

In pratica, abbiamo approssimato il valore dell'integrale con

Avremmo potuto scegliere anche f(xi-1) come altezza dei rettangoli. La formula ottenuta sarebbe stata simile ma avrebbe fornito un valore numerico differente dal precedente. Quale dei due valori è migliore? E soprattutto, cosa significa in questo caso la parola "migliore"?

Iniziamo l'analisi cercando una stima dell'errore che si commette sostituendo alla funzione f(x) una funzione costante a tratti (è un altro modo di descrivere quello che abbiamo fatto calcolando l'area dei rettangoli al posto dell'integrale di f(x)).
Qualunque risultato numerico, non solo il risultato del calcolo di un integrale ma anche, per esempio, la soluzione di un sistema lineare, o la soluzione approssimata di un equazione numerica o differenziale, non è utile se presentato da solo, senza altre informazioni. Ha significato solo se accompagnato da una stima dell'errore commesso o comunque da qualcosa che permetta di valutarne l'attendibilità.
Nel caso in cui si cerchino, le soluzioni di un equazione, se si usa ad esempio il metodo di bisezione, il metodo stesso fornisce una stima dell'errore. Se invece si usa il metodo delle corde, non siamo in grado di fornire una stima dell'errore ma, se x0 è il valore ottenuto per la soluzione, |f(x0)| può essere utile a valutare il risultato ottenuto.

Cerchiamo una stima dell'errore per la nostra formula di quadratura. Sia A il valore fornito dalla formula. Vogliamo trovare due quantità A1 e A2 tali che

A1< I[f]< A2, A1< A< A2.

In questo caso potremo affermare che A è un valore che approssima I[f] e che l'errore commesso non supera A2-A1. Il numero di cifre "esatte" di A è dato dal numero di ordini di grandezza di differenza tra A2-A1 ed A stesso (si osservi che se il valore esatto è 2.000 e come valore approssimato otteniamo 1.999, si dice che il valore ottenuto ha 3 cifre esatte anche se in effetti sono tutte differenti.

Le cifre esatte sono soltanto un modo per quantificare l'errore relativo ).
Nel nostro caso, se come altezza dei rettangolini scegliamo il min{f(xi-1), f(xi)}, otteniamo una approssimazione per difetto del valore di I[f]:

Analogamente, se invece del minimo scegliamo il massimo tra i 2 valori, otteniamo una approssimazione per eccesso:

Possiamo quindi scrivere una prima stima dell'errore:

Se supponiamo che la funzione sia derivabile con derivata continua, allora, per il teorema di Lagrange si ha:

f(xi)-f(xi-1)= f'(t)(xi-xi-1)

dove xi-1< t < xi. Quindi, se

M=max{|f'(t)| : a< t < b},
si ha:
|f(xi)-f(xi-1)|< M |xi-xi-1|.

Utilizziamo la maggiorazione appena ottenuta nella formula scritta sopra per l'errore:

Una scelta piuttosto comune è prendere gli xi equidistanti; in questo caso,

xi - xi-1 = (b-a)/n.

La stima che ricava è piuttosto semplice da interpretare:

Errore < M (b-a)2/n

Si vede chiaramente in che modo l'errore dipenda dal numero delle suddivisioni. Inoltre si può osservare che le funzioni più "difficili" da integrare sono quelle fortemente oscillanti (cioè con M grande), cosa che era da aspettarsi, data la maggiore difficoltà di approssimarle fedelmente con i rettangoli.

L'applicazione che segue utilizza il metodo dei rettangoli, modificato in vari modi, ma utilizza anche alcuni metodi più raffinati, come ad esempio il metodo dei trapezi. Anche per i trapezi è possibile, anche se è un po' più complicato, ottenere stime dell'errore:

Erroretrapezi = M2 (b-a)3/n3
dove
M2=max{|f"(t)| : a< t < b}.

Con f" si indica la derivata seconda di f.
Il metodo dei trapezi è più accurato del precedente, come si può osservare dal fatto che l'errore dipende dal cubo del numero di suddivisioni. Questo permette di raggiungere ottime approssimazioni con un numero relativamente piccolo di suddivisioni.
I metodi proposti fanno tutti parte di una stessa categoria, quella delle formule di Newton-Côtes. Viene calcolato esattamente, cioè senza errori dovuti al metodo, l'integrale del polinomio interpolante la funzione negli xi scelti (interpolante significa che passa per i punti (xi, f(xi)) ). A seconda del grado scelto per il polinomio si ottengono i vari metodi; per esempio se il grado è 0 si ha il metodo dei rettangoli, se è 1 si ha il metodo dei trapezi.

Questa applet calcola l'area sotto il grafico di una funzione data. La funzione di default è f(t) = 1/t. È possibile cambiare la funzione clickando sul pulsante "Function" e gli estremi di integrazione sia spostandoli con il mouse che scrivendo il valore nelle finestre. Si può anche scegliere il metodo di integrazione. È possibile soprattutto confrontare, anche visivamente, la differenza di accuratezza tra il metodo dei rettangoli e quello dei trapezi a parità del numero n di suddivisioni.

Ringraziamo il prof. Paul Caprioli caprioli@middlebury.edu, visiting professor al Middlebury College, che ci ha gentilmente concesso la possibilità di utilizzare l'applet, di cui è unico autore, in questo testo.


2. Soluzione di equazioni non lineari

Quando si vuole risolvere numericamente una equazione, di qualsiasi tipo essa sia, è conveniente prima di tutto riportarla nella forma

f(x) = 0

dove f(x) è una opportuna funzione. Prima di esaminare degli algoritmi per determinare le soluzioni di una equazione è necessario chiarire che cosa significa "risolvere" l'equazione. Infatti solo raramente è possibile determinarne le soluzioni esatte e questo può accadere per vari motivi.
Facciamo un paio di esempi:

A seconda del contesto, cioè a seconda del metodo usato, dire che un valore x* è una soluzione approssimata di una equazione f(x) = 0 può voler dire sia che |f(x*)| è "piccolo" a sufficienza sia che x* è "vicino" alla soluzione esatta dell'equazione.


Il metodo di bisezione

Consideriamo l'equazione
x3 - x - 1 = 0.

Osserviamo che

f(1) = -1e f(2) = 5,

allora, poiché la funzione f(x) è continua, essa assume tutti i valori compresi tra -1 e 5 (teorema di Weierstrass). In particolare deve esistere almeno un valore x* compreso tra 1 e 2 tale che f(x*) = 0. Studiando il segno della derivata f'(x) si potrebbe anche stabilire a priori che c'è esattamente un valore che annulla f(x) nell'intervallo [1, 2].
Potremo prendere il punto medio del segmento come prima approssimazione della soluzione:

x* = 1.5 con errore assoluto < 0.5

Per cercare di individuare la soluzione valutiamo f(x) nel punto medio dell'intervallo; si ha:

f(1.5) = 0.875 > 0 > -1 = f(1).

Quindi ora sappiamo che la soluzione x* è nell'intervallo [1, 1.5]. In altri termini:

x* = 1.25 con errore assoluto < 0.25.

Ripetendo lo stesso ragionamento, calcoliamo f(1.25):

f(1.25) = -0.296... < 0 < 0.875 = f(1.5),
quindi
x*=1.375 con errore assoluto < 0.125.

Proseguendo, si individua una sequenza di intervalli, di ampiezza sempre più piccola: ogni volta che si calcola la f in un nuovo punto l'errore dimezza. Questo procedimento prende il nome di metodo di bisezione:

Data una funzione f(x) continua nell'intervallo [a0, b0] e tale che f(a0)f(b0) < 0,
  1. Per n = 0, 1, 2, ..., finché non soddisfatti
    m = (an+bn)/2
    Se f(an)f(m) < 0 allora
    an+1 = an
    bn+1 = m
    altrimenti
    an+1 = m
    bn+1 = bn
  2. L'equazione f(x) = 0 ha una soluzione x = m con errore assoluto < bn+1 - an+1.
È necessario precisare cosa significa l'espressione finché non soddisfatti . Un algoritmo, per essere tale, ha bisogno di un termine. Il procedimento usato invece potrebbe essere portato avanti all'infinito: è sempre possibile trovare un nuovo punto medio, calcolarci la funzione, trarre le conclusioni e ripartire daccapo. È necessario quindi aggiungere un criterio di arresto al procedimento: se non lo facessimo, un programma che realizzasse il metodo, proseguirebbe a fare calcoli in eterno senza fornire mai una risposta. Per evitarlo è sufficiente prevedere che ci si fermi dopo un certo numero massimo di iterazioni, qualunque sia l'esito dell'elaborazione. D'altra parte, se ad un certo punto si avesse f(m) = 0, avremmo trovato la soluzione e non avrebbe senso proseguire oltre. Possiamo quindi riscrivere l'algoritmo nella maniera seguente:
Data una funzione f(x) continua nell'intervallo [a0, b0], tale che f(a0)f(b0) < 0 e dato un intero nmax
  1. Per n=0, 1, 2, ..., nmax
    m = (an+bn)/2
    Se f(m) = 0 salta al punto 2.
    Se f(an)f(m) < 0 allora
    an+1 = an
    bn+1 = m
    altrimenti
    an+1 = m
    bn+1 = bn
  2. L'equazione f(x) = 0 ha una soluzione x = m con errore assoluto < bn+1 - an+1.
Per i metodi iterativi è bene prevedere un criterio che arresti le iterazioni allorché il risultato raggiunto è giudicato buono a sufficienza - per esempio se l'errore relativo è minore di una tolleranza fissata. Per il metodo di bisezione, dato che l'errore relativo dimezza ad ogni ciclo, un criterio di questo tipo coincide esattamente col fissare un numero massimo di iterazioni. In generale però i due concetti sono distinti, come si può vedere nel caso del metodo delle corde.


Il metodo delle corde

Sia f(x) una funzione continua nell'intervallo [a0, b0], tale che f(a0)f(b0) < 0. Consideriamo la retta passante per i punti (a0, f(a0)) e (b0, f(b0)). Sia c0 il punto in cui la retta interseca l'asse delle x.
Se f(c0) = 0 abbiamo trovato la soluzione.
Se f(c0)f(b0) < 0, la soluzione si trova in [a1, b1] = [a0, c0];
Se f(c0)f(b0) > 0, la soluzione si trova in [a1, b1] = [c0, b0].
Possiamo iterare il procedimento ottenendo una successione di valori cn che convergono alla soluzione.

A differenza di quanto accade per il metodo di bisezione, non sappiamo come si riduce l'ampiezza dell'intervallo [an, bn]: addirittura non è neanche detto che l'ampiezza tenda a zero. Come criterio di arresto va previsto un numero massimo di iterazioni. Va inoltre introdotto un modo per valutare se il valore ottenuto è soddisfacente: possiamo decidere di ritenerci soddisfatti (cioè di interrompere il metodo prima del numero massimo di iterazioni previste) quando, ad esempio, |f(cn)| < tolleranza.

Concludiamo osservando che il metodo delle corde è più efficiente di quello di bisezione, nel senso che la convergenza alla soluzione è più rapida.


Bibliografia

  • S.D.Conte - C. De Boor, Elementary numerical analisys, McGraw-Hill, 1972
  • F. Fontanella - A. Pasquali, Calcolo numerico, Ed. Pitagora, 1982
  • G. Zwirner, Elementi di analisi matematica, C.E.D.A.M., 1975