• Sistemi di equazioni lineari
• Equazioni non lineari
• Integrazione numerica

Introduzione

Un sistema di n equazioni lineari in n incognite x₁, x₂, ......, x_n può scriversi in notazione algebrica:
 

a₁₁x₁ + a₁₂x₂ + a₁₃x₃ + ....... + a_1nx_n = b₁

a₂₁x₁ + a₂₂x₂ + a₂₃x₃ + ....... + a_2nx_n = b₂

...................................................................... (1.1)

.......................................................................

a_n1x₁ + a_n2x₂ + a_n3x₃ + ...... + a_nnx_n = b_n

ove i termini a_{i j} sono i coefficienti del sistema, mentre b₁, b₂, ..... , b_n sono i termini noti o in forma matriciale compatta:

A X = B (1.2)

dove:

A == matrice dei coefficienti del sistema (1.3)

Il rango della matrice  A (massimo numero di righe o di colonne di A linearmente indipendenti tra loro) si dice anche rango del sistema. Se esso è pari al numero di equazioni (ovvero, nel nostro caso, al numero delle incognite) , il sistema lineare si dice normale.

X =  = vettore colonna delle incognite (1.4)

B =  = vettore colonna dei termini noti (1.5)

Indichiamo inoltre con A’ la seguente matrice:

A’ == matrice dei coefficienti e termini noti o matrice completa (1.6)

Una equazione del sistema può essere scritta anche nella forma:

(i = 1, 2, .... , n) (1.7)

Relativamente alla soluzione di un sistema di questo tipo, si presentano quattro possibilità:
 

• il sistema ammette una unica soluzione, individuata dalla ennupla x₁, x₂, ... , x_n; questo si ha quando il rango della matrice A è pari al numero delle equazioni (ovvero delle incognite), ossia quando il sistema è normale, quindi deve aversi: det A ¹ 0;
• il sistema non ammette soluzioni, caso che si presenta quando le matrici A ed A’ non hanno lo stesso rango (sistema incompatibile);
• il sistema ammette infinite soluzioni, caso che si presenta quando A ed A’ hanno stesso rango, minore del numero delle incognite n;
• il sistema ammette soluzione banale: x₁ = x₂ = ....... = x_n = 0 (sistema di equazioni omogeneo, con rango uguale al numero delle incognite e con vettore colonna dei termini noti nullo).

 Definizioni relative alle matrici
 

• Si chiama matrice diagonale D una matrice quadrata con tutti gli elementi nulli eccetto quelli sulla diagonale principale; come esempio consideriamo la matrice diagonale di ordine 3:

D =  (1.8)

• Si chiama matrice identità I una matrice diagonale avente i termini della diagonale principale tutti uguali ad 1. Essa gode della proprietà: I M = M dove M è una matrice avente numero di righe pari all’ordine di I.
• Si chiama  matrice triangolare una matrice quadrata avente gli elementi da uno stesso lato della diagonale principale tutti nulli (i restanti possono essere nulli o diversi da 0).

 In particolare, una matrice triangolare superiore U ha gli elementi al di sotto della diagonale
principale tutti nulli, come la seguente di ordine 3:

U =  (1.9)

Viceversa per una matrice triangolare inferiore L.

Per la risoluzione del sistema, è opportuno ricordare alcune proprietà e definizioni relativi ai sistemi stessi, alle matrici e al calcolo matriciale.

• Ogni equazione (riga di A’) del sistema può essere moltiplicata per una costante senza modificarne la soluzione.
• Il determinante di una matrice A non cambia se si moltiplicano e dividono tutti gli elementi di una riga o di una colonna per una costante arbitraria diversa da 0.
• L’ordine delle equazioni (righe di A’) può essere cambiato senza modificare la soluzione del sistema .
• Scambiando tra loro due righe o due colonne di una matrice A, il suo determinante non cambia in valore assoluto, ma solo in segno.
• Ogni equazione del sistema può essere sostituita da una combinazione lineare della stessa con ogni altra equazione del sistema, senza modificarne la soluzione.
• Il determinante di una matrice A non cambia se si moltiplicano o dividono tutti gli elementi di una riga o di una colonna per una costante arbitraria.
• Se il determinante di una matrice quadrata A è nullo, la matrice A è detta singolare.
• Se gli elementi di una riga o di una colonna di A sono tutti nulli,il determinante di A è nullo, quindi A è singolare.

