2. CATEVA METODE NUMERICE CLASICE

 Metoda lui Gauss este una din cele mai populare metode în rezolvarea sistemelor de ecuatii liniare, calculul determinantilor si determinarea inversei matricilor. Aceasta metoda foloseste algoritmul eliminarii succesive a necunoscutelor, transformând matricea sistemului initial într-o matrice superior triunghiulara.

 2.1.1.1. Rezolvarea sistemelor de ecuatii liniare si calculul determinantilor

Fie un sistem de ecuatii de forma:

 sau scris matriceal sub forma:

 unde:

 Folosind metoda eliminarii a lui Gauss urmarim sa aducem sistemul de ecuatii (2.1) la forma echivalenta:

 unde:

 Determinarea elementelor matricilor A' si se face prin eliminarea succesiva a necunoscutelor sistemului initial, cu ajutorul relatiilor:

 respectiv:

 Solutiile sistemului initial se obtin pe calea substitutiei inverse, cu relatiile:

 Valoarea determinantului unei matrici superior triunghiulare de forma matricii A' se obtine simplu cu ajutorul relatiei:

 Pentru a putea folosi aceasta metoda este necesar ca elementele generatoare () sa nu fie nule. In plus, pot aparea erori numerice daca elementele generatoare au valori prea mici. Pentru a evita aceasta situatie se pot inversa ecuatiile sistemului între ele, fapt care duce însa la schimbarea semnului determinantului. Operatia de permutare a liniilor se numeste pivotare. Linia de pivotare este linia l corespunzatoare celui mai mare element din coloana k, adica elementul pivot va fi a_lk=max(a_ik), .

Daca eliminarea Gaussiana este folosita pentru determinarea inversei unei matrici, atunci este economicos ca detaliile transformarii matricii A, date de o relatie de tipul (2.5) sa fie stocate sub diagonala, astfel încât la eliminarea kpentru linia i vom avea:

iar ulterior pentru coloana j:

In plus, daca au aparut permutari de linii, permutarea trebuie memorata într-un vector întreg:

care la permutarea liniilor l si k devine:

 2.1.1.2. Inversa unei matrici

 Fie matricea patrata:

 Pentru a gasi inversa acestei matrici se utilizeaza relatia:

unde I este matricea unitate.

Daca matricea inversa A^-1 este:

 atunci pe baza relatiei (2.8) obtinem:

 unde d _ij este simbolul lui Kronecker.

Cele n sisteme de ecuatii obtinute pentru au aceeasi matrice A. Numai termenii constanti sunt diferiti si problema pate fi rezolvata prin metoda lui Gauss pentru a se gasi cele n² elemente ale matricii A^-1.

  Sa presupunem ca sistemul de ecuatii care trebuie rezolvat este scris sub forma:

 unde matricile A, si sunt date de (2.2).

Dupa cum este cunoscut, cu ajutorul metodei Gauss, solutiile sistemului de ecuatii se obtin simplu daca matricea sistemului este adusa la o forma triunghiulara. Pornind de la aceasta observatie vom scrie matricea sistemului sub forma unui produs de doua matrici L=[l_ij] si U=[u_ij], prima inferior triunghiulara, iar cea de-a doua superior triunghiulara cu toate elementele de pe prima diagonala egale cu unitatea:

A=LU (2.11)

unde:

 Sistemul initial se scrie acum sub forma:

 Daca notam , atunci rezolvarea sistemului precedent se reduce la rezolvarea urmatoarelor doua sisteme:

 cu avantajul ca matricile acestor sisteme sunt matrici triunghiulare.

