Date post: | 06-Feb-2016 |
Category: |
Documents |
Upload: | moisa-iulian |
View: | 109 times |
Download: | 11 times |
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
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
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
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
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
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
{ 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
}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
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
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
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
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
{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
î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
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
{ 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
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
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
#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
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
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
,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
{ 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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
.
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
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
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
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
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
{ 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
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
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
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
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
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
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
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
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
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
{ 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
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
(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
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
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
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
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
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
[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