• Il determinante di una matrice diagonale e di una matrice triangolare (sia essa inferiore che superiore) è dato dal prodotto degli elementi della diagonale principale.
• E’ possibile trasformare ogni matrice non singolare in una matrice triangolare, cosicché il suo determinante può essere facilmente determinato come detto in precedenza; seguiremo questa strada nell’algoritmo di calcolo del determinante di una matrice.
• Considerate due matrici quadrate A e B, se il prodotto tra le due è uguale all matrice identità, ossia:

A B = I allora la matrice B si dice matrice inversa di A, e si indica con B = A^-1.

Vale la proprietà commutativa rispetto al prodotto con la matrice inversa, ossia:

A^-1 A = A A^-1 = I

I  metodi per  risolvere sistemi di equazioni lineari possono suddividersi in due categorie:

• metodi diretti;
• metodi iterativi.

I metodi iterativi invece, non prevedono un numero preefissato di iterazioni, e portano ad una soluzione "approssimata", nel senso che essa sará comunque affetta da un errore, che teoricamente  puó essere reso piccolo a piacere aumentando il numero delle iterazioni.
Dal punto di vista pratico si fissa il  vettore soluzione X⁽⁰⁾  (prima iterazione)  arbitrariamente, si sostituiscono i suoi valori  nel sistema e se lo soddisfa vuol dire che esso rappresenta la soluzione e quindi si interrompe il calcolo,  altrimenti, in base all'errore commesso si determina un nuovo valore X⁽¹⁾ e si itera il procedimento fino a quando la successione dei valori generati converge verso la soluzione.
Alcuni esempi di metodi iterativi sono il metodo di Jacobi ed il metodo di Gauss - Seidel.
 
Anche se non esiste un criterio che possa indicare a priori quale metodo risolutivo adottare, risulta che i metodi diretti che certamente offrono maggiori garanzie, al crescere delle dimensioni del sistema vengono sostituiti dai i metodi iterativi.

Il metodo di eliminazione di Gauss-Jordan (una variante del metodo di Gauss) è basato sulla trasformazione della matrice

A' nella forma:
                                                                                            _{~                        ~}
cui corrisponde la seguente soluzione del sistema: x₁ = b₁  .......  x_n = b_n

Tale metodo si basa su n iterazioni identificate dall'indice k durante le quali vengono trasformati i coefficienti del sistema..
Durante la prima iterazione k=1 si vuole ridurre a 1 il coefficiente di x1 nella prima equazione, e a 0 in tutte le altre, pertanto:

• si divide la riga k-esima per  a_kk e ciascuna delle rimanenti righe per il coefficiente di x1;
• si sottrae la riga k da tutte le altre righe;

 
Come detto tale trasformazione è stata derivata eseguendo:
sulla riga  k la trasformazione dei coefficienti a_kj  in  a'_kj = a_kj /a_kk  ( j = i, i + 1, .... DIM)
e sulle  righe i la trasformazione dei coefficienti: a_ij in a'_ij = a_ij-a'_kj *a_ik  ( j = i, i + 1, .... DIM)

Si procede quindi con l’iterazione successiva.
Durante la seconda iterazione k=2 si vuole ridurre a 1 il coefficiente di x2 nella seconda equazione, e a 0 in tutte le altre, pertanto:

• si divide la riga k-esima per  a_kk e le rimanenti righe  per il coefficienti di x2;
• si sottrae la riga k da tutte le altre righe;

 
all'ultima iterazione (quella di indice n) la matrice si sarà trasformata nella seguente:

