1

ALCUNI ESEMPI DI ALCUNI ESEMPI DI ALGORITMI DI CALCOLO NUMERICOALGORITMI DI CALCOLO NUMERICO

La risoluzione di modelli matematici complessi richiede spesso l’uso di metodi di calcolo approssimato:

calcolo “numerico”

Riferimenti Bibliografici:

L.M. Barone, E. Marinari, G. Organtini, F. Ricci-Tersenghi, Programmazione Scientifica, Pearson Education, Milano, 2006G. Monegato, Fondamenti di Calcolo Numerico, Libreria Editrice Universitaria Levrotto&Bella, Torino

2

CALCOLO NUMERICOCALCOLO NUMERICOSi usa quando:

non esistono algoritmi esatti di risoluzione, oppuregli algoritmi esatti non sono efficienti (tempo, memoria): alto costo computazionale

Soluzione di problemi attraverso algoritmi euristici o approssimanti:

Determinano soluzioni "approssimate" del problema

Ad esempio, si basano su tecniche del tipo:discretizzazione di domini continui in domini discretiapprossimazioni successive => metodi iterativi

3

DISCRETIZZAZIONEDISCRETIZZAZIONEValori di un dominio continuo vengono rappresentati da elementi di un dominio discreto:

ξ viene rappresentato da V3

Con quale approssimazione? Con quali effetti sui risultati successivi (propagazione dell’approssimazione)?

V1V2V3V4

Dcontinuo

Ddiscreto

ξ

4

APPROSSIMAZIONI SUCCESSIVEAPPROSSIMAZIONI SUCCESSIVEvalore iniziale V0 (di primo tentativo)

formula iterativa => ottengo V1

valore di secondo tentativo V1 formula iterativa => ottengo V2

valore di terzo tentativo V2.......

L'algoritmo termina quando due valori consecutivi Vi e Vi+1sono "sufficientemente vicini":

| Vi+1 - Vi |<= ε scarto assolutoscarto assoluto| (Vi+1 - Vi) / Vi+1 | <= ε scartoscarto relativorelativo

Il valore ε rappresenta la precisione desiderata e viene stabilito arbitrariamente

5

Note:

Il valore di primo tentativo dovrebbe essere stimato nel modo migliore possibile (conoscenza del problema, ordine di grandezza, segno, …)Se ogni successiva approssimazione è più vicina alla risposta corretta e la differenza tende a calare ad ogni passo, l'iterazione converge (come capirlo su un numero finito di iterazioni?)Non tutte le iterazioni convergono (divergenti) o convergono in modo soddisfacente per essere utili

APPROSSIMAZIONI SUCCESSIVEAPPROSSIMAZIONI SUCCESSIVE

6

ESEMPIO: CALCOLO ESEMPIO: CALCOLO DELLA RADICE QUADRATADELLA RADICE QUADRATA

Dato un numero reale positivo R, calcolare:Q = SQRT(R)

1) Inizializzazione: Q1 = R/2; i=22) Formula Iterativa (passo i-esimo, i = 2, 3, ....)

Qi = (Qi-1 + R/Qi-1)/23) Condizione di arresto:

se |(Qi - Qi-1)/Qi)|<= ε or i >=MAXITERaltrimenti

i = i+1, torna al passo 2

Infatti:Supponiamo Q1 circa uguale (≈) a Q0, allora:(Q1 - Q0)2 ≈ 0; (Q1

2 - 2Q1*Q0 + Q02) ≈ 0

2Q1*Q0≈ Q02 + R (Q12 ≈ R)

Q1 ≈ (Q0 + R/Q0)/2

7

CALCOLO DELLA RADICE QUADRATACALCOLO DELLA RADICE QUADRATAEsempio:

R = 9Q1 = 9/2 = 4.5Q2 = (4.5 + 9/4.5)/2 = (6.5)/2 = 3.25Q3=(3.25+9/(3.25))=(3.25+ 2.769)/2=(6.019)/2=3.0095.......

8

Radice Quadrata: ALGORITMO IN CRadice Quadrata: ALGORITMO IN C#include <stdio.h>#include <math.h>#define MAX 1000#define EPS 0.00001

double radice_q(double R){ double Q, Qold;

int i, fine=0;Qold=R/2;for (i=0; i<MAX && !fine; i++)

{ Q=(Qold + R/Qold)/2;if (abs((Q-Qold)/Q)<EPS)

fine=1;else Qold=Q;

}if (i>=MAX) printf("Numero iterazioni eccessivo\n");return Q;

}

9

INVOCAZIONEINVOCAZIONE

main(){double V;printf("dammi il valore: ");scanf("%lf", &V);printf("\nRisultato:%lf (valore atteso:

%f)\n", radice_q(V), sqrt(V));

}

10

ESEMPIO: seno() DI UN ANGOLOESEMPIO: seno() DI UN ANGOLODato un angolo X espresso in radianti, calcolare:

S = seno(X)

Sviluppo in Serie di Taylor:

S = X − X3

3!+ X

5

5!− X

7

7!+ ... = (−1)i . X

2i+1

(2i +1)!i=0

∞∑

⇒ S = Tii=0

∞∑

Ti = (−1)i X2 i+1

(2 i + 1)!=

= −(−1)i−1X 2( i−1)+1

(2(i − 1) + 1)!⋅

X 2

(2i + 1)2i=

= −Ti−1⋅X 2

(2i + 1)2i

11

SENO DI UN ANGOLOSENO DI UN ANGOLOSi arresta lo sviluppo in serie dopo M termini,

con M (<= MAXITER) tale che:

Passo iterativo:

Termine i-esimo:Ti = (-1)i X2i+1 / (2i+1)! = -(-1)i-1 (X2(i-1)+1/(2(i-1)+1)! ) * X2/((2i+1)2i)= - X2/((2i+1)2i )* Ti-1

se S = (−1)i X 2 i +1

(2i + 1)!i= 0

M

∑X 2 M +1

(2M + 1)!≤ ε

12

UNA POSSIBILE IMPLEMENTAZIONEUNA POSSIBILE IMPLEMENTAZIONE#include <stdio.h>#include <math.h>#define EPS 0.0000001#define MAX 10000

double seno(double x){ double T, Told, sum;int i, fine=0;Told=x; sum=x;for (i=1; i<MAX && !fine; i++) {

T= -Told*x*x/((2*i+1)*2*i);sum = sum + T;if (abs(T-Told)<EPS) fine=1;else Told=T;

}if (i==MAX) printf("troppe iterazioni\n");return sum;

}

13

INVOCAZIONEINVOCAZIONE

Ad esempio:

main(){ double v;

printf("Dammi x: ");scanf("%lf", &v);printf("Seno di %lf: %lf (valore atteso %lf)\n",

v, seno(v), sin(v));}

14

EQUAZIONI: SOLUZIONI APPROSSIMATEEQUAZIONI: SOLUZIONI APPROSSIMATEIn problemi di ingegneria (e anche talora di ingegneria chimica…) è necessario spesso trovare i valori delle radici di un'equazione

Per equazioni lineari o quadratiche la soluzione èsemplice (esiste forma analitica chiusa):

aX + b = 0 ==> X = -b/a

aX2 + bX + c = 0 ==> X1= -( b + SQRT(b2-4ac))/2aX2= -( b - SQRT(b2-4ac))/2a

[ Hp: ∆ = (b2-4ac) >= 0 ]

15

Ricerca delle (eventuali) radici reali di una funzione che supporremo definita e continua in un certo intervallo dell'asse delle ascisse

La ricerca delle radici approssimate si compone di due fasi:

1) separazione delle radici, che consiste nel determinare degli intervalli [a, b] contenenti una sola radice

2) calcolo di un valore approssimato della radice e miglioramento di tale valore fini ad ottenere la precisione desiderata (tramite iterazioni successive)

CALCOLO DELLE RADICI DI UNA FUNZIONECALCOLO DELLE RADICI DI UNA FUNZIONE

16

FUNZIONI CONTINUE: PROPRIETFUNZIONI CONTINUE: PROPRIETÀÀProprietà

Se una funzione continua f(x) assume in due punti a e bvalori di segno opposto, esiste almeno un valore ξ (o un numero dispari di punti) compreso fra a e b in cui f(x)=0

ab

f

xξ

17

METODO DELLA BISEZIONEMETODO DELLA BISEZIONE

Dividiamo l'intervallo [a,b] a metà c = (a + b)/2-> calcoliamo f(c)

Se f(c) = 0, la radice è trovata (procedimento terminato)Se f(c) > 0, trascuriamo l'intervallo destro [c, b], poniamo b = c e ripetiamo la procedura nel nuovo intervallo [a, b]Se f(c) < 0, trascuriamo l'intervallo sinistro [a, c], poniamo a = c e ripetiamo la procedura nel nuovo intervallo [a, b]Proseguendo in questo modo o si trova la radice esatta o si raggiunge un intervallo [a, b] < 2ε

