1.5 1.3. Reprezentarea numerelor în calculator.............................................................................. 6 1.5.1 1.3.1. Sisteme de numeraţie.................................................................................................. 6 1.5.2 1.3.2. Reprezentarea internă a numerelor în calculator ........................................................ 7
1.5.2.1 Numere întregi 7 1.5.2.2 Numere Reale 7
1.6 1.4. Erori prin rotunjire................................................................................................................ 9 1.6.1 1.4.1. Erori prin rotunjire în reprezentarea internă a numerelor ............................................ 9 1.6.2 1.4.2. Efecte şi cauze ale erorilor prin rotunjire ................................................................... 10
3.1 3.1. Introducere.......................................................................................................................... 41 3.2 3.2. Formula de interpolare Lagrange ..................................................................................... 42 3.3 3.3. Formule de interpolare Newton prin noduri echidistante .............................................. 43
3.3.1 3.3.1. Formula de interpolare Newton ascendentă.............................................................. 43 3.3.2 3.3.2. Formula de interpolare Newton descendentă............................................................ 45
3.4 3.4. Formula de interpolare Newton prin noduri neuniform distribuite............................... 46
3.4.1 Algoritmul de calcul pentru tabela diferenţelor divizate ....................................................... 47 3.5 Funcţii spline cubice ................................................................................................................. 47
4.2.1 Regulile Simpson (1/3 şi 3/8) ............................................................................................... 58 4.2.2 Algoritmul regulii de cuadratură extinse a lui Simpson (1/3)................................................ 59
4.3 Formule de cuadratură Newton-Cotes .................................................................................... 60 4.4 Cuadratura Gauss ..................................................................................................................... 61 4.5 Cuadratura Gauss-Legendre .................................................................................................... 61 4.6 Cuadratura Gauss-Kronrod...................................................................................................... 65 4.7 Algoritmul de cuadratură automaţi şi adaptivi ....................................................................... 67
4.7.1 Algoritm de cuadratură adaptiv [ ] ....................................................................................... 69 4.8 Integrarea funcţiilor date în puncte ......................................................................................... 71 4.9 Integrale duble ........................................................................................................................... 72
4.9.1 Algoritm de calcul a integralelor multiple [ ] ......................................................................... 75 4.9.1.1 Exemplul Aproximarea integralelor duble 76
5.5.1 Exemplul............................................................................................................................... 86 6. METODE DE REZOLVARE A ECUAŢIILOR NELINIARE ................................... 88
6.2.1 Algoritmul ............................................................................................................................. 91 6.3 7.3 Metoda secantei (corzii)...................................................................................................... 91 6.4 7.4 Metoda aproximaţiilor succesive....................................................................................... 92
6.4.1 Algoritmul ............................................................................................................................. 94 6.5 7.5 Metoda lui Newton (tangentei) ........................................................................................... 95
7.1 Introducere................................................................................................................................. 98 7.2 Metoda dezvoltării în serie Taylor............................................................................................ 99 7.3 Metoda lui Euler (metoda liniilor poligonale)........................................................................ 100 7.4 12.4 Metodele Runge-Kutta .................................................................................................... 101
5.1 Erori absolute şi erori relative
1. 5. NUMERE APROXIMATIVE
1.1 5.1 Erori absolute şi erori relative
Dacă A este valoarea exactă a unei cantităţi şi a este o aproximaţie cunoscută a acesteia, atunci eroarea absolută a aproximaţiei a se consideră de obicei a fi mărimea Δa care satisface relaţia
|a - A| ≤ Δa (5.1) Pentru scopuri practice este convenabil să se ia pentru Δa valoarea cea mai mică pentru care este satisfăcută inegalitatea (5.1) sub circumstanţele date. Numărul exact A poate fi scris atunci
A = a ± Δa (5.2)
Eroarea relativă δa a unui număr aproximativ a este egală cu raportul dintre eroarea absolută Δa a numărului şi valoarea sa absolută:
δa = Δa|a| (a ≠ 0) (5.3)
Uneori, eroarea relativă se defineşte ca raportul Δa/|A|, cu A valoarea exactă însă necunoscută. Introducând Δa din (5.3) în (5.2), rezultă limite pentru numărul exact A:
A = a(1 ± δa) (5.4)
Fie, spre exemplu, numărul aproximativ a = 3.14, utilizat în locul numărului exact A = π. Având în vedere că 3.14<π<3.15, rezultă că |a - π| < 0.01 şi se poate lua pentru eroarea absolută Δa = 0.01. Corespunzător, avem eroarea relativă δa = 0.01/3.14 ≈ 0.003 şi deci putem scrie π = 3.14(1 ± 0.003). Observând că 3.14<π<3.142 obţinem estimări mai bune pentru eroarea absolută, Δa = 0.002, şi respectiv pentru eroarea relativă, δa = 0.002/3.14 = 0.0006. În consecinţă, putem scrie π = 3.14(1 ± 0.0006).
1.2 5.2 Sursele şi clasificarea erorilor
Erorile care apar în rezolvarea numerică a problemelor matematice se împart, în esenţă, în cinci categorii: erori de problemă, erori de metodă, erori iniţiale, erori de rotunjire şi erori de trunchiere.
1
5. NUMERE APROXIMATIVE
Erorile de problemă sunt cauzate de faptul că nu totdeauna formularea matematică descrie exact procesul modelat, deoarece de multe ori, pentru a reduce complexitatea formulării, suntem forţaţi să acceptăm condiţii simplificatoare. Erorile de metodă se datorează faptului că uneori este dificilă, dacă nu chiar imposibilă rezolvarea formulării exacte a problemei. În aceste cazuri, problema este înlocuită cu o problemă aproximativă, pentru care există tehnici adecvate de rezolvare şi care are un rezultat foarte apropiat. Metodele numerice sînt în majoritatea contextelor în care sînt utilizate metode de aproximare. Erorile iniţiale (inerente) sînt erori în valorile datelor, cauzate de incertitudini în măsurători sau de natura inerent aproximativă a reprezentării numerelor cu ajutorul unui număr finit de cifre. Erorile de rotunjire sînt cauzate de reprezentarea numerelor (date iniţiale sau rezultate ale unor calcule) cu un număr finit de cifre semnificative exacte. De exemplu, reprezentarea rezultatului operaţiei 1/3 sub forma 0.333 implică o eroare de trunchiere de aproximativ 3x10-4. Erorile de rotunjire depind de particularităţile hardware ale calculatorului, de modul de reprezentare internă al diferitelor tipuri de date utilizate în calcule. Erorile de rotunjire se acumulează prin creşterea numărului de calcule, mai ales al celor care implică scăderea unor valori aproximativ egale. Erorile de. trunchiere (reziduale) provin din natura infinită a unor procese utilizate în descrierea soluţiei problemelor matematice. Faptul că, practic, aceste procese trebuie întrerupte după un număr finit de paşi, introduce o eroare de trunchiere. Considerând exemplul dezvoltării funcţiilor în serie, însumarea trebuie întreruptă la un anumit termen, trunchiind seria şi introducând o anumită eroare de trunchiere. Metodele numerice permit, în principiu, controlul erorilor de rotunjire şi al celor de trunchiere, cu condiţia transcrierii adecvate a algoritmului. Astfel, în cazul erorilor de rotunjire, trebuie considerată acea formă a expresiilor matematice, în care sunt evitate operaţiile care introduc acest tip de erori. În ceea ce priveşte erorile de trunchiere, acestea se află sub controlul deplin al programatorului, care însă trebuie să fie preocupat de construirea unor algoritmi finiţi şi eficienţi. Dacă finitudinea algoritmului este obligatorie pentru convergenţa rezultatelor, eficienţa este o reflectare directă a experienţei şi stilului de programare. Se poate afirma (exagerând întrucâtva) că minimizarea erorilor de trunchiere constituie punctul central al domeniului analizei numerice.
Orice număr pozitiv poate fi reprezentat în baza 10 sub forma
(5.7)
unde αi sînt cifrele numărului a, cu αm≠ 0. De exemplu,
(m = 2 ⇒ n = 5) Situaţiile concrete implică utilizarea unor numere aproximative cu un număr finit de cifre. Toate cele n cifre zecimale reţinute αi, i = m, m-1, m - n + 1 , se numesc cifre semnificative. Unele dintre cifrele semnificative pot fi egale cu zero, cu excepţia primei cifre, αm. De exemplu, în numărul 0,00208 = 2⋅10-3+ 0⋅10-4 + 8⋅10⋅5, primele trei zerouri nu sînt cifre semnificative, deoarece ele servesc numai pentru fixarea poziţiei punctului zecimal în scrierea zecimală a numărului.
2
5.4 Erori ale operaţiilor elementare
Definiţie: Spunem că un număr aproximativ a are n cifre semnificative exacte, αm, αm-1, … αm-
n+1, dacă eroarea absolută a numărului nu depăşeşte jumătate de unitate în poziţia n, numărând de la stânga la dreapta, adică
(5.8) De exemplu, numărul a = 36.00 este în raport cu numărul exact A =35.97 o aproximaţie cu trei cifre semnificative exacte. într-adevăr, având în vedere că |a - A|=0.03<1/2-10-1, rezultă m-n+1=-1, de unde n =3 (m=1). Termenul "n cifre semnificative exacte" nu trebuie luat ad literam, deoarece nu este obligatoriu ca într-un număr aproximativ a având n cifre exacte, primele n cifre semnificative să coincidă cu cifrele corespunzătoare ale numărului exact A. De exemplu, numărul a=9.995 este o aproximare cu trei cifre corecte a numărului exact A=10 şi totuşi are toate cifrele diferite.
Teoremă. Dacă un număr pozitiv a are n cifre exacte, eroarea relativă δa a acestui număr satisface inegalitatea
(5.9) unde am este prima cifră semnificativă a numărului a. Într-adevăr,
Numărul de cifre semnificative exacte corespunzător unei erori relative δa date este conform relaţiei (5.9)
(5.10)
De exemplu, numărul aproximativ a=3.15, utilizat în locul numărului exact A=π, are n ≤ 2.72, adică două cifre exacte. Numărul aproximativ a = 3.14 va avea n ≤ 3.29, adică trei cifre exacte. Pentru numărul de cifre semnificative exacte, se poate obţine o estimare grosieră însă foarte utilă în practică, considerând că eroare relativă are forma δa"≈ 10-p şi αm ≈ 5: n ≈ p (5.11) Acest criteriu poate fi utilizat în situaţiile în care, cunoscându-se eroarea relativă a unui număr, se impune estimarea rapidă însă nu neapărat foarte precisă a numărului de cifre exacte ale acestuia.
1.4 5.4 Erori ale operaţiilor elementare
1.4.1 5.4.1 Eroarea sumei Teoremă. Eroarea absolută pentru suma algebrică a mai multor numere aproximative nu depăşeşte suma erorilor absolute ale numerelor. Într-adevăr, considerând suma algebrică a numerelor aproximative a1, a2,…,an
a = ± a1 ± a2 ± …± an
avem în mod evident
3
5. NUMERE APROXIMATIVE
Δa = ± Δa1 ± Δa2 ± …± Δan
şi deci
|Δa| ≤ |Δa1| + |Δa2| + …+ |Δan|
Pentru eroarea absolută a sumei avem, prin urmare,
Δa ≤ Δa1 + Δa2 + …+ Δan (5. 12)
iar pentru eroarea absolută limită,
Δa* = Δa1 + Δa2 + …+ Δan (5,13)
Din această relaţie rezultă că eroarea absolută limită nu poate fi mai mică decît eroarea absolută a celui mai puţin exact termen din sumă
Δa* ≥ max(Δa1, Δa2,… , Δan) (5.14)
în consecinţă, ceilalţi termeni, cu grad de precizie mai mare (cu erori absolute mai mici) nu pot ameliora precizia rezultatului. Teoremă. Dacă toţi termenii unei sume au acelaşi semn, eroarea relativă (limită) a sumei nu depăşeşte cea mai mare eroare relativă a termenilor. Într-adevăr, fie numerele aproximative pozitive a1, a2, …, an, ale căror valori exacte sînt respectiv A1, A2,…An. Notând cu a suma valorilor aproximative şi cu A suma valorilor exacte, are loc pentru eroarea relativă a sumei
şi deci de unde, eroarea relativă limită este
(5.16)
1.4.2 5.4.2 Eroarea diferenţei Considerăm diferenţa a două numere aproximative
a = a1 - a2 Conform relaţiei (5.13), eroarea absolută limită a diferenţei este egală cu suma erorilor absolute ale celor doi termeni
Δa* = Δa1 + Δa2
şi eroarea relativă limită a diferenţei va fi
unde A este valoarea exactă a modulului diferenţei. Dacă numerele aproximative a1 şi a2 sînt foarte apropiate ca valoare, atunci diferenţa exactă A este mică şi, chiar în condiţiile în care erorile relative δa1 şi δa2 sunt mici, eroarea relativă limită a diferenţei, δa*, poate fi foarte mare.
4
5.4 Erori ale operaţiilor elementare
Pentru a exemplifica cele arătate mai sus, considerăm numerele aproximative a1 = 47.132 şi a2= 47.111 , fiecare având cinci cifre semnificative exacte, adică o eroare absolută de cel mult 0.0005. Diferenţa a = 47.132 - 47.111 = 0.021 are doar două cifre semnificative exacte, iar eroarea absolută limită a diferenţei este
şi vom avea următoarele erori relative limită
După cum se observă, eroarea relativă limită a diferenţei este de aproximativ 5000 de ori mai mare decît erorile relative ale termenilor. Este de aceea de dorit, ca în calcule numerice să se rescrie expresiile care implică scăderea unor numere aproximativ egale.
1.4.3 5.4.3 Eroarea produsului Teoremă. Eroarea relativă a produsului mai multor numere aproximative nenule nu depăşeşte suma erorilor relative ale numerelor. Fie a = a1a2...an. Presupunând pentru simplitate că toţi factorii sînt pozitivi, avem
ln a = ln a1 + ln a2 +…+ ln an Utilizând formula de aproximare
obţinem Are loc inegalitatea
(5
Dacă Ai sînt valorile exacte ale factorilor şi dacă, aşa cum se întîmplă în mod obişnuit |Δai| sînt mici în comparaţie cu ai, putem considera aproximaţiile
şi relaţia (5.17) poate fi rescrisă
(5 Este evident că această inegalitate rămâne valabilă şi dacă factorii au semne diferite. Corolar. Eroarea relativă limită a produsului este egală cu suma erorilor relative ale factorilor
(5 De asemenea, avem
5
5. NUMERE APROXIMATIVE
Dacă toţi factorii produsului sînt exacţi cu excepţia unuia, atunci conform relaţiei (5.19) eroarea relativă a produsului coincide practic cu eroarea relativă a factorului aproximativ. În cazul particular al unui factor aproximativ şi al unuia exact (k ≠ 0)
a = ka1
avem
adică, atunci când se înmulţeşte un număr aproximativ printr-un factor exact, eroarea relativă rămâne neschimbată, în timp ce eroarea absolută creşte de |k| ori.
1.4.4 5.4.4 Eroarea împărţirii Dacă luăm în considerare cazul împărţirii a două numere aproximative a1 şi a2 (pentru simplitate, considerate pozitive), a = a1/a2, avem
de unde Această formulă arată că eroarea relativă a unei împărţiri nu depăşeşte suma erorilor relative ale deîmpărţitului şi împărţitorului.
1.4.5 5.4.5 Eroarea relativă a puterii Considerăm operaţia de ridicare a numărului aproximativ a1 la o putere reală p, a = a1
p. Prin logaritmare obţinem ln a = p ln a1 şi
⇒ adică, eroarea relativă a puterii p a unui număr aproximativ este de |p| ori mai mare decît eroarea relativă a numărului. Este evident că dacă |p|<1, rezultă, de fapt, o scădere a erorii relative.
1.5 1.3. Reprezentarea numerelor în calculator
Considerăm necesară cunoaşterea câtorva noţiuni fundamentale în acest domeniu pentru o înţelegere ulterioară mai bună a cauzelor care conduc la apariţia erorilor prin rotunjire. Cunoaşterea acestor cauze permite, de multe ori, evitarea efectelor nedorite a erorilor prin rotunjire printr-o tehnică corecta de elaborare, a algoritmilor.
1.5.1 1.3.1. Sisteme de numeraţie Sistemul de numeraţie în care lucram zilnic este sistemul zecimal având baza 10. Din motive tehnologice, calculatoarele electronice utilizează pentru reprezentarea internă a numerelor şi pentru calcule sistemul de numeraţie binar. Analizând limbajele interne ale calculatoarelor vom observa că sunt utilizate şi alte sisteme de numeraţie, cum ar fi sistemul octal care are baza 8 şi sistemul hexazecimal având baza 16. în sistemul binar fiecare cifră utilizată pentru reprezentarea unui număr are valoarea 0 sau 1, în sistemul octal cifrele sunt cuprinse în intervalul 0-7 iar în sistemul hexazecimal fiecare cifră este cuprinsă în intervalul 0-15, cifrele de la 10 la 15 fiind exprimate prin literele A,B,C,D,E şi F. Valorile zecimale corespunzătoare cifrelor din sistemul hexazecimal sunt: zecimal 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
6
1.3. Reprezentarea numerelor în calculator
hexazecimal 0 1 2 3 4 5 6 7 8 9 A B C D E F Baza unui sistem de numeraţie mai poartă şi denumirea de rădăcină şi poate fi specificată printr-un indice inferior. Utilizând scrierea poziţională, un număr oarecare, cu parte fracţionara, poate fi reprezentat intr-un sistem de numeraţie cu baza b, în forma:
±(anan-1…a1a0,a-1a-2…a-k)b unde ai, i = -k,...,-2,-1,0,1,...,n-1, sunt cifre. De exemplu (101,1101)2, (16C6,92)16 Numărul real din sistemul zecimal corespunzător unui număr exprimat într-un sistem de numeraţie oarecare este dat de următoarea relaţie polinomială:
Da exemplu, numărului (101,11)2 îi corespunde în sistemul zecimal numărul real 5,75, deoarece:
Prima cifră diferită de zero de la stânga numărului se numeşte prima cifră semnificativă.
1.5.2 1.3.2. Reprezentarea internă a numerelor în calculator Am arătat anterior că, din motive tehnologice, majoritatea calculatoarelor utilizează reprezentarea internă a numerelor în sistemul binar. Considerând un număr oarecare reprezentat în sistemul binar, cifrele 0 sau 1 poartă denumirea de bit. Bit este acronimul pentru cifră binară şi în ceea ce priveşte calculatorul este utilizat atunci când ne referim la un dement al memoriei având poziţia "deschis" sau "închis". Byte este un set de opt biţi şi este considerat ca unitate de referinţă. Există două categorii de numere care pot fi utilizate în calculele efectuate pe un calculator: numere întregi şi numere în virgulă mobilă (numere reale). Ne vom referi în acest capitol numai la modul de reprezentare a numerelor pe un calculator IBM - PC în limbajul de programare C.
1.5.2.1 Numere întregi Fie un întreg exprimat în sistemul binar: (anan-1…a1a0) unde ai, i = n, n-1,...,0 sunt cifre binare 0 sau 1. Într-un calculator valoarea maximă a numărului de biţi este limitaţi. Pe un IBM - PC (Turbo Pascal) se utilizează 2 bytes pentru reprezentarea unui întreg (variabile de tip integer). Primul bit din stânga înregistrează semnul: 0 → pozitiv şi 1 → negativ. Restul de 15 biţi sunt utilizaţi pentru reprezentarea cifrelor binare a14, a13,...,a1, a0. Astfel cel mai mare întreg pozitiv este:
Valoarea corespunzătoare în zecimal pentru acest număr este 214 + 213 + 212 + . . .+22+21+20= 32767 Cel mai mic întreg negativ ce poate fi reprezentat pe 2 bytes este (-32767)10.
1.5.2.2 Numere Reale Reprezentarea numerelor în virgulă mobilă, în simplă precizie, se realizează pe un IBM - PC (în Turbo Pascal) pe 6 bytes (48 biţi) în forma binară normalizată.
7
5. NUMERE APROXIMATIVE
Să considerăm un număr binar oarecare, având parte întreagă şi parte fracţionară, exprimat în scriere poziţională, în forma generală: (anan-1…a1a0,a-1a-2…a-k)2 În forma normalizată numărul este dat de relaţia: (0,anan-1…a1a0,a-1a-2…a-k)2(n+1) De exemplu numărul binar (101,11) se va scrie în forma normalizată (0,10111)x23. Numărul real corespunzător în zecimal numărului binar dat este: (1*2-1+0*2-2+1*2-3+1*2-4+1*2-5)*23 = 5,7510 În mod asemănător se vor reprezenta în forma normalizată şi numerele binare subunitare. De exemplu (0,0001011)2 se va scrie în forma normalizată (0,101l)x2J. Revenind la reprezentarea internă a unui număr în virgulă mobilă pe un IBM - PC (Turbo Pascal), cei 48 de biţi sunt utilizaţi după cum urmează: primul bit de la stânga pentru semn, următorii 8 biţi pentru exponent şi restul de 39 de biţi pentru mantisă:
0 12345678 9 10 11 ...... 45 46 47 s eeeeeeee m1, m2... m37 m38 m39, 1 bit pentru semn 8 biţi pentru exponent 39 biţi pentru mantisă s, e, m1...,m39 exprimă cifrele binare 0 sau 1. Pornind de la această reprezentare, un număr zecimal oarecare x, în virgulă mobilă este exprimat în următoarea forma normalizată: x = (semn)(0,1m1m2…m38m39)2(g-p) dacă g > 0, şi x = 0, dacă g = 0 unde am notat cu g exponentul translabil şi cu p constanta de translare. Exponentul g este întreg şi se reprezintă pe cei 8 biţi pentru exponent, ceea ca înseamnă că poate lua orice valoare cuprinsă în intervalul [0,255]. În particular, 0 şi 255 sunt valori rezervate. Constanta p are valoarea 128 (pentru un calculator IBM - PC în Turbo Pascal). Exprimând în acest mod exponentul lui 2 din forma normalizată se evită cerinţa unui bit suplimentar pentru semnul acestui exponent. În forma normalizată prima cifră semnificativă este întotdeauna 1 astfel încât numărul efectiv de cifre binare pentru mantisă este 40. Intervalul în care poate lua valori exponentul lui 2 în forma normalizată este -127 ≤ g-p ≤ +126, de unde 2-127 ≤ 2(g+p) ≤ 2-126. Cel mai mic număr care poate fi reprezentat pe o mantisă de 40 biţi este (0,1)2 = 0,510 iar cel mai mare este (0,11...11)2 = (1 - 2-40)10. Obţinem în acest mod cel mai mic număr în virgulă mobilă care poate fi reprezentat în calculator: 0,5 x 2-127 = 2,938... x 10-39, respectiv cel mai mare număr în virgulă mobilă reprezentabil: (1 - 2-40)* 2127 = 1,70141...x 1038. Deoarece numărul de biţi utilizat pentru reprezentarea mantisei şi. a exponentului este limitat, rezultă o primă concluzie importantă: pe un calculator IBM - PC (în Turbo Pascal), în simplă precizie, există un număr finit de numere reale, în virgulă mobilă, care permit o reprezentare exactă: 248 - 1. Limbajul de programare C oferă posibilitatea utilizării şi altor tipuri de date reale în funcţie de necesităţi (vezi tabelul 1.2.).
Tip Dimensiunea Domeniu unsigned char 8 biţi 0 la 255
char 8 biţi -128 la 127 enum 16 biţi -32,768 la 32,767
unsigned int 16 biţi 0 la 65,535 short int 16 biţi -32,768 la 32,767
8
1.4. Erori prin rotunjire
int 16 biţi -32,768 la 32,767 unsigned long 32 biţi 0 la 4,294,967,295
long 32 biţi -2,147,483,648 la 2,147,483,647float 32 biţi 3.4 x 10-38 la 3.4 x 10+38
double 64 biţi 1.7 x 10-308 la 1.7 x 10+308 long double 80 biţi 3.4 x 10-4932 la 1.1 x 10+4932
near (pointer) 16 biţi not applicable far (pointer) 32 biţi not applicable
1.6 1.4. Erori prin rotunjire
An demonstrat în capitolul anterior că există un număr limitat de numere reale care pot fi reprezentate exact în calculator, altfel spus într-un calculator nu există o continuitate a numerelor reale. De exemplu, dacă cel mai mic număr real care permite o reprezentare exactă în calculator este 2,938...x10-39, rezultă că nici un număr cuprins în intervalul (0, 2,9x10-39) nu poate fi reprezentat. Aceeaşi afirmaţie este valabilă şi pentru orice număr real aparţinând intervalului (-2,9x10-39, 0). Cum putem determina intervalul dintre două numere reale consecutive care pot fi reprezentate exact în calculator ş Pentru aceasta vom afla mai întâi o constantă importantă a oricărui calculator şi anume constanta epsilon care reprezintă diferenţa dintre numărul 1 şi cel mai mic număr care este mai mare decât 1 şi care este reprezentabil în calculator, altfel spus, cel mai mic număr real, pozitiv care adunat la 1 va da în calculator un număr diferit de 1. Reprezentarea binară în formă normalizată pentru 1 este (pe un calculator IBM - PC în Turbo Pascal): (1)10 = (0,1000 0000 0000 0000 0000 0000 0000 0000 0000 0000)2 x 21 Cel mai mic număr mai mare decât 1 care poate fi reprezentat pe mantisa de 40 biţi este: (0,1000 0000 0000 0000 0000 0000 0000 0000 0000 0001)2 x 21, Numărul real zecimal corespunzător acestui număr binar este (1*2-1 + 1*2-40)*21 = 1 + 2-39 = 1 + 1,81898... x10-12 = 1,000000000000181898...
Rezultă astfel pentru constanta epsilon: ε = 1,81898... x10-12 Cunoscând constanta epsilon a calculatorului putem afla diferenţa dintre orice număr real ri, care permite o reprezentare exactă în calculator şi următorul număr real ri+1, reprezentabil exact pe baza relaţiei: ri+1 - ri = epsilon x ri
1.6.1 1.4.1. Erori prin rotunjire în reprezentarea internă a numerelor Principala cauză a erorilor prin rotunjire este generată de reprezentarea numerelor reale în forma binară printr-un număr limitat de biţi.
Am arătat anterior că nici un număr real cuprins între 1 şi 1 + ε nu poate fi reprezentat în calculator. în acest caz, dacă în urma unor calcule obţinem un număr de forma 1 + α, unde α ∈ (0,ε), acest număr va fi aproximat, prin tăiere, cu 1 dacă: 0 < α < ε/2 sau va fi rotunjit la valoarea 1 + ε dacă ε > α > ε/2
9
5. NUMERE APROXIMATIVE
Acest proces de tăiere sau rotunjire este cunoscut în literatura de specialitate şi sub denumirea de rotunjire simetrică. Să generalizăm. Fie intervalul dintre două numere reale consecutive cu reprezentare exactă în calculator (ri, ri+1). Orice număr real de forma ri+α, cuprins în acest interval α ∈ (0, ri+1 - ri,)) va fi aproximat prin tăiere la valoarea r, dacă: 0 < α < ε⋅ ri/2 sau va fi rotunjit la valoarea ri+1, dacă α ≥ ε⋅ri/2 De exemplu, numărul 0,143 nu poate fi reprezentat exact în calculator şi este aproximat, prin tăiere, la valoarea 0,14299999986...; numărul 0,1431 nu poate fi reprezentat exact şi este aproximat prin rotunjire la valoarea 0,14310000025... Eroarea care se face la reprezentarea unui număr oarecare în calculator este aproximativ egală cu ε x r/2 dacă numărul este rotunjit şi ε x r dacă este tăiat.
1.6.2 1.4.2. Efecte şi cauze ale erorilor prin rotunjire Unul dintre efectele erorilor prin rotunjire se manifestă la adunarea şi scăderea numerelor. În general scăderea şi adunarea numerelor în virgulă mobilă, normalizate, ca urmare a modului în care se efectuează, necesită reprezentarea rezultatului cu un număr mai mare de cifre semnificative decât al numerelor însumate pentru a obţine o precizie satisfăcătoare. Să adunăm numerele (1)10 şi (0,0000001)10. Reprezentarea binară în formă normalizată a acestor numere este ( pe un IBM -PC în Turbo Pascal şi simplă precizie): (1)10 = (0,1000 0000 0000 0000 0000 0000 0000*0000 0000 0000 )2 x 21 (0,00000001)10 = (0,1010 0111 1100 0101 1010 1100 0110 1001 1100 1011)2 x 2-26
Pentru a aduna aceste două numere este necesar să le aducem la acelaşi exponent (exponentul mai mare): 1. Această operaţie conduce la denormalizarea numărului cu exponent mai mic, altfel spus se deplasează mantisa cu 27 de poziţii spre dreapta, aceste poziţii fiind ocupate de cifra 0. Următoarea etapă constă în adunarea mantiselor. Evident că rezultatul exact necesită o reprezentare pe o mantisă de 40 + 27 = 67 de biţi. Mantisa pe care o avem la dispoziţie fiind de numai 40 de biţi se vor pierde ultimele 27 de cifre binare semnificative. Rezultatul real este: (0,1000 0000 0000 0000 0000 0000 0001 0100 1111 1000)2 x 21
ceea ce în zecimal înseamnă (1,000000010008)10. Dacă veţi aduna (0,00000001)10 de 10000 de ori la 1, la fiecare adunare se va însuma şi o eroare de 8x10-12. Rezultatul va fi 1,000100008. O primă concluzie care se desprinde din acest exemplu este: evitaţi însumarea numerelor diferite considerabil ca ordin de mărime. În cazul în care aveţi o problemă în care trebuie să realizaţi o sumă de mai mulţi termeni diferiţi ca ordin de mărime, adunaţi mai întâi termenii de ordin mic. Pentru o mai bună înţelegere să considerăm două numere reale x şi y având următoarea reprezentare binară în forma normalizată pe o mantisă de 40 de biţi:
x10 = (0,1m1m2…mpmp+1mp+2…mk)⋅2n
y10 = (0,1m'1m'2…m'pm'p+1m'p+2…m'k)⋅2p Rezultatul diferenţei x - y este
(x - y)10 = (0,1m"1m"2…m"pm"p+1m"p+2…m"k00…0)⋅2n-p unde
10
1.5. Concluzii
(m"p+1m"p+2…m"k)2 = (mp+1mp+2…mk)2 - (m'pm'p+1m'p+2…m'k)2 Ce efect are asupra rezultatului modul în care s-a efectuat scăderea ş Mai întâi trebuie să ţinem cont de faptul ci un număr foarte mic de numere reale zecimale admit o reprezentare exactă în binar cu un număr finit de cifre binare. La scăderea a două numere binare în virgulă mobilă în formă normalizată şi având valori apropiate, mantisa rezultatului se deplasează spre stânga cu un număr da poziţii egal cu numărul cifrelor semnificative comune celor două numere care se scad iar exponentul se va micşora cu acest număr. În partea dreaptă a mantisei va rezulta un număr de cifre binare 0 egal cu numărul poziţiilor cu care s-a deplasat mantisa spre stânga. În situaţia în care două numere reale care se scad sunt aproximativ egale sa va pierde un număr de cifre binare semnificative egal cu numărul poziţiilor cu care s-a deplasat mantisa spre dreapta. Astfel, cu cât numerele care se scad sunt mai apropiate, numărul de cifre semnificative care se va pierde este mai mare. Mai apar erori prin rotunjire şi la scăderea numerelor mult diferite ca ordin de mărime. Rezultă în acest fel o a doua regulă de care trebuie să se ţină cont în elaborarea unui algoritm: evitaţi scăderea numerelor aproximativ egale sau mult diferite ca ordin de mărime.
1.7 1.5. Concluzii
O analiză superficială a erorilor prin rotunjire ne poate duce la concluzia greşită că efectul acestora este neglijabil. În realitate, într-o problemă complexă care implică milioane de calcule erorile prin rotunjire pot afecta considerabil rezultatul final. Acest efect poate fi evitat prin elaborarea unor algoritmi stabili numeric. Un algoritm este stabil numeric dacă modificări mici ale datelor de intrare vor genera modificări mici ale rezultatului final. Din păcate nu există o regulă unică care să ne permită elaborarea unor algoritmi stabili. Câteva strategii posibil de aplicat astfel încât să se reducă efectul erorilor prin rotunjire sunt:
1. Efectuarea calculelor în dublă sau multiplă precizie, opţiune existentă în Turbo Pascal. Dezavantajul acestei alegeri constă în creşterea considerabilă a timpului de calcul şi în faptul că erorile prin rotunjire nu sunt eliminate ci numai micşorate.
2. Metoda grupării. De exemplu am văzut în capitolul 1.4.2. că dacă adunăm de 10000 de ori 0,00000001 la 1 vom obţine 1,00100008. Dacă veţi efectua calculele adunând de cate o sută de ori 10-8, rezultatul adunându-l de o sută de ori la 1, veţi obţine 1,001.
3. Rescrierea ecuaţiilor pe baza cărora se elaborează algoritmul.
11
2.SISTEME DE ECUAŢII LINIARE
2. 2.SISTEME DE ECUAŢII LINIARE
2.1 2.1. Metode directe
2.1.1 2.1.1. Introducere Una dintre problemele frecvent întâlnite în calcule ştiinţifice constă în rezolvarea unor sisteme de ecuaţii liniare. În acest capitol vom descrie o serie de metode directe de rezolvare a sistemelor liniare, urmând ca în capitolul următor să abordăm şi tehnicile iterative. Tehnicile directe sunt metodele care generează soluţia finală într-un număr dat de paşi şi nu sunt afectate decât de erori prin rotunjire. Un sistem liniar de n ecuaţii cu n necunoscute este de forma:
Introducând notaţia matricial vectorială, sistemul (2,1) poate fi scris ca fiind: Ax = b unde A aste matricea coeficienţilor, x este vectorul coloană al necunoscutelor, şi b este vectorul coloană al termenilor neomogeni (liberi):
A = ⎝⎜⎛
⎠⎟⎞a11 a12 … a1n
a21 a22 … a2n…
an1 an2 … ann
x = ⎝⎜⎛
⎠⎟⎞x1
x2…xn
, b = ⎝⎜⎛
⎠⎟⎞b1
b2…bn
Se cunosc multe metode de rezolvare a sistemelor liniare. O altă posibilitate teoretică de rezolvare a sistemelor liniare Ax = b constă în scrierea soluţiei sub forma x = A-1b, unde A-1 este inversa matricei A. După cum vom vedea la capitolul 2.1.5. sau la PROGRAMUL 2.1., inversarea matricii A necesită de fapt rezolvarea unui sistem matricial de forma AX = I, unde A, X şi I sunt matrici de dimensiune n x n; altfel spus este necesară rezolvarea a n sisteme liniare de ordin n ceea ce implică un efort de calcul mai mare decât pentru sistemul Ax = B. Vom prezenta în continuare două dintre cele mai cunoscute tehnici de rezolvare, a sistemelor de ecuaţii liniare.
2.1.2 2.1.2. Eliminarea Gauss şi eliminarea Gauss-Jordan Vom considera pentru început metoda eliminării Gauss în rezolvarea unui sistem liniar de ecuaţii algebrice Ax = b , considerând matricea A nesingulară.(O matrice A este singulară dacă
12
2.1. Metode directe
nu admite inversă sau, dacă determinantul matricii A este egal cu zero sau, dacă există linii ale matricii liniar dependente sau dacă există un vector z diferit de zero astfel încât Az = 0.) Eliminarea Gauss se bazează pe ideea reducerii unui sistem de ecuaţii la o formă mai simplă care este mai uşor de rezolvat. Pentru o mai rapidă înţelegere a metodei să considerăm mai întâi o aplicaţie practică după care vom analiza algoritmul general al eliminării Gauss.
2.1.2.1 Exemplul 2.1.1. Eliminarea Gauss şi substituţia inversă Fie sistemul de ordinul 3:
La primul pas vom elimina necunoscuta x, din ecuaţiile doi şi trei. Pentru aceasta vom scădea prima ecuaţia înmulţită eu 0,5 din a doua ecuaţia şi vom scădea prima ecuaţia înmulţită cu -0,3 din a treia ecuaţie. Cu această transformare a ecuaţiilor doi şi trei sistemul devine
La următorul pas vom utiliza a doua ecuaţie pentru a elimina x2 din a treia ecuaţie. Vom scădea a doua ecuaţie înmulţită cu -0,04 = -0,1 /2,5 din a treia ecuaţie. Rezultă:
Am obţinut în acest mod un sistem a cărui matrice A este superior triunghiulară. Operaţiile aplicate până la acest pas cu scopul de a aduce sistemul la o formă superior triunghiulară fac parte dintr-o primă etapă a metodei Gauss, etapă cere poartă denumirea de eliminare directă. Următoarea etapă importantă, substituţia inversă, constă în determinarea necunoscutelor. Din a
treia ecuaţie obţinem
Această valoare este substituită în a doua ecuaţie din care rezultă x2:
În final, substituim valorile x3 şi x2 în prima ecuaţie şi aflăm x1: Să încercăm în continuare să generalizăm procedura Gauss. Pentru o mai bună raţionalizare, formal, vom construi matricea extinsă a sistemului prin concatenarea matricilor A şi b:
unde elementele de pe coloana n+1 sunt elementele vectorului b, ai,n+1 = b, pentru i = 1, 2, ...,n.
13
2.SISTEME DE ECUAŢII LINIARE
Prima etapă a metodei, eliminarea directă, constă în n-1 paşi. La pasul k, multipli ai ecuaţiei k se scad din restul de n-k ecuaţii pentru a elimina variabila k. La acest pas elementele matricii à se modifică conform formulei
aij ← aij - ⎝⎜⎛
⎠⎟⎞aik
akk ⋅akj
Termenii aik/akk poartă denumirea de multiplicatori. Aceste calcule sunt operate pentru i > k şi j > k. Elementul akk poartă denumirea de pivot. (în capitolul următor vom analiza ce este de făcut în cazul în care la un pas oarecare vom întâlni un pivot având valoare foarte mică). Mai precis, eliminarea directă constă în generarea unei serii de matrici extinse Ã(1), Ã(2),..., Ã(n) unde Ã(1) este identică cu à dată de relaţia (2.3.) şi Ã(k), pentru k = 2,3,...,n, conţine elementele aij
(k), unde:
Matricea extinsă Ã(k) obţinută la pasul k va fi de forma
şi va reprezenta sistemul liniar echivalent la pasul în care a fost eliminată variabila xk-1, din ecuaţiile k, k+1,...,n. În finalul procesului de eliminare Gauss obţinem următoarea matrice extinsă care exprimă sistemul original într-o formă echivalentă superior triunghiulară, mai simplă:
Elementele primei linii, a1,j
(1), j = 1,2,…, n+1 sunt identice cu elementele a1,j ale matricii extinse iniţiale, Ã Următoarea etapă, substituţia inversă, începe de la ultima ecuaţie din care se poate afla x:
În continuare se determină necunoscutele xn-1, …, x1, prin aplicarea relaţiei recursive
14
2.1. Metode directe
Sintetizând, metoda eliminării Gauss cu substituţie inversă poate fi exprimată în următorul algoritm:
2.1.3 Algoritmul eliminării Gauss cu substituţie inversă Rezolvă un sistem liniar de n ecuaţii cu n necunoscute DATE DE INTRARE:
• numărul n al necunoscutelor şi ecuaţiilor • matricea extinsă Ã = (aij) unde 1≤i≤n şi 1≤j≤n+1
DATE DE IEŞIRE: • soluţia x1,x2,...,xn, sau mesaj: "Sistemul nu are soluţie unică"
(Eliminarea directă:) Pasul 1. PENTRU k = 1, 2,...,n-1 EFECTUEAZĂ paşii 2-5:
Pasul 2. DETERMINAŢI cel mai mic întreg p, 1≤ p ≤n, astfel încât apk ≠ 0. Pasul 3. DACĂ nu există un astfel de întreg p atunci SCRIE "Sistemul nu are soluţie unică" şi EXIT.
Pasul 4. DACĂ p ≠ k ATUNCI SCHIMBAŢI între ele liniile p şi k: (Lp) ↔ (Lk) Pasul 5. PENTRU i = k+1,...,n (i - indice de linie), EFECTUAŢI pasul 6:
Pasul 6. MODIFICAŢI linia i: (Li) = (Li - mik Lk), unde mik = aikakk
Pasul 7. DACĂ ann = 0 atunci SCRIE mesaj: "Sistemul nu are soluţie unică" şi EXIT. (Substituţia inversă:)
Pasul 8. CALCULEAZĂ: xn = an n+1ann
Pasul 9. PENTRU i = n-1,...,1 SUMA = 0 Pasul 10 PENTRU j = i+1,…,n CALCULEAZĂ SUMA = SUMA + aijxj
Pasul 11. CALCULEAZĂ xi = ai n+1 - SUMAaii
Pasul 12. SCRIE (x1,x2, .. ,xn) STOP. La pasul 2 se schimbă între ele linia k cu prima dintre liniile k+1,...,n care conţine pe coloana p elementul apk, diferit de zero. Astfel evităm împărţirea cu zero în cazul în care akk = 0. Se poate demonstra că numărul total de operaţii aritmetice necesar la rezolvarea unui sistem liniar de n ecuaţii necunoscute, prin aplicarea algoritmului anterior, este
- Înmulţiri şi împărţiri: - adunări şi scăderi:
15
2.SISTEME DE ECUAŢII LINIARE
Vom analiza în continuare o altă variantă a metodei Gauss şi anume eliminarea Gauss - Jordan. Să începem cu o aplicaţie:
2.1.3.1 Exemplul 2.1.2. Eliminarea Gauss - Jordan Să se rezolve următorul sistem liniar:
Matricea extinsă a sistemului este: Prima etapă este, la fel ca fi la metoda Gauss, eliminarea directă. Aplicând acest proces obţinem următoarea matrice extinsă:
Următoarea etapă se desfăşoară după cum urmează: se împarte ultima linie cu a33 = 11/7
La următorul pas se scade linia a treia înmulţită cu a12 = 1/2 din linia a doua şi linia a treia înmulţită cu a13 = -3 din prima linie:
A doua linie se împarte la a22 = 7:
A doua linia înmulţită cu a12 = 1 se scade din prima linie:
în final sa împarte prima linie la a11 = 2: În urma acestui proces ultima coloană conţine soluţia sistemului (x1 = 1, x2 = 3, x3 = 2) iar matricea compusă din primele trei linii şi trei coloane are toate elementele egale cu zero cu excepţia diagonalei principale care are toate elementele egale cu 1. Procedeul descris în acest exemplu poate fi extins la sisteme liniare de orice dimensiune. Prima etapă constă în eliminarea directă şi a fost descrisă la metoda Gauss. A doua etapă începe de la matricea extinsă a sistemului triangularizat şi transformă în 1 toate elementele de
16
2.1. Metode directe
pe diagonala principală (diagonala pivoţilor) iar restul elementelor situate deasupra diagonalei principale le face egale cu 0. La primul pas ultima linie a matricii extinse triangularizate se împarte la an,n
(n), împărţire care modifică această linie: 0 0 0 … 1 an,n+1 unde
La următorul pas se elimină coeficienţii de pe coloana n a liniilor n-1, n-2,...,1 scăzând succesiv linia n înmulţită cu ai,n
(i), i = n-1, n-2,...,1, din liniile n-1, n-2,...,1. Obţinem:
unde
În continuare se împarte linia n-1 la an-1,n-1
(n-1), după care se elimină coeficienţii de pe coloana n-1 şi liniile n-2, n-3,...,1 scăzând din fiecare dintre aceste linii, linia n-1 înmulţită cu coeficientul ai,n-1
(n-1), i = n-2, n-3,...,1. Repetând în acest mod procedura de eliminare obţinem în final următoarea matrice extinsă:
Ceea ce este demn de observat în matricea extinsă (2.10.) este că toţi coeficienţii de pe primele n coloane sunt egali cu 0, cu excepţia pivoţilor care sunt egali cu 1 şi în al doilea rând că matricea (2.10.), exprimând un sistem de ecuaţii, fiecare linie poate fi interpretată sub forma:
altfel spus, coloana n+1 conţine soluţia finală a sistemului. În metoda Gauss - Jordan cea de-a doua etapă care aduce sistemul la forma (2.10.) poartă denumirea de eliminarea inversă. Spre deosebire de metoda Gauss, în metoda Gauss - Jordan eliminarea directă şi eliminarea inversă pot fi aplicate simultan, fiecare pivot putând fi utilizat pentru eliminarea atât a coeficienţilor situaţi sub el, cât şi a coeficienţilor situaţi deasupra lui în matricea extinsă a sistemului.
2.1.4 2.1.3. Pivotarea şi eliminarea Gauss standard În subcapitolul anterior am aplicat eliminarea Gauss în cazul unor sisteme liniare ideale. De ce ideale ş
17
2.SISTEME DE ECUAŢII LINIARE
În primul rând deoarece primul coeficient din linia întâi şi nici unul dintre coeficienţii care s-au obţinut pe diagonala principală în procesul de eliminare directă nu au fost egali cu zero, ceea ce ar fi condus la o împărţire la zero. A doua situaţie delicată poate să apară în cazul în care unul dintre pivoţi are o valoare foarte mică comparativ cu restul elementelor situate pe aceeaşi coloană, sub pivotul respectiv. Acest caz poate conduce la erori considerabile. O ilustrare a cestei dificultăţi se regăseşte în exemplul următor:
2.1.4.1 Exemplul 2.1.3. Pivoţi foarte mici Sistemul liniar
admite soluţia exactă x1 = 10 şi x2 = 10 Să rezolvăm acest sistem prin metoda eliminării Gauss. Vom lucra cu patru zecimale exacte. Primul pivot este a11
(1) = 0,001 şi multiplicatorul corespunzător este:
Scăzând prima linia înmulţită cu m21, din a doua linie obţinem:
Substituţia inversă generează soluţia Se observă în exemplul anterior că eroarea este de acelaşi ordin ca şi soluţia. În condiţiile în care nu a existat o acumulare de erori ca urmare a mii de operaţii aritmetice şi matricea A nu este cvasisingulară, care a fost cauza care a generat o eroare aşa de mare ? Dificultatea apare ca urmare a alegerii unui pivot foarte mic la primul pas al eliminării. În consecinţă, multiplicatorul corespunzător este 104 şi ecuaţia a doua rezultată în urma eliminării directe conţine coeficienţi de 103 ori mai mari decât cei din sistemul original. Eroarea obţinută în determinarea necunoscutei x2 este 0,0008. Această eroare s-a mărit foarte mult la calculul necunoscutei x1 ca urmare a ordinului de mărime în care s-a lucrat. Dacă în exemplul anterior schimbăm între ele cele două ecuaţii, soluţia care o vom obţine prin metoda Gauss va coincide cu soluţia exactă. Aceasta conduce la următoarea regulă: dacă multiplicatorii sunt mai mici sau cel mult egali cu 1. eliminarea Gauss devine stabilă şi va genera o soluţie la fel de precisă ca cea obţinută prin aplicarea oricărui alt algoritm. Păstrarea multiplicatorilor (în valoare absolută) mai mici decât 1 se face printr-un proces de pivotare: la pasul k al eliminării directe se alege ca pivot cel mai mare element (în valoare absolută) de pe coloana k şi liniile k,...,n. Linia conţinând acest element se schimbă cu linia k rezultând astfel pivotul akk
(k) maxim. Aceeaşi schimbare trebuie aplicată şi elementelor din vectorul termenilor liberi b. Necunoscutele x nu trebuie reordonate deoarece coloanele matricii A nu se schimbă între ele în cadrul procesului de pivotare. Să exemplificăm aceste afirmaţii:
18
2.1. Metode directe
2.1.4.2 Exemplul 2.1.4. Pivotarea în metoda Gauss Să se rezolve sistemul: Matricea extinsă a sistemului este:
Eliminarea elementelor de pe prima coloană şi liniile 2 şi 3 nu se poate face dacă luăm în considerare matricea extinsă scrisă în această formă deoarece primul element de pe prima linie este 0 ceea ce ar conduce la o împărţire imposibilă. La prima pivotare vom schimba intre ele liniile unu şi trei deoarece linia a treia conţine pe prima poziţia valoarea maximă pe prima coloană. După pivotare matricea extinsă devine:
În continuare eliminăm primul element de pe a doua linie scăzând din a doua linie prima linie înmulţita cu 1/2: primul element de pe linia a treia este egal cu zero aşa că vom trece la eliminarea celui de al doilea element de pe linia a treia. Acest element, 10, este mai mare decât al doilea element de pe linia a doua, motiv pentru care vom aplica a doua pivotare schimbând între ele liniile doi şi trei. Rezultă:
După ce eliminăm al doilea element de pe coloana a treia (prin înmulţirea celei de-a doua linii cu 1/10) am încheiat etapa eliminării directe: Prin substituţie inversă obţinem următoarea soluţie a sistemului:
Următorul algoritm implementează metoda eliminării Gauss cu pivotare. În general în rezolvarea sistemelor liniare complexe pivotarea nu rezolvă, toate problemele, fiind indicat să se lucreze şi în dublă precizie. 2.1.5 Algoritmul eliminării Gauss standard cu pivotare Rezolvă un sistem liniar de n ecuaţii cu n necunoscute DATE DE INTRARE:
• aij coeficienţii sistemului. i, j = 1,...,n; • bi termenii liberi
DATE DE IEŞIRE: • soluţia x1, x2,.., xn sau mesaj "Sistemul nu are soluţie unică"
(Eliminarea directă:) Pasul 1. PENTRU k = 1,...,n REPETĂ (k - pasul curent al eliminării directe) Pasul 2. CAUTĂ aik maxim PENTRU i = k,...,n
19
2.SISTEME DE ECUAŢII LINIARE
Pasul 3. DACĂ aik = 0 ATUNCI: SCRIE "Sistemul nu are soluţie unică" STOP.
Pasul 4. DACĂ akk ≠ max(aik)k = 1…n ATUNCI PERMUTĂ linia k cu linia ce conţine max(aik)
Pasul 5. DACĂ k ≠ 1 Pasul 6. PENTRU i = k+1,...,n REPETĂ (i - indice de linie) Pasul 7. PENTRU j = k,...,n+1 CALCULEAZĂ (j - indice coloană)
Pasul 9. PENTRU i = n-1,n-2, ..., 1 REPETĂ SUMA = 0 Pasul 10. PENTRU j = i+1, i+2, .. . ,n CALCULEAZĂ SUMA = SUMA+aijxi Pasul 11. CALCULEAZĂ xi = (ai,n+1 - SUMA)/aii Pasul 12. SCRIE (x1,x2, . . . ,xn) STOP
2.2 2.1.4. Operaţii matriciale
Am introdus în capitolele anterioare conc^pr.u.: at ;v. .. .., ea o metodă convenţională de exprimare şi opeis-f a:n!'Ti.t,i.i liniare. Ne vom ocupa în acest subcapitol de câteva noţiuni referitoare la operaţiile algebrice permise în calculul matricial şi vom descrie un algoritm care permite calculul inversei unei matrici. Să definim mai întâi dimensiunea unei matrici. Considerând o matrice cu m linii şi n coloane vom spune că este dreptunghiulară de dimensiune m x n. În cazul în care numărul liniilor este egal cu numărul coloanelor, de exemplu n, matricea respectivă este pătrată, de dimensiune n x n. Două matrici A şi B sunt egale dacă sunt de aceeaşi dimensiune, de exemplu m x n, şi dacă aij = bij, pentru orice: i = 1,...,m şi j = 1,...,n. Dintre operaţiile importante în calculul matricial ne vom ocupa pentru început de adunarea (scăderea) şi înmulţirea matricilor.
2.2.1 a) Adunarea şi scăderea. Două matrici de aceeaşi dimensiune pot fi adunate (scăzute) prin însumarea, respectiv scăderea, tuturor elementelor de pe aceleaşi poziţii din cele două matrici: A ± B = C unde C este o matrice cu fiecare dintre elemente dat de cij = aij + bij Introducând şi noţiunea de transpusă a unei matrici, care se obţine schimbând în matricea respectivă liniile cu coloanele, vom putea lua în considerare câteva dintre proprietăţile adunării şi scăderii matricilor. Fie matricile A,B şi C de dimensiune m x n. Adunarea şi scăderea matricilor îndeplinesc următoarele proprietăţi:
A + B = B + A (A ± B)T = AT ± BT (A + B) + C = A + (B + C) A+0 = 0+A = A (2.11.)
20
2.1.4. Operaţii matriciale
A + (-A) = -A + A = 0
2.2.2 b) Înmulţirea Orice matrice poate fi înmulţită cu un număr real, rezultatul obţinându-se prin înmulţirea fiecărui element al matricii cu numărul real dat. Făcând legătura dintre proprietăţile adunării şi scăderii matricilor şi înmulţirea unei matrici cu un scalar vom găsi următoarele proprietăţi: Fie matricile A,B şi C de dimensiune mxn şi numerele reale a şi b. Proprietăţile adunării şi scăderii matricilor pot fi suplimentate cu: a(A±B) = aA + aB
(a ± b)A = aA ± bA (2.12.) a(bA) = (ab)A În general înmulţirea a două matrici se efectuează după următoarea definiţie: Fie matricea A de dimensiune m x n şi matricea B de dimensiune n x p. Produsul matricial AB este o matrice C de dimensiune m x p ale cărei elemente sunt date de relaţia:
cij = ∑k=1
naik⋅bkj
pentru fiecare i = 1,2,...,m şi j = 1,2,...,p. Fiecare element cij se obţine prin înmulţirea elementelor de pe linia i a matricii A cu elementele corespondente de pe coloana j a matricii B urmată de însumarea acestor produse. Ţinând cont de relaţia de definiţie, rezulta că dacă matricea A este de dimensiune m x n şi matricea B este de dimensiune n x m, se pot defini produsele AB şi BA dar AB ≠ BA. Această observaţie rămâne valabilă chiar şi în situaţia în care cele două matrici sunt pătrate (m = n). Înmulţirea matricilor îndeplineşte următoarele proprietăţi:
(A + B)C = AC + BC; C(A*B) = CA + CB; A(BC) = (AB)C Un caz particular de matrici îl constituie vectorii coloană şi vectorii linie (îi vom nota cu litere mici).
Un vector coloană este o matrice având o singură coloană. De exemplu, x = ⎝⎜⎜⎛
⎠⎟⎟⎞x1
x2x3
, y = ⎝⎜⎜⎛
⎠⎟⎟⎞y1
y2y3
.
Ordinul unui vector coloană este dat de numărul de elemente pe care îl conţine. Vectorii x şi y daţi de (2.1.14.) sunt de ordinul trei. Un vector linie poate fi privit ca o matrice cu o singură linie. Ordinul său este dat de numărul de elemente pe care îl conţine. De exemplu, vectorul linie a, a = (a1, a2, a3, a4) este de ordinul patru. Vectorii linie şi coloană fiind cazuri particulare de matrici, toate regulile şi proprietăţile matricilor sunt valabile şi în cazul vectorilor.
Adunarea vectorilor este definită de relaţia: x ± y = z unde x, y şi z sunt vectori de acelaşi tip şi ordin iar elementul i al vectorului z este dat de xi ± yi = zi Prin înmulţirea unei matrici cu un vector coloană obţinem un vector coloană Ax = y unde elementul yi al vectorului coloană y poate fi scris explicit ca fiind yi = ∑
k=1aikxk
21
2.SISTEME DE ECUAŢII LINIARE
2.2.3 2.1.5. Inversa unei matrici După cum am amintit în introducerea acestui capitol putem afirma că o matrice A, de dimensiune n x m, este nesingulară dacă există o matrice A-1, cu proprietatea AA-1 = A-1A = I. Matricea A-1 poartă denumirea de inversa matricii A. O matrice care nu admite inversă poartă denumirea de matrice singulară. Matricea care am notat-o cu I este matricea unitate ale cărei elemente sunt zero, cu excepţia elementelor de pe diagonala principală care sunt egale cu unu. Putem defini matricea unitate de ordinul n, In ca fiind matricea ale cărei elemente sunt date de relaţia:
aij = ⎩⎪⎨⎪⎧1 pentru i=j0 pentru i≠j
Inversa unei matrici poate fi calculată aplicând eliminarea Gauss - Jordan. Să considerăm sistemul: Ax = y (2.15.) unde A este o matrice pătrată iar x şi y sunt vectori coloană. Presupunând că nu este necesară pivotarea dacă înmulţim ecuaţia (2.15.) cu matricea pătrată G, având aceeaşi dimensiune cu A, obţinem: GAx = Gy (2.16.) Dacă G este aleasă astfel încât să fie inversa matricii A, ecuaţia (2.16.) se reduce la forma: x=A-1y (2.17.) care este chiar soluţia sistemului. Cu alte cuvinte, eliminarea Gauss-Jordan este echivalentă cu o amplificare G = A-1. Din acest motiv dacă vom aplica matricii unitate I, aceleaşi operaţii ca şi în eliminarea Gauss-Jordan, matricea unitate se va transforma în A-1. Aceasta poate fi scrisă simbolic: GI = A-1 (2.18.) Pentru a calcula A-1, vom scrie matricea extinsă a matricilor A şi I sub forma:
⎝⎜⎜⎛
⎠⎟⎟⎞a11 a12 a13 1 0 0
a21 a22 a23 0 1 0a31 a32 a33 0 0 1
(Am considerat matricea A de dimensiuni 3x3)
în continuare vom aplica eliminarea Gauss-Jordan în acelaşi mod ca şi la rezolvarea unui sistem de ecuaţii liniare. Când jumătatea din stânga a matricii extinse se va reduce la matricea unitate, jumătatea din dreapta a matricii extinse va deveni A-1. Ex 2.1.5 Inversa unei matrici Să se calculeze inversa matricii A. Vom scrie matricile A şi I în forma extinsă:
În prima etapa se aplici eliminarea directă. Se scade prima linie înmulţită cu (-1/2) din a doua linie şi prima linie înmulţită cu 3/2 se scade din a treia linie:
Se scade linia a doua înmulţită cu -0,5/3,5 = -1/7 din a treia linie: Următoarea etapă aste eliminarea inversă. Se împarte ultima linia la 1,5714:
22
2.1.4. Operaţii matriciale
Continuând eliminarea inversă conform algoritmului Gauss - Jordam obţinem în final:
Cele trei coloane din jumătatea din dreapta a matricii extinse constituia inversa matricii A. Se poată verifica prin calcul că AA-1 = A-1A = I:
şi Inversa unei matrici poate fi calculata şl prin metoda Gauss. Programul 2.1. implementează acest algoritm pentru calculul inversei unei matrici.
2.2.4 2.1.6. Determinantul unei matrici În cazul unei matrici pătrate, A, este posibil să se calculeze determinantul det(A) al acestei matrici. Fie matricea A de dimensiuni n x n. Dacă n = 1, prin definiţie, det(A) = a11. Înainte de a descrie modul de calcul al unui determinant să dăm câteva definiţii utile. Dacă la o matrice pătrată A, se elimină linia i şi coloana j, se obţine tot o matrice pătrată, de dimensiune (n-1) x (n-1). Determinantul acestei matrici poartă denumirea de minorul matricii A şi se notează det(Mij),. Cofactorul matricii A se notează Aij
c şi este definit de relaţia: Aij
c = (-1)i+j det(Mij) (2.19.) Determinantul matricii A este dat de formula
det(A) = ∑j=1
naij⋅Aij
c pentru orice i = 1,2,…,n (2.20) sau det A = ∑i=1
naij⋅Aij
c pentru orice i = 1,2,…,n
(2.21)
Ca exemplu să considerăm matricea A = ⎝⎛
⎠⎞a11 a12
a21 a22
Minorii matricii A sunt: det(M11) = a22 ; det(M12) = a21 ; det(M21) = a12 ; det(M22) = a11 (2.22.) Să aplicăm formula (2.20.) pentru i = 1; conform relaţiilor (2.19.) şi (2.22.) rezultă:
det(A) = ∑j=1
2a1j⋅A1j
c = a11⋅A11c + a12⋅A12
c = a11(-1)1+1det(M11) + a12(-1)1+2det(M12) = a11a22 - a12a21
Acelaşi rezultat se obţine şi în cazul în care se aplică formula (2.20.) pentru i = 2. În mod asemănător, det(A) se poate calcula şi prin aplicarea relaţiei (2.21.), pentru j = 1 sau j = 2, rezultatul fiind acelaşi. Formulele (2.20.) şi (2,21) permit să se stabilească o serie de proprietăţi ale determinanţilor: Un determinant are valoarea zero dacă: o linie (o coloană) conţine toate elementele egale cu zero sau, două linii (coloane) sunt proporţionale.
23
2.SISTEME DE ECUAŢII LINIARE
Dacă o linie (coloană) este amplificată cu o constantă reală şi valoarea determinantului se va amplifica cu acea constantă. Este de asemenea important să amintim că: Dacă toate elementele unei linii (coloane) sunt înmulţite cu un factor oarecare şi adunate la altă linie (coloană), valoarea determinantului nu se schimbă. Dacă două linii (coloane) se schimbă între ele, determinantul îşi schimbă semnul.
det(AB) = det(A) det(B), în timp ce det(A+B) ≠ det(A) + det(B)
2.2.5 2.1.7. Matrici particulare Ne vom referi pentru început la matricile triunghiulare pe care le-am întâlnit deja la metoda eliminării Gauss. Există două tipuri de astfel de matrici: superior triunghiulare şi inferior triunghiulare. O matrice superior triunghiulară U (upper), de dimensiune n x n cu elementele
uij = 0 pentru i = j+1, j+2, ... , n şi j = 1, 2.... , n (adică 1≤j<i≤n) O matrice inferior triunghiulară L (lower), de dimensiune n x n
lij = 0 pentru i = 1, 2.... ,j-1 şi j = 1, 2.... , n (adică 1≤i<j≤n) (Abrevierile consacrate L, U provin din limba engleză L - lower şi U - upper) O proprietate a matricilor superior şi inferior triunghiulare este că determinantul acestor matrici este dat de produsul elementelor de pe diagonala principală:
det(L) = ∏i=1
nlii şi det(U) = ∏
i=1
nuii (2.25)
Ţinând cont de această proprietate rezultă că problema calculului determinantului unei matrici poate fi simplificată prin reducerea mai întâi a matricii la formă triunghiulară şi apoi prin aplicarea uneia dintre formulele (2.25.). De exemplu, matricea
prin aplicarea eliminării Gauss se reduce la următoarea formă superior triunghiulară:
Considerând a doua formulă (2.25.) obţinem: det(A = det(Ã(3) = 10⋅2.5⋅6.2 = 155 Dacă în procesul de eliminare are loc o interschimbare a liniilor pentru a se obţine la fiecare pas pivotul cu valoarea maxim posibilă, se va ţine cont de numărul interschimbărilor în stabilirea semnului determinantului: un număr par de interschimbări nu schimbă semnul, în timp ce un număr impar schimbă semnul determinantului matricii triangularizate. O teoremă importantă în algebra matricială, ale cărei consecinţe la vom detalia în capitolul următor, este următoarea: Teoremă: Dacă putem aplica sistemului liniar Ax = b procedura eliminării Gauss, fără interschimbarea liniilor, atunci matricea A poate fi descompusă într-un produs dintre o matrice inferior triunghiulară L şi o matrice superior triunghiulară U, A = LU, unde elementele lij şi uij ale matricilor L şi U sunt definite, pentru orice j, de relaţiile:
uij = ⎩⎨⎧aij
(i) pentru i≤j0 pentru i>j şi lij =
⎩⎪⎨⎪⎧0 pentru i<j1 pentru i=jmij pentru i>j
(2.26.)
24
2.1.4. Operaţii matriciale
unde aij(i) sunt elementele matricii triangularizate în urma procesului de eliminare, iar mij sunt
multiplicatorii. Dacă în procesul eliminării Gauss este necesară o interschimbare a liniilor, matricea A poate fi descompusă într-un produs LU, unde elementele matricii U sunt date de relaţia (2.26.) dar matricea L nu mai este inferior triunghiulară. Pentru a păstra matricea L într-o formă inferior triunghiulară se introduce o nouă matrice P, care poartă denumirea de matricea permutărilor. Matricea permutărilor se obţine plecând de la matricea unitate ale cărei elemente sunt egale cu 0, cu excepţia elementelor de pe diagonala principală care sunt egale cu unu. Dacă într-o matrice A de dimensiune n x n are loc interschimbarea liniilor i şi j, această interschimbare este echivalentă cu înmulţirea PA, unde P este matricea permutărilor formată prin interschimbarea liniilor i şi j din matricea In. De exemplu:
unde în matricea permutărilor s-a înregistrat interschimbarea liniilor 2 şi 3. Se poate demonstra uşor că P-1 = PT. Dacă interschimbarea liniilor operată în procesul eliminării Gauss se aplică şi matricii unitate I, rezultă matricea permutărilor, P, cu proprietatea LU = PA. Aceasta implică A = (P-1L)U, unde L este inferior triunghiulară. În subcapitolul 2.1.4. am introdus noţiunea de transpusa unei matrici. O matrice identică cu transpusa sa este o matrice simetrică. De exemplu,
Deoarece A = AT, rezultă că matricea A este simetrică. Un alt tip de matrici des întâlnite în practică îl constituie matricile bandă.
O matrice de dimensiune nxn este o matrice bandă dacă, există doi întregi p şi q (p, q ∈(1,n)) cu proprietatea că aij = 0 pentru i+p ≤ j şi j+q ≤ i. Lăţimea benzii pentru o astfel de matrice este l = p+q-1. Un caz special de matrice bandă care intervine în practică obţine pentru p = q = 2. Aceste matrici (numite tridiagonale) având lăţimea benzii l = 2+2-1 = 3 se întâlnesc în studiul interpolării spline şi a ecuaţiilor diferenţiale ordinare cu condiţii pe frontieră. Au următoarea formă:
O altă clasă importantă de matrici este dată de matricile strict dominante diagonal, care poartă această denumire deoarece au proprietatea:
25
2.SISTEME DE ECUAŢII LINIARE
|aij| > ∑j=1;j≠i
n |aij| pentru orice i = 1,2,...,n.
Dacă matricea A, de dimensiuni n x n este strict dominantă diagonal, atunci A este şi nesingulară şi se poate aplica eliminarea Gauss pentru rezolvarea oricărui sistem Ax = b, fără interschimbări de linii. Ultima clasă de matrici particulare care o vom analiza în acest capitol este clasa matricilor pozitiv definite. Să asociem unei matrici pătrate A de dimensiuni n x n, produsul xTAx unde x este un vector coloană de dimensiune n, x ≠ 0. Dacă xTAx > 0 pentru orice x ≠ 0 vom spune că matricea A este pozitiv definită.
Dacă o matrice pătrată A este pozitiv definită atunci det(A) ≠ 0, ceea ce înseamnă că matricea A este nesingulară. Reciproca nu este valabilă, altfel spus, dacă det(A) ≠ 0 nu putem concluziona că A este pozitiv definită.
2.2.6 2.1.8. Descompunerea LU An văzut în capitolul anterior ca o matrice pătrată A poate fi descompusă într-un produs de două matrici particulare: PA = LU, unde P este matricea permutărilor, L este o matrice inferior triunghiulară iar U superior triunghiulară. Să analizăm întâi poate fi utilizată această tehnică la rezolvarea sistemelor de ecuaţii liniare şi care sunt avantajele ei. Fie sistemul liniar Ax = b, pe care dorim să îl rezolvăm.
• La pasul 1 se descompune matricea A într-un produs de forma PA = LU, descompunere care se poate obţine fără a lua în considerare termenul liber b.
• La pasul 2 se reordonează elementele vectorului liber b în funcţie de interschimbările de linii care s-au efectuat la pasul 1. Această reordonare se face prin înmulţirea matricii permutărilor cu vectorul coloană b: Pb.
• La pasul 3 se rezolvă sistemul liniar simplificat Lz = b.
• La pasul 4 se rezolvă sistemul liniar Ux = z obţinând în acest fel soluţia finală, x. Care este avantajul acestei metode ? Să presupunem că avem un set de sisteme de ecuaţii liniare de forma Ax = bi, i = 1, 2, ..., m, în care matricea A este aceeaşi şi numai termenul liber se schimbă. Aplicând algoritmul descris anterior se observă că este suficient să aplicăm o singură dată descompunerea PA = LU, în continuare sistemele rezolvându-se prin parcurgerea paşilor 2-4 de câte ori este necesar. La prima vedere s-ar părea că avem de-a face cu un nou algoritm. În realitate această metodă este eliminarea Gauss scrisă sub formă matricială. Să detaliem fiecare pas în parte, considerând pentru simplificarea calculelor un sistem liniar de trei ecuaţii cu trei necunoscute. Descompunerea LU pentru matricea A, de dimensiune 3x3 poate fi scrisă sub forma (fără pivotare):
26
2.1.4. Operaţii matriciale
Pentru a determina elementele necunoscute lij, şi uij, este necesar să se impună trei condiţii suplimentare pe care trebuie să le îndeplinească elementele matricii L sau ale matricii U. Vom considera în acest caz metoda Doolittle care impune condiţia ca elementele de pe diagonala principală ale matricii L să fie egale cu unitatea, l11 = l12 = l13 = 1. Astfel, descompunerea (2.29.) devine
Pentru a afla elementele lij, şi uij, fără pivotare, vom înmulţi mal întâi prima linie din L cu cele trei coloane din U, de unde rezultă că prima linie din D este identică cu prima linie a matricii A:
u11 = a11, u12 = a12, u13 = a13
înmulţind succesiv liniile doi şi trei din L cu prima coloană din U obţinem a21 = l21⋅u11, a31 = l31⋅u11 de unde
l21 = a21u11
; l31 = a31u11
Înmulţind linia a doua din L cu coloana a doua respectiv a treia din U rezultă
a22 = l21u12 + u22 ⇒ u22 = a22 - l21u12
a23 = l21u13 + u23 ⇒ u23 = a23 - l21u13 înmulţim linia a treia din L cu coloana a doua din U:
de unde
În final, u33 se obţine din înmulţirea ultimei linii din L cu ultima coloană din U:
de unde Generalizând pentru o matrice A de dimensiune nxn, descompunerea LU se obţine prin aplicarea următoarelor relaţii: 1. a) prima linie din U este dată de relaţia: u1j = a1j, j = 1,2…,n (2.35)
1. b) prima coloană din L este dată de relaţia: li1 = ai1u11
, i = 1,2…,n (2.36)
2. a) elementele uii de pe diagonala principală din U sunt date de relaţia: uii = aii - ∑k=1
i-1lik⋅uki
2. b) linia i din U (i = 2,3,...,n-1) este data de relaţia: uij = aij - ∑k=1
j-1lik⋅ukj , j = i+1…n
2. c) coloana i din L (i = 2,3,...,n-1) este dată de relaţia: lji = ⎝⎜⎜⎛
⎠⎟⎟⎞aji - ∑
k=1
i-1ljk⋅uki
uii , j = i+1…n
3. elementul un,n este dat de relaţia: unn = ann - ∑k=1
n-1lnk⋅ukn
27
2.SISTEME DE ECUAŢII LINIARE
Următorul algoritm incorporează metoda descompunerii LU laolaltă cu o procedură de pivotare, substituţia directă şi substituţia inversă necesare în rezolvarea unui sistem de ecuaţii liniare.
2.2.7 Algoritmul descompunerii LU Rezolvă un sistem liniar de n ecuaţii cu n necunoscute, Ax = b, prin descompunerea matricii A în produsul LU urmată de rezolvarea sistemelor Lz = b şi Ux = z (se consideră lii = 1, i = 1,2,...,n) DATE DE INTRARE:
• n - dimensiunea sistemului • aij - elementele metricii A, i şi j = 1/2,.. .,n • bi - elementele vectorului coloană b, i = 1,2,...,n
DATE DE IEŞIRE: • x1, x2, ...,xn: soluţia sistemului sau mesaj "Sistemul nu are soluţie unică"
Pasul 1.(Caută pivotul) Fie p cel mal mic întreg astfel încât 1 ≤ p ≤ n şi |ap1| = max1≤i≤n|ai1| DACĂ |ap1| = 0 ATUNCI AFIŞEAZĂ mesaj: "Sistemul nu are soluţie unică" şi EXIT.
Pasul 2. DACĂ p ≠ 1, SCHIMBĂ între ele liniile p şi 1 din A şi P (P iniţial este In.) Pasul 3. l11 = 1, u11 = a11 (Condiţie rezultată din a11 = l11u11) Pasul 4. PENTRU j = 2,3,...,n: u1j = a1j, (Prima linie din U); lj1 = aj'1/u11 (Prima coloană din L) Pasul 5. PENTRU i = 2,..., n-1, EFECTUEAZĂ paşii 6-9: Pasul 6. CAUTĂ pivotul i Fie p cel mai mic întreg astfel încât i ≤ p ≤ n şi
⎪⎪⎪⎪
⎪⎪⎪⎪api - ∑
k=1
i-1lpk⋅uki = maxi≤j≤n
⎪⎪⎪⎪
⎪⎪⎪⎪aji - ∑
k=1
i-1ljk⋅uki
DACĂ maximum este zero ATUNCI SCRIE mesaj: "Sistemul nu are soluţie unică" STOP.
Pasul 7. DACĂ p ≠ i SCHIMBĂ între ele liniile p şi i din A, L şi P.
n-1lnk⋅ukn . DACĂ unn = 0 SCRIE mesaj: "Sistemul nu are soluţie unică".
STOP. (Rezolvă sistemul inferior triunghiular Lz = b) Pasul 11. REARANJAŢI elementele vectorului coloană liber b: b = Pb
28
2.1.4. Operaţii matriciale
Pasul 12. z1 = b1
Pasul 13. PENTRU i = 2,...,n CALCULEAZĂ zi = bi - ∑j=1
i-1lij⋅zj
(Rezolvă sistemul superior triunghiular Ux = z) Pasul 14. xn = zn/unn
Pasul 15. PENTRU i = n-1, n-2,...,1, CALCULEAZĂ xi = zi - ∑j=i+1
nuij⋅xj
Pasul 16. SCRIE (x1, x2,...,xn) STOP După cum se observă din algoritm, în cazul în care se schimbă sistemul prin modificarea numai a termenului liber b, descompunerea matricii A rămâne aceeaşi, calculele reluându-se de la Pasul 11. În cazul în care matricea A este pozitiv definită, se poate aduce o îmbunătăţire algoritmului descompunerii LU, astfel Încât să reducem numărul de operaţii aritmetice. ………
29
2.SISTEME DE ECUAŢII LINIARE
2.3 2.2. Metode iterative
2.3.1 2.2.1. Introducere În diverse aplicaţii rezultă sisteme de ecuaţii liniare la care matricea coeficienţilor, A, conţine multe elemente cu valoarea zero. Rezolvarea unor astfel de sisteme se poate efectua deseori aplicând tehnici iterative, avantajul acestor metode faţă de cele directe constând într-un efort de calcul redus şi eliminarea erorii rezultată la un pas oarecare la pasul următor (în cazul unor iteraţii convergente). Ne vom ocupa în acest capitol de cele mai cunoscute tehnici iterative. Pentru o corectă Înţelegere a convergenţei secvenţei de vectori care rezultă intr-o tehnică iterativă către soluţia exactă a unui sistem vom analiza pentru început noţiunile de normă vectorială şi normă matricială. Vom descrie în continuare tehnicile iterative de rezolvare a sistemelor liniare. În final vom aborda problema estimării erorii şi a condiţionării sistemului pe care dorim să îI rezolvăm indiferent de metoda pe care o alegem.
2.3.2 3.2.2. Norme vectoriale şi matriciale Ca metodă de comparaţie a numerelor suntem obişnuiţi să utilizăm noţiunea matematică de valoare absolută: |8|>|3| sau |-5|>|-2,1| Pentru a măsura mărimea generală a elementelor unui vector, altfel spus, pentru a defini o distanţă în Rn vom utiliza noţiunea de normă vectorială. O normă vectorială definită în Rn este orice funcţie definită în Rn cu valori în R, notată || ||, care satisface proprietăţile:
||x|| > 0 pentru orice x∈Rn, x≠0 ||x|| = 0 dacă şi numai dacă x = (0,0, .. .,0)T = 0 (2.43.)
||cx|| =|c|-||x||, pentru orice scalar c∈R şi vector x∈Rn
||x+y||≤||x||+||y||, pentru orice x,y∈Rn Orice vector definit în Rn este un vector coloană. Din motive practice, atunci când vom reprezenta un vector în funcţie de componentele sale, vom utiliza transpusa unui vector. De
exemplu, vectorul x ∈ Rn, x =
⎝⎜⎜⎛
⎠⎟⎟⎞x1
x2x3…xn
îl vom scrie x = (x1, x2, …, xn)T.
Normele vectoriale cel mai dea utilizate sunt norma 1 || ||1, norma a doua || ||2, cunoscută şi ca lungimea euclidiană şi norma infinită || ||∞ denumită şi max-norma. Pentru un vector oarecare, x∈Rn, x = (x1, x2, …, xn)T, acesta norme sunt definite de relaţiile:
||x||1 = ∑i=1
n|xi| ; ||x||2 =
⎝⎜⎛
⎠⎟⎞
∑i=1
n|xi|2
1/2 ; ||x||∞ = max1≤i≤n|xi|
……… Un alt concept utilizat în tehnicile iterative este cel de normă matricială. O normă matricială ||A|| , poate fi definită ca şi o normă vectorială şi satisface proprietăţile: ||A|| ≥ 0
30
2.2. Metode iterative
||A|| = 0 dacă şi numai dacă matricea A are toate elementele egale cu zero ||cA|| = |c||A||; ||A+B|| = ||A||+||B||; ||AB|| ≤ ||A|| ||B||
unde A,B ∈ Rnxn În acest curs ne vom referi numai la norme matriciale naturale care sunt definite după cum urmează: Teoremă: Dată fiind o normă vectorială oarecare || || în Rn, atunci ||A|| = max||x||=1||Ax|| (2.51.)
defineşte o normă matricială naturală pentru A ∈ Rnxn. Normele matriciale naturale pe care le vom lua în considerare sunt:
||A||∞ = max||x||∞=1||Ax|| şi ||A||2 = max||x||2=1||Ax||
Norma matriciala infinită are următoarea reprezentare raportată la elementele matricii A: Teoremă: Norma infinită a unei matrici A, de dimensiune n x n este dată de relaţia:
||A||∞ = max1≤i≤n ∑j=1
n|aij|
2.3.3 2.2.3. Metoda Jacobi şi metoda Gauss - Seidel Vom descrie în acest subcapitol două dintre cele nai cunoscute tehnici iterative de rezolvare a sistemelor liniare de forma Ax = b unde A aste o matrice reală de dimensiune n x n. Orice tehnică iterativă de rezolvare a unui asemenea sistem porneşte de la o aproximaţie iniţială x(0), pentru soluţia x a sistemului şi generează un şir de vectori (x(k))k=0
∞ care converge către soluţia exactă x. Majoritatea tehnicilor iterative implică un proces care transformă sistemul Ax = b în forma echivalentă x = Tx + c, unde matricea T este reală de dimensiune n x n iar c este un vector coloană de ordin n. După alegerea vectorului iniţial x(0) şirul de vectori care converge către soluţia x este generat de relaţia recursivă: x(k) = Tx(k-1) + c pentru k = 1,2,3,... Tehnicile iterative sunt recomandabile în cazul sistemelor liniare de mici dimensiuni sau al sistemelor de mari dimensiuni a căror matrice A conţine multe elemente egale cu zero. Astfel de sisteme rezultă de obicei în problemele cu condiţii iniţiale pentru ecuaţii diferenţiale ordinare sau la ecuaţiile diferenţiale cu derivate parţiale. Pentru simplificarea explicaţiilor să considerăm sistemul de 4 ecuaţii cu 4 necunoscute în forma generală:
Pentru a aduce sistemul Ax = b la forma echivalentă x = Tx+c vom rezolva fiecare ecuaţie i funcţie de xi, pentru i = 1,2,3,4. Obţinem în acest mod:
31
2.SISTEME DE ECUAŢII LINIARE
T şi c sunt:
Pentru aproximaţia iniţială x(0) = (x1
(0), x2(0), x3
(0), x4(0)) se va genera x(1) = (x1
(1), x2(1), x3
(1), x4(1))
prin relaţiile:
Următoarele iteraţii, x(k) = (x1
(k), x2(k), x3
(k), x4(k)) sunt generate într-un mod similar.
Criteriul de stopare a iteraţiilor succesive este dat de următoarea inegalitate care va opri iteraţiile la pasul k:
||x(k) - x(k-1)||∞||x(k)||∞
< ε, unde ε > 0 este eroarea maxim admisă şi se impune iniţial.
(unde ||x - y||∞ = max1≤i≤n|xi - yi| este "norma infinită a diferenţei vectoriale", iar ||x||∞ = max1≤i≤n|xi|) Metoda descrisă poartă denumirea de metoda iterativă Jacobi şi generalizând, constă în generarea fiecărei componente x(k), cu ajutorul componentelor x(k-1)
i, k ≤ 1 pe baza relaţiei:
xi(k) =
∑j=1; j≠i
n(-aijxj
(k-1)) + bi
aii , pentru i = 1…n, aii ≠ 0
Metoda poate fi scrisă sub forma x(k) = Tx(k-1) + c (2.60.)
prin descompunerea matricii A după cum urmează
32
2.2. Metode iterative
unde D este matricea diagonală care conţine elementele de pe diagonala principală a matricii A, -L este partea strict inferior triunghiulară a metricii A şi -V este partea strict superior triunghiulara a matricii A. Ecuaţia matricial vectorială Ax = b, echivalentă cu (D-L-U)x =b poate fi transformată în forma Dx = (L+U)x + b şi în final x = D-1(L+U)x + D-1b În acest mod, metoda Jacobi exprimată în formă matricial vectorială este:
x(k) = D-1(L+U)x(k-1) + D-1b, K = 1,2,… Practic, se utilizează în calcule relaţia (2.59.).
2.3.4 Algoritmul iterativ Jacobi Rezolvă un sistem liniar Ax = b cu o aproximaţie iniţială dată x(0) DATE DE INTRARE:
• numărul de ecuaţii şi necunoscute, n • elementele aij, i şi j = 1,2,...,n ale matricii A • elementele bi,, i = 1,2,...,n ale termenului liber b; • toleranţa TOL ; • numărul maxim de iteraţii M; • aproximaţia iniţială x1
(0), x2(0), …, xn
(0) DATE DE IEŞIRE:
• soluţia aproximativă x1, x2, …, xn, sau mesaj: "S-a depăşit numărul maxim de iteraţii" Pasul 1. k = 1 Pasul 2. CÂT TIMP k ≤ M EFECTUEAZĂ paşii 3-6: Pasul 3. PENTRU i = 1,2,...,n CALCULEAZĂ
xnoui =
- ∑j=1;j≠i
n(aij⋅xi) + bi
aii
Pasul 4. DACĂ ||xnou - x||∞||xnou||∞
≤ TOL ATUNCI SCRIE(xnou1,xnou2,..., xnoun,) (soluţia) STOP
(unde ||x - y||∞ = max1≤i≤n|xi - yi| este "norma infinită a diferenţei vectoriale", iar ||x||∞ = max1≤i≤n|xi|) Pasul 5. k = k+1 Pasul 6. PENTRU i = 1,2,...,n: xi = xnoui, Pasul 7. SCRIE("S-a depăşit numărul maxim de iteraţii") STOP (Procedura a eşuat)
33
2.SISTEME DE ECUAŢII LINIARE
La pasul 3 al algoritmului trebuie să fie îndeplinită condiţia aii ≠ 0. pentru orice i = 1,2,...,n. Dacă nu este îndeplinită această condiţie sistemul se va reordona. Pentru mărirea vitezei de convergenţă a soluţiei este indicat ca sistemul să fie ordonat astfel încât aii să fie cât mai mare. O analiză sumară a relaţiei (2.59.) ne permite o îmbunătăţire a metodei. În relaţia (2.59.) pentru a se calcula xi
(k) se utilizează xi (k-1). Deoarece, pentru i > 1, componentele x1
(k), x2(k), x3
(k), …xi-1(k)
sunt deja calculate şi se presupune că reprezintă o aproximaţie mai bună decât x1(k-1), x2
pentru i = 1,2,...,n. Aplicând relaţia (2.63.) pentru sistemul de 4 ecuaţii cu 4 necunoscute analizat la metoda Jacobi vom obţine următoarele relaţii care se vor aplica pentru k = 1,2,... şi aii ≠ 0:
Această tehnică poartă denumirea de metoda iterativă Gauss-Seidel. Pentru a o scrie sub formă matricial vectorială, să considerăm relaţiile (2.64.) în care vom înmulţi ambii membri cu aii, i = 1,2,3,4 şi vom trece în membrul din stânga componentele x(k):
Rezultă că în formă matricial vectorială metoda Gauss-Seidel este dată de ecuaţia: (D - L)x(k) = Ux(k-1) + b (2.66) sau x(k) = (D - L)-1Ux(k-1) + (D - L)-1b (2.67) pentru k = 1,2,... La fel ca şi în metoda Jacobi şi în metoda Gauss-Seidel se utilizează în aplicaţii relaţia (2.63.)
2.3.5 Algoritmul iterativ Gauss - Seidel Rezolvă un sistem liniar Ax = b cu o aproximaţie iniţială x(0) DATE DE INTRARE:
• numărul de ecuaţii şi necunoscute n • elementele aij, i şi j = 1,2,...,n ale matricii A • elementele bi ale vectorului liber • aproximaţia iniţială x1(0), x2(0), …, xn(0) • toleranţa impusă TOL • numărul maxim de iteraţii M.
34
2.2. Metode iterative
DATE DE IEŞIRE: • soluţia aproximativă x1, x2, …, xn, sau mesaj: "S-a depăşit numărul maxim de iteraţii"
Pasul 1. k = 1 Pasul 2. CÂT TIMP k ≤ M EFECTUEAZĂ paşii 3-6 Pasul 3. PENTRU i = 1,2,...,n CALCULEAZĂ
xnoui =
-∑j=1
i-1(aij⋅xnouj) - ∑
j=i+1
n(aij⋅xj) + bi
aii
Pasul 4. DACĂ ||xnou - x||∞||xnou||∞
≤ TOL ATUNCI
(unde ||x - y||∞ = max1≤i≤n|xi - yi| este "norma infinită a diferenţei vectoriale", iar ||x||∞ = max1≤i≤n|xi|) SCRIE (xnou1,..., xnoun) (soluţia) STOP. Pasul 5. k = k+1 Pasul 6. PENTRU i = 1,2,...,n: xi = xnoui. Pasul 7. SCRIE("S-a depăşit numărul maxim de iteraţii") STOP (Procedura a eşuat). Pentru a studia convergenţa tehnicilor iterative vom considera formula: x(k) = Tx(k-1) + c, k = 1, 2, … unde aproximaţia iniţială a vectorului soluţie, x(0) este arbitrară. O primă condiţie referitoare la convergenţa tehnicilor iterative este dată de următoarea teoremă:
Teoremă: şirul de vectori x(k) = Tx(k-1) + c, pentru k ≥ 1 şi c ≠ 0, converge către soluţia unică x = Tx + c,oricare ar fi aproximaţia iniţială x(0) ∈ Rn dacă şi numai dacă ρ(T) < 1; ρ(T) este raza spectrală a matricii T. O a doua condiţie de convergenţi a tehnicilor Jacobi şi Gauss - Seidel este dată de Teoremă: Dacă matricea A este strict dominantă diagonal atunci metodele Jacobi şi Gauss - Seidel generează" şirul (x(k))k
∞ = 0 convergent către soluţia unică a sistemului Ax = b, oricare ar fi aproximaţia iniţială x(0). în general nu putem afirma care dintre cele două metode, Jacobi şi Gauss - Seidel, este indicată pentru o problemă oarecare. În anumite condiţii răspunsul la această Întrebare este dat de următoarea teoremă:
Teoremă (Stein - Rosenberg): Dacă aij ≤ 0 pentru i ≠ j şi aij > 0 pentru i = 1,2,...,n, atunci una şi numai una dintre următoarele propoziţii este valabilă:
a) 0 < ρ(TG) < ρ(TJ) < 1
b) 1 < ρ(TJ) < ρ(TG)
c) ρ(TJ) = ρ(TG) = 0
d) ρ(TJ) = ρ(TG) = 1 unde TJ este matricea iteraţiilor asociată metodei Jacobi iar TG este matricea iteraţiilor asociată metodei Gauss - Seidel.
35
2.SISTEME DE ECUAŢII LINIARE
Pentru cazul particular dat de teorema Stein - Rosenberg observăm că dacă una dintre metode este convergentă, atunci ambele vor da convergenţă, metoda Gauss - Seidel fiind mai rapidă decât metoda Jacobi. Viteza de convergenţă a unei proceduri depinzând de raza spectrală a matricii iteraţiilor asociată metodei, o posibilitate de alegere a procedurii care va conduce la accelerarea convergenţei constă în alegerea metodei a cărei matrice asociată are raza spectrală minimă. În subcapitolul următor vom descrie o procedură de alegere optimă.
2.3.6 2.2.4. Metodele relaxării Înainte de a descrie procedura de selectare a metodei iterative optime pentru rezolvarea unui sistem liniar Ax = b să analizăm posibilitatea de a introduce o măsură pentru aprecierea diferenţei dintre soluţia aproximativă şi soluţia exactă.
Fie x~ ∈ Rn soluţia aproximativă a unui sistem liniar Ax = b. Există două posibilităţi de apreciere a preciziei soluţiei:
Eroarea: e = x - x~ (2.66.)
şi vectorul reziduu: r = b - Ax~ = A(x - x~ ) = Ae (2.69.) unde x este soluţia exactă a sistemului. Reziduul măsoară diferenţa dintre vectorul liber b şi ecuaţiile sistemului în care se substituie soluţia calculată. O analiză sumară a relaţiilor (2.68.) şi (2.69.) ne arată că, dacă matricea A nu este singulară, eroarea este egală cu zero dacă reziduul este zero şi reciproc. Această afirmaţie nu implică însă că eroarea şi reziduul vor avea acelaşi ordin de mărime. Vom analiza în detaliu această afirmaţie în subcapitolul următor. Ceea ce trebuie reţinut în această etapă este că ceea ce urmărim prin îmbunătăţirea tehnicilor iterative aste obţinerea unei soluţii aproximative în sensul unui reziduu cât mai mic. În tehnicile Jacobi şi Gauss - Seidel la fiecare pas al iteraţiilor rezultă un vector al soluţiei aproximative căruia îi putem asocia un vector reziduu. În ideea celor afirmate mai sus ne propunem să obţinem o metodă iterativă care să genereze o secvenţă de aproximări căruia să li corespundă un şir de vectori reziduu convergent către zero. Fie ri
(k) = (r1i(k), r2i
(k), …, rni(k))T vectorul reziduu corespunzător în metoda Gauss - Seidel
vectorului soluţiei aproximative (x1(k), x2
(k), …, xi-1(k), xi
(k-1),…, xn(k-1))T. Elementul m al vectorului
reziduu ri(k), este
sau
unde m= 1,2,...,n. în particular, componenta i a vectorului reziduu r(k) este
de unde în metoda Gauss - Seidel x(k)
i, este dat de formula
36
2.2. Metode iterative
astfel încât, ecuaţia (2.73.) poate fi rescrisă sub forma:
sau
relaţie care exprimi tehnica Gauss - Seidel într-o altă formă. Practic, modificând relaţia (2.74.) în forma
pentru anumite valori pozitive ale factorului ω se obţin creşteri semnificative ale convergenţei soluţiei. Metodele bazate pe ecuaţia (2.75.) poartă denumirea de metodele relaxării. Pentru valori ale factorului de accelerare ω cuprinse între 0 şi 1, 0 < ω < 1 procedurile poartă denumirea de metodele subrelaxării şi pot fi utilizate pentru a obţine convergenţa în cazul unor sisteme către nu sunt convergente prin metoda Gauss - Seidel.
în cazul unor valori ale factorului de accelerare ω > 1, procedurile poartă denumirea de metodele supra-relaxării şi se utilizează pentru accelerarea convergenţei sistemelor care sunt convergente prin metoda Gauss - Seidel. Aceste metode sunt cunoscute şi sub denumirea de suprarelaxări succesive cu abrevierea SRS sau SOR (Succesive - Over - Relaxation) şi sunt utile în soluţionarea unor sisteme generate în rezolvarea numerică a unor ecuaţii diferenţiale cu derivate parţiale. Ţinând cont de relaţiile (2.72.), ecuaţia (2.75.) poate fi reformulată sub forma:
Mu există o regulă generală care să permită alegerea unei valori optime pentru factorul de accelerare ω pentru un sistem liniar în forma sa generală. În anumite situaţii particulare răspunsul la această întrebare este dat de următoarele teoreme: Teoremă (Kahan)
Dacă aii ≠ 0 pentru i = 1,2,...,n, rezultă că
ceea ce implică (Tω) < 1 numai daci 0 < ω < 2. Tω este matricea iteraţiilor pentru metoda SRS şi este data de relaţia:
Teoremă (Ostrowakl - Reich): Dacă A este o matrice pozitiv definită şi 0 < ω < 2, metoda SRS este convergentă pentru orice aproximaţie iniţială a vectorului soluţie, x(0). Teoremă: Dacă A este o matrice pozitiv definită şi tridiagonală, atunci
şi valoarea optimă pentru factorul ω în metoda SRS este
37
2.SISTEME DE ECUAŢII LINIARE
pentru această valoare a factorului ω, (Tg) = -1.
2.3.7 Algoritmul SRS Rezolvă un sistem liniar Ax = b date fiind aproximaţia iniţială x(0) şi factorul de accelerare DATE DE INTRARE:
• numărul de ecuaţii şi necunoscute n • elementele aij, i şi j = 1,2,...,n ale matricii A • elementele bi ale vectorului liber b • aproximaţia iniţială (x1(0), x2(0), … xn(0)) • factorul de accelerare ω • toleranţa impusă TOL • numărul maxim de iteraţii M
DATE DE IEŞIRE:
• soluţia aproximativi (x1,...,xn) sau mesaj "S-a depăşit numărul maxim de iteraţii" Pasul 1. k = 1 Pasul 2. Cât timp k ≤ M efectuează paşii 3-6 Pasul 3. Pentru i = 1,2,...,n calculează
Pasul 4. Dacă ||xnou - x||∞ ||xnou||∞
≤ TOL atunci SCRIE (xnou,,... ,xnou.) (soluţia) STOP.
Pasul 5. k = k + 1 Pasul 6. Pentru i = 1,2,...,n: xi = xnoui Pasul 7. SCRIE ("S-a depăşit numărul maxim de iteraţii") STOP. (Procedura a eşuat)
38
2.3. Sisteme rău condiţionate
2.4 2.3. Sisteme rău condiţionate
Să considerăm un sistem liniar Ax = b şi x~ o aproximare a soluţiei exacte x a sistemului. La o analiză superficială suntem tentaţi să afirmăm că dacă vectorul reziduu r = b - Ax~ este mic, mai corect dacă norma ||r|| este mică, atunci şi norma erorii
||e|| = || x - x~ || este mică. În general aceasta afirmaţie este valabilă, dar există în practică şi sisteme care nu se supun acestei proprietăţi. Exemplul 3.3.1. Sisteme rău condiţionate Fie sistemul liniar
care admite soluţia unică x = (1,1)T.
Fie soluţia aproximativă x~ = (2,5 , 0)T. Vectorul reziduu este
Norma vectorială ||r||∞ > 0,00025. Cu toata acestea soluţia aproximativă x~ = (2,5 , 0)T este departe de cea exactă, ordinul de mărime al erorii fiind egal ea cal al soluţiei; practic
||e||∞ = ||x - x~ ||∞ = 1,5 Să modificăm sistemul (2.79.) după cum urmează:
Soluţia unică a sistemului astfel modificat este x = (0,5 , 4/3)T. Am remarcat două fenomene deosebite la sistemul analizat în Exemplul 2.3.1. Unei soluţii aproximative a sistemului depărtată de cea reală îi corespunde totuşi un reziduu mic şi o modificare aparent minoră a sistemului conduce la o schimbare majoră a soluţiei. O analiză grafică poate da explicaţia acestor dificultăţi. Soluţia sistemului (2.79.) reprezintă intersecţia dreptelor
39
2.SISTEME DE ECUAŢII LINIARE
(d1): 2x1 + 3x2 = 5 şi (d2): 2.0001x1 + 3x2 = 5.0001 Punctul (2,5 , 0 ) se afla pe dreapta d1, care este foarte apropiată de dreapta d2. Aceasta înseamnă că punctul (2,5 , 0) este foarte apropiat de dreapta (d2) chiar dacă diferă considerabil de punctul de intersecţie (1,1). Prin modificarea celei de a doua ecuaţii obţinem dreapta (d2'): 2,0002x1 + 3x2 = 5,0001 Dreapta (d'2) a rezultat printr-o modificare a pantei dreptei (d2). Cele două drepte fiind foarte apropiate, o deplasare oricât de mică a uneia dintre drepte conduce la o modificare considerabilă a punctului lor de intersecţie (a soluţiei sistemului). Sistemul descris în exemplu este un sistem rău condiţionat. Efecte similare pot apare şi în cazul sistemelor de mai multe ecuaţii rău condiţionate care sunt afectate de erori în generarea coeficienţilor. Grafic putem analiza cel mult un sistem de trei ecuaţii cu trei necunoscute. Să încercăm mai întâi să sintetizăm care sunt câteva dintre simptomele matricii A corespunzătoare unui sistem rău condiţionat [ 8 ]:
a. o schimbare minoră a elementelor matricii A generează modificări semnificative ale soluţiei sistemului ;
b. elementele de pe diagonala principală a matricii A sunt mai mici (ca ordin de mărime) decât restul elementelor;
c. produsul determinanţilor det(A) det(A-1) este semnificativ diferit de 1; d. inversa inversei matricii A, (A-1)-1 este diferită considerabil de A;
e. produsul metricilor A⋅A-1 este semnificativ diferit de matricea unitate I de acelaşi ordin eu A.
Punctele c) d) şi e) evidenţiază influenţa considerabilă a erorilor prin rotunjiră la calculul inversei unei matrici rău condiţionate. Cum putea aprecia condiţionarea unui sistem ? Dacă vom introduce la acest punct noţiunea de număr de condiţionare al unei metrici pătrate. Definiţie: Numărul de condiţionare C(A), al unei matrici pătrate nesingulare A, corespunzător unei norme naturale || ||, este
C(A) = ||A||⋅||A-1|| (2.80.) într-o primă etapă teoretică, importanţa numărului de condiţionare poate fi dedusă din următoarea teoremă: Teoremă: Dacă x este o aproximare a soluţiei sistemului Ax= b şi A este o matrice nesingulară, atunci, pentru orice normă naturală || || sunt satisfăcute inegalităţile:
Pentru x ≠ 0 şi b ≠ 0 (||x - x~ || = ||e|| şi ||r|| = ||b - Ax||). Înlocuind numărul de condiţionare dat de relaţia (2.80.) în (2.81.) şi (2.82.), obţinem:
||x - x~ || ≤ C(A) ||r||||b||
şi
40
3.1. Introducere
3. 3. INTERPOLAREA NUMERICĂ
3.1 3.1. Introducere
Interpolarea este operaţia prin care unei funcţii "fr(x)" determinată printr-un număr finit de puncte (xi, fr(xi)), i = 1,...,n, i se estimează valoarea corespunzătoare unui argument x ≠ xi, unde x ∈ (inf xi, sup xi), i = 1, . . . , n. Prin inf xi şi sup xi am notat elementul inferior, respectiv superior din mulţimea X = {xi / i = 1, . . .,n}. În viaţa de zi cu zi se întâlnesc nenumărate exemple de interpolare. Atunci când minutarul unui ceasornic se află între două diviziuni ale cadranului, valoarea diviziunii fiind de 5 minute aprecierea numărului de minute corespunzătoare orei respective se face prin interpolare. Evident, eroarea comisă nu peste depăşi numărul de minute aferent intervalului dintre două diviziuni consecutive ale cadranului. Aprecierea temperaturii dintr-o încăpere cu ajutorul unui termometru, citirea presiunii unui fluid cu ajutorul unui manometru sau al unui piezometru, măsurarea tensiunii şi a intensităţii unui curent electric au la bază operaţia de interpolare. Din punct de vedere matematic problema interpolării monovariabile se pune corect astfel încât să avem:
a) yi = f(xi), (∀) i = 1,...,n;
b) fr(α) - f(α)fr(α) < ε, (∀) α ∈ (inf xi, sup xi), α ≠ xi, i = 1 … n
Prin fr(α) am notat valoarea reală (exactă) a funcţiei cercetate în x = α, iar prin ε un număr pozitiv impus, numit eroare relativă maximă admisibilă. Funcţia y = f(x) care satisface numai condiţia a) se numeşte funcţie de interpolare sau - pe scurt - interpolant. Interpolantul care satisface şi condiţia b) se va numi interpolant admisibil. Fie d(xi) = f'r(xi) valoarea derivatei funcţiei exacte (supusă operaţiei de interpolare), în argumentul xi, (i = 1,...,n). Interpolantul y = f(x) cu proprietatea f'(xi) = d(xi), (∀) i = 1,...,n, ,se numeşte interpolant Hermite sau osculator. Un asemenea interpolant se va reprezenta grafic printr-o curbă tangentă graficului funcţiei exacte, y = fr(x), în punctele de abscise xi.
41
3. INTERPOLAREA NUMERICĂ
De regulă, un interpolant y = y(x) se construieşte de forma y = ∑j=1
nαj⋅fj(x) unde fj(x), (j = 1,...,n) se
numesc funcţii de bază iar coeficienţii αj (j = 2,...,n), urmează a se determină din condiţiile:
∑j=1
nαj⋅fj(xi) = yi, i = 1,…,n
Evident, aceste condiţii reprezintă un sistem liniar de ecuaţii în necunoscutele αj, (j = 1,...,n). Se cuvine să facem menţiunea că pot exista diferite funcţii de bază încât sistemul (3.2.) să fie compatibil. În aceste condiţii problema care se impune este de a alege cel mai potrivit ansamblu de funcţii de bază fj(x), j = 1,...,n, astfel încât să se obţină un interpolant optim sau, cel puţin, admisibil. (Prin interpolant optim se înţelege acel interpolant caracterizat prin cele mai mici erori relative). Determinarea acestui interpolant reprezintă obiectul acestui capitol. În ce priveşte funcţiile de bază utilizate, acestea vor fi de tip polinomial definite fie pe intervalul [inf xi, sup xi]i =
1…n, fie pe subintervale incluse în intervalul precedent.
3.2 3.2. Formula de interpolare Lagrange
Fie f(x) o funcţie definită pe [a,b] cunoscută în n+1 puncte arbitrare, distincte şi în ordine crescătoare x0, x1, . . . ,xn ∈ [a,b]
x x0 x1 … xn-1 xnf(x) f0 f1 … fn-1 fn
Pentru a înţelege ideea de bază a formulei lui Lagrange să considerăm produsul: L̄n,0(x) = (x - x1)(x - x2)...(x - xn)
Funcţia L̄n,0(x) este un polinom de grad n, egal cu zero pentru x = x1, x = x2, …, x = xn. Dacă împărţim L̄n,0(x) la L̄n,0(x0), funcţia rezultantă
Funcţia Ln,i(x) este un polinom de grad n care devine egal cu unitatea pentru x = xj, j ≠ i. Dacă înmulţim Ln,0 cu f0, Ln1(x) cu f1, … Lnn(x) cu fn şi însumăm expresiile astfel obţinute, suma va fi un polinom de grad n, egal cu fi, pentru i = 0,1,...,n. Formula de interpolare Lagrange dedusă în acest mod poate fi scrisă sub forma:
Ecuaţia (3.4.) este echivalentă cu o sumă de puteri ai cărei coeficienţi rezultă din rezolvarea unui sistem de ecuaţii liniare de forma (3.2.). Exemplul 3.1 - Interpolarea Lagrange în trei puncte. Să se afle interpolantul Lagrange care trece prin punctele (x0,y0) = (-1,2), (x1,y1) = (1,1) şi (x2,y2) =
42
3.3. Formule de interpolare Newton prin noduri echidistante
(2,1).
Vom determina Figura 3.1. ilustrează primele cinci funcţii Lagrange. Punctele xi au fost alese ca fiind cinci puncte egal depărtate pe intervalul [0,1]
3.3 3.3. Formule de interpolare Newton prin noduri echidistante
În capitolul anterior am studiat formula de interpolare Lagrange care poate fi utilizată atât pentru interpolarea în noduri echidistante cât şi pentru reţele la care distanţa dintre noduri este variabilă. Interpolarea Lagrange prezintă o serie de dezavantaje:
• un număr mare de calcule ;
• modificarea valorii x pentru care se face interpolarea implică de fiecare dată acelaşi număr de calcule deoarece nu se pot utiliza rezultate parţiale de la interpolările anterioare;
• modificarea numărului de noduri generează reluarea calculelor din acelaşi motiv ca la paragraful anterior ;
• evaluarea dificilă a erorii. Utilizarea interpolării Newton evită în bună măsură aceste dificultăţi. Pentru a opera o interpolare de tip Newton pentru un set de date cunoscut este necesară dezvoltarea unei tabele a diferenţelor. Odată obţinută tabela diferenţelor se pot obţine foarte uşor formule de interpolare pentru diferite seturi de noduri consecutive (de exemplu pentru i = 0,1,2,3 sau i = 3,4,5,6 sau i = 2,3,4 etc.). Vom descrie în continuare două versiuni pentru interpolarea Newton prin noduri echidistante: interpolarea Newton ascendentă şi interpolarea Newton descendentă. Matematic, cele două interpolări sunt echivalente chiar dacă au expresii diferite.
3.3.1 3.3.1. Formula de interpolare Newton ascendentă Să considerăm nodurile echidistante cu distanţa dintre două noduri succesive egală cu h. Vom nota setul de date cu (xi,fi), i = 0,1,...,n. Pentru a evalua formula ascendentă de interpolare Newton să definim mai întâi tabela diferenţelor ascendente şi coeficienţii binomiali. Diferenţele ascendente sunt definite de următoarele relaţii:
Δ0fi = fi (diferenţa ascendentă de ordinul zero)
Δ1fi = fi+1 - fi (diferenţa ascendentă de ordinul întâi)
Δkfi = Δk-1fi+1 - Δk-1fi (diferenţa ascendentă de ordinul k)
43
3. INTERPOLAREA NUMERICĂ
Tabelul 1 ilustrează o tabelă a diferenţelor până la ordinul cinci pentru un număr de şase noduri. Începând de la coloana a treia, fiecare coloană conţine diferenţele calculate pornind de la coloana anterioară. Fiecare linie include un set de diferenţe ascendente pentru punctele corespunzătoare în care s-au calculat.
Se observă uşor că în cazul particular în care fi = f(xi), dacă f(x) este un polinom de gradul k, diferenţele de ordinul k sunt egale şi diferenţele de ordin superior lui k sunt nule. Într-adevăr, diferenţele de ordinul întâi sunt polinoame de gradul k-1 ş.a.m.d. până la diferenţele de ordinul k, ce se reduc la un termen constant. Coeficienţii binomiali sunt daţi de relaţiile:
⎝⎛
⎠⎞s
0 = 1; ⎝⎛
⎠⎞s
1 = s; ⎝⎛
⎠⎞s
2 = 12! s(s - 1); ⎝⎛
⎠⎞s
3 = 13! s(s - 1)(s - 2); ⎝
⎛⎠⎞s
n = 1n! s(s - 1)(s - 2)…(s - n + 1);
unde s = (x - x0)/h iar h = xi+1 - xi. Polinomul de interpolare Newton ascendent care trece prin k+1 puncte f0,f1,..., fk poate fi scris sub forma:
Se poate arăta că ecuaţia (3.9.) este egală succesiv cu f0,f1,f2, ... ,fk pentru x = x0, x = x1, …, x = xk după cum urmează: s = 0: N(x0) = N(x0 + 0) = f0
s = k: N(xk) = f0 + kΔf0 + k(k - 1)/2⋅Δ2f0 …= fk Primii m+1 termeni din ecuaţia (3.9.) formează un polinom de interpolare de gradul m în m+1 puncte (x0, x1, …, xm). Similar, primii m+2 termeni formează un polinom de interpolare de gradul m+1 în m+2 puncte. Rezultă că gradul polinomului de interpolare se poate schimba uşor modificând numărul de diferenţe pe care le luăm în considerare de pe prima linie din tabelul 1.
44
3.3. Formule de interpolare Newton prin noduri echidistante
Dacă în ecuaţia (3.9.) înlocuim x0 şi f0 respectiv cu x2 şi f2 formula devine:
N(x2 + sh) = ∑n=0
k
⎝⎛
⎠⎞s
n Δnf2 (3.10.)
unde s = (x - x2)/h. Valoarea lui s devine egală cu 0 pentru x = x2 şi succesiv egală cu 1,2,3,... pentru x = x3, x = x4, x = x5... Ecuaţia (3.10.) este un polinom de gradul k ce va verifica punctele x2, x3, ...,xk+2 şi care utilizează diferenţele aflate de-a lungul celei de a treia linii din tabelul 1. Aceasta demonstrează că, odată întocmită tabela diferenţelor, se pot obţine uşor formule de interpolare pentru diferite seturi de puncte succesive.
3.3.2 3.3.2. Formula de interpolare Newton descendentă Polinomul de interpolare Newton descendent este funcţie de diferenţele descendente şi coeficienţii binomiali. Diferenţele descendente sunt definite de următoarele relaţii:
∇0fi = fi (diferenţa descendentă de ordinul zero)
∇1fi = fi - fi-1 (diferenţa descendentă de ordinul întâi)
∇2fi = ∇1fi - ∇1fi-1 = fi - 2fi-1 + fi-2 (diferenţa descendentă de ordinul doi)
∇kfi = ∇k-1fi-1 - ∇k-1fi (diferenţa descendentă de ordinul k) Coeficienţii binomiali utilizaţi în formula de interpolare Newton descendentă sunt:
⎝⎛
⎠⎞s-1
0 = 1; ⎝⎛
⎠⎞s
1 = s; ⎝⎛
⎠⎞s+1
2 = 12! s(s + 1); ⎝
⎛⎠⎞s+2
3 = 13! s(s + 1)(s + 2);
⎝⎛
⎠⎞s+n-1
n = 1n! s(s + 1)(s + n - 2)…(s + n - 1);
Polinomul de interpolare Newton descendent care trece prin punctele x = xj,x = xj-1, x = xj-2, … , x = xj-k este dat de relaţia:
N(x) = N(xj + sh) = ∑p=0
k
⎝⎛
⎠⎞s+p-1
n ∇pfj , -k ≤ s ≤ 0 (3.11.)
unde s = (x - xj)/h, ⎝⎛
⎠⎞s+p-1
p este coeficientul binomial şi ∇nfj este diferenţa descendentă. Echivalenţa dintre diferenţa ascendentă şi diferenţa descendentă este dată de relaţia:
Δpfj = ∇pfj-p (3.12.) Ţinând cont de (3.12.), ecuaţia (3.11.) poate fi scrisă în funcţie de diferenţele ascendente:
N(x) = ∑p=0
k
⎝⎛
⎠⎞s+p-1
n Δpfj-p , -k ≤ s ≤ 0
45
3. INTERPOLAREA NUMERICĂ
sau mai explicit
3.4 3.4. Formula de interpolare Newton prin noduri neuniform distribuite
Formulele de interpolare Newton descrise anterior sunt valabile numai în cazul nodurilor echidistante. În practică, de cele mai multe ori, nu întâlnim acest caz fericit. Vom determina în continuare o metodă de interpolare Newton pentru cazul nodurilor neuniform distribuite. Să presupunem că Pn este polinomul Lagrange de grad cel mult n care coincide cu o funcţie f în punctele x0,x1,...,xn,. Metoda descrisă în continuare, care se bazează pe noţiunea de diferenţe divizate, poate fi dedusă demonstrând că polinomul P, este de forma:
Pn(x) = a0 + a1(x - x0) + a2(x - x0) (x - x1) + … + an(x - x0)(x - x1) … (x - xn-1) (3.14.) pentru anumite valori ale constantelor a0,a1,...,an. Pentru a determina prima constantă, a0, vom ţine cont de faptul că Pn poate fi scris sub forma dată de ecuaţia (3.14.); evaluând Pn în punctul x0, obţinem a0 = Pn(x0) = f(x0). Similar, atunci când evaluăm Pn, în punctul x1, singurii termeni diferiţi de zero din expresia pentru Pn(x1) sunt constanta şi termenul liniar.
f(x0) + a1(x - x0) = Pn(x1) = f(x1)
de unde rezultă a1 = f(x1) - f(x0)x1 - x0
(3.15)
La această etapă vom introduce noţiunea de diferenţă divizată. Diferenţa divizată zero a funcţiei f(x) în punctul xi, se notează f[xi] şi este chiar valoarea funcţiei în punctul respectiv: f[xi] = f(xi). Următoarele diferenţe divizate se definesc prin inducţie; prima diferenţă divizată a funcţiei f în punctele xi şi xi+1, se notează f[xi,xi+1] şi este dată de relaţia
f[xi,xi+1] = f[xi+1] - f[xi]xi+1 - xi
(3.16)
Considerând că am determinat diferenţele divizate de ordinul (k-1), f[xi,xi+1,xi+2, … ,xi+k-1] şi f[xi+1,xi+2, … ,xi+k-1,xi+k], diferenţa divizată de ordinul k în punctele xi,xi+1,xi+2, … ,xi+k-1,xi+k este dată de relaţia :
Cu aceste notaţii coeficientul a1 dat de relaţia (3.15.) poate fi exprimat sub forma a1 = f[x0,x1] iar polinomul de interpolare Pn dat de (3.14.) devine:
Pn(x) = f[x0] + f[x0,x1](x - x0) + a2(x - x0)(x - x1) + … + an(x - x0)(x - x1) … (x - xn-1) Repetând raţionamentul anterior şi folosindu-ne de notaţiile introduse pentru diferenţele divizate, polinomul Pn(x), va fi exprimat de relaţia: Pn(x) = f[x0] + f[x0,x1](x - x0) + f[x0,x1,x2](x - x0)(x - x1) + … + f[x0,x1,…,xn](x - x0)(x - x1) … (x - xn-1) echivalentă cu
Ecuaţia (3.17.) este cunoscută sub denumirea de formula de interpolare Newton cu diferenţe divizate.
46
Funcţii spline cubice
O schemă de calcul a diferenţelor divizate este ilustrată in tabelul 2. Tabelul 2
z f(x) Prima diferenţă divizată A doua diferenţă divizată x0 f[x0] f[x0,x1] = f[x1] - f[x0]
x1 - x0 f[x0,x1,x2] = f[x1 x2] - f[x0 x1]
x2 - x0
x1 f[x1] f[x1,x2] = f[x2] - f[x1]x2 - x1
f[x1,x2,x3] = f[x2 x3] - f[x1 x2]x3 - x1
x2 f[x2] f[x2,x3] = f[x3] - f[x2]x3 - x2
f[x2,x3,x4] = f[x3 x4] - f[x2 x3]x4 - x2
x3 f[x3] f[x3,x4] = f[x4] - f[x3]x4 - x3
x4 f[x4]
3.4.1 Algoritmul de calcul pentru tabela diferenţelor divizate Permite calculul coeficienţilor polinomului de interpolare P în n+1 puncte. DATE DE INTRARE :
Pasul 2. SCRIE (Q0,0,Q1,1, Qn,n) STOP (Qi,i este f[x0,x1, … ,xi] )
3.5 Funcţii spline cubice
Vom descrie în acest capitol o metodă de interpolare care utilizează polinoame definite pe subintervale şi care, cu excepţia punctelor extreme ale domeniului de interpolare, nu necesită informaţii referitoare la derivatele funcţiei. Cea mai utilizată aproximare polinomială pe subintervale este interpolarea spline cubică care poate fi definită după cum urmează: Definiţie: Considerând o funcţie f definită pe intervalul [a,b] şi un set de noduri a = x0 < x1 < ….< xn = b, un interpolant spline cubic S, pentru funcţia f este o funcţie care satisface următoarele condiţii :
a. S este un polinom cubic, notat Si pe intervalul [xi,xi+1] pentru i = 0,1,...,n-1 ; b. S(xi) = f(xi) pentru i = 0,1,....n; c. Si+1(xi+1) = Si(xi+1) pentru i = 0,1,... ,n-1; d. S'i+1(xi+1) = S'i(xi+1) pentru i = 0,1,...,n-2 ; e. S"i+1(xi+1) = S"i(xi+1) pentru i = 0,1,...,n-2. f. este satisfăcută una din următoarele condiţii :
a. f1) S"(x0) = S"(xn) = 0 b. f2) S'(x0) = f(x0) şi S'(x.) = f'(xn)
Observăm ca un interpolant spline admite şi derivata de ordinul doi continuă, condiţie care constituie un avantaj în majoritatea aplicaţiilor. Condiţia din definiţie ca S"(x) să fie continuă in
47
3. INTERPOLAREA NUMERICĂ
cele n-2 puncte interioare nu este suficientă, fiind necesare informaţii suplimentare referitor şi la derivatele in punctele extreme x0 şi xn ale domeniului de interpolare. Vom deduce mai întâi expresia generală pentru un interpolant spline cubic, după care vom analiza mai multe posibilităţi de determinare ale unui spline unic în funcţie de condiţiile îndeplinite in punctele extreme x0 şi xn. Să considerăm intervalul [xi, xi+1] de lungime hi = xi+1 - xi (vezi figura 3.3.). Pe acest interval putem defini polinomul cubic
Si(x) = ai + bi (x - xi) + ci (x - xi)2 + di (x - xi)3 (3.53.) unde xi < x < xi+1
Fig. 3.5. Notaţii in interpolarea spline cubică Impunând mai întâi condiţia ca S(x) să fie egal cu valorile funcţiei f(x) în punctele x = xi şi x = xi+1 obţinem:
fi = ai (3.54.) fi+1 = ai + bihi + cihi
2 + dihi3 (3.55.)
în plus este necesar ca S' şi S" să fie continue în punctele x = xi şi x = xi+1 cu polinoamele cubice definite pe intervalele adiacente [xi-1, xi] şi [xi+1, xi+2]. Derivata de ordinul doi a polinomului cubic (3.53.),
S"(x) = 2ci + 6di (x - xi) (3.56.) este definită în punctele i şi i+1 de relaţiile:
S"i = 2ci (3.57.) S"i+1 = 2ci + 6di hi (3.58.)
Rezultă coeficienţii ci şi di, funcţie de S"i şi S"i+1 :
ci = S"i2 (3.59.); di = S"i+1 - S"i
6hi (3.60.)
Coeficientul ai este deja dat de expresia (3.54.). Înlocuind relaţiile (3.54.), (3.59.)/(3.60.) în (3.55.) vom elimina ai, ci şi di, şi vom obţine coeficientul bi:
Derivata de ordinul întâi a expresiei (3.62.) este dată in punctele x = xi şi x = xi+1 de relaţiile :
unde hi = xi+1 - xi. Pentru intervalul xi-1 < x < xi, (3.64.) devine
48
Funcţii spline cubice
Din ecuaţia de continuitate pentru derivata de ordinul întâi în nodul i, rezultă că expresiile date de relaţiile (3.63)şi (3.65.) pentru S'i trebuie să fie egale. Egalând aceste expresii rezultă relaţia
înlocuind ai şi ci, daţi de (3.54.) şi (3.59.), ecuaţia (3.66.) poată fi rescrisă sub forma
Aplicând relaţia (3.67.) pentru i = 1,2,...,n-1 rezultă un sistem liniar de n-1 ecuaţii. Necunoscute vor fi variabilele ci (i = 0,1,...,n), deoarece variabilele hi (i = 0,1,...,n-1) şi ai (i = 0,1,...,n) sunt cunoscute ca distanţele dintre noduri şi valorile funcţiei f date. Sistemul astfel obţinut având n-1 ecuaţii şi n+1 necunoscute, să analizăm posibilităţile de a elimina nedeterminarea. O primă posibilitate de a ridica nedeterminarea constă în impunerea condiţiei S"(a) = S"(b) = 0. Se obţine astfel un spline cubic natural. Fizic această condiţie consideră că graficul interpolantului este o dreaptă in afara intervalului [x0,xn]. Utilizând notaţiile existente, această condiţie conduce la cn = S"(x.) = 0 şi 0 = S"(x0) = 2c0 + 6d0(x0 - x0) de unde c0 = 0. Cele două ecuaţii c0 = 0 şi cn = 0 împreună cu ecuaţiile (3.67.) generează un sistem liniar descris de ecuaţia matricial vectorială Ax = b, în care A este matricea de dimensiuni (n+1)x(n+1):
iar b şi s sunt vectorii :
Matricea A este strict dominantă diagonal ceea ce înseamnă că sistemul admite o soluţie unică c0,c1,... ,cn. Soluţia unei probleme de interpolare cubică spline în care utilizăm condiţiile la limită S"(x0) = S"(xn) = 0 poate fi obţinută prin aplicarea următorului algoritm:
3.5.1 Algoritm spline cubic natural Construieşte un interpolant spline cubic S pentru funcţia f, definit în punctele x0 < x1 <...< xn şi îndeplinind condiţiile S"(x0) = S"(xn) = 0 DATE DE INTRARE:
49
3. INTERPOLAREA NUMERICĂ
• n, x0,x1,…,xn. Se generează ai = f(xi) pentru i = 0,1,...,n sau se introduc valorile ai = 0,1,…,n date tabelar
DATE DE IEŞIRE: • aj, bj, cj, dj pentru j = 0,1,...,n-1
(Obs : S(x) = Sj(x) = aj + bj(x - xj) + cj(x - xj)2 + dj(x - xj)3 pentru xj < x < xj+1) Pasul 1. PENTRU i = 0,1,...,n-1 CALCULEAZĂ hi = xi+1 - xi Pasul 2. PENTRU i = 1,2,...,n-1 CALCULEAZĂ (elementele vectorului termenilor liberi b)
αi = 3[ai+1hi-1 - ai(xi+1 - xi-1) + ai-1hi]hi-1hi
Pasul 3. (Paşii 3,4,5 şi parţial Pasul 6 rezolvă un sistem liniar tridiagonal)
l0 = 1; μ0 = 0; z0 = 0. Pasul 4. PENTRU i = 1,2,...,n-1 CALCULEAZĂ
li = 2(xi+1 - xi-1) - hi-1μi-1; μi = hi/li; zi = (αi - hi-1zi-1)/li Pasul 5. IMPUNE
cj = zj - μjcj+1; bj = (aj+1 - aj)/hj - hj(cj+1 + 2cj)/3 ; dj = (cj+1 - cj)/(3hj) Pasul 7. PENTRU j = 0,1,...,n-1SCRIE (aj, bj, cj, dj) STOP A doua condiţie care permite calcularea unui spline cubic este atunci când se cunosc valorile primei derivate ale funcţiei care urmează să fie aproximată Ia punctele x0 şi xn. În acest caz se impun condiţiile S"(x0) = f'(x0) şi S'(xn) = f'(xn). Să analizăm modul de rezolvare al problemei în această situaţie. Ţinând cont de (3.54.) ,şi (3.59.), relaţia (3.61.) care exprimă coeficientul bi poate fi rescrisă sub forma:
bi = 1hi (ai+1 - ai) -
hi3 (2ci + ci+1) (3.68.)
Utilizând egalitatea evidentă S'(a) = S'(x0) = b0, ecuaţia (3.61.) devine pentru i = 0
f'(a) = a1 - a0h0
- h03 (2c0 + c1) (3.68.)
În consecinţă,
Similar,
f"(b) = bn = bn-1 + hn-1(cn-1 + cn) de unde ecuaţia (3.61.) implică pentru i = n-1 :
50
Funcţii spline cubice
şi
Ecuaţiile (3.67.) împreună cu (3.69.) şi (3.70.) determină sistemul liniar Ax = b, unde A este următoarea matrice de dimensiune (n+1)x(n+1) :
iar b şi x sunt vectorii :
Matricea A este de asemenea strict dominantă diagonal, ceea ce înseamnă că sistemul admite o soluţie unică.
3.5.2 Algoritm spline cubic condiţionat la extremităţi Construieşte un interpolant spline cubic S pentru funcţia f, definit în punctele x0 < x1 <...< xn şi îndeplinind condiţiile S'(x0) = f'(x0) şi S'(xn) = f'(xn) DATE DE INTRARE:
• n, x0,x1,…,xn. Se generează ai = f(xi) pentru i = 0,1,...,n sau se introduc valorile ai = 0,1,…,n date tabelar.
• D0 = f'(x0), DN = f'(xn) DATE DE IEŞIRE:
• aj, bj, cj, dj pentru j = 0,1,...,n-1 (Obs : S(x) = Sj(x) = aj + bj(x - xj) + cj(x - xj)2 + dj(x - xj)3 pentru xj < x < xj+1) Pasul 1. PENTRU i = 0,1,...,n-1 CALCULEAZĂ hi = xi+1 - xi Pasul 2. CALCULEAZĂ
α0 = 3(a1 - a0)/h0 - 3D0; αn = 3DN - 3(an.- an-1)/hn-1 Pasul 3. PENTRU i = 1,2,...,n-1 CALCULEAZĂ
αi = 3[ai+1hi-1 - ai(xi+1 - xi-1) + ai-1hi]hi-1hi
Pasul 4. (Paşii 4,5,6 şi parţial Pasul 7 rezolvă un sistem liniar tridiagonal) l0 = 2h0; μ0 = 0,5; z0 = α0/l0 Pasul 5. PENTRU i = 1,2,...,n-1 CALCULEAZĂ
cj = zj - μjcj+1; bj = (aj+1 - aj)/hj - hj(cj+1 + 2cj)/3 ; dj = (cj+1 - cj)/(3hj) Pasul 8. PENTRU j = 0,1,...,n-1SCRIE (aj, bj, cj, dj) STOP Există şi alte posibilităţi care permit ridicarea nedeterminării pentru sistemul (3.67.). Amintim dintre acestea
• Determinarea derivatelor S" la extremităţi prin extrapolare dinspre interiorul domeniului de interpolare. S"0 poate fi extrapolat cu ajutorul relaţiei :
(3.71.) Similar poate fi exprimat şi S"n
(3.72.) Eliminând S"0 şi S"n din sistemul (3.67.) cu ajutorul ecuaţiilor (3.71.) şi (3.72.) va rezulta un sistem tridiagonal de (n-1) ecuaţii cu (n-1) necunoscute, necunoscute fiind S"i, i = 1,2,...,n-1. Problema mai permite rezolvarea şi în cazul în care funcţia cubică spline este o funcţie periodică cu perioada egală cu lungimea domeniului de interpolare [x0,xn] (de exemplu, o curbă închisă). Această condiţie implică
f0 = S(xn+1), S'(x0) = S'(xn+1), S"(x0) = S"(xn+1) ceea ce va permite generarea unui spline complet. În literatura de specialitate există descrise şi alte posibilităţi care fac posibilă generarea unui spline cubic complet.
3.6 3.11. Interpolarea bidimensională
Schemele de interpolare bidimensională sunt de două feluri. Primul tip utilizează de două ori interpolarea monodimensională şi poartă denumirea de interpolare dublă. Al doilea tip se bazează pe polinoame de interpolare bidimensionale definite pe subintervale şi este întâlnită in metoda elementului finit. Pentru a explica interpolarea dublă să presupunem că se dau valorile unei funcţii f(x,y) pe o reţea (xi,yj). Notăm cu fi,j = f(xi,yj) valoarea funcţiei intr-un punct (xi,yj). Interpolarea dublă implică doi paşi, în fiecare dintre paşi aplicându-se o metodă oarecare de interpolare monodimensională. Să presupunem că vrem să estimăm valorile funcţiei în punctul de coordonate (x,y), unde în xi-1 ≤ x ≤ xi şi yi-1 ≤ y ≤ yi (vezi Fig. 3.6.). Pentru simplificarea explicaţiilor vom aplica în cei doi paşi interpolarea liniară. La primul pas vom aplica interpolarea liniară pe direcţia y şi vom determina valorile funcţiei în punctele E şi F:
52
3.11. Interpolarea bidimensională
Fig. 3.6. Interpolare biliniară într-
un domeniu bidimensional
La al doilea pas aplicăm interpolarea liniară între fE si fF:
înlocuind (3.79.) în (3.80.) obţinem următoarea formulă:
Ordinea paşilor de interpolare poate fi schimbată fără să afecteze rezultatul final. Altfel spus, putem determina mai întâi valorile funcţiei prin interpolare pe direcţia x şi apoi aplicând interpolarea liniară pe direcţia y, obţinem g(x,y). Această schimbare nu afectează rezultatul. Interpolarea liniară care am utilizat-o în acest subcapitol poate fi înlocuită cu oricare dintre metodele de interpolare monovariabilă analizate anterior.
53
CUADRATURA NUMERICĂ pg 166
4. CUADRATURA NUMERICĂ pg 166
4.1 Introducere
În numeroase aplicaţii practice se întâlnesc integrale de tipul ⌡⌠a
bf(x) dx (4.1.) în situaţia când nu
se poate obţine primitiva funcţiei f sub forma unei sume cu număr finit de funcţii elementare. Evaluarea integralelor definite (4.1.) când nu se cunoaşte primitiva funcţiei f se face numai prin metode numerice şi se numeşte cuadratură. Cuadratura are deseori la bază operaţia de interpolare a funcţiei de sub integrala (4.1.). Aceasta înseamnă să înlocuim funcţia f printr-o funcţie de interpolare, definită pe domeniul [a,b], căreia i se poate determina primitiva. De regulă, funcţia de interpolare utilizată va fi polinomială, permiţând o evaluare simplă a integralei (4.1.). Uneori se defineşte pe întregul domeniu [a,b] un singur polinom de interpolare, iar alteori se definesc mai multe polinoame, fiecare aplicându-se pe câte un subdomeniu inclus în domeniul dat [a,b]. Se cuvine să menţionăm faptul că pe lângă evaluarea integralelor (4.1.) (utilizând o funcţie de interpolare oarecare pe domeniul [a,b]), trebuie să se precizeze şi ordinul fie mărime al erorii comise prin cuadratură. Consideraţii similare se pot face şi pentru integralele duble sau triple care implică evident, metode de studiu mai complexe. Revenind asupra integralei (4.1.), aceasta poate fi evaluată prin utilizarea funcţiei de interpolare f, obţinându-se o relaţie de cuadratură de forma :
I = ⌡⌠a
bf(x) dx = ∑
i=0
nwi⋅f(xi) dx + Rn
unde wi şi xi, (i = 1,n) se numesc ponderi, respectiv noduri, iar Rn - eroare sau rest. În paragrafele care urmează vom prezenta câteva reguli simple de cuadratură (regula dreptunghiurilor, a trapezelor şi a lui Simpson) precum şi regulile lui Gauss şi regula decuadratură Gauss - Konrod (considerată deosebit de eficientă).
4.2 Regula dreptunghiului şi regula trapezului
Aceste două reguli sunt dintre cele mai cunoscute în cuadratura numerică. Nu sunt utile aplicate în forma elementară dar există metode practice care la iau în considerare. Deoarece admit şi o interpretare geometrică vom face referire şi la graficul din Figura 4.1.
54
Regula dreptunghiului şi regula trapezului
Pentru a deduce formula care exprimă regula dreptunghiului vom dezvolta funcţia f(x) în serie Taylor, în punctul x̄ = (a+b)/2.
...+24
)m-)(x x(f+6
)m-)(xx("f+
2
)m-)(x x(f"+m)-)(x x(f+) xf( = f(x)
4IV3
2
′
′
(4.3)
Integrând expresia astfel obţinută între limitele a şi b obţinem următoarea formulă care exprimă regula dreptunghiului:
...+24
)a-(b) x(f"+a)-)(b xf( = x f(x)d = I3b
a∫ (4.4)
unde termenul care conţine f"(x) reprezintă eroarea dată de această formulă.
Am arătat în introducere că principiul de bază în rezolvarea integralei ⌡⌠a
bf(x) dx constă în
aproximarea acesteia printr-o sumă de forma ∑i=0
nwi⋅f(xi) dx .
Regula trapezului (şi regula lui Simpson descrisă în subcapitolul următor) se bazează pe polinoamele de interpolare descrise în capitolul 3. INTERPOLAREA NUMERICĂ. Pentru a determina formula trapezului de
aproximare a integralei: ⌡⌠a
bf(x) dx vom considera
nodurile x0 = a, x1 = b, distanţa dintre noduri h = b - a şi vom utiliza polinomul de interpolare Lagrange de gradul unu:
)xf()x-x()x-(x
101
0+)xf()x-x(
)x-(x = (x)L 010
11
Considerând şi eroarea interpolării Lagrange Figura 4.1 Regula dreptunghiului şl regula trapezului
⌡⌠a
bf(x) dx =
⌡⎮⌠
x0
x1
⎣⎢⎡
⎦⎥⎤(x0 - x1)
(x - x1) f(x0) + (x - x0)(x1 - x0) f(x1) dx +
12 ⌡⌠
x0
x1
f"(ξ(x)) (x - x0)(x - x1) ( )
unde ξ(x) se află în [a,b] pentru orice x. Deoarece produsul (x - x0)(x - x1) nu îşi schimbă semnul pe intervalul [x0,x1], se poate aplica teorema lui Rolle termenului care exprimă eroarea:
55
CUADRATURA NUMERICĂ pg 166
În consecinţă, ecuaţia (4.5.) devine:
sau
(4.6.) relaţie care exprimă regula trapezului. Motivul pentru care această formulă poartă denumirea de regula trapezului este că atunci când f este o funcţie care are numai valori pozitive, integrala
⌡⌠a
bf(x) dx este aproximată de aria trapezului reprezentat în graficul din Figura 4.1.
Să remarcăm faptul că regula trapezului generează un rezultat exact atunci când este aplicată oricărei funcţii a cărei derivată de ordinul doi este identic egală cu zero, altfel spus, oricărui polinom de grad cel mult unu. Există şi o a doua demonstraţie pentru regula trapezului. Dacă în dezvoltarea în serie Taylor a funcţiei f(x) înlocuim succesiv x = a şi x = b şi adunăm cele două serii care rezultă obţinem:
f(x̄) = [f(a) + f(b)]/2 - f"(x̄)(b - a)2
8 - fiv(x̄)(b - a)4
384
(Termenii care conţin derivatele de ordin impar s-au redus). Substituind relaţia astfel generată în formula dreptunghiului şi rearanjând termenii rezultă formula trapezului:
I = ⌡⌠a
bf(x) dx =
b - a2 [f(a) + f(b)] - f"(x̄)
(b - a)3
12
Să considerăm un interval de integrare [a,b], larg, pe care funcţia f(x) prezintă variaţii considerabile. Dacă ar fi să aplicăm regula trapezului în forma elementară dată de (4.6.),
evident că vom obţine o aproximaţie grosieră pentru ⌡⌠a
bf(x) dx
În consecinţă vom împărţi intervalul [a,b] în n subintervale egale (h = (b-a)/n) şi vom aplica regula trapezului pe fiecare subinterval [a + ih, a + (i+1)h], i = 0,1,...,n-1 (vezi Figura 4.2). Adunând relaţiile astfel obţinute rezultă regula extinsă a trapezului:
I = ⌡⌠a
bf(x) dx =
h2⎣⎢⎢⎡
⎦⎥⎥⎤
f(a) + 2 ⌡⌠j=1
n-1f(a+jh) + f(b) + E
care poate fi scrisă şi sub următoarea formă echivalentă:
56
Regula dreptunghiului şi regula trapezului
I = h2 (f0 + 2f1 + 2f2 + … + 2fn-1 + fn) + E
unde f0 = f(a), f1 = f(a+h) şi fi = f(a+ih) În regula extinsă a trapezului eroarea este definită ca suma erorilor rezultate pe fiecare dintre subintervale. Considerând că s-a aplicat regula extinsă a trapezului pe un interval [a,b] care este divizat de punctele a = x0, x1,...,xn = b în n subintervale, eroarea devine
) x(f"n
)a-(b121- i
n
1=i3
3
∑ E _
( )
unde xi = (xi-1 + xi)/2 şi h = (b - a)/n. Dacă definim f" ca media derivatelor de ordinul doi calculate în punctele xi,
Figura 4.2 Regula extinsă a trapezului
n n
)x(f" = "f
i
n
1=i∑
ecuaţia (4.8.) poate fi rescrisă sub forma
E = -112 (b - a)h2f"̄ (4.9.)
Se observă din (4.9) că eroarea în regula extinsă a trapezului este proporţională cu h2. În analiza anterioară a erorii am presupus că f(x) este o funcţie analitică pe [a,b]. În caz contrar eroarea nu mai este proporţională cu h2. O aplicaţie importantă a analizei erorii la regula extinsă a trapezului o constituie integrarea Romberg. Să presupunem că rezultatul obţinut prin regula extinsă a trapezului, pentru h = (b-a)/n (n-1 subintervale), este Tn şi pentru h'=h/2 (2n-2 subintervale), T2n-1. Deoarece eroarea în regula extinsă a trapezului este proporţională cu h2, erorile obţinute în cazul intervalelor de lungime h, respectiv h/2 pot fi scrise sub forma
unde C este o constantă. În acelaşi timp, valoarea exactă a integralei este
,E+T = E+T = I 1-2n1-2nnn de unde T-T=E-E n1-2n1-2nn
Înlocuind (4.10) în (4.11) obţinem )T-T(h34 = C n1-2n
2-
Astfel, valoarea aproximativă pentru En calculată cu prima relaţie (4.10) este:
Cunoscând valorile pentru Tn şi T2n-1 se poate calcula o valoare mai precisă pentru integrală aplicând relaţia:
57
CUADRATURA NUMERICĂ pg 166
Eroarea obţinută prin aplicarea relaţie (4.12') este proporţională cu h4 ceea ce înseamnă că rezultatul are o precizie mai bună decât Tn sau T2n-1. Această tehnică poartă denumirea de integrarea Romberg.
4.2.1 Regulile Simpson (1/3 şi 3/8) Prima regulă a lui Simpson (regula 1/3) se bazează pe un polinom de interpolare de gradul doi (vezi Figura 4.3). Integrând între limitele x0 = a şi x2 = b polinomul de interpolare Newton care verifică punctele x0, x1, x2, obţinem regula lui Simpson (1/3):
I = ⌡⌠a
bf(x) dx =
h3[(f(a) + 4f(x̄) + f(b)] + E Figura 4.3 Regula 1/3 a lui Simpson
unde h = (b - a)2 şi x̄ = (a + b)
2
Ecuaţia (4.13) poate fi scrisă într-o formă echivalentă:
I = h3 [(f0 + 4f1 + f2]+ E
unde fi = f(xi) = f(a+ih), i = 0,1,2 În ceea ce priveşte eroarea am putea încerca să integrăm expresia erorii interpolării pătratice, dar rezultatul integrării ar fi egal cu zero, ceea ce evident că nu ar exprima eroarea reală. Cauza acestei consecinţe înşelătoare provine ca urmare a aproximării ξ = x. Vom aborda în continuare o altă metodă de deducere a erorii. Pentru regula lui Simpson (1/3) vom utiliza dezvoltarea în serie Taylor. Dezvoltările în serie Taylor pentru f0 şi f2 în punctul x1 (echivalent cu x = (a+b)/2) sunt:
...+fh241+"fh6
1"+fh21+hf+f = f
...-fh241+"fh6
1"-fh21+hf-f = f
IV1
41
31
2112
IV1
41
31
2110
′
′
′
′
Înlocuind cele două dezvoltări în ecuaţia (4.14) obţinem:
E+...+fh361"+fh3
1+hf2 = I IV1
51
31 ( )
Concomitent, dezvoltarea în serie Taylor a funcţiei f(x) în punctul x1 este
...+fz241+"fz6
1"+fz21+zf+f = f(x) IV
14
13
12
11 ′′ ( )
unde x = x1 + z, echivalent cu z = x - x1. Integrarea analitică a acestei dezvoltări pe intervalul [a,b] ne conduce la următorul rezultat:
...+fh601"+fh3
1+hf2 = x d f(x) IV1
51
31
b
a∫ ( )
care este integrala exactă a funcţiei dezvoltate în serie Taylor. Scăzând ecuaţia (4.15) din ecuaţia (4.17) şi efectuând trunchierea după primul termen care nu se reduce, obţinem eroarea formulei de integrare (4.15.):
58
Regula dreptunghiului şi regula trapezului
E = -190 h5f1
iv (4.9.)
unde f1iv = fiv(x1) Eroarea dată de regula lui Simpson (1/3) este egală cu zero dacă f(x) este un polinom de grad cel mult trei. Se poate obţine şi regula extinsă lui Simpson (1/3) în modul următor: se divizează intervalul [a,b] în n subintervale şi se aplică regula elementară a lui Simpson pe fiecare pereche de subintervale consecutive. Deoarece aplicarea regulii lui Simpson necesită două subintervale, n trebuie să fie un număr par, adică n = 2m, m oarecare (vezi Fig. ) Prin însumare obţinem regula 1/3 extinsă a lui Simpson:
unde h = (b-a)/2m şi a = x0 < x1 < ….< x2m = b unde xj = x0 + jh pentru j = 0,1,2,…,2m
Eroarea este f180ha)-(b- = f
90h
2N- E IV
4IV
5
_ unde N
)( xf = f
iIV
n
1=iIV∑
.
Următorul algoritm implementează regula extinsă a lui Simpson (1/3) pentru un interval de integrare divizat în n= 2m subintervale. Este un algoritm de cuadratură des întâlnit în practică.
4.2.2 Algoritmul regulii de cuadratură extinse a lui Simpson (1/3)
Aproximează integrala I = ⌡⌠a
bf(x) dx
DATE DE INTRARE: limitele de integrare a şi b ; numărul natural m DATE IEŞIRE: XI: aproximarea integralei I Pasul 1. CALCULEAZĂ h = (b-a)/(2m); Pasul 2. INIŢIALIZEAZĂ XI0 = f(a)+f(b); XI1= 0 ; (suma termenilor f(x2j-1)) ; XI2= 0; (suma termenilor f(x2j)); Pasul 3. PENTRU i = 1,2,...,2m-1 EFECTUEAZĂ Paşii 4 şi 5: Pasul 4. CALCULEAZĂ X = a+ih ; Pasul 5. DACĂ i este par CALCULEAZĂ XI2 = XI2 + f(x) ALTFEL CALCULEAZĂ XI1 = XI1+f(x)
Pasul 6. CALCULEAZĂ XI = h(XI0 + 2⋅XI2 + 4⋅XI1)/3 ; Pasul 7 SCRIE (XI) STOP O altă regulă de cuadratură de tip Simpson este regula 3/8 care poate fi dedusă prin integrarea unei formule de interpolare polinomială de gradul trei. Pentru un domeniu de integrare [a,b] divizat în trei subintervale, această regula se scrie sub forma:
59
CUADRATURA NUMERICĂ pg 166
unde h = (b - a)3 şi fi = f(a + ih), i = 0,1,2,3
Termenul care exprimă eroarea E este dat de ) x(fh83- E IV5_ unde x̄ = a + b
2 . Această
expresie a erorii se poate deduce pornind de la dezvoltarea în serie Taylor şi urmărind aceeaşi paşi descrişi la regula lui Simpson (1/3). După cum am explicat deja regula extinsă 1/3 este aplicabilă pentru un număr par de subintervale, în timp ce regula extinsă 3/8 se poate aplica în cazul în care numărul de subintervale este multiplu de trei. Atunci când numărul de subintervale este impar şi nu este multiplu de trei se poate aplica regula 3/8 pentru primele trei sau ultimele trei subintervale iar pe restul par de subintervale regula extinsă 1/3. Ordinul de mărime al erorii este acelaşi în ambele reguli.
4.3 Formule de cuadratură Newton-Cotes
Regula trapezului şi regulile lui Simpson sunt exemple de reguli de cuadratură aparţinând clasei de formule Newton-Cotes care rezultă prin integrarea formulelor de interpolare Newton. Există două tipuri de formule Newton-Cotes: formule închise care utilizează limitele a şi b ale domeniului de integrare şi formule deschise care nu întrebuinţează aceste două limite. Formulele Newton-Cotes închise utilizează nodurile xi = x0 + ih pentru i = 0,1,...,n, unde x0 = a şi xn = b, h = (b - a)/n şi pot fi scrise în forma:
⌡⌠a
bf(x) dx = αh[w0f0 + w1f1 + … + wnfn] + E (4.22)
unde α şi w sunt constante (vezi Tabelul 4.1) şi fi = f(xi) pentru i = 0,1,...,n Tabelul 4.1 Constante pentru formule Newton-Cotes închise
n α wi, i = 0, 1,…n E 1 1/2 1 1 -1/12 h3 f" 2 1/3 1 4 1 -1/90 h5 fiv
5888 989 Se observă că pentru n = 1 avem regula trapezului iar pentru n = 2 şi n = 3 regulile Simpson 1/3 respectiv 3/8. Formula de integrare (4.22) poate fi extinsă şi în afara domeniului de integrare dat.
60
Cuadratura Gauss
Formulele Newton-Cotes deschise se obţin prin extinderea domeniului de integrare cu un subinterval la stânga primului punct dat şi cu un subinterval la dreapta ultimului punct al domeniului de integrare. Aceste formule sunt de forma:
unde h = (b-a)/(n+2). Dintre constantele α şi wi (vezi Tabelul 4.2) w0 şi wn+2 sunt întotdeauna egale cu zero deoarece corespund extremităţilor domeniului extins. Din acest motiv f0 şi fn+2 nu sunt necesare.
Tabelul 4.2 Constante pentru formule Newton-Cotes deschise
n α wi, i = 0,1 ,...,n+2 E 1 3/2 0 1 1 0 1/4 h3 f" 2 4/3 0 2 -1 2 0 28/90 h5 fiv
Comparând o formulă deschisă cu o formulă închisă, ambele utilizând aceleaşi număr de puncte n, eroarea generată prin aplicarea formulei deschise este semnificative mai mare decât cea dată de formula închisă. Formulele deschise prezintă totuşi avantajul că pot fi aplicate în cazul în care nu se cunosc valorile funcţionale în punctele extreme. Se observă în Tabelul 4.1 şi Tabelul 4.2 că pentru n mare valorile pentru wi devin mari şi prezintă schimbări de semn. Din acest motiv erorile prin rotunjire riscă să devină considerabile, motiv pentru care formulele Newton-Cotes de ordin mare nu sunt indicate. Să ne reamintim că în capitolul 3 am arătat că polinoamele de interpolare de grad mare nu dau rezultate mai precise decât cele de grad mic. Ţinând cont de faptul că formulele de cuadratură Newton-Cotes se deduc prin integrarea polinoamelor de interpolare, vom avea un motiv în plus să evităm formulele Newton-Cotes de ordin mare. Din acest punct de vedere este de preferat utilizarea unei reguli extinse (a trapezului sau Simpson) care are la bază o formulă de cuadratură de ordin mic.
4.4 Cuadratura Gauss
4.5 Cuadratura Gauss-Legendre
Înainte de a prezenta această regulă de cuadratură să descriem pe scurt polinoamele Legendre. Polinoamele Legendre sunt de forma:
P0(x) = 1; P1(x) = x; P2(x) = [3x2 - 1]/2; P3(x) = [5x3 - 3x]/2; Orice polinom Legendre de ordin egal sau mai mare ca doi poate fi dedus cu ajutorul formulei recursive:
nPn(x) - (2n - 1)xPn-1(x) + (n - 1)Pn-2(x) = 0
61
CUADRATURA NUMERICĂ pg 166
Un polinom Legendre de ordin n poate fi exprimat şi de relaţia:
Pn(x) = 12nn!
dn(x2 - 1)n
dxn
O proprietate importantă a polinoamelor Legendre este ortogonalitatea:
⌡⌠-1
1Pm(x)⋅Pn(x) dx = 0, pentru m ≠ n, şi ⌡⌠
-1
1Pm(x)⋅Pn(x) dx =
22n + 1 , pentru m = n
Proprietăţi:
• Orice polinom de grad n-1 sau mai mic este ortogonal cu polinomul Legendre de ordin n, aceasta deoarece orice polinom de grad n-1 sau mai mic poate fi exprimat sub forma unei combinaţii liniare de polinoame Legendre de ordin cel mult n-1.
• orice polinom Legendre, Pn(x) admite n rădăcini, toate situate în intervalul (-1,1). Cuadratura Gauss-Legendre este o metodă de integrare numerică care utilizează care utilizează punctele Legendre (rădăcinile polinoamelor Legendre). Regulile de cuadratură Gauss nu pot fi aplicate pentru integrarea funcţiilor date tabelar în puncte egal depărtate deoarece punctele Legendre nu sunt egal depărtate. Sunt indicate în special pentru integrarea funcţiilor analitice. Avantajul regulilor de cuadratură Gauss constă în precizia semnificativ mai mare faţă de formulele Newton-Cotes. Vom începe discuţia cuadraturii Gauss reluând ideile rezultate la analiza erorii formulelor Newton-Cotes. Ecuaţia (4.6) exprimă faptul că la regula trapezului eroarea este proporţională cu f". Dacă aplicăm regula trapezului pentru a integra funcţiile f = 1, f = x, f = x2, f = x3,..., eroarea va fi egală cu zero pentru f = 1 şi f = x şi diferită de zero pentru f = x2, f = x3,.... Aplicând regula lui Simpson pentru a integra aceleaşi funcţii, eroarea va fi zero pentru f = 1, x, x2, x3 deoarece la această regulă eroarea este proporţională cu fIV. Generalizând, o formulă Newton-Cotes închisă şi de ordin n impar este exactă dacă integrandul este un polinom de ordin cel mult n, în timp ce pentru n par formula este exactă dacă integrandul este un polinom de ordin cel mult n+1. Vom schiţa în continuare modul de generare a unei formule de cuadratură Gauss, exactă în cazul integrării oricărui polinom de grad cel mult trei.
Fie I = ⌡⌠-1
1f(x) dx (4.24). Orice formulă de cuadratură Gauss în două puncte este de forma:
I = w1⋅f(x1) + w2⋅f(x2) + E (4.25)
unde ponderile w1 şi w2 şi punctele x1 şi x2 sunt necunoscute. Deoarece formula de cuadratură pe care dorim să o obţinem trebuie să fie exactă (E = 0) în cazul integrării oricărui polinom de grad cel mult trei, va fi exactă şi în cazul particular al polinoamelor f(x) = 1, x, x2 şi x3. Înlocuind aceste polinoame în (4.25) obţinem următorul sistem neliniar în necunoscutele w1, w2, x1 şi x2:
2 = w1 + w2; 0 = w1x1 + w2x2; 2/3 = w1x12 + w2x2
2; 0 = w1x13 + w2x2
3
Membrii din stânga ai celor patru ecuaţii reprezintă valorile exacte ale integralei I = ⌡⌠-1
1f(x) dx .
Soluţia unică a acestui sistem este x1 = 13 = 0.577350269; x1 = - 1
3 ; w1 = w2 = 1.
62
Cuadratura Gauss-Legendre
Pentru simplificarea calculelor am considerat domeniul de integrare [-1,1]. Acesta poate fi însă transformat în oricare altul printr-o schimbare de variabilă. Formula care am dedus-o în acest mod este cea mai simplă dintre regulile de cuadratură Gauss. Deducerea unor formule de cuadratură Gauss în mai multe puncte prin extinderea tehnicii anterioare ridică probleme deosebite ca urmare a sistemelor neliniare complicate care pot rezulta. Vom descrie în continuare formula generală de cuadratură Gauss de ordin n, exactă atunci când este aplicată pentru integrarea oricărui polinom de grad cel mult 2n-1. Forma generală a regulii de cuadratură Gauss de ordin n pentru domeniul de integrare [-1,1] este:
⌡⌠-1
1f(x) dx = ∑
k=1
nwkf(xk) (4.26)
Să presupunem că f(x) este un polinom de grad cel mult 2n-1. f(x) poate fi scrisă în funcţie de polinomul Legendre Pn(x) de grad n sub următoarea formă: f(x) = c(x)Pn(x) + r(x) (4.27)
unde c(x) şi r(x) sunt polinoame de grad cel mult n-1. Următoarele două aspecte ale demonstraţiei sunt deosebit de importante. Mai întâi, integrând ecuaţia (4.27) pe domeniul [-1,1] obţinem:
⌡⌠-1
1f(x) dx = ⌡⌠
-1
1r(x) (4.28)
unde primul termen s-a redus ca urmare a proprietăţii de ortogonalitate dintre polinomul Legendre de grad n, Pn(x) şi orice alt polinom de grad cel mult n-1. În al doilea rând, dacă înlocuim pe x cu oricare dintre rădăcinile polinomului Legendre, Pn(x)= 0, primul termen din ecuaţia (4.27) devine zero şi ecuaţia (4.27) se reduce la f(xi) = r(xi). Integrandul r(x) din ecuaţia (4.28) fiind un polinom de grad cel mult n-1, poate fi exprimat printr-un polinom de interpolare Lagrange de grad n-1 (vezi ecuaţia (3,4)):
r(x) = ∑k=1
n
⎣⎢⎡
⎦⎥⎤
∏j=1 j≠k
nx - xjxk - xj
r(xk) (4.29)
Întrucât xk, k = 0,1,...,n sunt rădăcinile polinomului Legendre Pn(x) şi r(xi) = f(xi), ţinând cont de ecuaţia (4.27) ecuaţia (4.29) devine:
r(x) = ∑k=1
n
⎣⎢⎡
⎦⎥⎤
∏j=1 j≠k
nx - xjxk - xj
f(xk) (4.30)
Formula de cuadratură Gauss rezultă înlocuind (4.30) în (4.28):
63
CUADRATURA NUMERICĂ pg 166
⌡⌠-1
1f(x) dx) = ∑
k=1
n
f(xk)
⌡⎮⌠
-1
1
⎣⎢⎡
⎦⎥⎤
∏j=1 j≠k
nx - xjxk - xj
dx (4.31)
Definind
wk =
⌡⎮⌠
-1
1
⎣⎢⎡
⎦⎥⎤
∏j=1 j≠k
nx - xjxk - xj
dx (4.32)
rezultă echivalenţa dintre ecuaţiile (4.31) şi (4.26). În formula (4.26) n este numărul punctelor Gauss, wi sunt ponderile şi xi sunt punctele Gauss date în tabelul 4.3. Semnele ± din tabel indică simetria punctelor x care sunt date în pereche, unul cu semnul plus şi celălalt cu minus.
Tabelul 4.3 Punctele şi ponderile Gauss pentru ∫f(x)dx
+ 0,652145155⋅f(0,339981043)+ 0,347854845⋅f(0,861136312) De obicei regulile de cuadratură Gauss sunt date tabelar pentru un domeniu fix (în tabelul 4.3, [-1,1]). Să considerăm cunoscută formula de cuadratură:
⌡⌠α
βf(x) dx = ∑
k=1
nwkf(xk) + Rn (4.33)
(în tabelul 4.3: α = -1 şi β = 1). Dacă integrala pe care dorim să o evaluăm este definită pe intervalul [α,β], valorile xi şi wi date tabelar se vor prelua direct. În cazul în care dorim să
64
Cuadratura Gauss-Kronrod
calculăm o integrală definită pe un interval oarecare [a,b] diferit de [α,β], I = ⌡⌠a
bg(t) dt , valorile xi
şi wi date tabelar vor putea fi utilizate după aplicarea unei schimbări de variabilă t → x (care transformă intervalul [a,b] în [α,β]) de forma:
t = (b - a)x + aβ - bα(β - α)
Utilizând această transformare, integrala poate fi scrisă ca:
I = ⎝⎜⎛
⎠⎟⎞b - a
β - α ⌡⎮⌠
α
β
g⎝⎜⎛
⎠⎟⎞(b - a)x + aβ - bα
(β - α) dx = ⎝⎜⎛
⎠⎟⎞b - a
β - α ∑k=1
n
wkg⎝⎜⎛
⎠⎟⎞(b - a)xk + aβ - bα
(β - α) + ℜn (4.34)
pentru α = -1 şi β = 1,
I = ⎝⎜⎛
⎠⎟⎞b - a
2 ⌡⎮⌠
-1
1
g⎝⎜⎛
⎠⎟⎞(b - a)x + a + b
2 x = ⎝⎜⎛
⎠⎟⎞b - a
2 ∑k=1
n
wkg⎝⎜⎛
⎠⎟⎞(b - a)xk + a + b
2 + Rn (4.35)
Această transformare prezintă proprietatea de liniaritate: dacă P(t) este un polinom de grad s, în
variabila t, atunci P⎝⎜⎛
⎠⎟⎞(b - a)x + aβ - bα
(β - α) este un polinom de acelaşi grad s, în variabila x. În
acelaşi timp, dacă Rn este zero atunci când f(x) este un polinom de grad s, ℜn va fi de asemenea egal cu zero pentru g(x) polinom de grad s.
4.6 Cuadratura Gauss-Kronrod
La metoda extinsă a trapezului am văzut că aplicând această regulă pentru h= (b-a)/n, unde n este numărul de puncte, se obţine aproximarea Tn. Prin dublarea numărului de subintervale, altfel spus pentru h'= h/2, rezultă evaluarea T2n-1, aproximativ de patru ori mai precisă decât Tn. Calculele necesare pentru Tn implică n evaluări funcţionale iar calculele necesare pentru T2n-1 implică numai n-1 evaluări funcţionale şi nu 2n-1 a_a cum ne-am aştepta poate la prima vedere. Într-adevăr, dacă h este distanţa dintre nodurile adiacente pentru Tn şi h'= h/2 pentru T2n-1, atunci Aplicând metoda de cuadratură Romberg pentru evaluarea unei integrale vor fi necesare în total (2n-1) evaluări funcţionale. Am remarcat în subcapitolul . că regulile de cuadratură Gauss sunt mai precise. Se pune problema dacă este posibil să se determine o pereche de reguli de cuadratură Gauss care, la fel ca şi în integrarea Romberg să permită calculul atât a unei integrale oarecare cât şi a erorii, cu un volum minim de evaluări funcţionale. Aplicând o regulă de cuadratură Gauss în n puncte (n-1 subintervale) rezultă pentru integrala considerată aproximarea Gn pentru care sunt necesare n evaluări funcţionale. Dublând numărul de subintervale obţinem pentru aceeaşi integrală aproximarea G2n-1, cu o precizie mai bună decât Gn iar pentru estimarea erorii putem utiliza |Gn-G2n-1|. La metoda extinsă a trapezului am văzut că pentru T2n-1 sunt necesare numai n-1 evaluări funcţionale deoarece Tn şi T2n-1 au n puncte comune în care se calculează evaluările funcţionale. Regulile de cuadratură Gauss sunt
65
CUADRATURA NUMERICĂ pg 166
disjuncte deoarece nu există două reguli Gauss care să utilizeze noduri comune (cu excepţia cazului n impar când există un nod comun, xi= 0). Acest fapt implică pentru aproximarea G2n-1 2n-1 evaluări funcţionale. Volumul total de evaluări funcţionale pentru Gn şi G2n-1 va fi 3n-1 sau 3n-2 pentru n impar. Pornind de la o regulă de cuadratură Gauss în n puncte Gn,
,)f( xw = G kk
n
1=kn ∑ ( )
matematicianul rus Kronrod a determinat o altă regulă K2n+1 de forma
,)f( +)f( ybxa = K ii
1+n
1=ikk
n
1=k1+2n ∑∑ ( )
care prezintă n puncte comune cu Gn (xk, k= 1,2...,n), n+1 puncte suplimentare şi ponderile distincte ak şi bi. Gn este de grad polinomial 2n-1. Valorile ak, bi şi yi au fost determinate astfel încât K2n+1 să fie de grad polinomial 3n+1. Cele două reguli (Gn, K2n+1) poartă denumirea de pereche Gauss - Kronrod şi volumul de calcule pentru o astfel de pereche necesită numai 2n+1 evaluări funcţionale şi nu 3n-1 sau 3n-2 ca pentru o pereche de genul (Gn, G2n+1). Aplicând această metodă vom pute utiliza K2n+1 pentru estimarea integralei considerate şi |Gn-K2n+1| pentru eroare. Din experienţă s-a recunoscut că |Gn-K2n+1| supraestimează eroarea şi s-a ajuns în cele din urmă pentru evaluarea erorii la relaţia empirică
la prima vedere aceasta pare mai mare decât |Gn-K2n+1| dar în realitate |Gn-K2n+1| este mult mai mic decât 1 şi prin ridicare la puterea 1,5 devine şi mai mic. Tabelul . conţine perechea Gauss -Kronrod (G7, K15). Se observă că cele 7 puncte Gauss se reg_sesc şi printre cele 15 puncte Konrod, evident cu ponderi corespunzătoare diferite.
Tabelul Perechea Gauss-Konrod (7,15) pentru ∫f(x)dx ═════════════════════════════════╤════════════════════════════
Cele 7 puncte Gauss (x1-x7) Ponderile w1-w7 ═════════════════════════════════╪════════════════════════════
Perechea Gauss-Kronrod (G7, K15) s-a dovedit practică ca fiind una dintre cele mai eficiente metode de calcul pentru integrale, mai ales atunci când este implementată într-un algoritm de cuadratură adaptiv (vezi subcapitolul ). Programul implementează această regulă de cuadratură.
4.7 Algoritmul de cuadratură automaţi şi adaptivi
Un algoritm de cuadratură automat necestă ca date de intrare funcţia f, [a,b] şi toleranţa admisă ε iar ca date de ieşire va genera rezultatul integrării Q şi estimarea erorii E care trebuie să îndeplinească condiţia
ε<E|<I-Q|
Un astfel de algoritm poate fi dezvoltat plecând de la regula extinsă a trapezului. Va calcula Tn cu n= 2,3,5,9,17 până când
. |<T-T|34
n1-2n ε
Pentru un astfel de algoritm de cuadratură punctele de evaluare ale funcţiei f sunt uniform distribuite pe intervalul [a,b]. Evident că în cazul în care funcţia f prezintă variaţii funcţionale mici în anumite zone ale intervalului [a,b], un astfel de algoritm va genera multe calcule inutile. În astfel de situaţii se utilizează un algoritm de cuadratură adaptiv care are proprietatea că îşi adaptează distanţa h dintre nodurile adiacente în funcţie de modul de variaţie al funcţiei f: h mare în zonele cu variaţii funcţionale mici şi mic în zonele cu variaţii funcţionale mari, astfel încât să se evite calculele inutile şi eroarea să fie uniform distribuită. Ideea de bază a acestui algoritm constă în înlocuirea intervalului iniţial [a,b] printr-un set de subintervale de mărimi diferite. Printr-una din metodele cunoscute se estimează pe fiecare subinterval integrala funcţiei şi eroarea locală. Dacă suma tuturor erorilor locale este mai mare decât toleranţa iniţial impusă ε, subintervalul pe care a rezultat cea mai mare eroare locală este subdivizat şi pe fiecare dintre cele două subdiviziuni se recalculează integrala funcţiei şi erorile locale. Dacă suma tuturor erorilor locale rezultă mai mică decât toleranţa impusă ε se
67
CUADRATURA NUMERICĂ pg 166
returnează ca aproximare finală pentru integrală, suma estimărilor integralelor locale şi suma erorilor locale ca eroarea totală.
Să considerăm integrala pe care dorim să o aproximăm cu o eroare specificată ε > 0. x d f(x)b
a∫
Primul pas al procedurii constă în aplicarea, de exemplu, a regulii Simpson (1/3) cu mărimea pasului h= (b-a)/2. Rezultă
,b)(a, pentru ,)(f90h-b) S(a,= x d f(x) IV
5b
a
εξξ∫ ( )
unde S(a,b) este formula lui Simpson: ( )
Următorul pas constă în a găsi o modalitate de estimare a preciziei acestei aproximări fără să fie necesar să calculăm fIV(ξ). Pentru aceasta vom aplica regula extinsă a lui Simpson pe [a,b,], cu mărimea pasului h= (b-a)/4. Obţinem ( )
pentru ξ oarecare în (a,b). Ecuaţia ( ) poate fi scrisă într-o formă echivalentă mai simplă ( )
unde ţi Pentru a estima eroarea vom considera ξ= ξ, altfel spus fIV(ξ)= fIV(ξ). Rezultatul va depinde de corectitudinea acestei ipoteze. Dacă ipoteza este îndeplinită cu o bună precizie, atunci ecuaţiile ( ) şi ( ) implică de unde rezultă Înlocuind această relaţie în ecuaţia ( ) obţinem următoarea aproximaţie pentru eroare:
b),2
b+aS(-)2
b+aS(a,-b)S(a,151
b),2
b+aS(-)2
b+aS(a,-x d f(x) b
a
≈
≈∫ ( )
Aproximaţia astfel obţinută pentru eroare trebuie să fie mai mică decât toleranţa impusă ε. În consecinţă, dacă membrul al doilea al ecuaţiei ( ) este mai mic decât 15ε, atunci
b),2
b+aS(+)2
b+aS(a,
68
Algoritmul de cuadratură automaţi şi adaptivi
este o aproximaţie acceptabilă pentru
. x d f(x)b
a∫
Dacă nu este îndeplinită inegalitatea
,15 <b),2
b+aS(-)2
b+aS(a,-b)S(a, ε ( )
se repetă procedura de estimare a erorii succesiv pe fiecare dintre subintervalele [a,(a+b)/2] şi [(a+b)/2,b] pentru a afla dacă aproximările locale ale integralei se fac în limita unei toleranţe de
ε/2. Dacă da, suma acestor aproximări locale va reprezenta estimarea pentru cu o
toleranţă ε. Dacă aproximarea locală pe unul dintre subintervale nu rezultă în limita unei toleranţe de ε/2, acel subinterval se va diviza la rândul său şi se va analiza dacă aproximările locale pe fiecare dintre subdiviziuni sunt în limita unei toleranţe de ε/4. Această procedură de înjumătăţire se va repeta până când se vor obţine pe toate subintervalele aproximări locale ale integralei în limita toleranţei corespunzătoare.
x d f(x)b
a∫
Următorul algoritm implementează detaliat metoda de cuadratură adaptivă. La pasul 1 toleranţa luată în considerare este 30ε şi nu 15ε pentru a acoperi eroarea care se face prin ipoteza fIV(ξ)=fIV(ξ). În situaţiile în care se cunoaşte că fIV prezintă variaţii mari această toleranţă poate fi luată chiar şi mai mică. Algoritmul aproximează întotdeauna integralele locale începând cu primul subinterval de la stânga, ceea ce implică existenţa unei proceduri de stocare şi apelare a evaluărilor funcţionale în modurile incluse în subintervalele din dreapta. Pentru a evita posibilitatea de a se ajunge la subintervale de lungimi excesiv de mici se cere ca dată de intrare şi numărul de înjumătăţiri, N.
4.7.1 Algoritm de cuadratură adaptiv [ ]
Aproximează integrala cu o toleranţă dată TOL > 0 x d f(x) = Ib
a∫
DATE DE INTRARE: extremităţile a şi b ; toleranţa TOL ; numărul maxim de înjumătăţiri N. DATE DE IEŞIRE: aproximaţia AP pentru ∫f(x)dx sau mesaj "S-a depăşit numărul maxim de înjumătăţiri" Pasul 1. INIŢIALIZEAZĂ AP= 0 ; i= 1 TOLi= 10 TOL ; ai= a ; hi= (b-a)/2 ;
69
CUADRATURA NUMERICĂ pg 166
FAi= f(a) ; FCi= f(a+hi) ; FBi= f(b) ; Si= hi (FAi+4FCi+FBi)/3 (aproximare prin metoda Simpson pe întreg intervalul) INJi= 1. Pasul 2. CÂT TIMP i> 0 EFECTUEAZĂ paşii 3,4 şi 5: Pasul 3. CALCULEAZĂ FD= f(ai+hi/2) ; FE= f(ai+3hi/2) ; S1= hi(FAi+4FD+FCi)/6 ; (aproximaţiile prin metoda Simpson pentru jumătăţile subintervalului) S2= hi(FCi+4FE+Fbi)/6 ; v1= ai ; (salvează datele la acest nivel) v2= FAi ; v3= FCi ; v4= FBi ; v5= hi ; v6= TOLi ; v7= Si ; v8= Li . Pasul 4. i= i-1 (Şterge nivelul) Pasul 5. DACĂ |S1+S2-v7| < v6 ATUNCI AP= AP+(S1+S2) ALTFEL DACĂ v8 ≥ N ATUNCI SCRIE ("Numar maxim de înjumătăţiri depăşit") STOP. ALTFEL (Adaug_ un nivel) i= i+1 (Datele pentru jumătatea ai= v1+v5 din stânga a subintervalului) FAi= v3 ; FCi= FE ; FBi= v4 ; hi= v5/2 ; TOLi= v6/2 ; Si= S2 ; Li= v8+1 ; şi i= i+1; (Datele pentru jumătatea din ai= v1 ; dreapta a subintervalului) FAi= v2 ; FCi= FD ; FBi= v3 ; hi= hi-1 ; TOLi= TOLi-1 ; Si= S1 ; Li= Li-1 ; Pasul 6. SCRIE (AP) ; (Ap aproximează I în limita STOP . toleranţei TOL)
70
Integrarea funcţiilor date în puncte
4.8 Integrarea funcţiilor date în puncte
Fie o secvenţă de date (puncte) de forma (xi,yi), x1 < x2 < …< xn rezultată în urma unor măsurători experimentale sau generată printr-o metodă oarecare. Acest set de date poate reprezenta valorile exacte ale unei funcţii necunoscute f(x): yi = f(xi) sau, situaţie frecventă mai ales în cazul măsurătorilor experimentale, poate fi afectat de erori: yi = fr(xi) + ei
Considerând primul caz, vom încerca în acest capitol să analizăm problema rezolvării integralei
I = ⌡⌠x1
xn
f(x) dx .
Secvenţa de date nefiind afectată de erori, o abordare rezonabilă a problemei constă în determinarea unui interpolant care să verifice datele, altfel spus să determinăm o funcţie g(x) de
formă cunoscută, astfel încât g(xi) = yi, după care să calculăm I ≈ I' = ⌡⌠x1
xn
g(x) dx
Putem utiliza ca interpolanţi polinoame unice, polinoame cubice Hermite sau funcţii spline. (În cazul în care datele sunt afectate de erori interploarea este inadecvată şi se poate întrebuinţa o metodă de aproximare, de exemplu metoda celor mai mici pătrate) Vom aborda problema considerând pentru început că g(x) este un interpolant polinomial liniar pe subintervale, ceea ce înseamnă că graficul funcţiei g(x) uneşte perechile de puncte succesive prin segmente de dreaptă. Notând acest polinom cu L(x), va fi definit prin relaţia
yx-x
x-x+yx-xx-x = L(x) g(x) 1+i
i1+i
ii
1+ii
1+i≡ (4.43)
pentru x ∈ [xi,xi+1], i= 1,2,...,n-1. Integrala pe care dorim să o calculăm va fi
(4.44) x d L(x) = x d L(x) = I Ix
x
1-n
1=i
x
x
1+i
i
n
1
∫∑∫′≈
Cunoscând forma interpolantului L(x), prin inte grare obţinem în final
Se observă că în cazul în care punctele sunt uniform distribuite, formula anterioară exprimă regula extinsă a trapezului. O altă posibilitate de abordare a problemei, cu rezultate mai bune practic, o constituie utilizarea ca interpolant a unui polinom cubic Hermite de forma (vezi capitolul .):
(x)cd+(x)cy = C(x) g(x) iiii
n
1=iˆ∑≡
unde am notat cu di derivata în punctul i. Rezultă următoarea relaţie pentru integrală:
(4.46) x d (x)+x d cdcy = x d C(x) = I I i
x
x
i
n
1=ii
x
xi
n
1=i
x
x
n
1
n
1
n
1
ˆ∫∑∫∑∫′≈
71
CUADRATURA NUMERICĂ pg 166
Cunoscând relaţiile de definiţie pentru ci şi ci (vezi ecuaţia ( )) obţinem prin integrare
(4.47) )d+y( = I iiii
n
1=iβα∑′
unde
⎪⎪⎪⎪
⎩
⎪⎪⎪⎪
⎨
⎧
n =i pentru ,2x-x
1-n 2,3,..., =i pentru ,2
x-x
1 =i pentru ,2
x-x
=
1-nn
1-i1+i
12
iα şi
⎪⎪⎪⎪
⎩
⎪⎪⎪⎪
⎨
⎧
n =i pentru ,12
)x-x(-
1-n 2,3,..., =i pentru ,12
)x+x2-x)(x-x(
1 =i pentru ,12
)x-x(
=
21-nn
1-ii1+i1-i1+i
212
iβ
Practic nu cunoaştem valorile derivatei di dar le putem determina printr-o interpolare spline sau printr-una din metodele descrise la capitolul 5.
În cazul punctelor uniform distribuite în formula (4.47) termenii care conţin βi devin egali cu zero, cu excepţia β1 şi β2, (4.47) reducându-se în final la forma:
I' = h⎣⎢⎡
⎦⎥⎤y1
2 + y22 + … +
yn-12 +
yn2 +
h12(d1 - dn) (4.48)
Dacă modelul fizic sugerează periodicitatea setului de date cu perioada (xn - x1), atunci d1 = dn şi relaţia (4.48) devine identică cu formula extinsă a trapezului.
4.9 Integrale duble
Ne propunem în acest capitol să analizăm posibilităţilor de rezolvare numerică a integralelor duble. Fie
A d y)f(x, = I .
A
∫∫
Vom descrie tehnicile de rezolvare a integralelor duble definite pe un domeniu A rectangular
{ }dyc ,bxa |y)(x, = A ≤≤≤≤
de unde
x dy d y)f(x, = I d
c
b
a∫∫
sau a integralelor definite pe un domeniu la care frontierele din stânga şi dreapta sunt linii verticale iar frontierele superioară şi inferioară sunt definite de curbele y= d(x) respectiv y= c(x):
x dy d y)f(x, = I d(x)
c(x)
b
a∫∫ ( )
Aceste tehnici pot fi modificate pentru a fi utilizate şi la rezolvarea integralelor duble pe alte domenii sau a integralelor multiple. Fie integrala ( ). Dacă definim
72
Integrale duble
y , d y)f(x, = G(x) d
c∫ ( )
relaţia ( ) devine
x d G(x) = I b
a∫ ( )
Astfel problema s-a redus la o combinaţie de două integrale simple care pot fi rezolvate prin oricare din metodele descrise anterior. Soluţia integrării numerice a ecuaţiei ( ) poate fi scrisă sub forma:
)G( xw I ii
n
1=i∑_
unde wi şi xi sunt ponderile, respectiv punctele regulii de cuadratură considerate. Valorile pentru G(xi) se evaluează tot numeric. Pentru x= xi ecuaţia ( ) devine:
y , d y),xf( )xG( i
d
ci ∫≡ ( )
care este tot o integrală simplă deoarece integrandul este funcţie numai de y; şi această integrală poate fi rezolvată prin oricare din metodele de cuadratură cunoscute. Ca aplicaţie să rezolvăm integrala dublă: ( )
utilizând regula extinsă a lui Simpson (1/3). Domeniul de integrare [c,d] este divizat în 2m subintervale egale de mărime hy= (d-c)/2m iar domeniul de integrare [a,b] este divizat în 2n subintervale egale de mărime hx= (b-a)/2n (vezi fig în care n= m= 2) i=0 i=1 i=2 i=3 i=4
a=x0 x1 x2 x3 x4=b Fig Domeniul de integrare bidimensional (rectangular) Vom utiliza mai întâi regula extinsă a lui Simpson pentru a evalua integrala
y , d y)f(x, d
c∫
în care considerăm x= constant. Pentru yj= c+jhy, unde j= 0,1,...,2m obţinem (neglijând eroarea): Astfel
. x d )yf(x,3h+x d )yf(x,
3h4
+
+x d )yf(x, 3h2
+x d )yf(x,3h x dy d y)f(x,
2m
b
a
y1-2j
b
a
m
1=j
y
2j
b
a
1-m
1=j
y0
b
a
yd
c
b
a
∫∫∑
∫∑∫∫∫ _
( )
În continuare se va aplica regula extinsă a lui Simpson (1/3) pentru rezolvarea fiecăreia dintre integralele conţinute în membrul din dreapta al ecuaţiei anterioare. Pentru xi= a+ih, unde i= 0,1,...,2n, obţinem pentru fiecare j= 0,1,...,2m ( )
Înlocuind ( ) în ( ) rezultă: ( )
Numărul evaluărilor funcţionale este egal cu produsul dintre numărul evaluărilor funcţionale necesar la integrarea pe direcţia x şi cel necesar la integrarea pe direcţia y. Pentru a reduce numărul mare al evaluărilor funcţionale rezultat prin aplicarea formulelor de tip Newton-Cotes se pot utiliza metode mai eficiente cum ar fi cuadratura Gauss, Gauss-Kronrod sau un algoritm de cuadratură adaptiv. Să încercăm să rezolvăm integrala ( ) prin metoda Gauss. Primul pas constă în modificarea domeniului de definiţie:
} dycb,xa|y){(x, = A ≤≤≤≤
în
} 1v1,-1u-1|v){(u, = A ≤≤≤≤~
Transformarea liniară care realizează această translaţie este dată de relaţiile:
c-dd-c-2y = v si
a-bb-a-2x = u
Cu această schimbare de variabilă obţinem o integrală care poate fi rezolvată prin metoda de cuadratură Gauss:
74
Integrale duble
v d u d v)f(u, 4
c)-a)(d-(b = x dy d y)f(x, 1
1-
1
1-
d
c
b
a∫∫∫∫ ( )
Considerând o formulă de cuadratură Gauss în n puncte vom prelua dintr-un tabel de tipul punctele ui, vi şi ponderile wi, wj. Astfel
),f( vuww 4
c)-a)(d-(b I jiji
n
0=j
n
0=i
••∑∑_ ( )
Evident că putem alege două reguli de cuadratură diferite în n puncte pentru variabila u şi în m puncte pentru variabila v. Într-o asemenea situaţie cele duouă sume din ecuaţia ( ) vor avea limite superioare diferite. Metodele de aproximare a integralelor duble nu se limitează numai la integralele definite pe domenii rectangulare. Aceste tehnici pot fi modificate şi pentru aproximarea integralelor de forma:
x dy d y)f(x, d(x)
c(x)
b
a∫∫ ( )
Pentru a rezolva o astfel de integrală putem alege oricare dintre metodele de cuadratură cunoscute pentru rezolvarea integralelor simple. Ca aplicaţie vom considera formula elementară a lui Simpson (1/3). Mărimea pasului pentru variabila x este h = (b-a)/2, în timp ce mărimea pasului pentru variabila y variază funcţie de x ( vezi fig ):
2c(x)-d(x) = (x)hy
Aplicând aceleaşi raţionament ca şi la integralele duble diferite pe domenii rectangulare rezultă următoarea formulă de aproximare Următorul algoritm implementează regula extinsă a lui Simpson pentru aproximarea unei integrale de forma ( ) .
4.9.1 Algoritm de calcul a integralelor multiple [ ]
Aproximează integrala . x dy d y)f(x, = I d(x)
c(x)
b
a∫∫
DATE DE INTRARE: limitele a şi b ; întregii pozitivi m şi n . DATE DE IEŞiRE: aproximarea I a integralei Pasul 1. Calculează h= (b-a)/2n. Pasul 2. Iniţializează I1= 0 ; (Termenii de la extremităţi) I2= 0 ; (Termenii pari) I3= 0 ; (Termenii impari) Pasul 3. Pentru i= 0,1,...,2n
75
CUADRATURA NUMERICĂ pg 166
Iniţializează x= a+ih (Regula extinsă Simpson pentru x fixat) HX= (d(x)-c(x))/(2m) ; K1= f(x,c(x))+f(x,d(x)) ; (Termenii de cap_t pentru fiecare x) K2= 0 ; (Termenii pari pentru fiecare x) K3= 0 ; (Termenii impari pentru fiecare x) pentru j= 1,2,...,2m-1 calculează y= c(x)+jHX ; z= f(x,y) ; dacă j este impar atunci K2= K2+Z altfel K3= K3+Z ; Calculează L=(K1+2K2+4K3)HX/3 ; dacă i= 0 sau i= 2n atunci I1= I1+L altfel dacă i este par I2= I2+L altfel I3= I3+L Pasul 4. Calculează I= (I1+2I2+4I3)h/3 Pasul 5. SCRIE (I) ; STOP └─────────────────────────────────────────────────────────────
4.9.1.1 Exemplul Aproximarea integralelor duble Să se evalueze următoarea integrală dublă: aplicând formula elementară Simpson (1/3). Limitele de integrare sunt:
)5x(+3 = d(x)
,x = c(x)
3 = b ,1 =a
exp
ln
4.9.1.1.1 Rezolvare Pentru formula elementară Simpson (1/3), absciselex sunt: x0= 1, x1= 2, x2= 3
76
Integrale duble
Vezi fig în care s-au reprezentat domeniul de integrare şi punctele re_elei Fig Reţea pentru integare dublă Aplicând formula elementară Simpson (1/3) primei integrale obţinem: unde hx= (b-a)/2= 1 şi
y , d y)+(x = )xG( )
15x( +3
(x) i cos
exp
ln∫
de unde rezultă Aplicăm formula elementară Simpson (1/3) primei integrale din paranteza dreaptă: Calcule similare ne conduc la următoarele rezultate pentru celelalte două integrale:
0,81097 y d y)+(3
1,18350- y d y)+(2
4,8221
1,0986
4,44918
0,6931
≈
≈
∫
∫
cos
cos
Valoarea finală pentru integrală este:
77
Derivarea (diferenţierea) numerică
5. Derivarea (diferenţierea) numerică
5.1 Introducere
Diferenţierea numerică este utilizată pentru evaluarea derivatelor unei funcţii în cazul în care valorile funcţionale sunt date în puncte; altfel spus, nu se cunoaşte expresia analitică a funcţiei. Să considerăm o funcţie f(x) căreia vrem să îi calculăm valoarea primei derivate, f'(x) în punctul x = x0. În cazul în care se cunosc valorile funcţiei în punctele x = x0, x = x0 + h, x = x0 - h, f'(x0) poate fi aproximat printr-una din următoarele formule : - aproximarea directă prin diferenţe:
f'(x0) = f(x0+h) - f(x0)
h (5.1)
- aproximarea inversă prin diferenţe:
f'(x0) = f(x0) - f(x0-h)
h (5.2)
- aproximarea centrală prin diferenţe:
f'(x0) = f(x0+h) - f(x0-h)
2h (5.3)
Cele trei formule sunt trei modalităţi de aproximare ale tangentei la graficul funcţiei f(x) în punctul x0. Dintre modalităţile de determinare a formulelor de diferenţiere numerică vom analiza:
• diferenţierea polinoamelor de interpolare şi
• metoda bazată pe dezvoltarea în serie Taylor a unei funcţii. Vom descrie de asemenea în acest capitol un algoritm de generare a formulelor de diferenţiere.
5.2 Diferenţierea polinomului de interpolare Lagrange
Fie n+1 numere distincte x0,x1,...,xn aparţinând unui interval oarecare I şi o funcţie f ∈ Cn+1(I). Conform celor studiate la cap. 3,
f(x) = ∑k=0
nLn k(x) f(xk) +
(x - x0)…(x - xn)(n+1)! f(n+1)(ξ(x))
78
Diferenţierea polinomului de interpolare Lagrange
pentru ξ(x) ∈ I, unde Ln,k(x) este coeficientul Lagrange k. Derivând această expresie obţinem:
f'(x) = ∑k=0
nL'n k(x) f(xk) + Dx⎣⎢
⎡⎦⎥⎤(x - x0)…(x - xn)
(n+1)! f(n+1)(ξ(x)) + (x - x0)…(x - xn)
(n+1)! Dx[f(n+1)(ξ(x))](5.4)
În cazul în care x = xk, oricare ar fi k = 0,1,...,n, termenul Dx[f(n+1)(ξ(x))] este egal cu zero şi formula ( ) devine
f'(xk) = ∑k=0
nL'n k(x) f(xk) +
f(n+1)(ξ(xk))(n+1)! ∏
j=0 j≠i
n(xk - xj) (5.5)
Ecuaţia (5.5) poartă denumirea de formulă de aproximare pentru f'(xi) în n+1 puncte, deoarece s-a utilizat o combinaţie liniară de n+1 valori f(xj), j = 0,1,...,n. În general, utilizarea unui număr mare de puncte în ecuaţia ( ) va genera o precizie mai mare, chiar dacă numărul mare al evaluărilor funcţionale şi creşterea erorii prin rotunjire nu ne încurajează în această direcţie. Cele mai des utilizate formule sunt cele în trei sau cinci puncte de evaluare. Să deducem mai întâi câteva formule utile în trei puncte. În acest caz,
L2,0(x) = (x - x1)(x - x2)
(x0 - x1)(x0 - x2) ; L'2,0(x) = 2x - x1 - x2
(x0 - x1)(x0 - x2)
Similar
L'2,1(x) = 2x - x0 - x2
(x1 - x0)(x1 - x2) ; L'2,2(x) = 2x - x0 - x1
(x2 - x0)(x2 - x1)
Înlocuind coeficienţii Lagrange derivaţi în ecuaţia (5.5), obţinem:
f'(xj) = ⎣⎢⎡
⎦⎥⎤2xj - x1 - x2
(x0 - xn)(x0 - x2) f(x0) + ⎣⎢⎡
⎦⎥⎤2xj - x0 - x2
(x1 - x0)(x1 - x2) f(x1)
+ ⎣⎢⎡
⎦⎥⎤2xj - x0 - x1
(x2 - x0)(x2 - x1) f(x1) + 16 f(3)(ξ(xj)) ∏
i=0 i≠j
2(xj - xi) (5.6)
pentru j = 0,1,2. Să analizăm cele trei formule date de ecuaţiile (5.6) în cazul nodurilor egal depărtate, altfel spus, atunci când x1 = x0 + h şi x2 = x0 + 2h pentru h ≠ 0. Pentru xj = x0, x1 = x0 + h şi x2 = x0 + 2h, ecuaţia (5,6) devine
f'(x0) = 1h⎣⎢⎡
⎦⎥⎤-
32f(x0) + 2f(x1) -
12f(x2) +
h2
3 f(3)(ξ0)
În mod asemănător, pentru xj = x1 obţinem
f'(x1) = 1h⎣⎢⎡
⎦⎥⎤-
12f(x0) +
12f(x2) -
h2
6 f(3)(ξ1)
şi pentru xj = x2,
f'(x2) = 1h⎣⎢⎡
⎦⎥⎤1
2f(x0) - 2f(x1) + 32f(x2) +
h2
3 f(3)(ξ2)
Deoarece x1 = x0 + h şi x2 = x0 + 2h, aceste formule pot fi exprimate şi sub următoarea formă:
f'(x0) = 1h⎣⎢⎡
⎦⎥⎤-
32f(x0) + 2f(x0+h) -
12f(x0+2h) +
h2
3 f(3)(ξ0) (5.7)
79
Derivarea (diferenţierea) numerică
f'(x0+h) = 1h⎣⎢⎡
⎦⎥⎤1
2f(x0) + 12f(x0+2h) -
h2
6 f(3)(ξ1) (5.8)
f'(x0+2h) = 1h⎣⎢⎡
⎦⎥⎤1
2f(x0) - 2f(x0+h) + 32f(x0+2h) +
h2
3 f(3)(ξ2) (5.9)
Dacă în relaţiile (5,8) şi (5.9) facem substituţiile variabile x0 + h = x0, respectiv x0 + 2h = x0, rezultă următoarele trei formule de aproximare pentru f'(x0):
Generalizând, obţinem următoarele formule în trei puncte pentru calculul primei derivate a unei funcţii într-un punct oarecare xi:
f'(xi) = 1
2h[ ]-3f(xi) + 4f(xi+h) - f(xi+2h) + h2
3 f(3)(ξ) (5.13)
unde ξ este cuprins între xi şi xi + 2h;
f'(xi) = 1
2h[ ]-f(xi-h) + f(xi+h) - h2
6 f(3)(ξ) (5.14)
unde ξ este cuprins între xi - h şi xi + h;
f'(xi) = 1
2h[ ]f(xi-2h) - 4f(xi-h) + 3f(xi) + h2
3 f(3)(ξ) (5.15)
unde ξ este cuprins între xi - 2h şi xi. Se observă că eroarea dată de formula (5.14) este aproximativ de două ori mai mică decât eroarea dată de formulele (5.13) şi (5.15). Prezentăm în continuare, fără demonstraţie, formulele în cinci puncte pentru calculul primei derivate a unei funcţii într-un punct oarecare xi:
Cunoscând expresia analitică a funcţiei, putem calcula valoarea exactă a derivatei într-un punct: f'(x) = (x+1)ex, f'(2,0) = 22,167168. Aproximarea derivatei f'(2,0) utilizând diferite formule în trei şi cinci puncte conduce la următoarele rezultate: Formule în trei puncte:
Aplicând (5.13) cu h = 0,1: 1/0.2⋅[-3f(2.0) + 4f(2.1) - f(2.2)] = 22.032310
Aplicând (5.14) cu h = 0,1: 1/0.2⋅[f(2.1) - f(1.9)] = 22.228790
Aplicând (5.14) cu h = 0,2: 1/0.4⋅[f(2.2) - f(1.8)] = 22.414163 Formula în cinci puncte Aplicând ( ) cu h = 0,1:
Erorile generate prin aplicarea acestor formule sunt 1,35⋅10-1, -6,16⋅10-2, -2,47⋅10-1 şi respectiv 1,69⋅10-4. Evident că formula în cinci puncte generează rezultatele cele mai precise. Remarcaţi de asemenea că eroarea dată de formula (5.14) cu h = 0,1 est aproximativ jumătate, ca ordin de mărime, faţă de eroarea dată de formula (5.13) cu h = 0,1. Un subiect important în studiul diferenţierii numerice îl constituie efectul erorii prin rotunjire. Să examinăm formula (5.11):
f'(x0) = 1
2h[ ]-f(x0-h) + f(x0+h) - h2
6 f(3)(ξ(x1))
Să presupunem că la evaluarea valorilor funcţionale f(x0+h) şi f(x0-h) au intervenit erorile prin rotunjire e(x0+h), respectiv e(x0-h). Aceasta înseamnă că între valorile calculate f~(x0+h) şi f~(x0-h) şi valorile exacte f(x0+h) şi f(x0-h) există următoarele relaţii:
f(x0+h) = f~(x0+h) + e(x0+h) şi f(x0-h) = f~(x0-h) + e(x0-h) În acest caz eroarea rezultată în aproximarea:
)(f6h-
2hh)-xe(-h)+xe( =
2hh)-x(f-h)+x(f-)x(f 1
(3)2
00000 ξ
~~′
va conţine o parte datorată rotunjirii şi o parte datorată trunchierii. Dacă presupunem că erorile prin rotunjire e(x ± h) sunt mai mici decât un număr oarecare ε > 0, şi că a treia derivată este mai mică decât un număr oarecare M > 0, atunci
. M6h+
h2hh)-x(f-h)+x(f-)x(f
200
0ε
≤′~~
Dacă h este mic, eroarea prin rotunjire, ε/h poate lua valori mari. Rareori este avantajos să alegem valori foarte mici pentru h deoarece erorile prin rotunjire devin semnificative. Practic, este imposibil să determinăm o valoare optimă pentru h atât timp cât nu avem informaţii referitoare la derivata de ordinul trei a funcţiei.
5.3 Utilizarea seriei Taylor
În acest subcapitol vom încerca să descriem o metodă de generare a formulelor de diferenţiere numerică utilizând dezvoltarea în serie Taylor. Ne aşteptăm ca prin această metodă să obţinem aceleaşi rezultate ca şi prin diferenţierea formulelor de interpolare.
81
Derivarea (diferenţierea) numerică
Pentru o derivată de ordin p, numărul minim de puncte necesar pentru o formulă de diferenţiere numerică, este p+1. De exemplu, pentru aproximarea primei derivate a unei funcţii sunt necesare cel puţin două puncte. Vom începe prin a determina o formulă de aproximare pentru f'i = f'(xi) în funcţie de fi = f(xi) şi fi+1 = f(xi+1). Funcţia f este reprezentată numeric în puncte. Dezvoltarea în serie Taylor pentru fi+1 în punctul xi este
...+f24h+"f
6h"+f
2h+hf+f = f IV
i
4
i
3
i
2
ii1+i ′′ (5.18)
Rezolvând ecuaţia (5.18) în raport cu fi' obţinem
-..."f61"-fh
21-
hf-f = f ii
i1+ii ′′ (5.19)
Dacă trunchiem ecuaţia (5.19) după primul termen obţinem formula de aproximare directă prin diferenţe, descrisă deja în introducere. De obicei se consideră că eroarea prin trunchiere este reprezentată de primul termen la care s-a renunţat (în acest caz -(h/2)fi") deoarece ceilalţi termeni descresc rapid comparativ cu primul atunci când h descreşte. Aproximarea directă prin diferenţe, incluzând şi efectul erorii prin trunchiere poate fi exprimată sub forma
f'i = fi+1 - fi
h + O(h) (5.20)
unde O(h) = -12 hfi".
Termenul O(h) exprimă faptul că eroarea prin trunchiere este aproximativ proporţională cu distanţa dintre nodurile adiacente, h. Aproximarea inversă prin diferenţe pentru prima derivată (în funcţie de fi-1 şi fi) se obţine similar. Dezvoltarea în serie Taylor pentru fi-1 este:
...-f24h+"f
6h"-f
2h+hf-f = f IV
i
4
i
3
i
2
ii1-i ′′ (5.21)
Rezolvând ecuaţia (5.21) în raport cu fi' obţinem
f'i = fi - fi-1
h + O(h) (5.22)
unde O(h) = -12 hfi".
Formula de aproximare centrată prin diferenţe a primei derivate fi' (utilizând fi-1 şi fi+1) poate fi dedusă plecând de la dezvoltarea în serie Taylor pentru fi+1 şi fi-1. Scăzând ecuaţia (5.21) din ecuaţia (5.18) şi rezolvând ecuaţia astfel obţinută în raport cu fi' obţinem
f'i = fi+1 - fi-1
2h + O(h2) (5.23)
unde O(h) = -16 h2fi".
Ceea ce trebuie remarcat este că în această formulă eroarea este proporţională cu h2 şi nu cu h, ceea ce este datorat dispariţiei termenului în f" atunci când se scad ecuaţiile (5.21) şi (5.18). În consecinţă, atunci când h descreşte, eroarea va descreşte mai rapid decât în primele două formule.
82
Utilizarea seriei Taylor
După cum am menţionat deja, o aproximare prin diferenţe pentru fi(p) necesită cel puţin p+1 puncte. Dacă se utilizează mai multe puncte va rezulta o formulă de aproximare prin diferenţe mai precisă. Pentru un număr dat de puncte, ecuaţia prin diferenţe care va da cea mai bună precizie este ecuaţia în care termenul care exprimă eroarea este de cel mai mare ordin posibil. Pentru a ilustra această afirmaţie, să deducem o formulă de aproximare pentru fi' în funcţie fi,fi+1 şi fi+2. Numărul minim de puncte necesar acestei formule fiind doi, înseamnă că avem la dispoziţie un punct în plus faţă de numărul minim de puncte. Dezvoltările pentru fi+1 şi fi+2 sunt:
...+f24h+"f
6h"+f
2h+hf+f = f IV
i
4
i
3
i
2
ii1+i ′′ (5.24)
...+f24h+16"f
6h"+8f
2h4+hf2+f = f IV
i
4
i
3
i
2
ii2+i ′′ (5.25)
Pentru a elimina termenul care conţine derivata a doua fi" scădem ecuaţia (5.23) din ecuaţia (5.22) înmulţită cu 4,
+..."fh32-hf2+f3 = f-f4 i
3ii2+i1+i ′′ (5.26)
astfel încât termenul care exprimă eroarea prin trunchiere să fie termenul care conţine derivata de ordinul trei. Rezolvând ecuaţia (5,24) în raport cu fi' rezultă:
)hO(+2h
f3-f4+f- = f 2i1+i2+ii′ (5.27)
unde O(h2) = 13 h2fi"'
Ecuaţia (5.25) este aproximarea directă în trei puncte pentru fi' iar eroarea este de acelaşi ordin ca la aproximarea centrată. Similar poate fi dedusă formula de aproximare inversă în trei puncte în funcţie de fi, fi-1 şi fi-2
)hO(+2h
f+f4-f3 = f 22-i1-iii′ (5.28)
unde O(h2) = 13 h2fi"'.Să analizăm în continuare modul de generare a formulelor de aproximare
pentru derivata a doua. Principiul de bază constă în eliminarea primei derivate din dezvoltarea în serie Taylor şi dacă este posibil, a cât mai mulţi termeni de ordin mai mare ca 2. De exemplu, să deducem formula de aproximare pentru fi" în funcţie de fi+1, fi şi fi-1. Dezvoltarea în serie Taylor pentru fi+1 şi fi-1 este dată de ecuaţiile (5.22) şi (5.21). Adunând cele două relaţii obţinem
...+fh121"+fh+f2 = f+f IV
i4
i2
i1-i1+i
Scăzând 2fi din ambii membri ai acestei ecuaţii rezultă
...+fh121"+fh = f+f2-f IV
i4
i2
1-ii1+i
Rescriind ecuaţia astfel obţinută sub forma
)hO(+h
f+f2-f = "f 22
1-ii1+ii (5.29)
83
Derivarea (diferenţierea) numerică
unde O(h2) = - 112 h2fiIV
am obţinut în acest mod formula de aproximare centrată prin diferenţe pentru fi", în care termenul O(h2) reprezintă eroarea prin trunchiere. Se poate deduce şi o altă formulă de aproximare pentru fi" în funcţie de fi,fi-1,fi-2. Scăzând de două ori ecuaţia (5.21) din ecuaţia (5.22) şi rezolvând ecuaţia astfel obţinută în raport cu fi" obţinem formula inversă de aproximare prin diferenţe pentru fi":
O(h)+h
f+f2-f = "f 2i1-i2-i
i (5.30)
unde O(h) = hfi"' Prin diverse combinaţii liniare între dezvoltările în serie Taylor se pot obţine formule de aproximare şi pentru derivatele de ordin mai mare decât 2. Demonstraţiile devin din ce în ce mai greoaie cu cât creşte numărul de puncte sau ordinul derivatei, existând riscul să apară greşeli în special la termenii care exprimă eroarea. Din acest motiv vom descrie în subcapitolul următor un algoritm şi un program de calcul pentru generare formulelor de aproximare a derivatelor de orice ordin. Până acum am presupus că funcţiile sunt reprezentate numeric în puncte egal depărtate. Pornind de la dezvoltarea în serie Taylor se pot deduce şi formulele de derivare în cazul punctelor arbitrare. Listăm în tabelul cele mai frecvent utilizate formule de aproximare pentru derivatele de ordinul 1,2 şi 3. Tabelul Derivata de ordinul 1 a) Formule directe de aproximare:
"fh21- = O(h) ,O(h)+
hf-f = f i
i1+ii′
"fh31=)hO( ,)hO(+
2hf3-f4+f- = f i
222i1+i2+ii ′′
fh41- = )hO( ,)hO(+
6hf11-f18+f9-f2 = f IV
i333i1+i2+i3+i
i′
b) Formule inverse de aproximare:
hf"21 = O(h) ,O(h)+
hf-f = f 1-ii
i′
"fh31 = )hO( ,)hO(+
2hf+f4-f3 = f 2222-i1-ii
i ′′
fh41 = )hO( ,)hO(+
6hf2-f9+f18-f11 = f IV3333-i2-i1-ii
i′
c) Formule centrate de aproximare:
"fh61- = )hO( ,)hO(+
hf-f = f i
2221-i1+ii ′′
84
fh31 = )hO( ,)hO(+
12hf+f8-f8+f- = f V
i4442-i1-i1+i2+i
i′
Derivata de ordinul 2 d) Formule directe de aproximare:
"hf- = O(h) ,O(h)+h
f+f2-f = "f i2i1+i2+i
i ′
fh1211 = )hO( ,)hO(+
hf2+f5-f4+f- = "f IV
i222
2i1+i2+i3+i
i
e) Formule inverse de aproximare:
"hf = O(h) ,O(h)+h
f+f2-f = "f i22-i1-ii
i ′
fh1211 = )hO( ,)hO(+
hf-f4+f5-f2 = "f IV
i222
23-i2-i1-ii
i
f) Formule centrate de aproximare:
fh121- = )hO( ,)hO(+
hf+f2-f = "f IV
i222
21-ii1+i
i
fh901 = )hO( ,)hO(+
h12f-f16+f30-f16+f- = "f IV
i444
22-i1-ii1+i2+i
i
Derivata de ordinul 3 g) Aproximarea directă:
fh23- = O(h) ,O(h)+
hf-f3+f3-f = "f IV
i3i1+i2+i3+i
i ′
h) Aproximarea inversă:
hf23 = O(h) ,O(h)+
hf-f3+f3-f = "f IV
i33-i2-i1-ii
i ′
i) Aproximarea centrată:
fh41- = )hO( ,)hO(+
h2f-f2+f2-f = "f V
i222
32-i1-i1+i2+i
i ′
5.5 Aproximarea derivatelor parţiale
Formulele pentru aproximarea derivatelor parţiale ale funcţiilor de mai multe variabile sunt în esenţă asemănătoare cu cele de diferenţiere pentru funcţiile de o singură variabilă. Să considerăm o funcţie de două variabile, f(x,y). Formula de aproximare pentru derivata
parţială fx = ∂∂x f(x,y) pentru x = x0 şi y = y0 poate fi dedusă fixând y = y0 şi considerând f(x,y0) ca
85
Derivarea (diferenţierea) numerică
o funcţie de o singură variabilă. Astfel, formulele directă, inversă şi centrată de aproximare pentru derivata parţială precedentă vor putea fi scrise sub forma:
x
)y,x(-f)yx,+x(f f 0000x Δ
Δ_ (5.31)
x
)yx,-xf(-)y,xf( f 0000x Δ
Δ_ (5.32)
x2
)yx,-xf(-)yx,+xf( f 0000x Δ
ΔΔ_ (5.33)
Formulele centrate de aproximare a derivatelor parţiale de ordinul doi ale funcţiei f(x,y) în x = x0 şi y = y0 sunt
x
)yx,-xf(+)y,x2f(-)yx,+xf( fx
= f 2000000
2
2
xx ΔΔΔ
∂∂ _ (5.34)
y
y)-y,xf(+)y,x2f(-y)+y,xf( fy
= f 2000000
2
2
yy ΔΔΔ
∂∂ _ (5.35)
yx4y)-yx,-xf(+y)+yx,-xf(-+
+yx4
y)-yx,+xf(-y)+yx,+xf( fy x
= f
0000
00002
yx,
ΔΔΔΔΔΔ
ΔΔΔΔΔΔ
∂∂∂ _
(5.36)
5.5.1 Exemplul Se dă funcţia de două variabile f(x,y) reprezentată numeric în tabelul următor:
y\x x = 1,0 1,5 2,0 2,5 3,0
1,0 1,63 2,05 2,50 2,98 3,49
1,5 1,98 2,51 3,08 3,69 4,33
2,0 2,28 2,91 3,61 4,37 5,17
2,5 2,64 3,25 4,08 5,00 5,98
3,0 2,65 3,50 4,48 5,57 6,76
a) Utilizând formulele centrate de aproximare să se evalueze următoarele derivate parţiale: x(2,2), fy(2,2), fyy(2,2), fxy(2,2) b) Utilizând formulele directe de aproximare în trei puncte să se evalueze următoarele derivate parţiale: fx(2,2), fy(2,2) Conform modului de definire a funcţiei Δx= Δy= 0,5. Derivatele parţiale cerute se vor calcula după cum urmează:
86
Aproximarea derivatelor parţiale
0,57 = ](0,5)[2
2,51+3,25-3,69-5,0 =
= )(2h
h)-h,2-f(2+h)+h,2+f(2-h)+h,2+f(2 (2,2)f
0,24- = )(0,5
3,08+(3,61)2-4,08 =
= h
h)-f(2,2+2f(2,2)-h)+f(2,2 (2,2)f
1,00 = (0,5)(2)3,08-4,08 =
2hh)-f(2,2-h)+f(2,2 (2,2)f
1,46 = (0,5)(2)2,91-4,37 =
2hh,2)-f(2-h,2)+f(2 (2,2)f a)
2
2yx,
2
2yy
y
x
•
•
•
•
_
_
_
_
1.01 = (0,5)2
(3,61)3-(4,08)4+(4,48)-=
= 2h
3f(2,2)-h)+4f(2,2+2h)+f(2,2- (2,2)f
1,48 = (0,5)2
(3,61)3-(4,37)4+(5,17)-=
2h
3f(2,2)-h,2)+4f(2+2h,2)+f(2- (2,2)f b)
y
x
•••
•••
_
_
Bibliografie:
87
Metode de rezolvare a ecuaţiilor neliniare
6. Metode de rezolvare a ecuaţiilor neliniare
6.1 Separarea rădăcinilor [i] pg 95
O ecuaţie de forma: f(x) = 0 (6.1)
se numeşte algebrică dacă funcţia f(x) este un polinom sau poate fi adusă la formă polinomială [i]. O ecuaţie este transcendentă dacă nu este algebrică. În general pentru ecuaţiile algebrice sau transcendente nu se pot aplica metode exacte de rezolvare. Presupunând că:
f:[xmin, xmax]→R, [xmin, xmax]⊂ R (6.2)
orice valoare ξ∈[xmin, xmax] pentru care f(ξ) = 0 se numeşte rădăcină a ecuaţiei (6.1) sau zero al funcţiei f(x).
În general, prin rădăcină aproximativă a ecuaţiei (6.1) se înţelege o valoare ξ' apropiată de valoarea exactă ξ. O rădăcină aproximativă poate fi definită în două moduri:
• numărul ξ' cu proprietatea că |ξ' - ξ|< ε (ε∈R, ε>0) sau,
• numărul ξ' cu proprietatea că |f(ξ')|< ε Aceste două moduri de definire a rădăcinii aproximative nu sunt însă echivalente. Determinarea rădăcinilor ecuaţiei f(x) = 0 implică, de obicei două etape: 1. separarea (izolarea) rădăcinilor,
adică stabilirea unei partiţii xmin = x1, x2, …, xm, …, xM = xmax, astfel încât orice subinterval [xm, xm+1] să conţină cel mult o rădăcină a ecuaţiei f(x)=0
Figura 6.1
2. calculul rădăcinilor cu precizia dorită, presupunând că s-a efectuat separarea lor. Una dintre metodele de separare a rădăcinilor se bazează pe următoarea consecinţă a teoremei lui Rolle: între două rădăcini reale consecutive ale ecuaţiei f'(x) = 0 există cel mult o rădăcină reală a ecuaţiei f(x) = 0.
88
7.2 Metoda bisecţiei (Metoda înjumătăţirii)
Fie x1 < x2 < … < xM, rădăcinile ecuaţiei f(x)=0. Se formează şirul lui Rolle: f(-∞), f(x1), f(x2) ,…, f(xm), …, f(xM),f(+∞). Conform consecinţei enunţate, în fiecare interval: (-∞, x1), (x1, x2), …,(xm, xm+1), …, (xM, +∞) se află o rădăcină reală a ecuaţiei f(x)=0 numai dacă funcţia ia valori de semne contrare la capetele intervalului, adică dacă f(xm)⋅ f(xm+1)<0. Deci ecuaţia f(x}=0 are atâtea rădăcini cîte variaţii de semn prezintă şirul lui Rolle. Din punct de vedere practic, această metodă prezintă un important dezavantaj, deoarece implică rezolvarea ecuaţiei f'(x)=0 ceea ce uneori este la fel de dificilă ca şi rezolvarea ecuaţiei f(x)=0. Foarte utilă pentru separarea rădăcinilor este următoarea teoremă. Teoremă. Dacă o funcţie continuă f(x) admite valori de semn opus la capetele unui interval [a,b], adică f(a)⋅f(b)<0, atunci acel interval conţine cel puţin o rădăcină a ecuaţiei f(x)=0, cu alte cuvinte, va exista cel puţin un număr ξ∈[a,b] , astfel încât f(ξ)=0 . ,
Rădăcina ξ este unică în intervalul [a,b] dacă derivata f'(x) există şi păstrează semnul în acest interval. În caz contrar, este posibil ca în intervalul [a,b] să existe mai multe rădăcini, ca în figura 7.2(a), şi, mai mult, acest lucru se poate întâmpla chiar dacă , f(a)f(b)>0, după cum se poate observa din figura 7.2(b). Procesul de separare a rădăcinilor implică determinarea semnelor funcţiei, f(x) în punctele unei partiţii {xm}, alegerea punctelor intermediare fiind dictată de particularităţile funcţiei. Dacă subintervalele [xm, xm+1] sunt suficient de mici, atunci, conform teoremei enunţate, fiecare subinterval va conţine cel mult o rădăcină a ecuaţiei f(x)=0. Astfel, în intervalele pentru care f(xm)⋅ f(xm+1)>0 , nu va exista nici o rădăcină, iar în intervalele pentru care f(xm)⋅ f(xm+1)≤ 0, va exista o singură rădăcină a ecuaţiei.
Figura 6.2
6.2 7.2 Metoda bisecţiei (Metoda înjumătăţirii)
Fie ecuaţia f(x) = 0 (6.3)
unde funcţia f(x) este continuă pe intervalul (a,b). Presupunem că în urma unui proces de separare (izolare) - a rădăcinilor ecuaţia (7.2) are cel mult o rădăcină în intervalul [a,b]. Dacă notăm cu ξ această rădăcină, sunt posibile următoarele cazuri: a) dacă f(a) =0 , atunci ξ = a b) dacă f(b) =0 , atunci ξ = b, c) dacă f(a)⋅f(b)<0, atunci ξ∈(a,b) d) dacă f(a)⋅f(b)>0, atunci ξ∉[a,b].
89
Metode de rezolvare a ecuaţiilor neliniare
Presupunând că se realizează cazul c) (tratarea celorlalte cazuri este banală); pentru a găsi rădăcina ξ∈(a,b), împărţim intervalul [a,b] în două părţi egale prin punctul x0=(a+b)/2, după cum se observă în figura 7.3. Dacă f(x0)=0. atunci ξ = x0 este rădăcina căutată. În caz contrar alegem semiintervalul [a1,b1] la capetele căruia funcţia are semne opuse, deoarece acesta conţine rădăcina ξ. Prin urmare, avem
a1 = a0, b1 = x0, dacă f(a)⋅f(x0)<0,
a1 = x0, b1 = b0, dacă f(a)⋅f(x0)>0, Noul interval [a1,b1] se înjumătăţeşte din nou, făcându-se aceleaşi teste asupra semnului funcţiei ca şi mai sus. Continuându-se acest procedeu, după un număr de i paşi (iteraţii) se obţine intervalul [ai,bi] ca urmare a înjumătăţirii intervalului [ai-1,bi-1] prin punctul:
xi-1 = ai-1 + bi-1
2 (6.4) Figura 6.3
astfel încât
ai = ai-1, bi = xi-1, dacă f(ai-1)⋅f(x i-1)<0, respectiv ai = xi-1, bi = b i-1, dacă f(a i-1)⋅f(xi-1)>0, Este evident că extremităţile noului subinterval satisfac relaţia:
f(ai)⋅f(bi)<0 (6.5)
iar lungimea sa este
bi - ai = b - a
2i (6.6)
Ca urmare a înjumătăţirii acestui interval prin punctul
xi = ai + bi
2 (6.7)
se obţine fie rădăcina exactă ξ = xi, fie un nou interval [ai+1, bi+1]. Deoarece extremităţile stângi ale intervalelor, a0, a1, …, ai formează un şir mărginit nedescrescător, iar extremităţile drepte, b0, b1, …, bi; formează un şir mărginit necrescător, rezultă prin trecere la limită în relaţia (7.4) că şirurile {ai} şi {bi} au o limită comună
ξ = lim ai i→∞
= lim bi i→∞
(6.8)
De asemenea, în virtutea continuităţii funcţiei f(x), inegalitatea (7.3) conduce prin trecere la limită la [(f(ξ)]2<0, de unde, f(ξ)=0, adică ξ este rădăcină a ecuaţiei (7.2). Evident, în practică se poate efectua doar un număr finit de iteraţii, deci se poate determina doar o valoare aproximativă a rădăcinii. Procesul de iterare se poate încheia, considerând ca rădăcină aproximativă ξ' = xi, când este îndeplinită una dintre condiţiile
|f(xi)|≤ ε
bi - ai ≤ ε (6.9)
În cazul în care xi ≠ 0, un criteriu de eroare mai util îl oferă eroarea relativă şi cea de-a doua condiţie poate fi rescrisă
90
7.3 Metoda secantei (corzii)
bi - ai
xi ≤ ε (6.10)
Prezentăm în continuare algoritmul bazat pe metoda bisecţiei. ε este eroarea absolută admisibilă utilizată în criteriile (7.6) şi (7.7) (nu se verifică anularea funcţiei sau egalitatea capetelor intervalului, deoarece din punct de vedere numeric aceste condiţii pot să nu fie satisfăcute niciodată), iar zero reprezintă o variabilă căreia i se atribuie o valoare de adevăr (adevărat sau fals) după cum există sau nu există rădăcină în intervalul [a, b]. Valorile funcţiei pentru extremităţile intervalului obţinut în urma înjumătăţirii anterioare, f(ai) şi f(bi), nu sunt reactualizate; deoarece prezintă importanţă numai semnul lui f(ai) şi acesta rămâne neschimbat pe parcursul procesului de înjumătăţire. Dacă este posibil (x ≠ 0) în locul erorii absolute se testează convergenţa procesului cu ajutorul erorii relative.
6.2.1 Algoritmul DATE DE INTRARE:
• a, b {capetele intervalului} • ε {eroarea maximă admisibilă}
DATE DE IEŞIRE: • x {rădăcina} • zero {variabilă logică de control}
Pasul 1 : CITEŞTE a, b, ε Pasul 2. zero = adevărat {există o rădăcină} Pasul 3. x = a; fa = f(x);
DACĂ |fa| ≤ ε ATUNCI EXIT {rădăcina este x=a} Pasul 4. x = b; fb = f(x);
DACĂ |fb| ≤ ε ATUNCI EXIT {rădăcina este x=b} Pasul 5.
DACĂ fa ⋅ fb < 0 ATUNCI {există rădăcină în intervalul (a,b)} REPETĂ
x = (a+b)/2; fx = f(x) {f((a+b)/2)} DACĂ fa ⋅fx < 0 ATUNCI b = x ALTFEL a = x dx = b-a; DACĂ x ≠0 ATUNCI dx = dx/x {eroarea}
PÎNĂ CÂND |dx|≤ ε sau |fx|≤ ε {testează convergenţa} ALTFEL zero = fals {nu există rădăcină}
Pasul 6. STOP
6.3 7.3 Metoda secantei (corzii)
Metoda secantei este utilizată, ca şi metoda bisecţiei, pentru găsirea rădăcinii ξ a ecuaţiei f(x)=0 izolată într-un interval [a,b], în cazul în care f(a)⋅f(b)<0. Principala deosebire faţă de metoda bisecţiei constă în faptul că subintervalele succesive [a1, b1], [a2, b2],…, [ai, bi] nu se obţin prin înjumătăţire, ci prin împărţirea intervalului anterior în raportul f(ai-1)/f(bi-1) . Geometric, după cum se poate observa în figura 7.4, metoda secantei este echivalentă cu înlocuirea funcţiei f(x) prin coarda care trece prin punctele (ai, f(ai)) şi (bi, f(bi)).
Figura 6.4
91
Metode de rezolvare a ecuaţiilor neliniare
Din ecuaţia corzii
xi - aibi - ai
= - f(ai)
f(bi) - f(ai) (6.11)
se poate obţine coordonata punctului de intersecţie xi al corzii cu axa absciselor
xi = aif(bi) - bif(ai)
f(bi) - f(ai) (6.12)
După un anumit număr de paşi se obţine, ca şi în metoda bisecţiei, fie o rădăcină exactă (în sens numeric) ξ = xi, astfel încât f(xi) = 0 , fie o secvenţă de intervale succesive [a0,b0], [a1, b1], …, [ai, bi], …cu
ai+1 = ai, bi+1= xi, dacă f(ai)⋅f(xi) < 0 sau
ai+1 = xi, bi+1= bi, dacă f(ai)⋅f(xi) > 0
astfel încât f(ai+1)⋅f(b i+1)<0 Deşi în general mai eficientă decât metoda bisecţiei, metoda secantei este avantajoasă tot numai pentru determinarea grosieră a rădăcinilor reale ale ecuaţiilor algebrice şi transcendente, furnizând aproximaţii iniţiale pentru rădăcini, necesare altor metode iterative.
6.4 7.4 Metoda aproximaţiilor succesive
Una dintre cele mai importante metode numerice pentru rezolvarea ecuaţiilor algebrice şi transcendente este metoda aproximaţiilor succesive. Această metodă poate fi utilizată pentru rafinarea aproximaţiilor iniţiale ale rădăcinilor, furnizate de alte metode (cum ar fi, de exemplu, metoda înjumătăţirii sau metoda secantei). Presupunem că se cere rezolvarea ecuaţiei: f(x) = 0 (6.13)
unde f(x) este o funcţie continua pe un interval [a,b]. Metoda aproximaţiilor succesive implică punerea acestei ecuaţii sub forma echivalentă
x = ϕ(x) (6.14)
astfel încât ecuaţiile (7.8) şi (7.9) să aibă aceleaşi rădăcini. Punerea ecuaţiei (7.8) sub forma (7.9) este întotdeauna posibilă, chiar dacă funcţia f(x) nu conţine în mod explicit termenul x (eventual; se poate aduna acest termen în ambii membri ai ecuaţiei (7.8)).,
Fie x0 o aproximaţie iniţială pentru rădăcina ξ a ecuaţiei (7.9). Înlocuind această valoare în membrul drept al ecuaţiei, se obţine o aproximaţie îmbunătăţită a rădăcinii
x1 = ϕ(x0) (6.15)
Înlocuind x1 în membrul drept al ecuaţiei (7.9) se obţine o nouă aproximaţie
x2 = ϕ(x1) (6.16)
Repetând acest procedeu, se obţine pornind de la x0 secvenţa de numere
xi+1 = ϕ(xi), i=0,1,2,… (6.17)
Dacă acest şir este convergent, atunci există limita ξ = lim xi. Trecând la limită în (7.10) şi admiţând că ϕ(x) este continuă, găsim
lim xi+1 i→∞
= ϕ(lim xi i→∞
) (6.18)
92
7.4 Metoda aproximaţiilor succesive
ξ = ϕ(ξ) (6.19)
şi deci ξ este rădăcină a ecuaţiei (7.9) şi poate fi calculată pe baza relaţiei de recurenţă (7.10) cu orice precizie.
Din punct de vedere geometric, o rădăcină reală ξ a ecuaţiei x = ϕ(x) este abscisa punctului de intersecţie a curbei y = ϕ(x) cu dreapta y = x (prima bisectoare). Relativ la valoarea derivatei ϕ'(x), sunt posibile patru cazuri, iar modul cum rezultă şirul aproximaţiilor succesive x1, x2, x3, ... plecând de la aproximaţia iniţială x0, este ilustrat în figurile 7.5 şi 7.6. Astfel, în figurile 7.5 a şi b procesul iterativ este convergent, iar în figurile 7.6 a şi b, acest proces este divergent. Din studiul acestor figuri se poate trage concluzia că procesul este convergent, şi deci metoda este aplicabilă numai în intervalele unde
|ϕ'(x)|<1 (6.20)
Pentru a putea aplica practic metoda aproximaţiilor succesive trebuie îndeplinite condiţiile prescrise de următoarea teoremă, care furnizează o condiţie necesară şi suficientă de convergenţă a şirului de iterare definit de relaţia (7.10).
Fig. 7.6
Teoremă. Fie ecuaţia
x = ϕ(x) (6.21)
cu funcţia ϕ(x) definită şi derivabilă pe [a,b] . Dacă este satisfăcută inegalitatea
|ϕ'(x)|≤ λ <1 (6.22)
pentru orice x∈[a,b] , atunci şirul de iterare definit de relaţia
xi+1 = ϕ(xi), i=0,1,2,.. (6.23)
converge către rădăcina (unică dacă există) ξ∈[a,b] a ecuaţiei, indiferent de valoarea iniţială x0. Demonstraţie. Scădem egalitatea evidentă ξ = ϕ(ξ) din ecuaţia (7.14)
xi+1 - ξ = ϕ(xi) - ϕ(ξ) (6.24)
Din teorema lui Lagrange rezultă că există un punct ζ cuprins între xi şi ξ astfel încât
ϕ(xi) - ϕ(ξ) = ϕ'(ζ)(xi - ξ) (6.25)
Având în vedere inegalitatea (7.13) rezultă din relaţiile (7.15) şi (7.16)
|xi+1 - ξ|≤ λ|xi - ξ| (6.26)
Dând pe rând lui i valorile 0,1,2,... obţinem succesiv
|x1 - ξ|≤ λ|x0 - ξ|
|x2 - ξ|≤ λ|x1 - ξ|≤ λ2|x0 - ξ|
……………
|xi+1 - ξ|≤ λ|xi - ξ|≤ …≤ λi+1|x0 - ξ| (6.27)
93
Metode de rezolvare a ecuaţiilor neliniare
Trecând la limită în ultima inegalitate (0 ≤ λ < 1 ) rezultă
lim xi+1 i→∞
= ξ (6.28)
şi deci limita şirului de iterare este chiar rădăcina ξ a ecuaţiei.
Pentru a arăta că în condiţiile teoremei ξ este unica rădăcină din intervalul [a,b], presupunem contrariul, adică faptul că ar mai exista o rădăcină ξ∈[a,b], astfel încât ξ' = ϕ(ξ)'. Atunci
ξ - ξ' = ϕ(ξ) - ϕ(ξ') = ϕ'(ζ)(ξ - ξ') (6.29)
unde ζ se află între ξ şi ξ' . Mai putem scrie
(ξ - ξ')[1-ϕ'(ξ)] = 0 (6.30)
Cum prin ipoteza [1-ϕ'(ξ)] ≠ 0, rezultă ξ = ξ', adică rădăcina ξ este unică. Din punct de vedere practic, ecuaţia (7.8) poate fi pusă în mai multe moduri sub forma (7.9). Oricum, pentru ca metoda aproximaţiilor succesive să fie aplicabilă, trebuie să fie îndeplinită condiţia (7.13). Cu cât va fi mai mic număruI λ, cu atât va fi mai rapidă convergenţa procesului iterativ către rădăcina ξ. În general, ecuaţia (7.8) poate fi înlocuită prin ecuaţia echivalentă x = x - qf(x), q>0 (6.31)
astfel încât
ϕ(x) = x - qf(x), (6.32)
Parametrul q trebuie ales în aşa fel încât să fie satisfăcută condiţia,(7.13), adică
|ϕ'(x)| =|1 - qf'(x)|≤ λ <1 (6.33)
În cazul banal în care se poate alege q = 1 , ecuaţia echivalentă (7.17) devine x = x - f(x)
În aceste condiţii, relaţia de recurenţă (7.10) a metodei aproximaţiilor succesive devine xi+1 = xi - f(xi), i=0,1,2,.. (6.34)
Practic, se iterează pe baza acestei relaţii până când eroarea absolută corespunzătoare la două iteraţii consecutive este mai mică sau egală cu o anumită eroare maximă admisibilă ε,
Δi = |xi+1 - xi|≤ ε, i=0,1,2, …, imax (6.35)
ceea ce este echivalent cu
Δi = |f(xi)|≤ ε, i=0,1,2, …, imax (6.36)
sau până când numărul de iteraţii depăşeşte o anumită valoare limită imax. Dacă xi≠0 se poate apela în locul erori absolute la criteriul mai util al erorii relative
δi = |1 - xi+1 / xi|≤ ε, i=0,1,2, …, imax (6.37)
Limitarea numărului de iteraţii permite detectarea situaţiilor în care metoda este divergentă.
6.4.1 Algoritmul Date de intrare:
• x {aproximaţie iniţială pentru rădăcină} • ε {eroarea maximă admisibilă} • imax {număr maxim de iteraţii}
Date de ieşire:
94
7.5 Metoda lui Newton (tangentei)
• x {rădăcina} • zero {variabilă logică de control}
Pasul 1 : CITEŞTE x, ε, imax Pasul 2. zero = adevărat {se admite că există o rădăcină} Pasul 3. i = 0 {iniţializează contorul iteraţiilor} Pasul 4.
REPETĂ i = i+1 {incrementează contorul} dx = f(x) {corecţia} x = x - dx {noua aproximaţie pentru rădăcină} DACĂ x ≠0 ATUNCI dx = dx/x {eroarea relativă}
PÎNĂ CÂND |dx| < ε SAU i > imax {testează convergenţa} DACĂ i ≥ imax ATUNCI zero = fals {nu există rădăcină}
Pasul 5. STOP
6.5 7.5 Metoda lui Newton (tangentei)
Presupunem că ecuaţia algebrică sau transcendentă f(x) = 0 (6.38)
are o rădăcină reală ξ izolată în intervalul [a,b] şi că f'(x) şi f"(x) sunt continue şi păstrează semnul pentru x∈ [a b] .
Având o aproximaţie iniţială x0 a rădăcinii ξ, metoda lui Newton-permite găsirea unei aproximaţii îmbunătăţite x1 , ca intersecţie a tangentei dusă la curba y = f(x) în punctul (x0, f(x0)) cu axa Ox, aşa după cum se poate observa din figura 7.7a. Avem în acest caz, alegând x0 =b,
Figura 6.5
f'(x) = f(x0)
x0 - x1 (6.39)
de unde
x1 = x0 - f(x0) f'(x0) (6.40)
Ducând acum tangenta în punctul de coordonate (x1, f(x1)) , obţinem o nouă aproximaţie x2 etc. Repetând procedeul acesta pe baza relaţiei ,
xi+1 = xi - f(xi) f'(xi) (6.41)
se găseşte un şir de aproximaţii x1, x2, x3, ... ale soluţiei ξ a ecuaţiei (7.21). Condiţiile convergenţei acestui şir către soluţia ξ sînt specificate de teorema demonstrată în cazul
95
Metode de rezolvare a ecuaţiilor neliniare
metodei aproximaţiilor succesive. Pentru aceasta introducem funcţia de iterare a metodei lui Newton
ϕ(x) = x - f(x) f'(x) (6.42)
relaţia de iterare (7.22) putând fi pusă sub forma
xi+1 = ϕ(xi) , i=0,1,2,… (6.43)
Se remarcă faptul că în cazul considerat în figura 7.7a f(x0)⋅f"(x0)>0. Dacă se consideră ca aproximaţie iniţială x0 = a, ca în figura 7.7b, şi prin urmare , f(x0)⋅f"(x0) <0, se obţine aproximaţia de ordinul întâi x1' situată în afara intervalului [a;b]. Această aproximaţie iniţială este mai puţin-avantajoasă, în sensul că sunt necesare mai multe iteraţii pentru a ajunge la precizia dorită. Se poate demonstra, ca o regulă generală, că cea mai avantajoasă satisface condiţia
f(x0)⋅f"(x0) > 0 (6.44)
Practic, în cazul metodei lui Newton se iterează, ca şi în cazul metodei aproximaţiilor succesive, pînă când eroarea absolută corespunzătoare unor aproximaţii consecutive x0 şi x1 devine mai mică sau egală cu o eroare maximă admisibilă ε,
Δi = |xi+1 - xi|≤ ε, i=0,1,2, …, imax (6.45)
ceea ce este echivalent cu
Δi = |f(xi)/f'(xi)|≤ ε, i=0,1,2, …, imax (6.46)
sau până când numărul de iteraţii depăşeşte o anumită valoare maximă imax. Dacă xi ≠ 0 este recomandabil să se utilizeze eroarea relativă în testarea convergenţei în locul erorii absolute. Cu mici deosebiri, referitoare exclusiv la modul de calcul al corecţiei rădăcinii algoritmul metodei lui Newton este identic cu cel al metodei aproximaţiilor succesive.
6.5.1 Algoritmul DATE DE INTRARE:
• x {aproximaţie iniţială pentru rădăcină} • ε {eroarea maximă admisibilă} • imax {număr maxim de iteraţii}
DATE DE IEŞIRE: • x {rădăcina} • zero {variabilă logică de control}
Pasul 1. CITEŞTE x, ε, imax Pasul 2. zero = adevărat {se admite că există o rădăcină} Pasul 3. i = 0 {iniţializează contorul iteraţiilor} Pasul 4.
REPETĂ i = i+1 {incrementează contorul} dx = f(x); DACĂ f'(x) ≠ 0 ATUNCI dx = dx/f'(x) {corecţia} x = x - dx {noua aproximaţie pentru rădăcină} DACĂ x ≠0 ATUNCI dx = dx/x {eroarea relativă}
PÎNĂ CÂND |dx|≤ ε SAU i ≥ imax {testează convergenţa} DACĂ i ≥ imax ATUNCI zero = fals {nu există rădăcină}
Pasul 5. STOP Metoda lui Newton este mai rapid convergentă decât metoda aproximaţiilor succesive. Astfel, în cazul ecuaţiei
96
7.5 Metoda lui Newton (tangentei)
x - sin x - 0.25 = 0 (6.47)
pornind de la aproximaţia iniţială x0 = 0, metodei lui Newton îi sunt necesare 10 iteraţii pentru a obţine soluţia x = 1,17123, în timp ce metodei aproximaţiilor succesive îi sunt necesare 28 de iteraţii. Datorită faptului că metoda lui Newton face necesar calculul derivatei funcţiei, în cazul în care acest lucru este dificil sau chiar imposibil, se recurge frecvent la metoda aproximaţiilor succesive.
Ecuaţiile diferenţiale apar în diverse domenii ale matematicii aplicate, frecvent sub forma unor probleme de tip Cauchy sau Dirichlet. În acest capitol vom trata câteva metode numerice elementare pentru rezolvarea numerică a problemelor cu condiţii iniţiale pentru ecuaţii diferenţiale ordinare (probleme Cauchy). Fie sistemul de ecuaţii diferenţiale ordinare de ordinul întâi y'(x) = fi(x,y1(x),y2(x),…,yn(x)), i = 1,2,...,n (7.1)
unde x ∈ [a,b], iar fi, sunt funcţii cunoscute. Se poate arăta, de asemenea, că orice ecuaţie diferenţială de ordinul n: y(n) = f(x,y,y', …, y(n-1)) poate fi pusă sub forma echivalentă a unui sistem de n ecuaţii diferenţiale de ordinul întâi, considerând ca funcţii necunoscute alături de y şi derivatele yi = y(i-1), i = 1,2,...,n. Sistemul (12.1) poate fi scris mai simplu sub formă matricială y'(x) = f(x,y) (7.2)
unde
y = ⎣⎢⎡
⎦⎥⎤y1(x)
y2(x)…
yn(x)
, f(x,y) = ⎣⎢⎡
⎦⎥⎤f1(x y1(x) y2(x) … yn(x))
f2(x y1(x) y2(x) … yn(x)) …
fn(x y1(x) y2(x) … yn(x))
(7.3)
Ataşându-i acestui sistem condiţiile iniţiale în punctul x0 yi(x0) = yi0, i = 1,2,...,n (7.4)
sau, matricial, y(x0) = y0 (7.5)
obţinem o problemă cu condiţii iniţiale sau o problemă Cauchy. Orice metodă numerică pentru rezolvarea problemei Cauchy (12.1)-(12.4) sau (12.2)-(12.5) se bazează pe partiţionarea intervalului [a,b] prin intermediul unei diviziuni a = x0,x1,x2,…,xM = b şi urmăreşte calcularea progresivă a soluţiei problemei pentru fiecare punct al reţelei. Metodele numerice destinate rezolvării problemei Cauchy pot fi împărţite, în esenţă, în două mari categorii. Din prima categorie, cea a metodelor directe (sau pas cu pas, sau unipas), fac parte algoritmii care furnizează soluţia la un pas xm numai pe baza informaţiilor din nodul imediat anterior, xm-1, şi, eventual, din intervalul [xm-1,xm]. Dintre aceste metode amintim metoda liniilor poligonale a lui Euler şi metodele Runge-Kutta. Din cea de a doua categorie, cea a metodelor indirecte (sau multipas), fac parte algoritmi care permit calculul soluţiei ecuaţiei la un pas xm pe baza informaţiilor de la mai mulţi paşi anteriori … xm-2,xm-1. Dintre aceste metode amintim
98
Metoda dezvoltării în serie Taylor
metodele Adams-Moulton, Milne, Fox-Goodwin şi Numerov. Fiecare dintre cele două categorii de metode prezintă avantaje şi dezavantaje specifice.
7.2 Metoda dezvoltării în serie Taylor
Una dintre cele mai vechi metode de rezolvare a ecuaţiilor diferenţiale se bazează pe dezvoltarea soluţiei în serie Taylor. Această metodă permite, teoretic, găsirea soluţiei oricărei ecuaţii diferenţiale, dar aplicarea ei practică se dovedeşte de multe ori ineficientă. Importanţa acestei metode constă, mai degrabă, în faptul că ea furnizează criterii pentru evaluarea preciziei metodelor de interes practic. Să presupunem că se cere găsirea soluţiei ecuaţiei diferenţiale y'(x) = f(x,y) (7.6)
cu condiţia iniţială y(x0) = y0 (7.7)
Admiţând că soluţia este dezvoltabilă în serie Taylor în vecinătatea punctului x0, avem
y(x) = y(x0) + x - x0
1! y'(x0) + (x - x0)2
2! y"(x0) + … (7.8)
şi, în particular, valoarea soluţiei pentru x = x1, este:
y1 = y0 + hy0' + h2
2 y0" + … (7.9)
unde h = x1 - x0 şi y0(i) = y(i)(x0), j = 0,1,2,…
Presupunând acum că se cere determinarea soluţiei pe reţeaua de puncte x0,x1,x2,…,xm,xm+1,… formula (12.9) oferă un procedeu de calcul a acesteia din aproape în aproape, ca o succesiune de soluţii ale unor probleme Cauchy. Astfel, propagarea soluţiei de la punctul de reţea xm la punctul xm+1 este echivalentă cu rezolvarea ecuaţiei diferenţiale (12.6) cu condiţia iniţială y(xm) = ym (7.10)
şi se realizează pe baza unei formule analoage relaţiei (12.9)
ym+1 = ym + hym' + h2
2 ym" + … (7.11)
cu h = xm+1 - xm şi ym
(j) = y(j)(xm), j = 0,1,2,… (7.12)
Cu cât pasul h este mai mic (decât raza de convergenţă a seriei) şi cu cât este mai mare numărul termenilor luaţi în considerare în seria (12.11), cu atât valoarea ym+1 este mai apropiată de valoarea exactă y(xm+1) a soluţiei. Se observă că pentru calculul valorii ym+1 este necesară determinarea valorilor y'm, y'm, ... ale derivatelor soluţiei în punctul xm. Din ecuaţia diferenţială (12.6) avem y'm = f(xm,ym) (7.13)
Derivând acum (12.6) în raport cu x găsim
(7.14)
99
12. INTEGRAREA ECUAŢIILOR DIFERENŢIALE ORDINARE
înlocuind (12.13) şi (12.14) în (12.11), obţinem
(7.15) unde notaţia O(h3) indică faptul că restul acestei formule este de ordinul h3, adică toţi termenii neglijaţi conţin h la puterea a treia sau mai mare. Altfel spus, utilizarea formulei (12.15) fără termenul O(h3) prezintă o eroare de trunchiere de forma Kh3, unde K este o constantă. Metoda dezvoltării în serie Taylor este considerată ca fiind o metodă directă (sau pas cu pas), deoarece pentru determinarea lui yM+1 sunt necesare numai informaţii de la punctul anterior, adică xm. Dacă este necesară o aproximare superioară a soluţiei în nodurile reţelei, trebuie luaţi în considerare mai mulţi termeni în seria (12.11). Derivatele de ordin superior implicate devin din ce în ce mai complicate, ducând la formule ineficiente pentru propagarea soluţiei. În unele situaţii, evaluarea derivatelor soluţiei este chiar imposibilă. În ciuda inconvenientelor de ordin practic pe care le prezintă, metoda dezvoltării în serie Taylor furnizează totuşi criterii de evaluare a preciziei altor metode destinate rezolvării numerice a problemelor Cauchy. Aceste criterii se referă la gradul în care soluţia obţinută cu ajutorul metodelor numerice coincide cu dezvoltarea în serie Taylor a soluţiei exacte, grad exprimabil ca puterea lui h din ultimul termen luat în considerare în dezvoltarea cu care este echivalentă soluţia aproximativă. Acest mod de apreciere este aplicabil chiar dacă metodele de interes practic nu implică calculul explicit al vreunei derivate a funcţiei f(x,y).
7.3 Metoda lui Euler (metoda liniilor poligonale)
Această metodă reprezintă cel mai simplu caz particular al metodei dezvoltării în serie Taylor, anume cazul în care se reţin primii doi termeni, obţinându-se aproximaţia liniară
ym+1 = ym + hf(xm,ym), m = 0,1,2,… (7.16)
cu h = xm+1 - xm. Pe baza formulei de recurenţă (12.16) se poate stabili, plecând de la valoarea iniţială y0 = y(x0), şirul aproximaţiilor y1,y2,... ale valorilor adevărate y(x1),y(x2),... ale soluţiei. Aşa cum se observă din Figura 7.1, din punct de vedere geometric metoda lui Euler revine la înlocuirea curbei y = y(x) cu o linie poligonală, (x0,y0), (x1,y1),… (x2,y2),… de unde şi numele alternativ al metodei, segmentul iniţial reprezentând tangenta la curbă în punctul x0.
Figura 7.1
Se remarcă faptul că metoda lui Euler este o metodă directă (pas cu pas) ca şi metoda dezvoltării în serie Taylor, iar algoritmul este de tip explicit. Eroarea comisă în calculul soluţiei pe baza formulei (12.16) se compune din eroarea de trunchiere, de ordinul h2, şi din eroarea de rotunjire. Aceste erori se propagă de la un pas la altul, mărimea erorii totale obţinută prin acumulare depinzând de numărul paşilor de integrare efectuaţi. Precizia metodei lui Euler fiind destul de mică, pentru obţinerea unor rezultate satisfăcătoare, se impune utilizarea unui pas de integrare h foarte mic.
100
12.4 Metodele Runge-Kutta
Pentru obţinerea unei precizii sporite, se poate utiliza algoritmul predictor-corector al lui Euler. În cadrul acestui algoritm se determină în locul aproximaţiei ym+1 a soluţiei în punctul de reţea xm+1 pe baza formulei (12.16), două aproximaţii. Prima dintre aceste aproximaţii, numită valoare de predicţie, este dată de formula ym+1
p = ym + hf(xm,ym), m = 0,1,2,… (7.17)
numită formulă predictoare. Valoarea obţinută din această relaţie este utilizată apoi pentru obţinerea unei aproximaţii corectate cu ajutorul formulei
ym+1c = ym +
h2 [f(xm,ym) + f(xm+1,ym+1
p)] (7.18)
care se numeşte formulă corectoare. Algoritmul predictor-corector al lui Euler implică la fiecare pas de integrare următoarele etape de calcul:
1. Se calculează ym+1p din formula predictor (12.17).
2. Se calculează recurent ym+1c din formula corector (12.18), înlocuind de fiecare dată
valoarea "predictor", ym+1p, cu valoarea "corector" anterioară, ym+1
c, până când se obţine precizia dorită, potrivit criteriului
|ym+1c - ym+1
p| ≤ ε, ε>0
7.4 12.4 Metodele Runge-Kutta
Metodele Runge-Kutta constituie o clasă largă de metode numerice destinate rezolvării problemelor de tip Cauchy. Diferitele metode din această categorie implică un volum mai mic sau mai mare de calcule în funcţie de precizia pe care trebuie să o asigure în rezultate. Metodele Runge-Kutta au următoarele proprietăţi definitorii:
• Sunt metode directe (pas cu pas).
• Soluţia furnizată de metoda Runge-Kutta de ordinul p coincide cu dezvoltarea în serie Taylor până la termeni de ordinul hp.
• Necesită numai evaluarea funcţiei f(x,y), nu şi a derivatelor sale. Fie ecuaţia diferenţială y' = f(x,y). Propagarea soluţiei acestei ecuaţii de la punctul de reţea xm la punctul xm+1 se poate realiza, aşa după cum s-a arătat în paragraful 12.2, pe baza formulei
(7.19) Reţinerea unui număr mare de termeni în această serie este imposibilă datorită creşterii rapide a complexităţii derivatelor implicate în coeficienţi. Metodele Runge-Kutta îşi propun construirea unei funcţii de h, care să depindă de f(x,y), dar nu şi de derivatele sale şi pentru care dezvoltarea în serie după puterile lui h să coincidă cu dezvoltarea (12.19) a soluţiei exacte până la termeni de ordin cât mai mare, indiferent de funcţia f{x,y). Căutăm ym+1 sub forma unei combinaţii liniare
unde (7.20)
101
12. INTEGRAREA ECUAŢIILOR DIFERENŢIALE ORDINARE
(7.21)
αi, βij şi wi fiind constante care trebuie determinate. În toate cazurile se consideră
Din relaţiile de mai sus reiese caracterul algoritmic al metodelor Runge-Kutta: diversele mărimi se calculează numai cu ajutorul mărimilor deja cunoscute. Vom determina în continuare constantele αi, βij şi wi astfel încât dezvoltarea în serie după puterile lui h a formulei (12.20) să coincidă până la o anumită putere a lui h cu dezvoltarea în serie (12.19), indiferent de funcţia f(x,y). Dând diverse valori lui p în (12.20), se obţin diverse formule Runge-Kutta. Pentru p = 1 rezultă ym+1 = ym + hf(xm,ym) (7.23)
şi se regăseşte metoda lui Euler. Pentru p = 2, din formulele (12.20) şi (12.22) rezultă ym+1 = ym + w1k1 + w2k2 (7.24)
Dezvoltând această expresie în serie după puterile lui h se obţine
Comparând această formulă cu (12.19) şi identificând coeficienţii aceloraşi puteri ale lui h, găsim
⎩⎪⎨⎪⎧
w1 + w2 = 1
w2α2 = 12
w2β21 = 12
Acesta este un sistem de trei ecuaţii cu patru necunoscute, care nu are soluţie unică. Se observă că trebuie să avem w2 ≠ 0. Alegând w2 = 1/2, obţinem α2 =β21 = 1; w1 = w2 = ½ şi, cu aceasta, formulele (12.24) şi (12.25) devin
102
12.4 Metodele Runge-Kutta
(7.26) (7.27) Regăsim astfel algoritmul predictor-corector al lui Euler, cu formula predictor înlocuită în formula corector. Varianta cea mai importantă a metodelor Runge-Kutta se obţine pentru p = 4. Urmând raţionamentul pentru determinarea constantelor αi, βij şi wi din cazul p = 2, se ajunge la un sistem de 10 ecuaţii cu 13 necunoscute, care nu are soluţie unică. S-au determinat diferite sisteme de constante care verifică acest sistem, căutând să se obţină constante cu structură simplă pentru uşurarea calculelor (anumite sisteme de constante permit chiar şi reducerea necesarului de memorie). Alegerea cea mai răspândită şi care duce la cea mai utilizată metodă Runge-Kutta, reprezentând un optim simplitate-precizie şi permiţând o codificare robustă şi eficientă, este definită de relaţiile
ym+1 = ym + 16(k1 + 2k2 + 3k3 + k4), m = 0,1,2,… (7.28)
unde
(7.29)
Eroarea de trunchiere pentru metoda Runge-Kutta de ordinul patru este de forma εT = Kh5 şi remarcăm faptul că funcţia f(x,y) este evaluată de patru ori pentru fiecare pas de integrare. Fără o măsură a erorii de trunchiere, este dificilă alegerea pasului de integrare potrivit. O regulă utilă
afirmă că dacă expresia ⎪⎪⎪
⎪⎪⎪k2 - k3
k1 - k2 devine mai mică decât câteva sutimi, pasul h trebuie micşorat.
Principalul avantaj al metodei Runge-Kutta îl constituie posibilitatea autopornirii datorită faptului că este o metodă de tip pas cu pas. Precizia destul de bună a metodei pentru p = 4, alături de o relativă simplitate a implementării, face ca metoda Runge-Kutta să fie utilizată pe scară largă, mai ales pentru determinarea valorilor iniţiale necesare metodelor multipas, care, sub rezerva faptului că nu se pot auto-porni, sunt în general mai eficiente, oferind precizii mai mari în obţinerea soluţiei. Considerând că în relaţiile (12.28) şi (12.29) notaţiile sunt matriciale, putem utilii aceste relaţii pentru rezolvarea problemelor Cauchy pentru sisteme de ecuaţii diferenţiale de ordinul întâi. Prezentăm în continuare, ca exemplu de implementare în limbajul Pascal procedura RK4, destinată rezolvării acestui tip de probleme. Procedura RK4 are rolul de a rezolva o problemă Cauchy "elementară", propagând soluţia între două puncte de reţea vecine.
procedure RK4(x,h: real; var y0,y: vector; n: integer); var i : integer; hh : real; f1,f2,f3,yy : vector; begin
103
12. INTEGRAREA ECUAŢIILOR DIFERENŢIALE ORDINARE
hh := 0.5*h; func(x,y0,f1); for i := 1 to n do yy[i] := y0[i]+hh*f1[i]; func(x+hh,yy,f2); for i := 1 to n do yy[i] := y0[i]+hh*f2[i]; func(x+hh,yy,f3); for i := 1 to n do begin yy[i] := y0[i]+h*f3[i]; f2[i] := f2[i]+f3[ij; end; func(x+h,yy,f3); hh := h/6.0; for i := 1 to n do y[i] := y0[i]+hh*(f1[i]+2*f2[i]+f3[i]); end;
rocedura RK4 are ca date de intrare x (xm în formularea matematică), h (pasul de integrare),
1 2 3 4
onstruire a procedurii func, destinată generării membrilor
'(x0) = y0'Această problemă pentru o ecuaţ ul doi poate fi pusă sub forma
⎩ 0
Procedura func corespunzătoare va a
Py0 (tabloul condiţiilor iniţiale ale sistemului de ecuaţii diferenţiale) şi r (numărul de ecuaţii diferenţiale din care se compune sistemul) şi are ca scop propagare.' soluţiei pe distanţa h. Componentele soluţiei în punctul xm+1 = xm + h sunt depuse în tablou y. Utilizatorul trebuie să furnizeze o procedură numită func, având trei parametri: x, y (având semnificaţiile de mai sus) şi f, care este un tablou în care se returnează vectorul valorilor funcţiilor f(x,y) din membrii drepţi ai ecuaţiilor diferenţiale care alcătuiesc sistemul. După cum se poate observa, în loc să se calculeze componentele vectorilor k , k , k şi k din motive de eficienţă se utilizează tablourile f1, f2 şi f3 pentru memorarea valorilor funcţiilor din componenţa acestor coeficienţi. De fapt, în final, tabloul f2 conţine suma componentelor vectorilor k2 şi k3, iar tabloul f3 conţine componentei vectorului k4. În tabloul yy sunt generate argumente pentru procedura func.
Pentru exemplificarea modului de cdrepţi ai ecuaţiilor diferenţiale care alcătuiesc sistemul, considerăm problema, Cauchy
⎨⎧y" + y = 0 ⎩y(x0) = y0; yie diferenţială de ordin
echivalentă a unui sistem de două ecuaţii diferenţiale de ordinul întâi
⎨⎧y1' = y2; y1(x0) = y0
' y2' = -y1; y2(x0) = yvea forma:
procedure func(x:real; var y,f:vector); begin f[1] := y[2]; f[2] := -y[1]; end;
[i] Beu, T. - Analiză numerică în Turbo Pascal, Ed. Microinformatica, Cluj-Napoca, 1992