ove nella colonna di indice n+1 troviamo le soluzioni del sistema.

 Al procedimento descritto può essere associato il seguente algoritmo:
 

 
in cui per aumentare l'efficienza sono state eseguite solo le operazioni che contribuiscono al calcolo della soluzione e non quelle che contribuiscono alla riduzione a 0 e 1 dei termini della matrice.

Programma per la risoluzione di sistemi di equazioni lineari con il metodo di  Gauss-Jordan

Si noti che un programma completo dovrebbe prevedere dei controlli sia sul determinante A (diverso da =)  sia sugli elementi a_kk della matrice che dovrebbero essere diversi da 0 (ciò può comunque essere ottenuto scambiando di posto alle equazioni).

                              EQUAZIONI NON LINEARI

Numerosi problemi di matematica applicata corrispondono ad equazioni non lineari del tipo ¦(x) = 0 di cui non si conosce la soluzione in forma chiusa (e quindi l’espressione analitica delle radici) o di cui è complesso il calcolo, pertanto sono indispensabili dei metodi numerici che conducano alla individuazione delle radici.

Quando  l’equazione non lineare ¦(x) ad una incognita soddisfa il teorema di esistenza degli zeri, ossia è reale e continua in un intervallo chiuso e limitato [a,b] ed assumere valori di segno opposto agli estremi dell’intervallo, possono essere usati alcuni metodi iterativi per individuare la soluzione.

Infatti se una funzione ¦(x) è reale e continua in un intervallo chiuso e limitato [a,b] e assume valori di segno opposto agli estremi dell’intervallo (ovvero ¦(a) ^.¦(b) < 0), allora, per il teorema di esistenza degli zeri esiste sicuramente un punto x, detto zero della funzione, interno all’intervallo [a,b] in cui ¦(x) assume il valore 0, che rappresenta la radice dell’equazione ¦(x) = 0.

La ricerca di tutti gli zeri della funzione deve essere preceduta quindi dalla tabulazione della funzione f(x) con passo adeguato in modo da individuare tutti gli intervalli chiusi [a,b] che contengono al loro interno le radici e nell’applicare a ogni singolo intervallo il metodo risolutivo scelto.

La valutazione numerica di ciascuna radice avviene tramite la costruzione di una opportuna successione di valori {x_i}, tale che risulti:

lim {x_n} = x     ove     x Î [a,b]     e        f (x ) = 0
^{n® ¥}

La successione {x_i} dei valori ottenuti dipende dal procedimento iterativo applicato.

Data una successione di valori {x_i} si definiscono:

errore (assoluto) di troncamento e_n: e_n = x - x_n

fattore di convergenza c: c = lim e _n+1/ e _n
                                                            ^{n® ¥}

Nei successivi paragrafi saranno analizzati i metodi di bisezione e delle corde (iterativi a due punti) ed il metodo di Newton-Raphson (iterativo a un punto).

#include <iostream.h>
main()
{
double prec=1.0;
while ((1.+ prec)!=1.0) prec/=2.;
prec*=2;
cout<<prec;
}
 

timeelapsed% a.out
2.22045e-16timeelapsed%
 
 

#include <iostream.h>
main()
{
float prec=1.0;
float a=1.0;
float x=2.0;
while ((a+ prec)!=a) prec/=x;
prec*=2;
cout<<prec;
}
 

timeelapsed% g++ b.c
timeelapsed% a.out
1.19209e-07timeelapsed%
 

Il metodo di bisezione è il metodo iterativo più semplice e sicuramente convergente per calcolare le radici di un’equazione non lineare ¦(x) reale e continua in un intervallo chiuso e limitato [a,b] che assuma valori di segno opposto agli estremi dell’intervallo.

Il metodo di bisezione consiste nel considerare il punto medio x_m =(a + b)/2 dell’intevallo [a,b] e verificare se x_m è una radice.
Per verificare questa condizione, non si controlla se ¦ (x_m ) = 0, ma (poiché si lavora con i numeri reali) se |¦ (x_m) | < p, ove p rappresenta la precisione scelta.
Se x_m non è una radice si valuta il segno del prodotto ¦(a) ^.¦(x_m) e si procede allo stesso modo nel sottointervallo [a, x_m] se il segno è negativo, o nel sottointervallo [x_m,b] se è positivo.

