Metode Numerice - Seminar

77
Universitatea din Bacău Facultatea de Ştiinţe METODE NUMERICE Note de seminar pentru studenţii Facultăţii de Inginerie

description

met num sem

Transcript of Metode Numerice - Seminar

Page 1: Metode Numerice - Seminar

Universitatea din BacăuFacultatea de Ştiinţe

METODE NUMERICENote de seminar

pentru studenţii Facultăţii de Inginerie

Conf. univ.dr. Mihai TALMACIUAsist. univ. drd. Alina – Mihaela PATRICIU

Page 2: Metode Numerice - Seminar

Cuprins

Introducere 3

Capitolul 1Rezolvarea numerică a ecuaţiilor algebrice neliniare 7

Capitolul 2Rezolvarea numerică a sistemelor de ecuaţii algebrice liniare prin metode directe 15

Capitolul 3Rezolvarea numerică a sistemelor de ecuaţii algebrice liniare prin metode iterative 23

Capitolul 4Determinarea rădăcinilor unui polinom 35

Capitolul 5Rezolvarea numerică a problemelor algebrice de valori şi vectori proprii 41

Capitolul 6Polinomul de interpolare Lagrange 51

Capitolul 7Diferenţe divizate 57

Capitolul 8Polinomul de interpolare Newton 61

Capitolul 9Interpolarea prin polinoame trigonometrice 67

Capitolul 10Diferenţe finite 73

Capitolul 11Formule de interpolare pe noduri echidistante 77

Capitolul 12Formule de derivare numerică pe noduri echidistante 87

Bibliografie 97

1

Page 3: Metode Numerice - Seminar

Introducere

Drumul parcurs pentru rezolvarea unei probleme dintr-un domeniu oarecare cu ajutorul calculatorului constă în: stabilirea unui model matematic al problemei concrete (model ce se poate încadra într-o categorie cum ar fi: o ecuaţie neliniară, un sistem de ecuaţii liniare sau neliniare), care fiind de multe ori de natură continuă trebuie discretizat; soluţia problemei discretizate trebuie să fie consistentă şi stabilă (robustă); modelul discretizat trebuie transpus într-un algoritm, descris de obicei într-un limbaj de programare.

Pentru parcurgerea şi utilizarea unui asemenea material, cititorul are nevoie de cunoştinţe de matematică şi de programare în limbajele de programare Pascal şi C la îndemâna studenţilor care au promovat primul an de studiu al oricărei facultăţi cu profil tehnic, matematico-informatic sau economic.

Materialul cuprins în aceastǎ carte îşi are originea deopotrivǎ în experienţa de programare (la C.T.C.E. Bacǎu) şi în seminariile/laboratoarele de analizǎ numericǎ, calcul numeric şi metode numerice ţinute la secţiile Facultǎţii de Inginerie şi la secţia de Matematicǎ-Fizicǎ a Facultǎţii de Ştiinţe a Universitǎţii din Bacǎu.

Ca limbaje de programare au fost alese limbajul Pascal şi C, datoritǎ, pe de o parte, faptului cǎ prin modernitatea structurilor implementate permit abordarea sistematicǎ a noţiunilor de programare, iar pe de altǎ parte, au fost alese ambele limbaje deoarece se predau în anii anteriori în facultǎţile mai sus menţionate.

Implementarea algoritmilor a urmǎrit evidenţierea specificului metodelor numerice şi a implicat un compromis între claritatea programelor (în ideea formǎrii unui stil de programare) şi performanţele (complexitatea) acestora. Din acest motiv, soluţiile de implementare nu sunt totdeauna cele optime.

Prezentarea metodelor, algoritmilor, programelor este însoţitǎ de exemple numerice, menite pe de o parte sǎ clarifice modul de comunicare a datelor şi de interpretare a rezultatelor sau sǎ evidenţeize o anumitǎ comportare specificǎ a soluţiei, iar pe de altǎ parte sǎ faciliteze înţelegerea algoritmilor.

Lucrarea are douǎsprezece capitole, fiecare cuprinzând tema, metoda, exemplu, algoritmul (în pseudocod), programul sursǎ Pascal, programul sursǎ C şi exerciţii.

2

Page 4: Metode Numerice - Seminar

CAPITOLUL 1Rezolvarea numerică a ecuaţiilor

algebrice neliniare

1.1. TemaRezolvarea numerică a ecuaţiilor algebrice neliniare sau transcendente prin

metoda de părţi proporţionale şi metoda tangentei.Fie o funcţie dată. Forma generală a unei ecuaţii algebrice neliniare sau

transcendente este .

1.2. MetodaFie un interval unde ecuaţia are o soluţie. Relaţia de recurenţă în metoda de părţi proporţionale este:

,

în timp ce pentru metoda tangentei, relaţia de recurenţă este:

.

Utilizăm verificarea condiţiei pentru convergenţa şirului.

1.3. Exemplu Să se parcurgă primii trei paşi din metoda lui Newton pentru rezolvarea ecuaţiei

, . Rezolvare.

Cum şi rezultă că , adică ecuaţia dată are soluţie în intervalul [2,3]. Alegem . Cum rezultă

adică soluţia căutată este

1.4. AlgoritmulAlgoritmul asociat metodei de părţi proporţionale este

Intrări: f = funcţia din metodăa, b = capetele intervalului unde ecuaţia are o soluţie

= precizia y = aproximaţia iniţială a soluţiei

Ieşiri: x = soluţia aproximativă

{

3

Page 5: Metode Numerice - Seminar

d + 1k 0

cât timp d {

x

d | x – y |y xk k + 1

} }

Algoritmul asociat metodei tangentei este Intrări: f, fd = funcţia, derivata, din metodă

a, b = capetele intervalului unde ecuaţia are o soluţie = precizia y = aproximaţia iniţială a soluţiei

Ieşiri: x = soluţia aproximativă { d + 1 cât timp d {

x

d | x – y |y x

} }1.5. Programul sursă PascalPROGRAM METODA_DE_PARTI_PROPORTIONALE;VAR A,B,D,EPS,X,Y:REAL; K:INTEGER;FUNCTION F(T:REAL):REAL;BEGIN F:=EXP(T)-3;END;BEGIN WRITE('DATI EPS: '); READLN(EPS); WRITE('DATI A: '); READLN(A); WRITE('DATI B: '); READLN(B); WRITE('DATI X: '); READLN(Y); D:=EPS+1; K:=0; WHILE (D>=EPS) DO

4

Page 6: Metode Numerice - Seminar

BEGIN X:=Y-(F(Y))/((F(A)-F(Y))/(A-Y)); D:=ABS(X-Y); Y:=X; K:=K+1; END; WRITELN; WRITELN('DUPA ',K,' ETAPE ',' X= ',X:12:10); READLN;END.

PROGRAM NEWTON_ECUATIE;VAR A,B,EPS,X,Y,D:REAL;FUNCTION F(T:REAL):REAL;BEGINF:=EXP(T)-2;END;FUNCTION FD(U:REAL):REAL;BEGINFD:=EXP(U);END;BEGINWRITE('A=');READLN(A);WRITE('B=');READLN(B);WRITE('EROARE=');READLN(EPS);WRITE('X=');READLN(Y);D:=EPS+1;WHILE D>=EPS DOBEGINX:=Y-F(Y)/FD(Y);D:=ABS(X-Y);Y:=X;END;WRITELN('SOLUTIA ESTE:',X:15:12);END.

1.6. Programul sursă C// Metoda de parti proportionale#include <stdio.h>#include <conio.h>#include <iostream.h>#include <math.h>#define max 10float funct(float t)

5

Page 7: Metode Numerice - Seminar

{ float f; f=t*t-3*t+1; return f;}void main (void){ float a,b,eps,x,y,d; int k; clrscr(); cout << "dati a " << endl; cin >> a; cout << "dati b " << endl; cin >> b; cout << "dati eroarea " << endl; cin >> eps; cout << "dati x" << endl; cin >> y; d=eps+1; k=0; while (d>=eps) { x=y-(funct(y))/((funct(a)-funct(y))/(a-y)); d=abs(x-y); y=x; k=k+1; } cout << "Dupa " << k << " etape solutia este " << x << endl; getche();}

// Metoda tangentei#include <stdio.h>#include <conio.h>#include <iostream.h>#include <math.h>#define max 10float funct(float t){ float f; f=t*t-3*t+1;