dove ε è la precisione desiderataIn pratica, si approssima la funzione con la retta che passa per i punti (a, sign(f(a))), (b, sign(f(b)))

ξa bc

f

x

18

METODO DELLA BISEZIONEMETODO DELLA BISEZIONE#include <stdio.h>#include <math.h>#define JMAX 1000#define EPS 0.00001

double myfun(double z) { // ad es. per myfun= z3/2return z*z*z/2; }

double bisez(double a, double b, double xacc)

{ int i, fine=0;double Xa, Xb, Fa, Fb, Xm, Fm, tmp;

if (a>b) // estremi non ordinati{tmp=a; a=b; b=tmp;} // scambio

Fa=myfun(a); Fb=myfun(b);if (Fa*Fb>=0) { printf("intervallo non valido\n");

exit(); }else { Xa=a; Xb=b; }

…

19

…for (i=0; i<JMAX && !fine; i++){

Fa=myfun(Xa); Fb=myfun(Xb);Xm=(Xa+Xb)*0.5;Fm=myfun(Xm);if (Fm*Fa<0) ;

else Xa=Xm;if (Xa!=Xm) Xb=Xm;if ((abs(Xb-Xa)< xacc) || (Fm==0)) fine=1;}

if (fine) return Xm;else

{ printf("troppe iterazioni\n");exit(); }

}

METODO DELLA BISEZIONEMETODO DELLA BISEZIONE

20

CALCOLO DELLE RADICI:CALCOLO DELLE RADICI:CORDE, SECANTI, NEWTONCORDE, SECANTI, NEWTON--RAPHSONRAPHSON

Idea di Base:

Detto ξ lo zero di f (appartenente all’intervallo [a,b]), sia x0 una arbitraria approssimazione dello zero ξnell’intervallo [a,b]:Si approssima la funzione con una retta passante per (x0,,f(x0)) la cui equazione è

y=K0(x-x0) +f(x0)L’intersezione tra la retta e l’asse delle ascisse dà origine alla nuova approssimazione x1 della radice ξ

f

ξx0

21

Quindi, si ripete iterativamente:Ad ogni iterazione i si ottiene una nuova approssimazione xi+1dello zero, risolvendo l’equazione:

Ki(x-xi) +f(xi)=0dove xi, soluzione approssimata dell’iterazione precedente

CALCOLO DELLE RADICI:CALCOLO DELLE RADICI:CORDE, SECANTI, NEWTONCORDE, SECANTI, NEWTON--RAPHSONRAPHSON

Questi metodi si basano tutti su approssimazioni successive della funzione f con rette del tipo:

y=Ki(x-xi) +f(xi)

Ogni metodo differisce dagli altri per la scelta del coefficiente angolare Ki

xia b x

f

22

METODO DELLE CORDE (Regula Falsi)METODO DELLE CORDE (Regula Falsi)

Dati xa, xb, due valori reali appartenenti ad [a, b] tali che sign(f(xa))!=sign(f(xb))

Retta passante per f(xa) , f(xb)

intersezione con l’asse x: x1

Ad ogni iterazione, il coefficiente angolare KKii della retta che approssima f è calcolato come:

dove:xlim=xa se sign(f(xa))!=sign(f(xi))xlim=xb se sign(f(xa))=sign(f(xi ))

Ki =f (xi ) − f (x lim )

xi − xlim

23

Coefficiente angolare: Ki =f (xi )− f (xlim )xi − xlim

x0

xa x1x2

xb

METODO DELLE CORDE (Regula Falsi)METODO DELLE CORDE (Regula Falsi)

24

METODO DELLE CORDE (Regula Falsi)METODO DELLE CORDE (Regula Falsi)…#define JMAX 1000#define EPS 0.00001

double myfun(double z) { // ad es. per myfun= z3/2return z*z*z/2; }

double K (double xi, double xa, double xb){double lim, cond;

cond=myfun(xi)*myfun(xa);if (cond<0) lim=xa; else lim=xb;return (myfun(xi)-myfun(lim))/(xi-lim);

}

double corde(double a, double b, double xacc){ int i, fine=0;

double Fa, Fb, Xi, Fi, Xold, tmp;if (a>b) /* estremi non ordinati */

{tmp=a; a=b; b=tmp; }

25

METODO DELLE CORDE (Regula Falsi)METODO DELLE CORDE (Regula Falsi)…

Fa=myfun(a);Fb=myfun(b);

if (Fa*Fb>=0){printf("intervallo non valido\n");exit();}

Xold=(a+b)*0.5; /*inizializzazione*/

for (i=0; i<JMAX && !fine; i++) {Xi=-myfun(Xold)/K(Xold,a,b) + Xold;Fi=myfun(Xi);if ((abs(Xi-Xold)< xacc) || (Fi==0)) fine=1;else Xold=Xi; }

if (fine) {printf("eseguite %d iterazioni\n",i); return Xi;}else { printf("troppe iterazioni\n");

exit(); }}

26

METODO DELLE SECANTIMETODO DELLE SECANTISia x0 , un valore reale appartenente ad [a, b]

Si traccia la retta passante per f(x0 ) e f(xa ) se sign(f(x0 )) != sign(f(xa ))f(x0 ) e f(xb ) se sign(f(x0 )) != sign(f(xb ))

individuo l’intersezione x1 con l’asse x

Ad ogni iterazione i, il coefficiente angolare KKii viene calcolato come:

Ki = f(xi)− f(xi−1)xi −xi−1

x0xa x1

x2

x3

27

Metodo delle SecantiMetodo delle SecantiDiversamente dal metodo di bisezione, il metodo delle secanti può non convergere:

in caso di convergenza, il metodo delle secanti è piùefficiente di quello delle corde

x1

x0x2

28

METODO DELLE SECANTIMETODO DELLE SECANTI…#define JMAX 1000#define EPS 0.00001

double myfun(double z) { // ad es. per myfun= z3/2return z*z*z/2; }

double Ks(double Xnew, double Xold) {return (myfun(Xnew)-myfun(Xold))/(Xnew-Xold);

}

double secanti(double a, double b, double xacc) {int i, fine=0;double Fa, Fb, Xi, Fi, Xnew,Xold, tmp;if (a>b) /* estremi non ordinati */

{tmp=a; a=b; b=tmp;}Fa=myfun(a);Fb=myfun(b);if (Fa*Fb>=0) {

printf("intervallo non valido\n");exit(); }

29

METODO DELLE SECANTIMETODO DELLE SECANTI…

Xold=a;Xnew=b;

for (i=0; i<JMAX && !fine; i++) {Xi=-myfun(Xold)/Ks(Xnew,Xold) + Xold;Fi=myfun(Xi);if ((abs(Xi-Xold)< xacc) || (Fi==0)) fine=1;else {Xold=Xnew;Xnew=Xi;}}

if (fine) {printf("eseguite %d iterazioni\n",i); return Xi;}else {

printf("troppe iterazioni\n");exit(); }

}

30

METODO DELLE TANGENTI METODO DELLE TANGENTI (Newton(Newton--Raphson)Raphson)

Ad ogni iterazione i, il coefficiente angolare KKii della retta che approssima f è la derivata della funzione nel punto xi

Ovviamente è necessario che la funzione sia derivabilenell’intervallo considerato

Ki = f ' (xi )

X0X1Xξ

31

Anche il metodo delle tangenti può non convergere:

imponiamo anche la condizione che derivata seconda sia continua e mantenga segno costantenell'intervallo: f(X)*f''(X) > 0

condizione sufficiente per la convergenza del metodo

X0X1

X2

X3

METODO DELLE TANGENTI METODO DELLE TANGENTI (Newton(Newton--Raphson)Raphson)

32

…#define JMAX 1000#define EPS 0.00001#define DELTA 0.001

double myfun(double z) { // ad es. per myfun= z3/2return z*z*z/2; }

double Kt(double xi){ return (myfun(xi+DELTA)-myfun(xi))/DELTA;}

double tangenti(double a, double b, double xacc) {int i, fine=0;double Fa, Fb, Xi, Fi, Xold, tmp;if (a>b) /* estremi non ordinati */

{tmp=a; a=b; b=tmp;}Fa=myfun(a);Fb=myfun(b);if (Fa*Fb>=0) {

printf("intervallo non valido\n");exit(); }

METODO DELLE TANGENTI METODO DELLE TANGENTI (Newton(Newton--Raphson)Raphson)

33

…Xold=(a+b)*0.5; /*inizializz. */for (i=0; i<JMAX && !fine; i++) {

Xi=-myfun(Xold)/Kt(Xold) + Xold;Fi=myfun(Xi);if ((abs(Xi-Xold)< xacc) || (Fi==0)) fine=1;else Xold=Xi;}

if (fine) {printf("eseguite %d iterazioni\n",i); return Xi;}

else { printf("troppe iterazioni\n");exit(); }

}

METODO DELLE TANGENTI METODO DELLE TANGENTI (Newton(Newton--Raphson)Raphson)