Naturalmente iterando il procedimento si genera una successione di termini {x_i}, ciascuno dei quali rappresenta l’intersezione con l’asse X della retta individuata dai punti di coordinate (a, segno f(a)) e (b, segno f(b)) che solo sotto le ipotesi fatte converge sicuramente alla radice richiesta.

Il metodo di bisezione, come risulta osservando la Fig.1.1 ad ogni iterazione dimezza l’intervallo contenente la radice che viene calcolata con la precisione desiderata.

Il punto di forza di questo metodo è nella sua semplicità e sicura convergenza, il suo punto debole è invece nella scarsa efficienza rispetto agli altri metodi infatti poiché ad ogni passo viene dimezzato l’intervallo contenente la radice, pertanto richiede un numero elevato di iterazioni.

PROGRAM EQUAZNLIN

Il metodo delle corde è un metodo iterativo più efficiente per calcolare le radici di un’equazione non lineare ¦(x) reale e continua in un intervallo chiuso e limitato [a,b] che assuma valori di segno opposto agli estremi dell’intervallo, condizione che, come detto in precedenza, assicura la presenza di almeno una radice.

Il metodo delle corde (o della falsa posizione), la cui interpretazione grafica è riportata in Fig.1.2, consiste nel calcolare se l’intersezione x_i della retta (corda) passante per i punti A,B di coordinate (a,ƒ(a)) e (b, ƒ(b)) rispettivamente con l’asse delle ascisse è una radice.

Poiché l’equazione di una retta passante per i punti A,B ha equazione:

x  - a    =   y- ƒ(a)
b - a            ƒ(b) -ƒ(a)

 
ponendo y=0 si ottiene:

x= a -   (b - a)  * f(a)
                ƒ(b) -ƒ(a)
o ancora meglio:

x=   a*f(b) - b* f(a)
            ƒ(b) -ƒ(a)
Per verificare questa condizione, non si controlla se ¦ (x_i  ) = 0, ma (poiché si lavora con i numeri reali) se |¦ (x_i ) | < p, ove p rappresenta la precisione scelta.
Se x_i  non è una radice si valuta il segno del prodotto ¦(a) ^.¦(x_i ) e si procede allo stesso modo nel sottointervallo [a, x_i ] se il segno è negativo, o nel sottointervallo [x_i ,b] se è positivo.

Pertanto iterando il procedimento si genera una successione di termini {x_i}, il cui termine (i+1)-esimo è:

x_i+1 =   x_i-1 *f( x_i ) -  x_i * f( x_i-1 )        i= 2,3,4…..
                    ƒ(x_i) - ƒ(x_i-1)

Poiché il procedimento su cui è basato il metodo delle corde approssima la funzione con la retta individuata dai punti di coordinate (a, f(a)) e (b, f(b)) si può dimostrare che l’errore commesso durante ogni iterazione è minore e pertanto è necessario un minor numero di iterazioni, come può essere provato facendo eseguire il programma.

METODO DELLE CORDE
Al procedimento descritto (sostituendo ad a x1 e a b x2) può essere associato il seguente algoritmo:
 

Si noti che nell'algoritmo abbiamo inserito anche il controllo su |f2-f1|<p  in quanto nel caso che fossimo vicini ad un punto singolare (massimo, minimo, flesso) si tenterebbe di dividere per 0 e ciò causerebbe un errore.
In un programma completo andrebbe previsto che in tale ipotesi si ricominci l'esecuzione con due diversi valori iniziali x1, x2.

METODO DI NEWTON-RAPHSON

Il metodo di Newton-Raphson è tra i metodi iterativi descritti il più efficiente.

Esso consiste nel calcolare se l’intersezione x₁ della retta tangente al grafico della funzione passante per un punto A di coordinate (x₀,ƒ(x₀)) con l’asse delle ascisse è una radice.

Per verificare questa condizione, non si controlla se ¦ (x₁ ) = 0, ma (poiché si lavora con i numeri reali) se |¦ (x₁) | < p, ove p rappresenta la precisione scelta.

Ricordando che:

• l’equazione del fascio di rette di centro P (x₀,y₀) è:        y – y₀= m (x– x₀)     ove m è il coefficiente angolare e che
• ƒ' (x₀), ossia la derivata prima di una funzione ƒ (x) calcolata nel punto x₀, coincide con il coefficiente angolare m della retta tangente al grafico della funzione nel punto di coordinate (x₀,ƒ(x₀)),

y – ƒ(x₀)= ƒ' (x₀) (x– x₀)

pertanto ponendo y=0 si ottiene:
 – ƒ(x₀)= ƒ' (x₀) x – ƒ' (x₀) x₀

quindi l'ascissa x₁in cui la tangente incontra l'asse delle X:

x₁= x₀ – ƒ(x₀) / ƒ' (x₀)

 
Se  se |¦ (x₁) | < p, ovvero se x₁  è una radice si interrompe l'elaborazione, altrimenti si procede con un altro punto x₂ intersezione della retta tangente al grafico della funzione passante per un punto di coordinate (x₁,ƒ(x₁)), applicando iterativamente il metodo si genera una successione di termini {x_i}, il cui termine (i+1)-esimo è:

x_i+1 = x_i – ƒ(x_i)/ ƒ' (x_i) i = 0,1,2….

Poiché il procedimento su cui è basato il metodo di Newton approssima la funzione con la retta tangente al grafico della funzione si può dimostrare che l’errore commesso durante ogni iterazione è minore rispetto a quello commesso nei metodi precedenti, infatti si dice che il metodo di Newton è del secondo ordine o che ha convergenza quadratica, cioè l’errore al passo n-esimo è proporzionale al quadrato dell’errore al passo precedente e pertanto richiede un minor numero di iterazioni, come può essere provato facendo eseguire il programma.
Si noti però che, a differenza degli altri metodi, la convergenza della successione { x_i } non è garantita per ogni x₀ e per ogni funzione ƒ.

Per poter codificare l’algoritmo di Newton occorre introdurre anche il calcolo delle derivate delle funzioni oltre che di ƒ(x).

Per quanto attiene il calcolo della derivata di una funzione ricordandone la definizione:

ƒ' (x_i) = lim       f(x+D x) - f(x)
            ^{D x® 0}           D x

basta fissare un valore molto piccolo per D x e calcolare numericamente il valore che la funzione assume in corrispondenza della x.

METODO DI NEWTON_RAPHSON

Data la funzione f(x) di cui si riporta il grafico, i punti P,  Q  e A individuano il triangolo rettangolo PQA.
 

Possiamo scrivere che :
 

Da cui deriva la corrispondenza tra il coefficiente angolare m ()

con il valore della derivata prima della funzione calcolata nel punto di coordinate x f(x).
 

La valutazione di integrali definiti quando non è nota la primitiva della funzione integranda o quando il procedimento analitico risulta complesso richiede l’applicazione di metodi numerici.
Le tecniche numeriche si basano sull’applicazione di una formula che usa i valori che la funzione integranda assume all’interno dell’intervallo di integrazione.
Ogni formula di integrazione numerica ha un suo grado di precisione infatti:

 _{b                      n}
ò f(x) =     å c_i f(x) + R_n (f)
^{a                     i=1}

il primo termine è detto parte approssimante e R_n (f) errore di troncamento.

La valutazione numerica dell’integrale avviene tramite la costruzione di una opportuna successione di valori {x_i}, tale che risulti:

                 _{n                         b}
lim     å c_i f(x) = ò f(x) e R_n (f) = 0
^{n® ¥    i=1                      a}

La successione {x_i} dei valori ottenuti dipende dal procedimento iterativo applicato.