return f;}float funct_der(float t){ float fd; fd=2*t-3;

return fd;

6

Page 8: Metode Numerice - Seminar

}void main (void){

float a,b,eps,x,y,d; int k;

clrscr(); cout << "dati a " << endl; cin >> a; cout << "dati b " << endl; cin >> b;

cout << "dati eroarea " << endl; cin >> eps; cout << "dati x" << endl; cin >> y; d=eps+1; k=0; while (d>=eps) { x=y-(funct(y))/funct_der(y); d=abs(x-y); y=x; k=k+1; } cout << "Dupa " << k << " etape solutia este " << x << endl; getche();}

1.7. Exerciţii1. Pentru ecuaţia să se determine cel puţin trei aproximaţii utilizând metodele: falsei poziţii, înjumătăţirii intervalului.

2. Să se elaboreze programele în limbajele de programare Pascal şi C pentru metodele Aitken, Steffensen, falsei poziţii şi înjumătăţirii intervalului.

7

Page 9: Metode Numerice - Seminar

CAPITOLUL 2Rezolvarea numerică a sistemelor de ecuaţii algebrice liniare prin

metode directe

2.1. TemaRezolvarea numerică a sistemelor algebrice liniare prin metoda de eliminare a

lui Gauss.Un sistem de n ecuaţii liniare cu n necunoscute are forma , unde A este o

matrice pătrată, nesingulară, de dimensiune , iar x şi b sunt vectori coloană de dimensiune n.

2.2. Metoda

Notăm, iniţial, , indicele superior reprezentând etapa. Relaţiile de recurenţă în metoda de eliminare a lui Gauss sunt:

. Sistemul se rezolvă prin metoda substituţiei inverse conform relaţiilor

2.3. Exemplu Să se aplice metoda de eliminare a lui Gauss pentru rezolvarea sistemul

Rezolvare.

Matricea asociată sistemului (în etapa 1) este , iar vectorul

termenilor liberi este .

k = 1

k = 2

8

Page 10: Metode Numerice - Seminar

Soluţiile sunt:

,

,

,

adică soluţia sistemului este

.

2.4. AlgoritmulAlgoritmul asociat metodei Gauss este

Intrări: n = numărul de ecuaţii şi necunoscute ale sistemuluiA = matricea sistemului

b = vectorul termenilor liberiIeşiri: x = vectorul soluţie

{ pentru k 1 : n – 1 { pentru i 1 : n pentru j 1 : n dacă atunci altfel

dacă atunci altfel pentru i 1 : n dacă atunci altfel } pentru i n – 1 : 1 {

9

Page 11: Metode Numerice - Seminar

pentru j i + 1 : n

} }

2.5. Programul sursă PascalPROGRAM GAUSS;VAR A:ARRAY[1..10,1..10,1..10] OF REAL; B:ARRAY[1..10,1..10] OF REAL; X:ARRAY[1..10] OF REAL; S:REAL; N,I,J,K:INTEGER;BEGINWRITE ('N=');READ(N);FOR I:=1 TO N DOFOR J:=1 TO N DOBEGINWRITE('A[',I,',',J,']=');READ(A[I,J,1]);END;FOR I:=1 TO N DOBEGINWRITE('B[',I,']=');READ (B[I,1]);END;FOR K:=1 TO N-1 DOBEGINFOR I:=1 TO N DOFOR J:=1 TO N DOIF(I<=K) THEN A[I,J,K+1]:=A[I,J,K]ELSE IF (J<=K) THEN A[I,J,K+1]:=0ELSE A[I,J,K+1]:=A[I,J,K]-A[I,K,K]*A[K,J,K]/A[K,K,K];FOR I:=1 TO N DOIF (I<=K) THEN B[I,K+1]:=B[I,K]ELSE B[I,K+1]:=B[I,K]-A[I,K,K]*B[K,K]/A[K,K,K];END;X[N]:=B[N,N]/A[N,N,N];FOR I:=N-1 DOWNTO 1 DOBEGINS:=0;FOR J:=I+1 TO N DOS:=S+A[I,J,I]*X[J];X[I]:=(B[I,I]-S)/A[I,I,I];END;FOR I:=1 TO N DOWRITELN(X[I]);

10

Page 12: Metode Numerice - Seminar

READLN;END.

2.6. Programul sursă C// Metoda Gauss de eliminare#include <stdio.h>#include <conio.h>#include <iostream.h>#include <math.h>#define max 15void main (void){

float s; float a[max+1][max+1][max+1],b[max+1][max+1],x[max+1]; int n,i,j,k;

clrscr();cout << "dati dimensiunea matricii" << endl;

cin >> n; cout << "dati matricea sistemului " << endl; for (i=1;i<=n;i++) for (j=1;j<=n;j++) { cout << "a[" << i << j << "]=";

cin >> a[i][j][1]; } cout << endl; cout << "dati vectorul termenilor liberi " << endl; for (i=1;i<=n;i++) { cout << "b[" << i << "]=";

cin >> b[i][1]; } cout << endl; for (k=1;k<=n-1;k++) {

for (i=1;i<=n;i++) for (j=1;j<=n;j++) { if (i<=k)a[i][j][k+1]=a[i][j][k]; else if (j<=k)a[i][j][k+1]=0; else a[i][j][k+1]=a[i][j][k]-a[i][k][k]*a[k][j][k]/a[k][k][k]; } for (i=1;i<=n;i++) if (i<=k)b[i][k+1]=b[i][k]; else b[i][k+1]=b[i][k]-a[i][k][k]*b[k][k]/a[k][k][k]; } x[n]=b[n][n]/a[n][n][n]; for (i=n-1;i>=1;i--)

11

Page 13: Metode Numerice - Seminar

{s=0; for (j=i+1;j<=n;j++) s=s+a[i][j][i]*x[j]; x[i]=(b[i][i]-s)/a[i][i][i]; } cout << "Solutia aproximativa este:" << endl; for (i=1;i<=n;i++) cout << x[i] << ' '; cout << endl; getche();}

2.7. Exerciţii1. Să se aplice metoda Cholesky pentru a rezolva sistemul liniar

.

2. Să se elaboreze programele în limbajele de programare Pascal şi C pentru metoda eliminării Gauss cu pivotare parţială.

CAPITOLUL 3Rezolvarea numerică a sistemelor de ecuaţii algebrice liniare prin

metode iterative

3.1. TemaRezolvarea numerică a sistemelor algebrice liniare prin metodele Jacobi şi

Gauss-Seidel.Un sistem de n ecuaţii liniare cu n necunoscute are forma , unde A este o

matrice pătrată, nesingulară, de dimensiune , iar x şi b sunt vectori coloană de dimensiune n.

3.2. Metoda

Considerăm aproximaţia iniţială un element din R , unde indicele superior reprezintă etapa.

Construim şirul definit prin formulele de recurenţă:

12

Page 14: Metode Numerice - Seminar

în metoda Jacobi:

, ;

în metoda Gauss-Seidel:

, .

Utilizăm verificarea condiţiei pentru convergenţa şirului.

3.3. Exemplua) Să se rezolve, folosind metoda iterativă Jacobi, sistemul

. Rezolvare.

Matricea asociată sistemului (în etapa 1) este iar vectorul

termenilor liberi este .

Şirul de recurenţe pentru soluţie este

Iniţializăm şi alcătuim tabelul

Etapa0 0 0 0

1

2

3

4

Se obţine astfel soluţia aproximativă .

b) Să se rezolve, folosind metoda Gauss-Seidel, sistemul .

Rezolvare.

13

Page 15: Metode Numerice - Seminar

Tabelul soluţiilor este

Etapa0 0 0 0

1

2

3

Se obţine astfel soluţia aproximativă .

3.4. AlgoritmulAlgoritmul asociat metodei Jacobi este

Intrări: n = numărul de ecuaţii şi necunoscute ale sistemuluiA = matricea sistemului

b = vectorul termenilor liberi y = aproximaţia anterioară a soluţiei = precizia

Ieşiri: x = aproximaţia soluţiei { repetă

14

Page 16: Metode Numerice - Seminar

{ pentru i 1 : n { pentru j 1 : n dacă atunci }

pentru i 1 : n pentru i 1 : n } până când ( )

}

Algoritmul asociat metodei Gauss-Seidel este Intrări: n = numărul de ecuaţii şi necunoscute ale sistemului