Mai întâi trebuie rezolvat sistemul (2.14) pentru a se determina solutiile, iar apoi sistemul (2.13) pentru a se obtine solutiile finale . Desigur ca în prealabil trebuie determinate elementele matricilor L si U. In acest scop, din egalitatea (2.11) se obtin urmatoarele relatii de recurenta:

 Cunoscându-se elementele matricii L, cu ajutorul ecuatiei (2.14) se obtin solutiile y_i:

 Cunoscând elementele y_i si cu ajutorul ecuatiei (2.13) se obtin în final solutiile sistemului initial:

 Daca matricea A este simetrica, atunci:

 Matricea A va fi folosita pentru memorarea matricilor L si U, iar vectorulpentru memorarea vectorilor sau , dupa necesitati.

  Fie sistemul liniar simetric:

 Urmarim sa factorizam matricea A sub forma particulara pentru matrici simetrice:

A=T^TT

unde:

 Solutiile sistemului initial se pot obtine prin rezolvarea a doua sisteme de ecuatii a caror matrici sunt T si respectiv T^T:

 si ale caror elemente se determina cu relatiile:

 Construirea solutiilor sistemului se face pe baza relatiilor:

 Pentru aceasta metoda, matricea A este stocata sub forma unui vector de dimensiune n(n+1)/2 care poate fi utilizat si pentru stocarea elementelor t_ij, vectorii si apoi putând fi stocati în vectorul . Se obtine astfel o economie maxima de memorie la rezolvarea sistemelor liniare simetrice.

  Fie un set de valori y_i coprespunzatoare unui set de puncte discrete x_i aflate într-un interval [a,b]. Nefiind cunoscuta forma analitica a functiei y=f(x) se pune problema determinarii unei functii, numita functie de interpolare, ale carei valori în punctele x_i (punctele de interpolare) sa fie tocmai valorile corespunzatoare y_i. Din punct de vedere geometric problema consta în gasirea curbei y=F(x) care sa treaca neaparat prin punctele (x_i,y_i) date.

Daca sunt date n+1 puncte de interpolare (x₀, x₁,...,x_n), atunci functia de interpolare trebuie cautata sub forma unui polinom de grad n (polinomul de interpolare Lagrange), L_n(x), cu conditia:

L_n(x_i)=y_i i=0,1,...,n (2.25)

In scopul determinarii polinomului Lagrange se construieste un polinom p_i(x), cu conditia:

 Deoarece acest polinom se anuleaza în n puncte el va avea forma:

 Cum p_i(x_i)=1, se obtine:

 si deci:

 Polinomul Lagrange care satisface conditia (2.25) este:

 cu conditia:

L_n(x_j)=y_j

Polinomul de interpolare Lagrange poate fi folosit în locul unei functii tabelate y(x_i) în scopul evaluarii numerice a integralei functiei respective. Daca distanta dintre punctele de interpolare x_i este constanta, atunci pasul de integrare va fi h=x_i-x_i-1, iar numarul de intervale pâna la un punct oarecare x este: q=(x-x₀)/h.
Cu aceste notatii, polinomul Lagrange poate fi scris:

 si tinând cont ca x₂-x₀=2h se obtine formula lui Simpson pentru acest caz simplu (n=2):

 Fie acum un numar par de puncte (n=2m) echidistante. Aplicând formula lui Simpson fiecarui subinterval dublu (x₀,x₂), (x₂,x₄),...,(x_2m-2,x_2m), se obtine formula lui Simpson generalizata:

 unde:

Integralele improprii contin functii care [31]:

Una dintre solutiile cele mai simple consta în schimbarea de variabila t=1/x, adica:

 Pentru integralele care au o singularitate de tip putere la limita inferioara,

adica (x - a), 0 1, pentru x=0, putem folosi:

 Pentru singularitatile frecvente de tip inversul radacinii patrate folosim:

Pentru integralele improprii convergente, fara singularitati, putem formula un procedeu general bazat pe utilizarea preciziei sistemului de calcul. Deoarece integrala este o suma, ultima valoare calculata a functiei care nu mai modifica valoarea sumei este data de variabila REPS din capitolul I. Folosind metoda lui Newton pentru ecuatia neliniara f(x)-REPS=0 aflam valoarea abscisei pentru care putem închia integrarea. Apoi, aplicam formula lui Simpson pentru calculul integralei între limitele de integrare care au astfel valori bine determinate. Numarul de pasi poate fi eventual ajustat tinând sema de eroarea calculata cu formula (2.35).