Nei successivi paragrafi saranno analizzati alcuni tra i metodi di integrazione numerica più semplici da realizzare: il metodo dei rettangoli, dei trapezi e di Cavalieri-Simpson.

Il metodo dei rettangoli è il metodo iterativo più semplice per calcolare l’integrale definito di una funzione ¦(x) reale e continua in un intervallo chiuso e limitato [a,b].

Esso consiste nel suddividere l’intervallo di integrazione [a,b] in un certo numero di sottointervalli n e nell’approssimare la funzione all’interno di ciascun intervallo con il valore costante che questa assume nel suo punto medio, pertanto il calcolo dell’integrale numerico viene quindi a corrispondere con la somma dei rettangoli aventi come base l’ampiezza costante dei sottointervalli h = (b-a)/n e come altezza il valore che la funzione assume nel punto medio di ciascun rettangolo.

La formula risultante, nel caso in cui si scelga n=1 risulta:
   _b
ò f(x) » h *f (a+h/2)      ove h= (b-a)
^a

Nel caso in cui n=2, h= (b-a)/n

_b
ò f(x) » h* f (a+h/2) + h f(a+h/2+h)
^a

mentre per n generico e h= (b-a)/2n si ha:
_{b                  n-1}
ò f(x) »  å h f(a+(1/2+i) h)
^{a                  i=0}
Il procedimento di calcolo dell’integrale, può essere semplicemente realizzato prevedendo n iterazioni che conviene numerare con i da 0 a n-1 e incrementare ad ogni passo, come è riportato nel programma, la somma integrale con il valore: F(A+(i+.5)*h) .
infatti la x dovrebbe assumere i valori:
A+h/2, A+h/2+h,A+ h/2+h+h,....ovvero ad ogni iterazione x= A+(i+0.5)h
Per una migliore approssimazione del risultato e per non incorrere in problemi di convergenza il calcolo dell’integrale viene ripetuto dimezzando a ogni passo l’ampiezza dei sottointevalli, fino a quando la differenza tra due valori successivamente calcolati è minore della precisione che si vuole raggiungere nel calcolo.

Per verificare questa condizione (poiché si lavora con i numeri reali) si controlla se |I-I1 | < p, ove p rappresenta la precisione scelta.

Nel metodo dei rettangoli, la cui interpretazione grafica è riportata in Fig.2.1 poiché si approssima la funzione, in ogni sottointervallo, con la retta parallela all’asse x di equazione y = f(x_i) viene commesso un errore di approssimazione che richiede un elevato numero di iterazioni prima di giungere al risultato.

Esso consiste nel suddividere l’intervallo di integrazione [a,b] in un certo numero di sottointervalli n e nell’approssimare la funzione all’interno di ciascun intervallo con la spezzata che congiunge i punti per cui passa il grafico della funzione agli estremi di ogni intervallo, pertanto il calcolo dell’integrale corrisponde alla somma dei trapezi aventi come altezza l’ampiezza dei sottointervalli SPEP = (b-a)/n e come base maggiore e minore il valore che la funzione assume agli estremi dei singoli sottointervalli.
La formula risultante, nel caso in cui si scelga n=1 risulta:

_b
ò f(x) » h/2 ( f (a+h) + f(a))       ove    h= (b-a)
^a
Nel caso in cui n=2, h= (b-a)/n

_b
ò f(x) » h/2 ( f (a+h) + f(a)) + h/2 ( f (a+2h) + f(a+h))
^a
mentre per n generico e h= (b-a)/2n si ha:

 _{b                  n-1}