A = matricea sistemului b = vectorul termenilor liberi y = aproximaţia anterioară a soluţiei = precizia

Ieşiri: x = aproximaţia soluţiei

{ pentru i 1 : n cât timp

{ pentru i 1 : n { pentru j 1 : n dacă atunci }

pentru i 1 : n pentru i 1 : n }

15

Page 17: Metode Numerice - Seminar

3.5. Programul sursă PascalPROGRAM JACOBI_SISTEM;VAR A:ARRAY[1..10,1..10] OF REAL; X,Y,B:ARRAY[1..10] OF REAL; EPS,S,S1:REAL; N,I,J:INTEGER; BEGIN WRITE('DATI N:=');READLN(N); WRITE('DATI EPS:=');READLN(EPS); FOR I:=1 TO N DO FOR J:=1 TO N DO BEGIN WRITE('A[',I,',',J,']='); READLN(A[I,J]); END; FOR I:=1 TO N DO BEGIN WRITE('B[',I,']='); READLN(B[I]); END; FOR I:= 1 TO N DO BEGIN WRITE('Y[',I,']='); READLN(Y[I]); END; REPEAT FOR I:=1 TO N DO BEGIN S1:=0; FOR J:=1 TO N DO IF (I<>J) THEN S1:=S1+A[I,J]*Y[J]; X[I]:=(B[I]-S1)/A[I,I]; END; S:=0; FOR I:=1 TO N DO S:=S+ABS(X[I]-Y[I]); FOR I:=1 TO N DO Y[I]:=X[I]; UNTIL (S<EPS); FOR I:=1 TO N DO WRITELN(X[I]:14:12); READLN; READLN; END.

PROGRAM GAUSS_SEIDEL;VAR A:ARRAY[1..10,1..10] OF REAL;

16

Page 18: Metode Numerice - Seminar

B,X,Y:ARRAY[1..10] OF REAL;EPS,S:REAL;N,I,J:INTEGER;BEGINWRITE('N='); READLN(N);FOR I:=1 TO N DOFOR J:=1 TO N DOBEGINWRITE('A[',I,',',J,']=');READLN(A[I,J]);END;FOR I:=1 TO N DOBEGINWRITE('B[',I,']=');READLN(B[I]);END;FOR I:=1 TO N DOBEGINWRITE('X[',I,']=');READLN(X[I]);END;WRITE('EPS=');READLN(EPS);FOR I:=1 TO N DOY[I]:=X[I];S:=EPS+1;WHILE S>=EPS DOBEGINFOR I:=1 TO N DOBEGINX[I]:=B[I]/A[I,I];FOR J:=1 TO N DOIF J<>I THEN X[I]:=X[I]-(A[I,J]/A[I,I])*X[J];END;S:=0;FOR I:=1 TO N DOS:=S+ABS(X[I]-Y[I]);FOR I:=1 TO N DOY[I]:=X[I];END;FOR I:=1 TO N DOWRITELN(X[I]);END.

3.6. Programul sursă C// Metoda Jacobi pentru rezolvarea sistemelor#include <stdio.h>#include <conio.h>#include <iostream.h>

17

Page 19: Metode Numerice - Seminar

#include <math.h>#define max 10void main (void){

float s,s1,eps; float a[max+1][max+1],b[max+1],x[max+1],y[max+1]; int n,i,j;

clrscr(); cout << "dati numarul de ecuatii si necunoscute " << endl; cin >> n; cout << "dati matricea sistemului " << endl; for (i=1;i<=n;i++) for (j=1;j<=n;j++) { cout << "a[" << i << j << "]=";

cin >> a[i][j]; } cout << endl; cout << "dati vectorul termenilor liberi " << endl; for (i=1;i<=n;i++) { cout << "b[" << i << "]=";

cin >> b[i]; } cout << endl; cout << "dati aproximatia initiala a solutiei " << endl; for (i=1;i<=n;i++) { cout << "x[" << i << "]=";

cin >> y[i]; } cout << endl; cout << "dati eroarea " << endl; cin >> eps; cout << endl; s=eps+1; while (s>=eps) { for (i=1;i<=n;i++) { s1=0; for (j=1;j<=n;j++) if (j!=i)s1=s1+a[i][j]*y[j]; x[i]=(b[i]-s1)/a[i][i]; } s=0; for (i=1;i<=n;i++) s=s+abs(x[i]-y[i]); for (i=1;i<=n;i++) y[i]=x[i]; }

18

Page 20: Metode Numerice - Seminar

cout << "Solutia aproximativa este " << endl; for (i=1;i<=n;i++) cout << x[i] << endl; getche();}

// Metoda Gauss_Seidel pentru rezolvarea sistemelor#include <stdio.h>#include <conio.h>#include <iostream.h>#include <math.h>#define max 10void main (void){

float s,s1,eps; float a[max+1][max+1],b[max+1],x[max+1],y[max+1]; int n,i,j;

clrscr(); cout << "dati numarul de ecuatii si necunoscute " << endl; cin >> n; cout << "dati matricea sistemului " << endl; for (i=1;i<=n;i++) for (j=1;j<=n;j++) { cout << "a[" << i << j << "]="; cin >> a[i][j]; } cout << endl; cout << "dati vectorul termenilor liberi " << endl; for (i=1;i<=n;i++) { cout << "b[" << i << "]="; cin >> b[i]; } cout << endl; cout << "dati aproximatia initiala a solutiei " << endl; for (i=1;i<=n;i++) { cout << "x[" << i << "]="; cin >> x[i]; y[i]=x[i]; } cout << endl; cout << "dati eroarea " << endl; cin >> eps; cout << endl; s=eps+1; while (s>=eps) { for (i=1;i<=n;i++) { s1=0; for (j=1;j<=n;j++)

19

Page 21: Metode Numerice - Seminar

if (j!=i)s1=s1+a[i][j]*x[j]; x[i]=(b[i]-s1)/a[i][i]; } s=0; for (i=1;i<=n;i++) s=s+abs(x[i]-y[i]); for (i=1;i<=n;i++) y[i]=x[i]; } cout << "Solutia aproximativa este " << endl; for (i=1;i<=n;i++) cout << x[i] << endl; getche();}

3.7. ExerciţiiSă se elaboreze programele în limbajele de programare Pascal şi C pentru

metoda relaxării.CAPITOLUL 4

Determinarea rădăcinilor unui polinom

4.1. TemaDeterminarea rădăcinilor unui polinom prin metoda Lobacevski. Presupunem că

rădăcinile ecuaţiei,

sunt distincte în modul şi le notăm cu .

4.2. Metoda

Notăm , Relaţiile de recurenţă sunt ( ):

cu dacă ori .Soluţiile sunt date de:

, .

4.3. Exemplu

Să se aproximeze rădăcinile ecuaţiei

utilizând metoda Lobacevski cu o dispersare (p = 1). Rezolvare.

Cum

20

Page 22: Metode Numerice - Seminar

,folosind relaţiile de recurenţă avem:

Soluţiile sunt date de:

Dacǎ atunci componenta k a soluţiei este altfel , unde este componenta k nenegativǎ iar cea negativǎ, .

4.4. AlgoritmulAlgoritmul asociat metodei Lobacevski este

Intrări: n = gradul polinomuluia = vectorul coeficienţilor polinomului

Ieşiri: x = vectorul soluţie cu componentele nenegative y = vectorul soluţie cu componentele negative { pentru k 1 : n { suma 0 pentru i 1 : n dacă ( ) şi ( ) atunci dacă i este par atunci altfel } pentru { ; } pentru

21

Page 23: Metode Numerice - Seminar

{ pentru ; ; pentru ; dacă atunci altfel } }4.5. Programul sursă PascalPROGRAM LOBACEVSKI;VAR A,B,X,Y,Z:ARRAY[0..100] OF REAL; SUMA,PX,PY:REAL; I,K,N:INTEGER; BEGIN WRITE('DATI N=');READLN(N); FOR I:=0 TO N DO BEGIN WRITE('A[',I,']='); READLN(A[I]); END; B[0] := A[0] *A[0]; FOR K:=1 TO N DO BEGIN B[K]:=A[K]*A[K]; SUMA:=0; FOR I:=1 TO K DO IF (K-I>=0) AND (K+I<=N) THEN IF (I MOD 2=0) THEN SUMA:=SUMA+A[K-I]*A[K+I] ELSE SUMA:=SUMA-A[K-I]*A[K+I]; B[K]:=B[K]+2*SUMA; END; FOR K:=1 TO N DO BEGIN X[K]:=SQRT(B[K]/B[K-1]); Y[K]:=-X[K]; END; FOR K:=1 TO N DO BEGIN PX:=A[0]; FOR I:=1 TO N DO PX:=PX*X[K]+A[I]; PX:=ABS(PX); PY:=A[0]; FOR I:=1 TO N DO PY:=PY*Y[K]+A[I];

22

Page 24: Metode Numerice - Seminar

PY:=ABS(PY); IF (PX<PY) THEN Z[K]:= X[K] ELSE Z[K]:= Y[K] WRITELN(Z[K]:10:10); END; READLN; END.

4.6. Programul sursă C// Metoda Lobacevscki#include <stdio.h>#include <conio.h>#include <iostream.h>#include <math.h>#define max 10void main (void){

float suma,px,py; float a[max+1],b[max+1],x[max+1],y[max+1]; int n,i,k;

clrscr();cout << "dati n " << endl;

cin >> n; for (i=0;i<=n;i++) { cout << "a[" << i << "]="; cin >> a[i]; } cout << endl; b[0]=a[0]*a[0]; for (k=1;k<=n;k++) { b[k]=a[k]*a[k]; suma=0; for (i=1;i<=k;i++) if ((k-i>=0) && (k+i<=n)) if (i % 2==0) suma=suma+a[k-i]*a[k+i]; suma=suma-a[k-i]*a[k+i]; b[k]=b[k]+2*suma; } for (k=1;k<=n;k++) { x[k]=sqrt(b[k]/b[k-1]); y[k]=-x[k]; } for (k=1;k<=n;k++) {

23

Page 25: Metode Numerice - Seminar

px=a[0]; for (i=1;i<=n;i++) px=px*x[k]+a[i]; px=abs(px); py=a[0]; for (i=1;i<=n;i++) py=py*y[k]+a[i]; py=abs(py); if (px<py) cout << x[k] << endl; else cout << y[k] << endl; } getche();}4.7. Exerciţii

Să se aproximeze rădăcinile ecuaţiei

utilizând metoda Lobacevski cu cel puţin două dispersări ( p = 2).

24

Page 26: Metode Numerice - Seminar

CAPITOLUL 5Rezolvarea numerică a problemelor algebrice de valori şi vectori

proprii

5.1. TemaDeterminarea valorilor şi vectorilor proprii ale unei matrici reale şi simetrice

prin metoda Jacobi. O problemă algebrică de valori şi vectori proprii se exprimă prin relaţia:

,unde A este o matrice pătratică de ordin n, x un vector (nenul) care se numeşte vector propriu (la dreapta) iar un număr care se numeşte valoare proprie.

5.2. MetodaFie , , indicele superior reprezentând etapa.Relaţiile de recurenţă sunt:

unde şi cu .

Utilizăm verificarea condiţiei pentru convergenţa şirului.

5.3. Exemplu Să se determine valorile şi vectorii proprii pentru matricea

folosind metoda lui Jacobi. Rezolvare. Luăm, de exemplu,

şi .Atunci:

.

25

Page 27: Metode Numerice - Seminar

Cum rezultă

,

adică .

:

. Astfel,

.

În momentul îndeplinirii condiţiei de oprire în algoritm, valorile proprii sunt elementele de pe diagonala principală şi vectorii proprii sunt coloanele matricei.

5.4. AlgoritmulAlgoritmul asociat metodei Jacobi este

Intrări: n = ordinul matricei sistemuluiA = matricea sistemului

= precizia Ieşiri: B = matricea ce are pe diagonala principală valorile proprii şi ale cărei coloane sunt vectorii proprii { cât timp {

pentru i 1 : n pentru j 1 : n dacă ( ) atunci dacă atunci {

} dacă atunci { }

26

Page 28: Metode Numerice - Seminar

altfel {

}

pentru i 1 : n dacă ( ) şi ( ) atunci {

pentru j 1 : n dacă ( ) şi ( ) atunci {

} } pentru i 1 : n pentru j 1 : n pentru i 1 : n pentru j 1 : n dacă atunci }

}5.5. Programul sursă PascalPROGRAM JACOBI_V_V_P;VAR A,B:ARRAY[1..10,1..10] OF REAL; EPS,S1,MAX,S,C:REAL; I,J,N,P,Q:INTEGER;BEGIN WRITE('DATI N:');READLN(N); FOR I:=1 TO N DO FOR J:=1 TO N DO BEGIN WRITE('A[',I,',',J,']='); READLN(A[I,J]);

27

Page 29: Metode Numerice - Seminar

END; WRITE('EROARE=');READLN(EPS); S1:=EPS+1; WHILE S1>=EPS DO BEGIN MAX:=ABS(A[1,2]); P:=1;Q:=2; FOR I:=1 TO N DO FOR J:=1 TO N DO IF J<>I THEN IF ABS(A[I,J])>MAX THEN BEGIN MAX:=ABS(A[I,J]); P:=I; Q:=J; END; IF A[P,P]=A[Q,Q] THEN BEGIN C:=SQRT(2)/2; S:=C; END ELSE BEGIN C:=COS(ARCTAN(2*A[P,Q]/(A[P,P]- A[Q,Q]))/2); S:=SIN(ARCTAN(2*A[P,Q]/(A[P,P]-A[Q,Q]))/2); END; B[P,P]:=A[P,P]*C*C+A[Q,Q]*S*S+2*A[P,Q]*S*C; B[Q,Q]:=A[P,P]*S*S+A[Q,Q]*C*C-2*A[P,Q]*S*C; B[P,Q]:=(A[Q,Q]-A[P,P])*S*C+(C*C-S*S)*A[P,Q]; B[Q,P]:=B[P,Q]; FOR I:=1 TO N DO IF (I<>P) AND (I<>Q) THEN BEGIN B[I,P]:=A[I,P]*C+A[I,Q]*S; B[P,I]:=B[I,P]; B[I,Q]:=-A[I,P]*S+A[I,Q]*C; B[Q,I]:=B[I,Q]; FOR J:=1 TO N DO IF (J<>P) AND (J<>Q) THEN BEGIN B[I,J]:=A[I,J]; B[J,I]:=B[I,J]; END; END; FOR I:=1 TO N DO FOR J:=1 TO N DO A[I,J]:=B[I,J]; S1:=0;

28

Page 30: Metode Numerice - Seminar

FOR I:=1 TO N DO FOR J:=1 TO N DO IF J<>I THEN S1:=S1+A[I,J]*A[I,J]; END; FOR I:=1 TO N DO WRITELN('VALOAREA PROPRIE',I,' ESTE',A[I,I]:14:12); FOR J:=1 TO N DO BEGIN WRITELN('VECTORUL PROPRIU',J,' ESTE'); FOR I:=1 TO N DO WRITELN(A[I,J]:14:12, ' '); END; READLN; END.

5.6. Programul sursă C// Valori si vectori proprii - Jacobi#include <stdio.h>#include <conio.h>#include <iostream.h>#include <math.h>#define max 10void main (void){

float s,max1,t,eps; float a[max+1][max+1],b[max+1][max+1]; int n,i,j,p,q;

clrscr();cout << "dati dimensiunea matricii"<<endl;cin>> n; cout <<"dati matricea"<< endl;

for (i=1;i<=n;i++) for (j=1;j<=n;j++) { cout << "a[" << i << j << "]="; cin >> a[i][j]; } cout << endl; cout << "dati eroarea in forma: 1.0e-k, unde k este natural fara 0"; cin >> eps; cout << endl; s=eps; while (s>=eps) { for (i=1;i<=n;i++) for (j=1;j<=n;j++) b[i][j]=a[i][j]; max1=abs(b[1][2]); p=1; q=2;

29

Page 31: Metode Numerice - Seminar

for (i=1;i<=n;i++) for (j=1;j<=n;j++) { if (i != j) if (abs(b[i][j])>max1) { max1=abs(b[i][j]); p=i; q=j; } } t=0.5* atan(2*b[p][q]/(b[p][p]-b[q][q])); a[p][p]=b[p][p]*cos(t)*cos(t) + b[q][q]*sin(t)*sin(t) + 2*b[p][q]*sin(t)*cos(t); a[q][q]=b[p][p]*sin(t)*sin(t) + b[q][q]*cos(t)*cos(t) - 2*b[p][q]*sin(t)*cos(t); a[p][q]=0; a[q][p]=0;

for (i=1;i<=n;i++) { if ((i!= p) && (i != q)) { a[i][p]=b[i][p]*cos(t) + b[i][q]*sin(t); a[i][q]=-b[i][p]*sin(t) + b[i][q]*cos(t); a[p][i]=a[i][p]; a[q][i]=a[i][q]; for (j=1;j<=n;j++) { if ((j != p) && (j != q)) { a[i][j]=b[i][j]; a[j][i]=a[i][j]; } } } } s=0; for (i=1;i<=n;i++) for (j=1;j<=n;j++) if (j != i) s = s + a[i][j]*a[i][j];

} cout << "Val. proprii sunt:" << endl; for (i=1;i<=n;i++) cout << a[i][i] << ' '; cout << endl; cout << "Vect. proprii sunt:" << endl; for (j=1;j<=n;j++) { cout << "Vect. propriu " << j << " este ";

30

Page 32: Metode Numerice - Seminar

for (i=1;i<=n;i++) cout << a[i][j] << ' ';

cout << endl; } getche();}

5.7. ExerciţiiSă se determine valorile şi vectorii proprii pentru matricea

.

CAPITOLUL 6Polinomul de interpolare Lagrange

6.1. TemaDeterminarea valorii polinomului Lagrange într-un punct care aproximează o

funcţie pe noduri de interpolare date.

6.2. Metoda Polinomul de interpolare Lagrange pentru funcţia pe nodurile este

.

31

Page 33: Metode Numerice - Seminar

6.3. Exemplu Să se găsească polinomul de interpolare Lagrange pentru funcţia şi nodurile . Rezolvare.

Astfel,

,

de unde

.

6.4. Algoritmul Algoritmul pentru determinarea valorii polinomului Lagrange este Intrări: n = numărul nodurilor de interpolare

a = punctul în care se calculează valoarea polinomului x = vectorul cu nodurile de interpolare f = funcţia din metodă

Ieşiri: = valoarea polinomului Lagrange în a { 0 pentru i 1 : n { p 1 pentru j 1 : n dacă atunci } }

6.5. Programul sursă PascalPROGRAM POLINOM_LAGRANGE;VAR X:ARRAY[1..10] OF REAL; P,L,A:REAL; I,J,N:INTEGER; FUNCTION F(T:REAL):REAL; BEGIN F:=EXP(T); END; BEGIN WRITE('DATI N:');READLN(N); FOR I:=1 TO N DO

32

Page 34: Metode Numerice - Seminar

BEGIN WRITE('X[',I,']='); READ(X[I]); END; WRITE('A=');READLN(A); L:=0; FOR I:=1 TO N DO BEGIN P:=1; FOR J:=1 TO N DO IF (J<>I)THEN P:=P*(A-X[J])/(X[I]-X[J]); L:=L+P*F(X[I]); END; WRITELN('L=',L:14:12); READLN; END.

6.6. Programul sursă C// Polinomul de interpolare Lagrange#include <stdio.h>#include <conio.h>#include <iostream.h>#define max 10float funct(float t){ float f; f=t*t+t+1;

return f;}void main (void){

float l,p,a; float x[max+1]; int n,i,j;

clrscr(); cout << "dati valoarea in care se doreste aproximarea functiei" << endl; cin >> a; cout << "dati numarul de noduri" << endl; cin >> n; cout << "dati nodurile" << endl; for (i=1;i<=n;i++) { cout << "x[" << i << "]=";

cin >> x[i]; } cout << endl; l=0; for (i=1;i<=n;i++) {

33

Page 35: Metode Numerice - Seminar

p=1; for (j=1;j<=n;j++) if (j !=i) p=p*(a-x[j])/(x[i]-x[j]); l=l+funct(x[i])*p; } cout << "Val. pol. in " << a << " este " << l << endl; getche();}

6.7. Exerciţii

Scrieţi polinomul de interpolare Lagrange pentru funcţia pe nodurile

, , , 0, , 1, .

CAPITOLUL 7Diferenţe divizate

7.1. TemaDeterminarea diferenţei divizate pentru o funcţie şi noduri de interpolare date.

7.2. MetodaDiferenţa divizată pentru funcţia pe nodurile este:

.

7.3. Exemplu

Să se determine diferenţa divizată de ordinul 3 pentru funcţia şi nodurile de interpolare 1, 3, 4, 9.

Rezolvare.

de unde rezultă

34

Page 36: Metode Numerice - Seminar

7.4. Algoritmul Algoritmul pentru determinarea diferenţelor divizate este Intrări: k = numărul nodurilor de interpolare x = vectorul cu nodurile de interpolare f = funcţia din metodă

Ieşiri: dd = diferenţa divizată { dd 0 pentru i 1 : k { p 1 pentru j 1 : k dacă atunci ; } } 7.5. Programul sursă PascalPROGRAM DIF;VAR X: ARRAY[1..10] OF REAL; K,I,J:INTEGER; DD,P:REAL;FUNCTION F(T:REAL):REAL;BEGIN F:=SIN(T);END;BEGIN WRITE('DATI K= '); READLN(K); FOR I:=1 TO K DO BEGIN WRITE('DATI X[',I,']= '); READLN(X[I]); END; DD:=0; FOR I:=1 TO K DO BEGIN P:=1; FOR J:=1 TO K DO IF (J<>I) THEN P:=P*(X[I]-X[J]);

35

Page 37: Metode Numerice - Seminar

DD:=DD+F(X[I])/P; END; WRITELN; WRITELN('DIFERENTA DIVIZATA: ',DD:14:12); READLN;END.

7.6. Programul sursă C// Diferenta divizata#include <stdio.h>#include <conio.h>#include <iostream.h>#define max 10float funct(float t){ float f; f=t*t+t+1;

return f;}void main (void){

float dd,p; float x[max+1]; int k,i,j;

clrscr();cout << "dati ordinul diferentei divizate" << endl; cin >> k;

cout << "dati nodurile" << endl; for (i=1;i<=k;i++) { cout << "x[" << i << "]="; cin >> x[i]; } cout << endl; dd=0; for (i=1;i<=k;i++) { p=1; for (j=1;j<=k;j++) if (j !=i) p=p*(x[i]-x[j]); dd=dd+funct(x[i])/p; } cout << "Diferente divizata de ordinul " << k << " este " << dd << endl; getche();}

7.7. ExerciţiiSă se elaboreze programele în limbajele de programare Pascal şi C pentru

determinarea diferenţei divizate prin recurenţă.

36

Page 38: Metode Numerice - Seminar

CAPITOLUL 8Polinomul de interpolare Newton

8.1. TemaDeterminarea valorii polinomului Newton într-un punct care aproximează o

funcţie pe noduri de interpolare date.

8.2. Metoda Polinomul de interpolare Newton pentru funcţia pe nodurile este

.

8.3. Exemplu

Să se scrie polinomul de interpolare Newton pentru funcţia şi nodurile de interpolare 1, 2, 5.

Rezolvare. Polinomul de interpolare Newton este

.

Astfel

de unde rezultă

37

Page 39: Metode Numerice - Seminar

.

8.4. Algoritmul Algoritmul asociat metodei Newton este Intrări: n = numărul nodurilor de interpolare a = numărul în care se calculează valoarea polinomului x = vectorul cu nodurile de interpolare f = funcţia din metodă

Ieşiri: = valoarea polinomului Newton în a { 0 pentru k 2 : n { p1 1 pentru 1 : k – 1 s1 0 pentru i 1 : k { p2 1 pentru j 1 : k dacă atunci } } }

8.5. Programul sursă PascalPROGRAM INTERPOLARE_NEWTON;VAR X:ARRAY[1..20] OF REAL; L1,S1,S2,P1,P2,A:REAL; I,J,K,L,N:INTEGER; FUNCTION F(T:REAL):REAL; BEGIN F:=T-1; END; BEGIN WRITE('DATI N=');READLN(N); WRITE('DATI A=');READLN(A); FOR I:=1 TO N DO

38

Page 40: Metode Numerice - Seminar

BEGIN WRITE('X[',I,']=');READLN(X[I]); END; S2:=0; FOR K:=2 TO N DO BEGIN P1:=1; FOR L:=1 TO K-1 DO P1:=P1*(A-X[L]); S1:=0; FOR I:=1 TO K DO BEGIN P2:=1; FOR J:=1 TO K DO IF J<>I THEN P2:=P2*(X[I]-X[J]); S1:=S1+F(X[I])/P2; END; S2:=S2+S1*P1; END; L1:=F(X[1])+S2; WRITELN('L=',L1:14:12); READLN; END.

8.6. Programul sursă C// Polinomul de interpolare Newton#include <stdio.h>#include <conio.h>#include <iostream.h>#define max 10float funct(float t){ float f; f=t*t+t+1;

return f;}void main (void){

float l,s,p,p1,a; float x[max+1]; int n,i,j,k,m;

clrscr();cout << "dati numarul de noduri" << endl;

cin >> n; cout << "dati nodurile" << endl; for (i=1;i<=n;i++) { cout << "x[" << i << "]=";

cin >> x[i]; }

39

Page 41: Metode Numerice - Seminar

cout << endl; cout << "dati valoarea in care se aproximeaza functia" << endl; cin >> a; l=funct(x[1]);for (k=2;k<=n;k++){

s=0; for (i=1;i<=k;i++) { p=1; for (j=1;j<=k;j++) if (j !=i) p=p*(x[i]-x[j]); s=s+funct(x[i])/p; } p1=1; for (m=1;m<=k-1;m++) p1=p1*(a-x[m]); l=l+s*p1;}

cout << "Val. pol. in " << a << " este " << l << endl; getche();}

8.7. Exerciţii Să se elaboreze programele pentru calcularea valorii polinomului Newton într-un punct pentru o funcţie pe noduri de interpolare date folosind proceduri/funcţii în Pascal, funcţii în C.

40

Page 42: Metode Numerice - Seminar

CAPITOLUL 9Interpolarea prin polinoame trigonometrice

9.1. TemaDeterminarea valorii polinomului trigonometric într-un punct care aproximează

o funcţie pe noduri de interpolare date.

9.2. Metoda Polinomul trigonometric de interpolare pentru funcţia pe nodurile

este

.

9.3. Exemplu

Să se scrie polinomul trigonometric de interpolare pentru funcţia în

nodurile de interpolare 2, 3, 4, 5 şi să se calculeze valoarea polinomului în 4.5. Precizare:

x 2 3 4 5Valoarea aproximativă a lui 23 74 185 388

Rezolvare.

Polinomul trigonometric de interpolare pentru funcţia pe nodurile 2,

3, 4, 5 este

41

Page 43: Metode Numerice - Seminar

de unde rezultă expresia

Valoarea polinomului în 4.5 este

adică, .

9.4. AlgoritmulAlgoritmul pentru determinarea valorii polinomului trigonometric de

interpolare este Intrări: 2n + 1 = numărul nodurilor de interpolare a = numărul în care se calculează valoarea polinomului x = vectorul cu nodurile de interpolare f = funcţia din metodă

Ieşiri: t = valoarea polinomului trigonometric de interpolare în a

{ t 0 pentru i 0 : 2n

42

Page 44: Metode Numerice - Seminar

{ p 1 pentru j 0 : 2n dacă atunci } }

9.5. Programul sursă PascalPROGRAM POLINOMUL_TRIGONOMETRIC_DE_INTERPOLARE;VAR X:ARRAY[0..20] OF REAL; P,T,A:REAL; I,J,N:INTEGER; FUNCTION F(U:REAL):REAL; BEGIN F:=EXP(U); END; BEGIN WRITE('DATI N:');READLN(N); FOR I:=0 TO 2*N DO BEGIN WRITE('X[',I,']='); READ(X[I]); END; WRITE('A=');READLN(A); T:=0; FOR I:=0 TO 2*N DO BEGIN P:=1; FOR J:=0 TO 2*N DO IF (J<>I)THEN P:=P*SIN((A-X[J])/2)/SIN((X[I]-X[J])/2); T:=T+P*F(X[I]); END; WRITELN('T=',T:14:12); READLN; END.

9.6. Programul sursă C// Polinomul trigonometric de interpolare#include <stdio.h>#include <conio.h>#include <iostream.h>#include <math.h>#define max 10float funct(float u){

43

Page 45: Metode Numerice - Seminar

float f; f=u*u+u+1;

return f;}void main (void){

float t,p,a; float x[max+1]; int n,m,i,j;

clrscr(); cout << "dati valoarea in care se doreste aproximarea functiei" << endl; cin >> a;

cout << "dati numarul de noduri (impar)" << endl; cin >> m; n=(m-1)/2; cout << "n= " << n << endl; cout << "dati nodurile" << endl; for (i=0;i<=2*n;i++) { cout << "x[" << i << "]=";

cin >> x[i]; } cout << endl; t=0; for (i=0;i<=2*n;i++) { p=1; for (j=0;j<=2*n;j++) if (j !=i) p=p*sin((a-x[j])/2)/sin((x[i]-x[j])/2); t=t+funct(x[i])*p; } cout << "Val. pol. in " << a << " este " << t << endl; getche();}

9.7. ExerciţiiSă se testeze programele în limbajele de programare Pascal şi C pentru:

a) funcţia în nodurile: 1; 2; 3; 4; 5;

b) funcţia în nodurile: 0; 0.1; 0.2; 0.3; 0.4; 0.5.

44

Page 46: Metode Numerice - Seminar

CAPITOLUL 10Diferenţe finite

10.1. TemaDeterminarea diferenţei finite pentru o funcţie şi noduri de interpolare date.

10.2. MetodaDiferenţa finită pentru funcţia pe nodurile echidistante (

) este

.

10.3. Exemplu Să se determine diferenţa finită de ordinul 3 pentru funcţia şi nodurile de interpolare 1, 2, 3, 4 în primul nod. Rezolvare.

de unde .

10.4. Algoritmul Algoritmul pentru determinarea diferenţelor finite este Intrări: n = numărul nodurilor de interpolare k = ordinul diferenţei finite x = vectorul cu nodurile de interpolare h = pasul

Ieşiri: df = diferenţa finită { dacă k este par atunci df altfel df pentru j 1 : k { c 1 pentru 1 : j dacă este par atunci df altfel df }

} 10.5. Programul sursă Pascal

45

Page 47: Metode Numerice - Seminar

PROGRAM DIF_FINITA; FUNCTION F(T:REAL):REAL; BEGIN F:=EXP(-(T*T)); END;VAR X:ARRAY[1..20] OF REAL; H,C,D,DF:REAL; L,I,J,K,N:INTEGER; BEGIN WRITE('DATI N=');READLN(N); WRITE('DATI H=');READLN(H); WRITE('DATI X[1]=');READLN(X[1]); FOR I:=2 TO N DO X[I]:=X[1]+(I-1)*H; WRITE('K=');READLN(K); IF (K MOD 2=0) THEN DF:=F(X[1]) ELSE DF:=-F(X[1]); FOR J:=1 TO K DO BEGIN C:=1; FOR L:=1 TO J DO C:=C*(K-L+1)/L; IF (K-J) MOD 2=0 THEN DF:=DF+C*F(X[1]+J*H) ELSE DF:=DF-C*F(X[1]+J*H);

END; WRITE('DF=',DF:14:12); READLN; END.

10.6. Programul sursă C// Diferente finite#include <stdio.h>#include <conio.h>#include <iostream.h>#include <math.h>#define max 10float funct(float t){ float f; f=exp(-(t*t));

return f;}void main (void){

float df,c,h; float x[max+1]; int n,k,i,j,l;

clrscr(); cout << "dati ordinul diferentei finite" << endl;

46

Page 48: Metode Numerice - Seminar

cin >> k; cout << "dati numarul de noduri" << endl; cin >> n; cout << "dati pasul" << endl; cin >> h; cout << "dati primul nod" << endl; cin >> x[1]; for (i=2;i<=n;i++) { x[i]=x[1]+(i-1)*h; } cout << endl;

if (k%2==0)df=funct(x[1]); else df=-funct(x[1]); for (j=1;i<=k;j++) { c=1; for (l=1;l<=j;l++) c=c*(k-l+1)/l; if ((k-j)%2==0)df=df+c*funct(x[1]+j*h); else df=df-c*funct(x[1]+j*h); } cout << "Diferenta finita de ordinul " << k << " este " << df << endl; getche();}

10.7. ExerciţiiSă se scrie programele Pascal şi C pentru calculul diferenţei finite utilizând

funcţii în limbajele de programare amintite.

CAPITOLUL 11Formule de interpolare pe noduri echidistante

47

Page 49: Metode Numerice - Seminar

11.1. TemaDeterminarea valorii polinomului Newton pe noduri echidistante în „primul” şi

respectiv „ultimul” nod.

11.2. MetodaFie funcţia şi nodurile de interpolare echidistante (

). Formula lui Newton progresivă este:

în timp ce formula lui Newton regresivă este:

.

11.3. Exemplu

Să se scrie polinomul Newton (varianta progresivă şi regresivă) pentru funcţia şi nodurile 1, 3, 5, 7.

Rezolvare. Polinomul lui Newton progresiv este:

de unde se obţine

în timp ce polinomul lui Newton regresiv este:

de unde

48

Page 50: Metode Numerice - Seminar

11.4. AlgoritmulAlgoritmul pentru determinarea valorii polinomului Newton pe noduri

echidistante în „primul” nod este Intrări: n = numărul nodurilor de interpolare x = vectorul cu nodurile de interpolare h = pasul a = numărul în care se calculează valoarea polinomului

Ieşiri: = valoarea polinomului în a { pentru k 1 : n – 1 { p 1 pentru 1 : k pentru 1 : k dacă k este par atunci altfel pentru j 1 : k { pentru dacă este par atunci altfel } }

} Algoritmul pentru determinarea valorii polinomului Newton pe noduri echidistante în „ultimul” nod este Intrări: n = numărul nodurilor de interpolare x = vectorul cu nodurile de interpolare h = pasul a = numărul în care se calculează valoarea polinomului

Ieşiri: = valoarea polinomului în a { pentru k 1 : n – 1 { p 1

49

Page 51: Metode Numerice - Seminar

pentru 1 : k pentru 1 : k dacă k este par atunci altfel pentru j 1 : k { pentru dacă este par atunci altfel } }

}

11.5. Programul sursă PascalPROGRAM NEWTON_PROGRESIV;VAR X:ARRAY[1..10] OF REAL;H,A,L1,P,FC,S1,C:REAL;N,J,K,L,S:INTEGER;FUNCTION F(T:REAL):REAL;BEGINF:=T+2;END;BEGINWRITE('N=');READLN(N);WRITE('X=');READLN(X[1]);WRITE('H=');READLN(H);WRITE('A=');READLN(A);FOR J:=1 TO N DOX[J]:=X[1]+(J-1)*H;L1:=F(X[1]);FOR K:=1 TO N-1 DOBEGINP:=1;FOR L:=1 TO K DOP:=P*(A-L+1);FC:=1;FOR L:=1 TO K DOFC:=FC*L;

50

Page 52: Metode Numerice - Seminar

IF K MOD 2=0 THENS1:=F(X[1])ELSES1:=-F(X[1]);FOR J:=1 TO K DOBEGINC:=1;FOR S:=1 TO J DOC:=C*(K-S+1)/S;IF (K-J) MOD 2=0 THENS1:=S1+C*F(X[1]+J*H)ELSES1:=S1-C*F(X[1]+J*H);END;L1:=L1+S1*P/FC;END;WRITELN('L1',L1:14:12);END.

PROGRAM NEWTON_REGRESIV;VAR X:ARRAY[1..10] OF REAL;H,A,LN,P,FC,S1,C:REAL;N,J,K,L,S:INTEGER;FUNCTION F(T:REAL):REAL;BEGINF:=T-2;END;BEGINWRITE('N=');READLN(N);WRITE('X=');READLN(X[1]);WRITE('H=');READLN(H);WRITE('A=');READLN(A);FOR J:=1 TO N DOX[J]:=X[1]+(J-1)*H;LN:=F(X[N]);FOR K:=1 TO N-1 DOBEGINP:=1;FOR L:=1 TO K DOP:=P*(A+L-1);FC:=1;FOR L:=1 TO K DOFC:=FC*L;IF K MOD 2=0 THENS1:=F(X[N-K])ELSES1:=-F(X[N-K]);FOR J:=1 TO K DO

51

Page 53: Metode Numerice - Seminar

BEGINC:=1;FOR S:=1 TO J DOC:=C*(K-S+1)/S;IF (K-J) MOD 2=0 THENS1:=S1+C*F(X[N-K]+J*H)ELSES1:=S1-C*F(X[N-K]+J*H);END;LN:=LN+S1*P/FC;END;WRITELN('LN',LN:14:12);END.

11.6. Programul sursă C// Polinomul de interpolare Newton-progresiv#include <stdio.h>#include <conio.h>#include <iostream.h>#include <math.h>#define max 10float funct(float t){ float f; f=exp(-(t*t));

return f;}void main (void){

float np,df,c,h,p,fc,a; float x[max+1]; int n,k,i,j,l;

clrscr(); cout << "dati valoarea unde se aproximeaza functia" << endl; cin >> a; cout << "dati numarul de noduri" << endl; cin >> n; cout << "dati pasul" << endl; cin >> h; cout << "dati primul nod" << endl; cin >> x[1]; for (i=2;i<=n;i++) { x[i]=x[1]+(i-1)*h; } cout << endl; np=funct(x[1]); for (k=1;k<=n-1;k++)

52

Page 54: Metode Numerice - Seminar

{ fc=1; for (l=1;l<=k;l++) fc=fc*l; p=1; for (l=1;l<=k;l++) p=p*(a-l+1); if (k%2==0)df=funct(x[1]); else df=-funct(x[1]); for (j=1;i<=k;j++) { c=1; for (l=1;l<=k;l++) c=c*(k-l+1)/l; if ((k-j)%2==0)df=df+c*funct(x[1]+j*h); else df=df-c*funct(x[1]+j*h); } np=np+p*df/fc; } cout << "Valoarea polinomului in " << a << " este " << np << endl; getche();}

// Polinomul de interpolare Newton-regresiv#include <stdio.h>#include <conio.h>#include <iostream.h>#include <math.h>#define max 10float funct(float t){ float f; f=exp(-(t*t));

return f;}void main (void){

float np,df,c,h,p,fc,a; float x[max+1]; int n,k,i,j,l;

clrscr(); cout << "dati valoarea unde se aproximeaza functia" << endl; cin >> a; cout << "dati numarul de noduri" << endl; cin >> n; cout << "dati pasul" << endl; cin >> h;

53

Page 55: Metode Numerice - Seminar

cout << "dati primul nod" << endl; cin >> x[1]; for (i=2;i<=n;i++) { x[i]=x[1]+(i-1)*h; } cout << endl; np=funct(x[n]); for (k=1;k<=n-1;k++) { fc=1; for (l=1;l<=k;l++) fc=fc*l; p=1; for (l=1;l<=k;l++) p=p*(a+l-1); if (k%2==0)df=funct(x[n-k]); else df=-funct(x[n-k]); for (j=1;i<=k;j++) { c=1; for (l=1;l<=k;l++) c=c*(k-l+1)/l; if ((k-j)%2==0)df=df+c*funct(x[n-k]+j*h); else df=df-c*funct(x[n-k]+j*h); } np=np+p*df/fc; } cout << "Valoarea polinomului in " << a << " este " << np << endl; getche();}

11.7. ExerciţiiSă se scrie programele Pascal şi C pentru calculul valorii polinoamelor Newton

progresiv şi regresiv într-un punct utilizând funcţii în limbajele de programare amintite.CAPITOLUL 12

Formule de derivare numerică pe noduri echidistante

12.1. TemaDeterminarea derivatei întâi a unei funcţii pe noduri de interpolare date derivând

formula lui Newton progresivă şi regresivă.

12.2. MetodaDerivata întâi pentru funcţia pe nodurile echidistante (

) este:

(varianta progresivă),

54

Page 56: Metode Numerice - Seminar

(varianta regresivă).

12.3. Exemplu

Să se calculeze derivata polinomului Newton în primul şi ultimul nod, polinom ce aproximează funcţia pe nodurile , , 0 , 1, 2. Rezolvare. Polinomul Newton progresiv este:

de unde se obţine

în timp ce polinomul lui Newton regresiv este:

de unde, prin înlocuire, se obţine