Valoarea integralelor unidimensionale, poate fi estimata cu ajutorul relatiei (Fig.2.1.):

 Daca alegem b-a=1, adica daca se transforma intervalul [a,b] în intervalul [0,1], atunci o evaluare a lui f(e ) se poate obtine dupa un numar mare de încercari n:

 Generarea valorilor x se poate face folosind generatorul intern de numere aleatoare caracteristic limbajului de programare folosit. Utilizarea numerelor aleatoare este justificata de faptul ca ordinea termenilor din suma (2.48) este indiferenta. Metoda Monte Carlo este utilizata în special pentru calculul comod al integralelor multiple.

Fie cazul test [6]:

 Daca generam C numere aleatoare, doar n dintre ele se vor afla în intervalul de definitie, astfel încât:

 Asadar vom avea:

 Metoda Monte-Carlo asigura obtinerea unor valori aproximative ale integralelor multiple, foarte apropiate de cele exacte dupa un numar mare de evaluari.

In cazul în care intervalul de definitie al variabilelor de integrare nu este [0,1] se face transformarea de variabila x_iy_i:

x_i=x_iminim+(x_imaxim-x_iminim)*y_i pentru 0y_i1.

 în care:

 este Jacobianul transformarii.

Daca notam:

 atunci avem:

 Dispersia si eroarea patratica se determina conform formulelor din Cap.3.

Desi, în principiu, derivarea poate fi efectuata analitic, în cazul functiilor complicate exista riscul unor greseli inerente si în plus în unele programe este necesar un calcul numeric al derivatelor.

Una dintre metodele folosite în derivarea numerica porneste de la folosirea polinomului Lagrange care poate sa treaca prin mai multe puncte echidistante. Fie polinomul de interpolare care trece prin trei puncte echidistante:

 Derivata acestui polinom este:

 unde:

x_j=x₀+jh, j=1,2,...,n (2.48)

Tinând cont de ultima relatie, derivatele de ordin 1 si 2 în punctele x₀,x₁ si respectiv x₂ sunt:

 In cazul unui polinom de interpolare care trece prin cinci puncte:

 Evident ca în acest caz, pentru primele si ultimele doua puncte, derivata trebuie calculata cu formulele anterioare (2.56).

O solutie mult mai expeditiva este data de formula de definitie a derivatei partiale:

2.4. METODA LUI NEWTON PENTRU REZOLVAREA ECUATIILOR SI SISTEMELOR DE ECUATII NELINIARE

  Pentru ecuatiile neliniare, metoda lui Newton foloseste procesul iterativ:

 dedus din dezvoltarea în serie Taylor a functiei f(x).

Fie cazul test: f(x)=x⁴+2x³-x-1 cu radacina 0.867.

Metoda descrisa mai sus pentru rezolvarea unei ecuatii neliniare poate fi generalizata pentru un sistem de ecuatii neliniare. Fie pentru început un sistem de doua ecuatii neliniare cu doua necunoscute:

f(x,y)=x³-3xy²-2x+2=0

g(x,y)=3x²y-y³-2y=0 (2.52)

cu solutiile: x=0.884646 si y=0.589743.

Dezvoltam în serie Taylor cele doua functii f si g:

 Daca ne oprim la primii doi termeni ai dezvoltarii obtinem urmatorul proces iterativ:

 La fiecare iteratie se vor actualiza valorile x^n-1 si y^n-1. Pentru cazul test calculul se continua pâna la precizia calculatorului.

Programele exemplificative folosesc o procedura mai generala, utilizând derivarea numerica.