ò f(x) » å h/2 [f(a+(i+1) h) + f(a+i h)]           oppure
^{a                  i=0}
_{b                                               n-1}
ò f(x) » h/2 [f(a+f(b)] + å h [f(a+i h)] (detta formula di Bezout)
^{a                                               i=0}

 
Il calcolo dell’integrale, di cui in Fig.2.2.ne è riportata l’ interpretazione grafica, può essere semplicemente realizzato prevedendo n iterazioni che conviene numerare con un indice i da 0 a n-1 e incrementare ad ogni passo la somma integrale con il valore:
h*.5* ((F(a+i*h)+F(a+(i+1)*h)))

Infatti durante il procedimento di calcolo si deve generare una successione di termini ( ciascuno rappresentante l’area di un trapezio), ove si moltiplica l’altezza costante dei trapezi, pari ad h, per la somma delle due basi diviso 2, pertanto, come è riportato nel programma:

I=I+(h*.5*(F(a+i*h)+F(a+(i+1)*h)))

Per rendere minimo l’errore e per non incorrere in problemi di convergenza il procedimento viene ripetuto dimezzando a ogni passo l’ampiezza dei sottointevalli, fino a quando la differenza tra due valori successivamente calcolati diventa trascurabile ovvero (poiché si lavora con i numeri reali) fino a quando |I2-I1 | < p, ove p rappresenta la tolleranza scelta.

S può dimostrare che l’errore commesso usando la formula trapezoidale, che risulta :

R(t) < 1/12 (x_i - x_i+1)^3• f¢ ¢ (x )

è molto minore rispetto a quello commesso usando la formula rettangolare, pertanto il risultato ottenuto a parità di sottointervalli risulta più approssimato.

 
SUBROUTINE TRAPEZI

METODO DI CAVALIERI-SIMPSON

Il metodo di Cavalieri-Simpson consiste nel suddividere l’intevallo [a,b] in un certo numero di sottointervalli n e nell’approssimare la funzione, all’interno di ogni sottointervallo, con una parabola individuata dai punti di coordinate:
(x_i, f(x_i)), (x_i+1, f(x_i+1)) e dal punto intermedio (x_i + (b-a)n) f(x_i + (b-a)n)).

La formula risultante, nel caso in cui si scelga n=1 risulta:
_b
ò f(x) » (h/3)[f(a)+4f (a+h)+f( b)]                  ove h= (b-a)/2
^a

Nel caso in cui n=2, h= (b-a)/2n
_b
ò f(x) » (h/3)[f(a)+4f (a+h)+f( a+2h)] + (h/3)[f(a+2h)+4f (a+3h)+f( b)]
^a

mentre per n generico e h= (b-a)/2n si ha:

_b
ò f(x) »   (h/3)[f(a)+4f (x₁)+2f(x₂)+....+ 4f(x_2n-1)+f(x_2n)]
^a

Il procedimento di calcolo dell’integrale, può essere semplicemente realizzato inizializzando la somma integrale con la sommatoria dei primi due termini e dell’ultimo termine:

I=(F(a)+F(b)+4*F(a+h*(2*n-1)))*h/3.0

tramite l’iterazione di n passi che conviene numerare tramite un indice i da 0 a n-1, ove la somma integrale verrà incrementata con la coppia di termini successivi di coefficiente 4 e 2:

I=I+h/3.0*(4*F(a+i*h)+2*F(a+(i+1)*h))

Si noti che il calcolo dell’integrale secondo questo metodo, come risulta osservando la Fig.2.3 prevede la ulteriore divisione a metà di ogni intervallo, pertanto dopo aver fissato il valore iniziale di n questo sarà subito raddoppiato per poter iniziare correttamente il calcolo.
Come nei casi precedenti per rendere minimo l’errore e per non incorrere in problemi di convergenza il procedimento viene ripetuto raddoppiando a ogni iterazione il numero dei sottointevalli, fino a quando la differenza tra due valori successivamente calcolati è minore della precisione che si vuole raggiungere nel calcolo, ovvero si controlla se |I2-I1 | < p, ove p rappresenta la tolleranza scelta.

Il metodo di Cavalieri-Simpson è uno tra i metodi più utilizzati per la notevole precisione che è in grado di raggiungere con pochi intervalli di calcolo, infatti è possibile dimostrare che l’errore di calcolo commesso risulta :

R(t) < 1/90 (x_i - x_i+1)^5• f¢ ¢ (x )

ovvero molto minore rispetto alle formule presentate in precedenza.