12.4. AlgoritmulAlgoritmul pentru aproximarea derivatei întâi în „primul” nod este

Intrări: n = numărul nodurilor de interpolare x = vectorul cu nodurile de interpolare h = pasul f = funcţia din metodă

Ieşiri: d1 = derivata (aproximativă) în primul nod { d1 0 pentru k 1 : n – 1 { dacă k este par atunci df altfel

55

Page 57: Metode Numerice - Seminar

df pentru j 1 : k { pentru dacă este par atunci altfel } dacă ( ) este par atunci altfel } }

Algoritmul pentru aproximarea derivatei întâi în „ultimul” nod este Intrări: n = numărul nodurilor de interpolare x = vectorul cu nodurile de interpolare h = pasul f = funcţia din metodă

Ieşiri: dn = derivata (aproximativă) în ultimul nod { dn 0 pentru k 1 : n – 1 { dacă k este par atunci d altfel d pentru j 1 : k { pentru dacă este par atunci altfel }

} }

12.5. Programul sursă PascalPROGRAM DERIVATA_NEWTON_PROGRESIV ; VAR X:ARRAY[1..100] OF REAL; H,D1,DF,C:REAL; I,J,K,L,N:INTEGER;

56

Page 58: Metode Numerice - Seminar

FUNCTION F(T:REAL):REAL; BEGIN F:=T*T+1; END; BEGIN WRITE('N=');READLN(N); WRITE('X1=');READLN(X[1]); WRITE('H=');READLN(H); FOR I:=1 TO N DO X[I]:=X[1]+(I-1)*H; D1:=0; FOR K:=1 TO N-1 DO BEGIN IF K MOD 2=0 THEN DF:=F(X[1]) ELSE DF:=-F(X[1]); FOR J:=1 TO K DO BEGIN C:=1; FOR L:=1 TO J DO C:=C*(K-L+1)/L; IF (K-J) MOD 2=0 THEN DF:=DF+C*F(X[1]+J*H) ELSE DF:=DF-C*F(X[1]+J*H); END; IF (K-1) MOD 2=0 THEN D1:=D1+DF/K ELSE D1:=D1-DF/K; END; D1:=D1/H; WRITELN('D1=',D1:14:12); READLN; END.

PROGRAM DERIV_NEWTON_REGRESIV;VAR X:ARRAY[1..10] OF REAL;H,DN,D,C:REAL;N,J,K,L:INTEGER;FUNCTION F (T:REAL):REAL;BEGINF:= T*T+1;END;BEGINWRITE('N='); READLN(N);WRITE('H='); READLN(H);WRITE ('X[1]='); READLN (X[1]);FOR J:= 1 TO N DOX[J]:= X[1]+ (J-1)*H;DN:= 0;FOR K:= 1 TO N-1 DO

57

Page 59: Metode Numerice - Seminar

BEGINIF K MOD 2 = 0 THEND:= F(X[N-K])ELSED:= -F(X[N-K]);FOR J:= 1 TO K DOBEGINC:= 1;FOR L:= 1 TO J DOC:= C*(K-L+1)/L;IF (K-J) MOD 2= 0 THEND:= D + C*F(X[N-K]+J*H)ELSED:= D-C*F(X[N-K]+J*H);END;DN:=DN+D/K;END;DN:=DN/H;WRITELN('DERIVATA ESTE:',DN:14:12);END.

12.6. Programul sursă C// Derivata polinomul Newton-progresiv (in primul nod)#include <stdio.h>#include <conio.h>#include <iostream.h>#include <math.h>#define max 10float funct(float t){ float f; f=exp(-(t*t));

return f;}void main (void){

float dnp,df,c,h; float x[max+1]; int n,k,i,j,l;

clrscr(); cout << "dati numarul de noduri" << endl; cin >> n; cout << "dati pasul" << endl; cin >> h; cout << "dati primul nod" << endl; cin >> x[1]; for (i=2;i<=n;i++) {

58

Page 60: Metode Numerice - Seminar

x[i]=x[1]+(i-1)*h; } cout << endl; dnp=0; for (k=1;k<=n-1;k++) { if (k%2==0)df=funct(x[1]); else df=-funct(x[1]); for (j=1;i<=k;j++) { c=1; for (l=1;l<=j;l++) c=c*(k-l+1)/l; if ((k-j)%2==0)df=df+c*funct(x[1]+j*h); else df=df-c*funct(x[1]+j*h); } if ((k-1)%2==0)dnp=dnp+df/k; else dnp=dnp-df/k; } dnp=dnp/h; cout << "Derivata polinomului in primul nod este " << dnp << endl; getche();}

// Derivata polinomul Newton-regresiv (in ultimul nod)#include <stdio.h>#include <conio.h>#include <iostream.h>#include <math.h>#define max 10float funct(float t){ float f; f=exp(-(t*t));

return f;}void main (void){

float dnp,df,c,h; float x[max+1]; int n,k,i,j,l;

clrscr(); cout << "dati numarul de noduri" << endl; cin >> n; cout << "dati pasul" << endl; cin >> h; cout << "dati primul nod" << endl; cin >> x[1];

59

Page 61: Metode Numerice - Seminar

for (i=2;i<=n;i++) { x[i]=x[1]+(i-1)*h; } cout << endl; dnp=0; for (k=1;k<=n-1;k++) { if (k%2==0)df=funct(x[n-k]); else df=-funct(x[n-k]); for (j=1;i<=k;j++) { c=1; for (l=1;l<=j;l++) c=c*(k-l+1)/l; if ((k-j)%2==0)df=df+c*funct(x[n-k]+j*h); else df=df-c*funct(x[n-k]+j*h); } dnp=dnp+df/k; } dnp=dnp/h; cout << "Derivata polinomului in ultimul nod este " << dnp << endl; getche();}

12.7. ExerciţiiSă se scrie programele Pascal şi C pentru calculul derivatei utilizând funcţii în

limbajele de programare specificate.

BIBLIOGRAFIE

[1] Bucur, C.M., Metode numerice, Ed. Facla, Timişoara, 1973.[2] Coman, G., Analiză numerică, Ed. Libris, Cluj, 1995.[3] Demidovici, B.P., Maron,I., Elements de calcul numerique, Ed. Mir de Mosou,

1973.[4] Dodescu, Gh., Toma, M., Metode de calcul numeric, E. D. P., Bucureşti, 1976.[5] Dodescu, Gh., Metode numerice în algebră, Ed. tehnică, Bucureşti, 1979.[6] Ichim, I., Marinescu, G., Metode de aproximare numerică, Ed. Academiei R.

S. R., Bucureşti, 1986.[7] Ignat, C., Ilioi, C., Jucan, T., Elemente de informatică şi calcul numeric,

Universitatea „Al. I. Cuza”, Iaşi, Facultatea de Matematică, 1989.[8] Juan Antonio Infante del Rio, Jose Maria Rey Cabezas, Metodos

Numericas, Teoria, problemas y practicas con MATLAB, Ed. Piramide, 2002.[9] Knut, D.E., Sortare şi căutare, vol. 3, Tratat de programarea

calculatoarelor, Ed. Tehnică, Bucureşti, 1976.

60

Page 62: Metode Numerice - Seminar

[10] Lucanu, D., Bazele proiectǎrii programelor şi algoritmilor, volume 1,2,3. Curs multiplicat la Editura Universitǎţii „Al. I. Cuza”, Iaşi, România, 1996.

[11] Mihu, C., Metode numerice în algebra liniară, Ed. Tehnică, Bucureşti, 1977.

[12] Press, W.H., Teuklosky, S.A., Vetterling, W.T., Flannery, B.P., Numerical Recipes in C: The Art of scientific Computing, Second Edition (Cambridge University Press, Cambridge, 1992).

[13] Scheiber, E., Metode numerice, Universitatea Transilvania din Braşov, Facultatea de Matematică – Informatică (electronic).

[14] Talmaciu, M., Metode numerice Volumul I, Editura Tehnica Info, Chişinǎu, 2005.

[15] Toma, M., Odăgescu, I., Metode numerice şi subrutine, Ed. Tehnică, Bucureşti, 1980.

[16] Vraciu, G., Popa, A., Metode numerice cu aplicaţii în tehnica de calcul, Scrisul românesc, Craiova, 1982.

[17] ***, Borland International, Inc. Borland C++, Programming Guide (Borland International, Scotts wally, CA, 1992).

61