Pentru un numar oarecare de ecuatii, sistemul neliniar se pune sub forma:

 In acest caz, pentru sistemul liniar utilizat iterativ în metoda lui Newton se obtine:

 unde se numesc corectii ale solutiei, iar procesul iterativ se scrie sub forma:

 sau scalar:

 Metoda lui Newton modificata se bazeaza pe observatia lui Kantorovici ca matricea A^-1 se modifica putin în vecinatatea solutiei si de aceea poate fi folosita la mai multe iteratii.

  Ecuatiile diferentiale descriu fenomene tranzitorii, cum ar fi încarcarea unui condensator sau ecuatia lui Schrödinger unidimensionala dependenta de timp.

Cel mai simplu caz, al fenomenelor care depind de o singura variabila, are forma:

 cu conditia initiala (la limita) y(x₀)=a,unde x0 si a sunt constante, provenite, de exemplu, din masuratori experimentale.

Ecuatia diferentiala este folosita pentru a afla valoarea lui y pentru un x diferit de x₀. Deci rezolvarea (integrarea) ecuatiei diferentiale este o metoda de extrapolare.

Cea mai simpla abordare este metoda lui Euler, sau a liniilor poligonale, care foloseste din nou dezvoltarea în serie a functiei y:

 Aceasta relatie permite aproximarea functiei în punctul x din aproape în aproape prin alegerea judicioasa a unui pas, h=x-x₀ suficient de mic.

O aproximare mult mai buna este data de metodele Runge-Kutta, în care se foloseste observatia ca dezvoltarea în serie Taylor produce un polinom de aproximare al functiei, iar derivatele de ordin superior pot fi calculate pornind de la f(x,y). Dupa numarul de termeni din dezvoltare, sunt descrise metodele Runge-Kutta de ordinele I, II, III, IV si V.

Sa scriem noua aproximatie a lui y în forma:

 unde variabilele si desemneaza vecinatatea lui x respectiv y fata de care se face aproximarea. In aceasta formula impunem conditia ca dezvoltarea în serie dupa puterile lui h sa coincida cu dezvoltarea în serie a functiei exacte y(x), pentru o putere cât mai mare, s.

 unde

 Fie s=1 si atunci

Deci, metoda lui Euler este echivalenta cu metoda Runge-Kutta de ordinul I. Functia trebuie calculata o singura data la fiecare crestere a variabilei cu h.

Fie s=2. Dezvoltarea în serie Taylor are forma:

 Cu notatiile Runge-Kutta avem:

 Identificând termenii dupa puterile lui h obtinem:

 Se alege P₂₂=1/2 si atunci P₂₁=1/2, a ₂=1, b ₂₁=1. Rezulta:

 La fiecare etapa de calcul, functia trebuie evaluata de doua ori.

Pentru alegerea judicioasa a pasului se verifica marimea:

Eroarea poate fi evaluata de asemenea prin repetarea calculului pentru un pas înjumatatit:

 Generalizarea pentru sisteme de ecuatii diferentiale ordinare se obtine simplu. Fie cazul unui sistem de doua ecuatii diferentiale ordinare:

 unde:

Ecuatia lui Scroedinger, unidimensionala, dependenta de timp are forma:

Pentru partea reala a funciei de unda avem reprezentarea grafica urmatoare:

Functia de unda se normeaza la unitate în modul cunoscut,

Pachetul de unde a fost astfel ales pentru a prezenta localizarea asteptata a unei particule în spatiu. Trebuie însa sa observam ca ecuatiile diferentiale admit întotdeauna o suprapunere de solutii, pe care le putem scrie în forma unei suprapuneri de armonici:

De aceea, pentru reprezentarea fazei se imprumuta compozitia spectrala a luminii:

Fie acum cazul unei bariere de potential:

Pentru acest caz obtinem urmatoarea reprezentare a undei incidente, reflectate, transmise si a alterarii fazei: