Universitatea “Dunărea de Jos”

METODE NUMERICE Gabriel FRUMUŞANU

Galaţi - 2008

Departamentul pentru Învăţământ la Distanţă şi cu Frecvenţă Redusă Facultatea de Mecanica Specializarea Inginerie economica si industriala Anul de studii / Forma de învăţământ II/IFR

REFERENŢI ŞTIINŢIFICI: Prof. dr. ing. Nicolae OANCEA, Prof. dr. ing. Cătălina MAIER,

Universitatea „Dunărea de Jos” din Galaţi.

3

PREFAŢĂ

Apariţia şi perfecţionarea continuă a calculatoarelor electronice a

adus în prim-plan un domeniu al matematicii care exista de multă vreme şi care, iniţial, era prezentat ca o anexă a matematicii clasice: metodele numerice. În prezent acestea s-au constituit într-un domeniu de sine stătător, la graniţa cu programarea calculatoarelor şi cu ingineria.

În condiţiile în care viteza de lucru a calculatoarelor a atins limite de neimaginat, chiar în urmă cu un deceniu, calculul numeric oferă soluţii pentru orice problemă matematică a cărei rezolvare e necesară în tehnică; mai mult, soluţiile „aproximative” obţinute cu ajutorul metodelor numerice tind să se transforme, pe zi ce trece, în soluţii a căror precizie o depăşeşte cu mult pe ce a sistemelor tehnice.

Această carte se adresează studenţilor de la specializarea „Inginerie şi management”, forma de studiu cu frecvenţă redusă şi a fost, prin urmare, concepută modular, într-un număr de capitole egal cu numărul de cursuri din programă; fiecare capitol include noţiuni teoretice, de regulă şi câte un algoritm, exemple de aplicare şi chestiuni de verificare propuse.

Autorul

CUPRINS

5

CUPRINS

Prefaţă ......................................................................................................... 3 Cuprins ........................................................................................................ 5 CAPITOLUL 1 - Reprezentarea algoritmilor în pseudo-cod ……………. 7 CAPITOLUL 2 - Rezolvarea ecuaţiilor algebrice prin metoda generală … 17 CAPITOLUL 3 - Rezolvarea ecuaţiilor algebrice prin metoda bisecţiei … 23 CAPITOLUL 4 - Rezolvarea ecuaţiilor algebrice prin metoda

Newton – Raphson ……….. 27 CAPITOLUL 5 - Rezolvarea ecuaţiilor algebrice prin metoda

punctului fix ……………… 33 CAPITOLUL 6 - Rezolvarea sistemelor de ecuaţii liniare prin

metoda Gauss …………….. 37 CAPITOLUL 7 - Rezolvarea sistemelor de ecuaţii liniare prin

metoda Jacobi …………….. 43 CAPITOLUL 8 - Rezolvarea sistemelor de ecuaţii liniare prin

metoda Gauss - Seidel ……. 49 CAPITOLUL 9 - Polinomul de interpolare a lui Lagrange ……………… 53 CAPITOLUL 10 - Polinomul de interpolare a lui Newton ………………. 59 CAPITOLUL 11 - Aproximarea cu abatere medie pătratică minimă ……. 67 CAPITOLUL 12 - Derivarea cu ajutorul polinomului de interpolare

Lagrange …………………. 71 CAPITOLUL 13 - Metode numerice pentru calculul integralei definite … 75 CAPITOLUL 14 - Exemple de aplicare a metodelor numerice

în ingineria tehnologică ….. 81 Bibliografie ................................................................................................. 90

CAPITOLUL 1 – Reprezentarea algoritmilor în pseudo-cod

7

Capitolul 1

REPREZENTAREA ALGORITMILOR ÎN PSEUDO-COD.

1. Caracterizarea lucrării

Rezolvarea problemelor ştiinţifice şi tehnice cu ajutorul calculatorului

presupune, mai întâi, găsirea unui algoritm specific şi, apoi, implementarea acestuia pe un sistem de calcul.

În afara algoritmului de rezolvare propriu-zis, un rol important în găsirea eficientă a soluţiei îl joacă alegerea corectă a structurilor de date utilizate. În esenţă, programul de calcul se bazează pe algoritm (descrierea operaţiilor ce vor fi efectuate pentru găsirea soluţiei) şi pe structura de date (modul în care se reprezintă datele de intrare, variabilele intermediare şi datele de ieşire).

Pentru descrierea structurilor de date şi a algoritmilor se poate utiliza un limbaj de programare (BASIC, PASCAL, FORTRAN, C++ etc.) sau un pseudo-limbaj cu o sintaxă mai permisivă (mai puţin rigidă). Pentru a evidenţia invarianţa algoritmilor la limbajul utilizat, de regulă se recurge la a doua soluţie.

Scopul acestei lucrări este de a familiariza utilizatorii cu gândirea algoritmică structurată (prin utilizarea pseudo-codului), de a evidenţia metodele de implementare a pseudo-codului în diferitele limbaje de programare şi, nu în ultimul rând, de a evidenţia importanţa tipurilor abstracte de date, cu caracter matematic (vectori, matrice, numere complexe, etc.).

METODE NUMERICE

2. Algoritmi de calcul. Stabilitatea şi condiţionarea algoritmilor.

Rezolvarea unei probleme practice, pornind de la un anumit set de date iniţiale, presupune utilizarea unui algoritm de calcul, pentru obţinerea datelor de ieşire căutate, potrivit schemei

8

Algoritmul de calcul reprezintă un sistem de reguli care transformă datele de intrare în date de ieşire, cu ajutorul unor operaţii succesive, unic determinate. Un algoritm de calcul trebuie să îndeplinească următoarele cerinţe:

Generalitate; se înţelege prin aceasta că nu este suficient ca un algoritm să permită rezolvarea doar a unei probleme oarecare, ci trebuie să permită rezolvarea tuturor problemelor din categoria respectivă.

Finitudine; numărul de transformări intermediare (iteraţii) aplicate datelor de intrare pentru a obţine datele de ieşire trebuie să fie finit.

Unicitate; transformările intermediare trebuie să fie unic determinate (să conducă de fiecare dată la obţinerea aceloraşi rezultate).

Un algoritm de calcul este stabil dacă aplicat unei probleme cu date

iniţiale „uşor perturbate” conduce la o soluţie apropiată, într-un anume sens, de soluţia problemei cu datele iniţiale.

O altă noţiune, înrudită cu cea de stabilitate a unui algoritm, este cea de condiţionare a acestuia; un algoritm este bine condiţionat dacă mici erori relative în datele iniţiale se transmit în mici erori relative ale datelor de ieşire.

Condiţionarea unui algoritm se poate ilustra simplu, spre exemplificare, pentru calculul unei funcţii z = f(x, y) într-un punct (x, y). Aplicând formula lui Taylor se poate scrie

( ) ( ) ...yfxfy,xfyy,xxff yx +Δ⋅+Δ⋅=−Δ+Δ+=Δ (1)

În cazul în care y,x ΔΔ sunt suficient de mici, termenii de rang mai mare sau egal cu doi din (1.19) se pot neglija; rezultă

yfxfyfxfdff yxyx Δ⋅+Δ⋅≤Δ⋅+Δ⋅=≅Δ . (2)

Δf reprezintă, în această situaţie, eroarea absolută care afectează rezultatele ca efect al erorilor Δx şi Δy apărute în datele iniţiale; mai semnificativă este, însă, eroarea relativă a rezultatelor, care, pentru f 0 este ≠

Date de intrare Algoritm de calcul Date de ieşire

CAPITOLUL 1 – Reprezentarea algoritmilor în pseudo-cod

yf

fx

ff

ff yx Δ⋅+Δ⋅≤

Δ, (3)

relaţie care, pentru x, y 0, se poate rescrie ca ≠

yy

f

fy

xx

ffx

ff yx Δ

⋅⋅

+Δ

⋅⋅

≤Δ

. (4)

Kx Ky

Coeficienţii Kx şi Ky de mai sus sunt denumiţi indici de condiţionare a problemei; se poate observa că algoritmul poate să nu fie bine condiţionat, în cazul în care aceşti indici au valori mari, chiar atunci când Δx şi Δy au valori relativ mici. 3. Descrierea pseudo-limbajului

Pseudo-codul este o metodă simplă şi eficientă pentru reprezentarea unui algoritm şi a structurilor de date asociate. Un algoritm redactat în pseudo-cod este, de fapt, un text alcătuit din linii (rânduri), fiecare dintre acestea conţinând, de regulă, o declaraţie (al cărui scop principal constă în descrierea datelor) sau o instrucţiune (care descrie o operaţie care urmează a fi efectuată).

Fiecare linie dintr-un text în pseudo-cod este alcătuită din cuvinte şi simboluri (caractere speciale nealfanumerice). Anumite cuvinte, cu o semnificaţie bine determinată, independentă de aplicaţie, se numesc „cuvinte cheie” şi, pentru a fi distinse faţă de celelalte cuvinte (specifice aplicaţiei), sunt, de obicei, subliniate (îngroşate).

Orice linie poate conţine precizări suplimentare, numite comentarii, care ajută la o mai bună înţelegere a algoritmului redactat în pseudo-cod, fără a face parte din descrierea propriu-zisă a acestuia sau din structura de date. Comentariile sunt amplasate la sfârşitul liniilor la care se referă şi încep cu caracterul „;”.

3.1. Structuri de date

Declaraţiile se referă la datele cu care se operează şi pot fi de tip simplu (fundamental) sau structurate (agregate).

Se consideră următoarele tipuri de date fundamentale:

logic - date cu două valori (0 = fals, 1 = adevărat);

întreg - date care pot lua doar valori întregi;

real - date ale căror valori aproximează numere reale;

caracter - literă, cifră sau semn special (aritmetic sau de punctuaţie). 9

METODE NUMERICE

10

Exemple de declaraţii ale unor variabile de tip fundamental:

logic l1, l2, l3 întreg i, j, s real a, m, x, j caracter c

După cum se poate constata, o declaraţie conţine cuvântul-cheie care specifică tipul, urmat de lista numelor variabilelor de tipul respectiv, separate între ele prin virgule. Numele variabilelor pot să fie alcătuite din litere şi cifre, primul caracter fiind, obligatoriu, o literă şi au semnificaţie limitată doar la aplicaţia respectivă. Se recomandă ca, în limita posibilităţilor, numele variabilelor să fie alese cât mai sugestiv iar, pentru eliminarea ambiguităţilor, fiecare mărime cu care se operează într-o problemă să aibă propria declaraţie, eventual însoţită de un comentariu lămuritor:

real m ; masa corpului real v ; viteza de aşchiere

Pentru rezolvarea unor probleme mai complexe, tipurile fundamentale de date nu sunt suficiente, recurgându-se la tipuri de date structurate, apelate prin cuvintele-cheie:

tablou - structură de date care conţine un număr cunoscut de elemente de acelaşi tip;

înregistrare - structură de date ce poate conţine elemente de tipuri diferite.

Exemple de declaraţii pentru tablouri:

tablou real V(3) ; V este un vector cu trei elemente reale tablou întreg A(5, 5); A este o matrice 5 x 5 cu elemente

întregi

Cu toate că tabloul are o singură denumire pentru întreaga structură de date, elementele acestuia se identifică prin index (număr de ordine), ca de exemplu V(2); A(1, 2); A(5, 4) etc.

3.2. Structuri de control

Instrucţiunile unui algoritm redactat în pseudo-cod descriu operaţiile pe care trebuie să le efectueze sistemul de calcul cu datele anterior descrise prin declaraţii.

Instrucţiunile sunt de două tipuri: simple şi structurate; cele simple sunt, la rândul lor, de trei feluri:

de atribuire;

CAPITOLUL 1 – Reprezentarea algoritmilor în pseudo-cod

de intrare;

de ieşire.

Instrucţiunea de atribuire are sintaxa:

variabilă = expresie

în care „variabilă” este numele unei variabile a cărei valoare urmează să fie modificată în urma instrucţiunii iar „expresie” este o construcţie sintactică alcătuită din constante, variabile, operatori şi paranteze, după regulile din algebră. Efectul execuţiei instrucţiunii de atribuire constă în evaluarea expresiei şi modificarea, în concordanţă, a variabilei al cărui nume se află în stânga semnului egal.

Se consideră că operanzii care intervin în expresii au valori corespunzătoare unuia dintre tipurile fundamentale. Dacă operanzii sunt de tip logic, atunci se admit operatori logici precum sau, şi, nu, ca în exemplele:

l1 = nu(l2) l3 = l1 sau l2 l4 = l1 şi l3

Dacă operanzii sunt numerici (real sau întreg), se admit operatori aritmetici: +, -, *, / sau de relaţie ( ≥≤><≠= ,,,,, ), ca în exemplele:

bx2bxay

−⋅+⋅

=

( ) ( )axsibx1l >≤=

În primul caz rezultatul este de tip numeric iar, în al doilea, de tip logic.

Instrucţiunile de intrare – ieşire au sintaxa:

citeşte variabilă (listă variabile) scrie variabilă (listă variabile)

Prima instrucţiune are ca efect transferul pe canalul de intrare a unei valori (de exemplu introducerea de la tastatură), care modifică valoarea variabilei specificate, iar a doua are ca urmare transferul valorii variabilei pe canalul de ieşire (de exemplu, afişarea pe ecran sau tipărirea la imprimantă).

De exemplu, pseudo-codul:

real x, y, s, p citeşte x, y

s = x + y p = x * y

scrie s, p

11stop

METODE NUMERICE

reprezintă algoritmul unui program care calculează suma şi produsul a două numere reale, ale căror valori se citesc de la tastatură.

Programele se încheie printr-o altă instrucţiune simplă, cu sintaxa

stop

care are ca efect terminarea execuţiei programului respectiv. Pentru realizarea unor operaţii mai complicate se utilizează, în afara

instrucţiunilor simple, instrucţiunile structurate, care sunt:

secvenţa;

decizia (cu sau fără alternativă);

ciclul (cu test iniţial, cu test final sau cu contor);

rutina (procedură sau funcţie).

Secvenţa sau blocul de instrucţiuni reprezintă un şir de instrucţiuni simple sau structurate (redactate câte una pe linie, în linii succesive), care se execută una după alta, în ordinea în care au fost scrise; programul de mai sus este un exemplu de secvenţă.

Decizia este o instrucţiune care permite controlul execuţiei şi are una dintre următoarele variante sintactice:

decizia simplă:

dacă condiţie atunci secvenţă

decizia cu alternativă:

dacă condiţie atunci secvenţa 1 altfel secvenţa 2

în care „condiţie” este o expresie de tip logic, iar „secvenţa” este o secvenţă de una sau mai multe instrucţiuni. Pentru a uşura înţelegerea acestei instrucţiuni, se constată că secvenţele sunt redactate indentat – decalate faţă de cuvântul-cheie dacă. În urma execuţiei acestei instrucţiuni, se evaluează expresia logică „condiţie”; dacă valoarea rezultată este adevărată, atunci se execută „secvenţa” (respectiv „secvenţa 1”), altfel se continuă cu instrucţiunea următoare (respectiv se execută „secvenţa 2”).

Exemple de instrucţiuni de decizie:

; calculul funcţiei modul, xy = dacă atunci 0x ≥ y = x

12

CAPITOLUL 1 – Reprezentarea algoritmilor în pseudo-cod

altfel y = -x

; calculul funcţiei ⎪⎩

⎪⎨

⎧

>

≤≤−−<

=

1x,x

1x1,x21x,0

y2

dacă x<-1 atunci y = 0 altfel dacă (x -1) şi (x≥ ≤1) atunci y = 2x altfel y = x2

Ciclul reprezintă o instrucţiune care permite repetarea unei secvenţe; se pot utiliza trei tipuri de cicluri: cu test iniţial, cu test final şi cu contor.

Ciclul cu test iniţial are sintaxa:

cât timp condiţie repetă secvenţă

Efectul este evaluarea expresiei logice „condiţie”; dacă rezultatul este afirmativ (adevărat), atunci se execută secvenţa şi ciclul se reia până când „condiţia” devine falsă, după care se sare peste „secvenţă” şi se continuă cu următoarea instrucţiune. Se constată că există riscul repetării la infinit a ciclului dacă valoarea logică a condiţiei rămâne mereu adevărată; este deci responsabilitatea programatorului să asigure caracterul finit al ciclului.

În exemplul următor:

k = 1 s = 0 cât timp ak > 0 repetă s = s + ak k = k + 1 scrie s

poate fi utilizat pentru a aduna elementele unui tablou (de tip şir) până la întâlnirea primului element negativ.

Ciclul cu test final are sintaxa:

repetă secvenţa până când condiţie

„secvenţa” fiind executată repetat până când „condiţia” devine adevărată. 13

METODE NUMERICE

Deosebirea faţă de ciclul cu test iniţial constă în aceea că, în acest caz, corpul ciclului este executat cel puţin o dată. În exemplul următor, se calculează cu o eroare maximă impusă, suma seriei cu termenul general (-1)k/k!.

t = 1 k = 0 s = 0 repetă t = -t / k s = s + t k = k + 1 până când t < eps

Ciclul cu contor permite repetarea unei secvenţe de un număr determinat de ori. Sintaxa ciclului cu contor este:

pentru contor = val_ini, val_fin, pas repetă secvenţă

în care „contor” este numele unei variabile de tip întreg, „val_ini”, „val_fin” şi „pas” sunt constante sau expresii de tip întreg. Secvenţa de instrucţiuni din corpul ciclului este repetată pentru valori succesive ale contorului, pornind de la valoarea iniţială „val_ini”, incrementat cu pasul „pas”, până când acesta depăşeşte valoarea finală „val_fin”. Dacă în instrucţiune lipseşte valoarea „pas”, se presupune că aceasta are valoarea implicită 1.

În exemplul următor:

s = 0 pentru k = 1, n s = s + ak scrie s

se calculează suma primelor n elemente ale unui tablou.

Se întâlnesc, frecvent, situaţii în care o anumită secvenţă de instrucţiuni trebuie executată de mai multe ori, în momente diferite ale execuţiei unui program. Pentru a se evita rescrierea de mai multe ori a acestei secvenţe se recurge la conceptul de rutină.

Rutina reprezintă o secvenţă de declaraţii şi instrucţiuni căreia i se atribuie un nume. Dacă într-un program se face apel la o rutină, controlul execuţiei se transferă rutinei, iar după încheierea acesteia se revine în programul apelant.

14

În consecinţă, o rutină presupune, pe de o parte, definiţia acesteia (prin specificarea declaraţiilor şi instrucţiunilor care o alcătuiesc) iar, pe de altă parte, apelarea ei. Utilizarea rutinelor este justificată chiar şi în cazul în care ele sunt apelate o singură dată într-un program, deoarece ele permit structurarea

CAPITOLUL 1 – Reprezentarea algoritmilor în pseudo-cod

15

modulară a unui algoritm. Pentru a realiza acest deziderat, o rutină trebuie să aibă o anumită consistenţă, să îndeplinească o funcţie bine definită, care să permită reutilizarea ei şi în cadrul altor programe.

După modul de apelare, rutinele se împart în două categorii: proceduri şi funcţii.

Procedura este definită printr-o construcţie sintactică de forma:

procedură nume(listă de parametri formali) secvenţă retur

în care „nume” este numele procedurii, alcătuit din caractere alfanumerice, iar lista parametrilor formali conţine nume de parametri, separate prin virgule. O parte din parametrii formali sunt parametri de intrare (ale căror valori provin din programul apelant şi se transferă procedurii) iar alţii sunt parametri de ieşire (ale căror valori se determină în procedură şi se transferă, apoi, programului apelant). Apelarea unei proceduri astfel definite se face prin invocarea numelui ei, urmat de lista parametrilor actuali:

nume(listă de parametri actuali)

Parametrii actuali trebuie să concorde, ca număr, tip şi ordine, cu cei formali (dar nu neapărat şi ca nume). În exemplul următor,

citeşte a, b sumprod(a, b, s, p) scrie s, p sumprod(s, p, s1, p1) scrie s1, p1

se face apel la procedura definită în continuare:

procedură sumprod(x, y, sumă, prod) real x, y, sumă, prod sumă = x + y prod = x * y retur

care are pe x şi y ca parametri formali de intrare şi sumă, prod – parametri formali de ieşire.

În urma primului apel al acestei proceduri (cu parametri actuali de intrare a şi b) se calculează şi afişează s = a + b, p = a * b iar în urma celui de-al doilea apel se calculează şi se afişează s1 = s + p + a * b, şi p1 = (a + b) * a * b.

În cazul unei proceduri, numărul parametrilor de intrare sau de ieşire este arbitrar (poate fi, inclusiv, nul).

Funcţia este o variantă de rutină la care toţi parametrii formali sunt parametri de intrare, dar rutina întoarce o valoare. Definiţia funcţiei se realizează prin construcţia sintactică

METODE NUMERICE

16

funcţia nume(listă de parametri formali) secvenţă întoarce valoare

Funcţia poate fi apelată ca operand într-o expresie, în particular în atribuirea

valoare = nume(listă de parametri actuali)

Funcţia se aseamănă cu o procedură cu un singur parametru de ieşire, dar, faţă de aceasta, are o flexibilitate suplimentară în apelare. Exemplul următor

real a, b citeşte a, b x = min(a, b) / max(a, b) scrie a, b

recurge la utilizarea a două funcţii, min şi max, definite prin

funcţia max(x, y) ; valoarea maximă real x, y, m dacă x > y atunci m = x altfel m = y întoarce m

funcţia min(x, y) ; valoarea minimă real x, y, m dacă x > y atunci m = y altfel m = x întoarce m

Funcţiile elementare (modul, radical, putere, exponenţială, logaritm, sinus, cosinus etc.) se consideră predefinite în pseudo-cod, deoarece majoritatea limbajelor de programare de nivel înalt le au implementate.

4. Chestiuni de verificare

4.1. Să se scrie o procedură care să calculeze produsul scalar a doi vectori. 4.2. Să se scrie o procedură de înmulţire a matricelor dreptunghiulare. 4.3. Să se scrie algoritmul unui program care să calculeze

⎪⎩

⎪⎨

⎧

>≤≤−−

−<=

3x,xln3x5,3x5

5x,xsiny

CAPITOLUL 2 – Rezolvarea ecuaţiilor algebrice prin metoda generală

Capitolul 2

REZOLVAREA ECUAŢIILOR ALGEBRICE PRIN METODA GENERALĂ.

1. Principiul metodei

Fie funcţia şi ecuaţia f(x) = 0. Dacă f este un polinom

sau poate fi adusă la formă polinomială, în urma unor transformări, ecuaţia se numeşte algebrică; în caz contrar este vorba de o ecuaţie transcendentă.

RRD:f →⊂

Algebra clasică oferă formule şi algoritmi care permit rezolvarea exactă doar pentru o mică parte din ecuaţiile care intervin în soluţionarea diverselor probleme practice din domeniul ingineriei; chiar şi atunci când astfel de posibilităţi există (dar aplicarea lor este complicată), ca şi în toate celelalte cazuri, rezolvarea numerică a ecuaţiilor este calea de urmat.

Rezolvarea oricărei ecuaţii, pe cale numerică, presupune parcurgerea a două etape:

Separarea rădăcinilor ecuaţiei, care presupune determinarea unui interval de existenţă pentru fiecare dintre rădăcinile reale ale ecuaţiei.

Calculul aproximativ pentru fiecare rădăcină reală a ecuaţiei şi evaluarea erorii.

În ceea ce priveşte separarea rădăcinilor, aceasta se poate realiza, spre exemplu prin metoda clasică a şirului lui Rolle; alte posibilităţi sunt analiza datelor problemei practice ce conduce la ecuaţia respectivă, sau aproximarea

17

METODE NUMERICE

funcţiei f printr-un polinom de interpolare, rădăcinile polinomului de interpolare servind ca valori aproximative iniţiale ale rădăcinilor ecuaţiei.

Cu toate că metoda generală de rezolvare numerică a unei ecuaţii poate fi apreciată ca fiind mai rudimentară decât metodele „clasice”, aplicabilitatea metodei este extrem de largă, ea ilustrând foarte bine mecanismul general al algoritmilor de căutare a soluţiei unei ecuaţii.

Fig. 1 - Metodă generală de rezolvare numerică a ecuaţiilor

Fie , despre care se cunoaşte că [ ] Rb,a:f → ( ) ( )b,a∈ξ∃ , unic, astfel

încât f(ξ) = 0; considerându-se, potrivit primului criteriu enunţat mai sus, o valoare maximă admisibilă a erorii er1 (renotată, pentru simplitate, cu er), se va considera soluţie a ecuaţiei orice număr xk pentru care ( ) erx kf ≤ .

Metoda presupune, mai întâi, generarea unui şir de (n+1) puncte echidistante între a şi b, potrivit relaţiei

( ) ( )1n...2,1k,n

ab1kax k −=−

⋅−+= , (1)

după care se testează condiţia ( ) erxf k ≤ ; în cazul în care răspunsul este afirmativ, punctul în cauză este inclus în mulţimea soluţiilor.

După cum se poate observa şi din figura 1, în general există mai multe puncte din şirul (1) care pot fi considerate soluţii (în cazul din figură - xk-2 , xk-1, xk ...). Numărul acestor puncte este egal cu numărul termenilor din şir cuprinşi în intervalul [α, β], unde α şi β sunt soluţiile din intervalul [a, b] ale ecuaţiei ( ) erxf = , α < β.

Pentru a alege din mulţimea soluţiilor aproximative pe cea mai apropiată de soluţia exactă, ξ, se poate recurge la două strategii:

18

CAPITOLUL 2 – Rezolvarea ecuaţiilor algebrice prin metoda generală

1) se calculează eroarea εk = f(xk), pentru fiecare soluţie potenţială, xk şi se alege ca soluţie finală cea pentru care eroarea este minimă, în valoare absolută;

2) se reduce valoarea erorii maxime admisibile, er, reducându-se astfel numărul de soluţii aproximative găsite; în cazul unui program de calcul construit pe baza acestei metode, prin rularea repetată a acestuia, pentru valori din ce în ce mai mici ale parametrului er, se poate ajunge, la un moment dat, la situaţia în care mulţimea soluţiilor potenţiale conţine un singur element, care poate fi declarat soluţie a ecuaţiei.

2. Algoritmul metodei program Metoda_generală real a,b ; domeniul de definiţie al funcţiei f real er ; eroarea maximă admisibilă a soluţiei real eps ; eroarea curentă întreg n ; numărul de iteraţii de efectuat întreg k ; contorul iteraţiilor

citeşte a, b, n, er pentru k = 1, n-1 repetă x = a + k * (b - a) / n dacă ( ) erxf < atunci eps = ( )xf scrie x ; soluţia aproximativă a ecuaţiei scrie eps ; eroarea între f(x) şi 0

stop 3. Exemplu de aplicare

Fie ecuaţia x = cos(x); dacă se consideră funcţia , f(x) = x - cos(x), se poate observa că f este continuă şi derivabilă că f(0) = -1, în

timp ce

RR:f →

22f π

=⎟⎠⎞

⎜⎝⎛ π .

Conform teoremei lui Rolle, rezultă că pe intervalul ⎟⎠⎞

⎜⎝⎛ π

2,0 funcţia f

se anulează cel puţin o dată; punctul ( )b,a∈ξ pentru care f(ξ) = 0 reprezintă soluţia exactă a ecuaţiei considerate.

Dacă se utilizează un program de calcul bazat pe această metodă pentru rezolvarea ecuaţiei, cu diverse valori ale parametrului n (din relaţia (1)) şi pentru diverse valori ale erorii maxime admisibile, er, se obţin rezultatele din tabelul 1.

19

METODE NUMERICE

După cum s-a menţionat deja, sunt situaţii în care se obţine o mulţime de soluţii potenţiale, dintre care trebuie selectată cea mai bună; dacă se utilizează a doua strategie de selecţie, pentru n = 100, în cazul în care er = 0.1 se găsesc 7 soluţii potenţiale, dacă er = 0.05 - 3 soluţii, în timp ce pentru er = 0.01 - o singură soluţie.

Evident, valoarea pentru er trebuie aleasă judicios, în corelaţie şi cu numărul de puncte, n; o reducere peste o anumită limită a erorii maxime admisibile conduce la a nu găsi nici o soluţie.

Dacă se recurge la prima strategie de selecţie a celei mai bune soluţii, programul trebuie să furnizeze şi valorile funcţiei f în fiecare din punctele corespunzătoare soluţiilor potenţiale (coloana a patra din tabelul 1); pe această bază se poate determina, cu uşurinţă, soluţia aproximativă a ecuaţiei.

Tab. 1

n

er

xk

f(xk)

0.691 -0.079 0.707 -0.054 0.722 -0.028 0.738 -0.001 0.754 0.025 0.770 0.052

0.1

0.785 0.078 0.722 -0.028 0.738 -0.001

0.05

0.754 0.025

100

0.02 0.738 -0.001 0.737 -0.004 0.738 -0.001 0.740 0.001

0.005

0.741 0.004 0.738 -0.001

1000

0.002 0.740 0.001

10000 0.0001 0.73906 -0.000043

Trebuie observat că metoda poate rezolva, la nivelul actual al calculatoarelor, o varietate largă de ecuaţii, rapid şi cu o precizie mai mult decât satisfăcătoare.

20

CAPITOLUL 2 – Rezolvarea ecuaţiilor algebrice prin metoda generală

21

4. Chestiuni de verificare

Să se rezolve prin metoda generală următoarele ecuaţii de forma f(x) = 0, cu precizia specificată, ştiind că admit câte o soluţie pe intervalele indicate:

4.1. f(x) = x2 – 5x + 4, er = 10-3, ( )3,0ξ∈ . 4.2. f(x) = x3 - 8x2 + 19x -12, er = 10-5, ( )2,1ξ −∈ . 4.3. f(x) = xlnxsin − , er = 10-6, ( )π,0ξ∈ .

CAPITOLUL 3 – Rezolvarea ecuaţiilor algebrice prin metoda bisecţiei

Capitolul 3

REZOLVAREA ECUAŢIILOR ALGEBRICE PRIN METODA BISECŢIEI.

1. Principiul metodei

Această metodă poate fi considerată o variantă mai evoluată a metodei

generale, atât deoarece viteza cu care se determină soluţia aproximativă este mai mare, cât şi pentru că permite cunoaşterea, de la bun început a numărului de iteraţii care trebuie parcurse pentru rezolvarea ecuaţiei, cu o precizie impusă.

Condiţia ca prin această metodă să se poată rezolva o ecuaţie de forma f(x) = 0 este ca f să existe un interval compact [a, b] astfel încât

[ ] [ ]( ) ( ) ( ) 0bfaf,b,aCf,Rb,a:f 0 <⋅∈→ . (1)

Condiţia (1) asigură faptul că, pe intervalul [a, b], f are un număr

impar de rădăcini (cel puţin una). Aplicarea metodei mai presupune şi impunerea iniţială a unei erori maxime admisibile, ε, pentru soluţia aproximativă.

Metoda constă, în principiu, într-o serie de înjumătăţiri succesive ale intervalului [a, b] şi derivatelor acestuia, prin care se generează un şir de intervale [a1 , b1], [a2 , b2] . . . [an, bn] . . ., lungimea fiecăruia fiind egală cu jumătatea intervalului precedent (a se vedea şi figura 1). Mijloacele acestor intervale - c1 , c2 . . . cn . . . se constituie, astfel, într-un şir de aproximări succesive ale soluţiei exacte, ξ, care se va demonstra că este convergent către ξ.

23

METODE NUMERICE

Fig.1 - Metoda înjumătăţirii intervalului

Algoritmul de aplicare a metodei este prezentat în cele ce urmează, în

ipoteza în care f(a) > 0 şi f(b) < 0 (cealaltă variantă de aplicare a metodei şi anume când f(a) < 0 şi f(b) > 0 derulându-se analog).

Se calculează 2

bac0+

= , după care se evaluează funcţia f în punctul c0;

în cazul în care f(c0) = 0, înseamnă că, absolut întâmplător, s-a găsit soluţia exactă a ecuaţiei, ξ = c0 .

Dacă f(c0) > 0, se renotează c0 cu a1 şi b cu b1 (cazul din fig.2.3) iar, în caz contrar, a devine a1 şi c0 - b1 .

Se calculează 2

bac 11

1+

= , după care se determină f(c1); dacă f(c1) = 0

rezultă ξ = c1 .

Se înjumătăţeşte intervalul [a1 , b1], rezultând [a2, b2], fie renotându-se c1 cu a2 şi b1 cu b2 (cazul figurat), fie c1 cu b2 şi a1 cu a2 , în aşa fel încât să avem f(a2) · f(b2) < 0.

Se determină c2 , apoi a3 şi b3 , c3 , a4 şi b4 . . . şi aşa mai departe, până când, la un moment dat cn se constituie într-o aproximaţie suficient de bună a soluţiei exacte şi se adoptă ξ = cn .

Lungimea lk a intervalului oarecare din şir, [ak ,bk], rezultat din

intervalul [a, b] după k înjumătăţiri succesive este

24

CAPITOLUL 3 – Rezolvarea ecuaţiilor algebrice prin metoda bisecţiei

kk

2abl −

= . (2)

Prin procedeul enunţat mai sus, se generează în paralel două şiruri - şi - primul dintre ele monoton crescător iar celălalt, monoton

descrescător, ambele având ca limită soluţia exactă a ecuaţiei: ( ) 0nna ≥ ( ) 0nnb ≥

ξ==∞→∞→

nnnnblimalim . (3)

Această afirmaţie poate fi susţinută dacă se observă că eroarea εn făcută prin aproximarea lui ξ cu cn este mărginită superior de un şir cu limita 0:

0lc nnn →≤−ξ=ε , atunci când . (4) ∞→n

Mai mult, dacă se impune, apriori, o limită maximă a erorii admisibile, er, din impunerea condiţiei ercn ≤ξ− se poate determina numărul minim de iteraţii (înjumătăţiri) care trebuie efectuate pentru a găsi soluţia aproximativă a ecuaţiei, Nmin:

12lner

abln1

erablogN 2min +

⎥⎥⎥⎥

⎦

⎤

⎢⎢⎢⎢

⎣

⎡ −

=+⎥⎦⎤

⎢⎣⎡ −

≥ , (5)

unde cu paranteze drepte s-a notat funcţia parte întreagă.

2. Algoritmul metodei program Metoda_bisecţiei real a,b ; domeniul de definiţie al funcţiei f real er ; eroarea maximă admisibilă întreg n ; numărul de iteraţii de efectuat întreg k ; contorul iteraţiilor

citeşte a, b, er n = [ln((b - a) / er) / ln(2)] + 1 ; calculul numărului de iteraţii k = 0 ; iniţializarea contorului

repetă k = k + 1 c = (a + b) / 2 dacă f(a) * f(c) > 0 atunci a = c altfel b = c

25

METODE NUMERICE

26

până când k > n

scrie x ; soluţia aproximativă a ecuaţiei top

. Exemplu de aplicare

x = c s 3

os(x), se poate observa că f este continuă şi derivabilă că

timp ce

Fie ecuaţia x = cos(x); dacă se consideră funcţia RR:f → , f(x) = x - c f(0) = -1, în

22⎟⎠

Considerând această ecuaţie şi pornind, de asemenea, de la intervalul

delimitat de a = 0 şi

f π=⎞⎜

⎝⎛ π .

2b = , s-a utilizat un program deπ calcul bazat pe algoritmul

meto

10-3 trebuie efectuate Nmin = 11 iteraţii, soluţia aproximativă fiind

-6 trebuie efectuate Nmin = 21 iteraţii, soluţia aproximativă fiind

ebuie efectuate Nmin = 31 iteraţii, soluţia aproximativă fiind = 0

fiind relativ ic, chiar în cazul unor erori maxime admisibile extrem de reduse.

. Chestiuni de verificare

dei da faţă. Rezultatele sunt următoarele:

- pentru er =x = 0.7386;

- pentru er = 10x = 0.7390859;

- pentru er = 10-9 trx .7390851328.

După cum se poate observa, soluţia ecuaţiei este găsită mult mai repede faţă de cazul metodei generale, numărul de iteraţii necesarem 4

teraţii necesare cu cel de la rezolvarea aceloraşi ecuaţii prin metoda generală.

.1. f(x) = x2 – 5x + 4, er = 10-3,

Să se rezolve prin metoda bisecţiei ecuaţiile de mai jos, de forma f(x) = 0, cu precizia specificată, ştiind că admit câte o soluţie pe intervalele indicate. Să se compare numărul de i

( )3,0ξ∈ .

.2. f(x) = x3 - 8x2 + 19x -12, er = 10-5,

4

( )2,1ξ −∈4 .

4.3. f(x) = , er = 10-6,

xlnxsin − ( )π,0ξ∈ .

CAPITOLUL 4 – Rezolvarea ecuaţiilor algebrice prin metoda Newton - Raphson

Capitolul 4

REZOLVAREA ECUAŢIILOR ALGEBRICE PRIN METODA NEWTON - RAPHSON.

1. Principiul metodei

Funcţionarea acestei metode se bazează, ca principiu, pe generarea

unui şir de aproximaţii succesive ale soluţiei exacte, ξ, pentru o ecuaţie de forma f(x) = 0. Şirul de aproximaţii succesive, ( ) Nnnx ∈ , se construieşte pornind de la o valoare de start, x0 , convenabil aleasă, utilizând o formulă de tip recurent

( )

0n,hxf

xxn

nn1n ≥−=+ , (1)

unde hn = h(xn) este o funcţie nenulă (eventual o constantă) aleasă, de asemenea, convenabil.

Potrivit modului de definire a funcţiei hn , există mai multe metode, aşa-numite „clasice”: metoda Newton - Raphson (metoda tangentelor variabile), metoda tangentelor fixe, metoda secantelor variabile (regula Falsi) şi metoda secantelor fixe.

În cazul metodei Newton - Raphson, funcţia hn din relaţia (1) se alege de forma hn = f ’(xn) (derivata funcţiei f). Rezultă o formulă de recurenţă de tipul

( )( ) 0n,xfxf

xxn

nn1n ≥

′−=+ . (2)

27

METODE NUMERICE

Formula (2) mai poartă şi numele de formula iterativă a lui Newton. Pentru rezolvarea numerică a unei ecuaţii de forma f(x) = 0 prin acest

procedeu, funcţia f trebuie să îndeplinească următoarele condiţii: , (adică f derivabilă de două ori pe [a, b], cu derivatele continui), iar

şi f păstrează semnul constant pe tot intervalul; în plus, mai trebuie ca f să aibă o singură rădăcină în intervalul considerat. Punctul de start, x0 , trebuie ales astfel încât

[ ] Rb,a:f →[ ]( b,aCf 2∈

f ′ ′′)

( ) ( ) 0xfx 00 >′′⋅f . Metoda are o interpretare geometrică sugestivă (figura 1).

A0

A1

A2

Fig.1 - Metoda Newton-Raphson - interpretare geometrică

Potrivit formulelor din Geometria analitică, ecuaţia tangentei dusă la curba y = f(x) într-un punct oarecare, de abscisă xk , este

( ) ( )( )kkk xxxfxfy −′=− . (3)

Dacă se notează cu xk+1 abscisa punctului în care această tangentă

intersectează axa Ox (obţinut pentru y = 0 în ecuaţia (3)) şi se exprimă xk+1 , se obţine formula lui Newton de mai sus.

În aceste condiţii, se poate observa că fiecare valoare, xk+1, din şirul aproximaţiilor succesive ale soluţiei ecuaţiei rezultă ca abscisa punctului de intersecţie, cu axa Ox, a tangentei la graficul funcţiei duse în punctul precedent, Ak , de abscisă xk . Metoda se mai numeşte şi metoda tangentelor variabile, deoarece panta fiecărei tangente, pe baza căreia se generează punctele din şir, are, în general, panta diferită de a celorlalte tangente (precedente sau următoare).

Funcţia φ, definită ca 28

CAPITOLUL 4 – Rezolvarea ecuaţiilor algebrice prin metoda Newton - Raphson

( ) ( )( )xfxfxx′

−=ϕ (4)

poartă numele de funcţie de iterare a metodei Newton - Raphson. Se poate verifica imediat că

( ) ξ=ξϕ , ( ) 0=ξϕ′ , ( ) 0≠ξϕ ′′ . (5)

Cu această notaţie, şirul aproximaţiilor succesive ale soluţiei ecuaţiei date se poate genera după formula xk+1 = φ(xk), k = 0, 1, 2, . . . (6)

Pentru estimarea modului în care se propagă eroarea în cazul metodei de faţă, precum şi pentru determinarea ordinului de convergenţă, se consideră dezvoltarea în serie Taylor, într-o vecinătate a punctului ξ,

( ) ( ) ( ) ( ) ( ) ( ) +ξϕ ′′⋅ξ−

+ξϕ′⋅ξ−

+ξϕ=ϕ!2

x!1

xx2

. . . (7)

Dacă în această dezvoltare se înlocuieşte x cu xn, se reţin doar termenii

de rang 0, 1 şi 2 şi se au în vedere relaţiile (5), rezultă

( ) ( ) ( )ξϕ ′′⋅ξ−

+ξ≅ϕ2

xx

2n

n . (8)

Ţinând cont şi de (6), se poate scrie egalitatea

( ) ( )2n1n x2

x ξ−⋅δϕ ′′

≅ξ−+ , (9)

ceea ce este echivalent cu

( ) 2n1n 2ε

ξϕ ′′≅ε + . (10)

Din relaţia (10) rezultă, pe de o parte că metoda Newton - Raphson

este caracterizată de o convergenţă pătratică, fiind una dintre metodele cele mai rapide de rezolvare a unei ecuaţii de forma f(x) = 0; pe de altă parte, însă, metoda reclamă evaluarea funcţiei f şi a derivatei sale la fiecare iteraţie, ceea ce

29

METODE NUMERICE

poate fi, uneori, dificil (sau chiar imposibil, pentru funcţii care nu sunt cunoscute analitic, ci doar tabelar).

2. Algoritmul metodei program Metoda Newton_Raphson

real x0 ; punctul de pornire a căutării iteraţiilor real er ; eroarea maximă admisibilă real xvechi, xnou ; valoarea precedentă / curentă a soluţiei întreg n ; numărul de iteraţii efectuate întreg k ; contorul iteraţiilor

funcţia f(x) ; funcţia care defineşte ecuaţia real x f = …………… ; expresia funcţiei f întoarce f

funcţia df(x) ; derivata funcţiei f real x df = …………… ; expresia funcţiei df întoarce df

citeşte x0, er k = 0 ; iniţializarea contorului xnou = x0 ; iniţializarea şirului aproximaţiilor succesive repetă xvechi = xnou k = k + 1 xnou = xvechi - f(xvechi) / df(xvechi) până când erxvechixnou <− scrie xnou ; soluţia aproximativă a ecuaţiei scrie k ; numărul de iteraţii parcurse

stop 3. Exemplu de aplicare

Pentru a da posibilitatea unei comparaţii reale între diferitele metode numerice de rezolvare a ecuaţiilor, se consideră aceeaşi ecuaţie ca la lucrările 2 şi 3: f(x) = x - cos(x) = 0.

Evident, f’(x) = 1 + sin(x). Dacă se alege ca punct de start 2

bx 0π

==

şi se utilizează un program de calcul scris pe baza metodei de faţă, rezultă următoarele:

30

CAPITOLUL 4 – Rezolvarea ecuaţiilor algebrice prin metoda Newton - Raphson

31

pentru er2 = 10-3 sunt necesari n = 3 iteraţii şi soluţia aproximativă este x = 0.7391;

pentru er2 = 10-6 sunt necesari n = 4 iteraţii şi soluţia aproximativă este x = 0.7390851;

pentru er2 = 10-9 sunt necesari n = 5 iteraţii şi soluţia aproximativă este x = 0.7390851332.

Rezultatele confirmă că metoda Newton - Raphson are o convergenţă mult mai bună decât a celorlalte prezentate până acum. 4. Chestiuni de verificare

Să se rezolve prin metoda Newton - Raphson ecuaţiile de mai jos, de forma f(x) = 0, cu precizia specificată, ştiind că admit câte o soluţie pe intervalele indicate. Să se compare numărul de iteraţii necesare cu cel de la rezolvarea aceloraşi ecuaţii prin metoda generală şi metoda bisecţiei.

4.1. f(x) = x2 – 5x + 4, er = 10-3, ( )3,0ξ∈ . 4.2. f(x) = x3 - 8x2 + 19x -12, er = 10-5, ( )2,1ξ −∈ . 4.3. f(x) = xlnxsin − , er = 10-6, ( )π,0ξ∈ .

CAPITOLUL 5 – Rezolvarea ecuaţiilor algebrice prin metoda punctului fix

Capitolul 5

REZOLVAREA ECUAŢIILOR ALGEBRICE PRIN METODA PUNCTULUI FIX.

1. Principiul metodei

Se spune că o funcţie are un punct fix , dacă g(ξ) = ξ

(fig.1). Altfel spus, o funcţie are unul sau mai multe puncte fixe dacă graficul său intersectează odată - sau de mai multe ori - prima bisectoare.

R:Dg → D∈ξ

Fig.1 - Interpretarea geometrică a punctelor fixe ale unei funcţii

33

METODE NUMERICE

Fie o ecuaţie de forma f(x) = 0, unde [ ] R:b,af → şi să presupunem că ecuaţia considerată are o unică soluţie reală [ ]b,a∈ξ .

Ecuaţia de mai sus poate fi scrisă, adunând x în ambii membri, în forma echivalentă g(x) = x (unde s-a notat f(x) + x = g(x)); după cum se va arăta într-un caz concret, această transformare a ecuaţiei iniţiale se poate realiza, în general, în mai multe moduri. Este uşor de observat că rădăcina ecuaţiei iniţiale satisface egalitatea g(ξ) = ξ, deci este un punct fix pentru funcţia g.

Cu ajutorul funcţiei g astfel obţinute se poate genera, pornind de la un număr arbitrar, x0, un şir de numere reale, utilizându-se formula de recurenţă

( ) 1n,xgx n1n ≥=+ . (1)

Dacă funcţia g este derivabilă pe intervalul [a, b], iar derivata acesteia

îndeplineşte condiţia ( ) ( ) [ ]b,ax,1xg ∈∀<λ≤′ , (2)

atunci se poate demonstra că şirul xn constituie un şir de aproximaţii succesive ale rădăcinii exacte a ecuaţiei, convergent către aceasta.

Într-adevăr, se observă că funcţia g îndeplineşte condiţiile de aplicare a teoremei lui Lagrange pe orice interval de forma [xn, ξ]; de aici rezultă că

astfel încât ( ) ( )ξ∈α∃ ,x nn ( ) ( ) ( )( )ξ−α′=ξ− nnn xggxg . (3)

Dacă se ţine cont de relaţia (1) şi de faptul că g(ξ) = ξ, egalitatea (3)

devine ( )( )ξ−α′=ξ−+ nn1n xgx . (4)

Utilizând în continuare şi relaţia (2) rezultă că ξ−⋅λ≤ξ−+ n1n xx , (5)

ceea ce este echivalent cu n1n ε⋅λ≤ε + . (6)

Dacă se dau valori de la 0 la n în relaţia (6) şi se înmulţesc cele

(n + 1) relaţii astfel obţinute, rezultă că ; cum λ < 1, şi

deci . Acesta este acelaşi lucru cu a afirma că . 0

n1n ε⋅λ≤ε + 0lim n

n=λ

∞→

ξ0lim nn=ε

∞→=

∞→nn

xlim

Tot din relaţia (6) rezultă şi că metoda punctului fix are convergenţă liniară.

34

CAPITOLUL 5 – Rezolvarea ecuaţiilor algebrice prin metoda punctului fix

2. Algoritmul metodei program Metoda punctului fix

real x0 ; punctul de pornire a căutării iteraţiilor real er ; eroarea maximă admisibilă real xvechi, xnou ; valoarea precedentă / curentă a soluţiei întreg n ; numărul de iteraţii efectuate întreg k ; contorul iteraţiilor

funcţia g(x) ; funcţia de punct fix real x g = …………… ; expresia funcţiei g întoarce g

citeşte x0, er k = 0 ; iniţializarea contorului

xnou = x0 ; iniţializarea şirului aproximaţiilor succesive

repetă xvechi = xnou k = k + 1 xnou = g(xvechi) până când erxvechixnou <− scrie xnou ; soluţia aproximativă a ecuaţiei scrie k ; numărul de iteraţii parcurse

stop 3. Exemplu de aplicare

Pentru a da posibilitatea unei comparaţii reale între diferitele metode numerice de rezolvare a ecuaţiilor, se consideră aceeaşi ecuaţie ca la lucrările 2, 3 şi 4: f(x) = x - cos(x) = 0.

Ecuaţia dată se poate scrie sub forma x = cos(x), deci se va adopta funcţia g de forma g(x) = cos(x). Cum g’(x) = -sin(x) se poate observa că

( ) ( ) ( ) ⎟⎠⎞

⎜⎝⎛ π

∈∀−∈′2

,0x,0,1xg , de unde rezultă că pentru această situaţie metoda

punctului fix este convergentă oscilant. Dacă se adoptă drept criteriu de stopare a căutării soluţiei Criteriul 2

(relaţia (2.6)) şi se utilizează un program de calcul bazat pe metoda punctului fix (a se vedea § 2.6.5), rezultă următoarele:

◊ pentru er = 10-3 şi x0 = 0 ⇒x1 = 1.0000, x2 = 0.5403, x3 = 0.8575, x4 = 0.6543, x5 = 0.7935, x6 = 0.7014 . . . x18 = 0.7388;

35

METODE NUMERICE

36

◊ pentru er = 10-3 şi x0 = 2π ⇒x1 = 0.0008, x2 = 0.9999, x3 = 0.5403,

x4 =0.8575, x5 = 0.6543, x6 = 0.7935, x7 = 0.7014 . . . x19 = 0.7388;

◊ pentru er = 10-6 şi x0 = 0 ⇒x35 = 0.7390855;

◊ pentru er = 10-9 şi x0 = 0 ⇒x53 = 0.7390851335.

Se poate observa cum oscilează valorile şirului aproximaţiilor succesive, de o parte şi de alta a soluţiei exacte; de asemenea, se poate remarca şi convergenţa mai lentă a metodei, faţă de metoda Newton – Raphson.

4. Chestiuni de verificare

A. Fie ecuaţia f(x) = x2 - 2x - 3 = 0, având rădăcinile ξ1 = -1 şi ξ2 = 3. Să se studieze posibilitatea rezolvării acestei ecuaţii, prin metoda punctului fix, în următoarele trei cazuri:

4.1. ( ) 3x2xg += , x0 = 4;

4.2. ( )2x

3xg−

= , x0 = 0;

4.3. ( )2

3xxg

2n −

= , x0 = 4.

B. Să se rezolve prin metoda Newton - Raphson ecuaţiile de mai jos,

de forma f(x) = 0, cu precizia specificată, ştiind că admit câte o soluţie pe intervalele indicate. Să se compare numărul de iteraţii necesare cu cel de la rezolvarea aceloraşi ecuaţii prin metoda generală şi metoda bisecţiei.

4.4. f(x) = x2 – 5x + 4, er = 10-3, ( )3,0ξ∈ . 4.5. f(x) = x3 - 8x2 + 19x -12, er = 10-5, ( )2,1ξ −∈ . 4.6. f(x) = , er = 10-6, xlnxsin − ( )π,0ξ∈ .

CAPITOLUL 6 – Rezolvarea sistemelor de ecuaţii liniare prin metoda Gauss

Capitolul 6

REZOLVAREA SISTEMELOR DE ECUAŢII LINIARE PRIN METODA GAUSS.

1. Principiul metodei

În primul rând, trebuie făcută precizarea că metodele numerice pentru

rezolvarea sistemelor de ecuaţii, care urmează a fi prezentate în continuare, sunt destinate abordării sistemelor de n ecuaţii cu n necunoscute, cu coeficienţi reali, despre care se cunoaşte, apriori, faptul că sunt compatibile.

Un sistem de forma

, (1)

⎪⎪⎩

⎪⎪⎨

⎧

=⋅++⋅+⋅

=⋅++⋅+⋅=⋅++⋅+⋅

nnnn22n11n

2nn2222121

1nn1212111

bxa...xaxa.............................................................bxa...xaxabxa...xaxa

în care pentru i = 1, 2, . . . n, j = 1, 2, . . . n, poartă numele de sistem de ecuaţii liniare.

Rb,a iij ∈

Detaliind, în continuare, trebuie arătat că metodele numerice utilizate pentru rezolvarea sistemelor de ecuaţii liniare se împart în două categorii: metode directe, care determină soluţia exactă a sistemului şi metode iterative, care - pornind de la o soluţie iniţială arbitrară - determină un şir de aproximaţii succesive ale vectorului - soluţie a sistemului. 37

METODE NUMERICE

Cea mai veche, cea mai bine cunoscută şi, totodată, cea mai utilizată metodă de rezolvare numerică a sistemelor de ecuaţii algebrice liniare este metoda Gauss, care face parte din categoria metodelor directe; principiul metodei este prezentat în continuare.

Pornindu-se de la un sistem de forma (1), se construieşte o matrice de forma

⎟⎟⎟⎟⎟

⎠

⎞

⎜⎜⎜⎜⎜

⎝

⎛

=

nnn2n1n

2n22221

1n11211

ba...aa...............ba...aaba...aa

A , (2)

rezultată prin alipirea la matricea coeficienţilor sistemului, A, a vectorului coloană a termenului liber, B. În continuare, algoritmul metodei urmăreşte transformarea matricei A într-o matrice ale cărei elemente situate sub diagonala a11 - a22 - . . . -ann să fie nule; trebuie parcurşi următorii paşi:

se calculează o serie de multiplicatori de forma 11

1i1i a

al = , după care din

fiecare linie „i” se scade linia întâi înmulţită cu li1, i = 2, 3, . . . n; în cazul în care a11 = 0, se schimbă prima linie cu o alta, care să aibă primul element nenul, după care se procedează la cele de mai sus. Rezultă, astfel, după renotare, o matrice de forma

⎟⎟⎟⎟⎟

⎠

⎞

⎜⎜⎜⎜⎜

⎝

⎛

=

)1(n

)1(nn

)1(2n

)1(2

)1(n2

)1(22

)1(1

)1(n1

)1(12

)1(11

)1(

ba...a0...............

ba...a0ba...aa

A ; (3)

se calculează o nouă serie de multiplicatori, de forma )1(

22

)1(2i

2ia

al = , după

care din fiecare linie „i” se scade linia întâi înmulţită cu li2, i = 3, 4, . . . n; în cazul în care a22

(1) = 0, se schimbă a doua linie cu o alta, care să aibă al doilea element nenul, după care se procedează la cele de mai sus. Rezultă, astfel, după o nouă renotare, o altă matrice de forma

⎟⎟⎟⎟⎟

⎠

⎞

⎜⎜⎜⎜⎜

⎝

⎛

=

)2(n

)2(nn

)2(2

)2(n2

)2(22

)2(1

)2(n1

)2(12

)2(11

)2(

ba...00...............

ba...a0ba...aa

A ; (4)

38

CAPITOLUL 6 – Rezolvarea sistemelor de ecuaţii liniare prin metoda Gauss

se continuă aplicarea procedurii până când, după n-1 paşi, se obţine matricea

⎟⎟⎟⎟⎟

⎠

⎞

⎜⎜⎜⎜⎜

⎝

⎛

=

−−

−−−

−−−−

−

)1n(n

)1n(nn

)1n(2

)1n(n2

)1n(22

)1n(1

)1n(n1

)1n(12

)1n(11

)1n(

ba...00...............

ba...a0ba...aa

A , (5)

în care toate elementele de sub diagonala a11

(n-1), a22(n-1), . . . ann

(n-1) sunt nule. Dacă se renunţă, în notaţiile din această ultimă matrice, la exponenţii care indică numărul pasului şi se reface, pe baza ei, sistemul de ecuaţii, rezultă:

(6)

⎪⎪⎪

⎩

⎪⎪⎪

⎨

⎧

=⋅=⋅+⋅

=⋅+⋅++⋅=⋅+⋅++⋅+⋅

−−−−−

−−

−−

nnnn

1nnn1n1n1n1n

2nn21n1n2222

1nn11n1n1212111

bxabxaxa

............................................................................bxaxa.........xabxaxa...........xaxa

din ultima ecuaţie a sistemului (6) se determina xn, după care, înlocuind în

penultima ecuaţie pe xn, se determină xn-1 şi aşa mai departe, până când, în final, se înlocuiesc în prima ecuaţie valorile, deja aflate, pentru x2, x3, . . . xn şi se calculează x1.

Conform [11], numărul N de operaţii aritmetice simple necesar a fi

efectuate pentru rezolvarea, prin metoda Gauss, a unui sistem de ecuaţii liniare de rang n este

3

nn4N3 −

= (7)

Metoda Gauss, fiind o metodă directă, conduce la găsirea soluţiei

exacte a sistemului, în măsura în care acesta este compatibil. 2. Algoritmul metodei program Metoda Gauss

real l ; valoarea curentă a multiplicatorului de linie real s ; sumă auxiliară întreg n ; dimensiunea sistemului

39 întreg i, j, k ; contori

METODE NUMERICE

tablou real A(n,n) ; matricea sistemului tablou real B(n) ; vectorul-coloană a termenului liber tablou real X(n) ; vectorul-soluţie curentă

; introducere date de intrare

citeşte n pentru i = 1, n repetă citeşte b(i) pentru j = 1, n repetă citeşte a(i,j)

; eliminare

pentru k = 1, n-1 repetă pentru i = k+1, n repetă l = a(i,k) / a(k,k) pentru j = k+1, n repetă a(i,j) = a(i,j) - l·a(k,j) b(i) = b(i) - l·b(k)

; retro-substituţie

x(n) = b(n) / a(n,n) scrie x(n) pentru i = n-1, 1, -1 repetă s = b(i) pentru j = n, i+1, -1 repetă s = s - a(i,j)·x(j) x(i) = s / a(i,i) scrie x(i) stop 3. Exemplu de aplicare

Să se rezolve sistemul ⎪⎩

⎪⎨

⎧

=++=++

=++

31x17x4x228x10x3x4

11xxx2

321

321

321

.

Se construieşte matricea extinsă a sistemului, ⎟⎟⎟

⎠

⎞

⎜⎜⎜

⎝

⎛=

31174228103411312

A .

Se calculează multiplicatorii 224

aa

l11

2121 === , 1

22

aa

l11

3131 === şi se

efectuează operaţiile L2 - l21 · L1 şi L3 - l31 · L1 (unde Li înseamnă linia „i”).

40

CAPITOLUL 6 – Rezolvarea sistemelor de ecuaţii liniare prin metoda Gauss

41

Rezultă ⎟⎟⎟

⎠

⎞

⎜⎜⎜

⎝

⎛=

2014306410

11312A

)1(; se calculează, apoi, multiplicatorul

313

a

al

)1(22

)1(32

32 === şi se efectuează L3 - l32 · L2, după care se obţine

⎟⎟⎟

⎠

⎞

⎜⎜⎜

⎝

⎛=

22006410

11312A

)2(. Sistemul iniţial poate fi rescris, acum, sub forma

. Urmează imediat, din ultima ecuaţie, x3 = 1, apoi, din a

doua, x2 = 2 şi, în fine, din prima ecuaţie, x1 = 3.

⎪⎩

⎪⎨

⎧

==+=++

2x26x4x

11x3xx2

3

32

321

4. Chestiuni de verificare

Să se rezolve, prin metoda Gauss, sistemele de ecuaţii:

4.1. ⎪⎩

⎪⎨

⎧

−=−+=++−−=+−

.8z6yx2;3z2y4x;5z2y3x7

4.2.

⎪⎪⎩

⎪⎪⎨

⎧

−=−++−=+−+

−=++−=−+−

.1x4xxx;1x2x6x2x

;2xx3x7x;3x2xxx5

4321

4321

4321

4321

CAPITOLUL 7 – Rezolvarea sistemelor de ecuaţii liniare prin metoda Jacobi

Capitolul 7

REZOLVAREA SISTEMELOR DE ECUAŢII LINIARE PRIN METODA JACOBI.

1. Principiul metodei

Metoda lui Jacobi este o metodă iterativă, utilizabilă pentru rezolvarea

numerică a sistemelor de ecuaţii algebrice liniare. Fie sistemul

, (1)

⎪⎪⎩

⎪⎪⎨

⎧

=⋅++⋅+⋅

=⋅++⋅+⋅=⋅++⋅+⋅

nnnn22n11n

2nn2222121

1nn1212111

bxa...xaxa.............................................................bxa...xaxabxa...xaxa

în care pentru i = 1, 2, . . . n, j = 1, 2, . . . n. Rb,a iij ∈

Algoritmul metodei Jacobi începe prin a pune în evidenţă, în partea stângă a semnului egal, în ecuaţia numărul k a variabilei xk, k = 1, 2, . . . n. Evident, pentru a putea rescrie astfel sistemul, trebuie ca toţi coeficienţii de forma akk să fie nenuli; în cazul în care această condiţie nu este satisfăcută, se schimbă ecuaţiile între ele, în aşa fel încât pe diagonala principală a matricei sistemului să nu existe nici un element nul şi se renumerotează necunoscutele în concordanţă cu schimbările efectuate.

Rezultă astfel:

43

METODE NUMERICE

( )

( )

( )⎪⎪⎪

⎩

⎪⎪⎪

⎨

⎧

⋅−−⋅−⋅−⋅=

⋅−−⋅−⋅−⋅=

⋅−−⋅−⋅−⋅=

−− .xa...xaxaba1x

.................................................................................

;xa...xaxaba1x

;xa...xaxaba1x

1n1nn22n11nnnn

n

nn2323121222

2

nn1313212111

1

(2)

Pentru uşurinţa scrierii, se trece la transcrierea sistemului sub forma

unei ecuaţii matriceale x = A · x + u , (3)

în care s-au făcut notaţiile

⎟⎟⎟⎟⎟

⎠

⎞

⎜⎜⎜⎜⎜

⎝

⎛

=

n

2

1

x

xx

xM

,

⎟⎟⎟⎟⎟⎟⎟⎟

⎠

⎞

⎜⎜⎜⎜⎜⎜⎜⎜

⎝

⎛

=

nn

n

22

2

11

1

ab

abab

uM

,

⎟⎟⎟⎟⎟⎟⎟⎟

⎠

⎞

⎜⎜⎜⎜⎜⎜⎜⎜

⎝

⎛

−=

0aa

aa

aa

0aa

aa

aa

0

A

nn

2n

nn

1n

22

n2

22

21

11

n1

11

12

L

MMMM

L

L

. (4)

Principiul metodei lui Jacobi se bazează pe relaţia (3), unde în

membrul stâng se scrie vectorul soluţiilor la pasul k + 1, iar în membrul drept acelaşi vector, la pasul k. Rezultă, astfel

uxAx )k()1k( +⋅=+ . (5)

Evident, demararea calculelor impune adoptarea unui vector iniţial, x0,

care poate fi ales arbitrar. Egalitatea (5) mai poate fi scrisă şi sub forma

⎟⎟⎟

⎠

⎞

⎜⎜⎜

⎝

⎛⋅−= ∑

≠=

+ n

ij1j

)k(jiji

ii

)1k(i xab

a1x , i = 1, 2, . . . n. (6)

Cu ajutorul relaţiei (5) se poate genera un şir de aproximaţii succesive

ale vectorului-soluţie a sistemului, care - în anumite condiţii - converge către soluţia exactă.

În această situaţie, se pot formula condiţiile de convergenţă ale metodei lui Jacobi; dacă 44

CAPITOLUL 7 – Rezolvarea sistemelor de ecuaţii liniare prin metoda Jacobi

1) ; 0Alim kk

=∞→

2) toate valorile proprii ale matricei A sunt, în modul, mai mici decât unitatea,

atunci şirul de aproximaţii succesive ale vectorului-soluţie, generat cu relaţia (5), converge către soluţia exactă.

Condiţia 2) de mai sus este echivalentă, [13], cu relaţia

,1aan

ij1j ii

ij <∑≠=

i = 1, 2, . . . n, (7)

formă sub care este mai uşor de verificat, din punct de vedere practic.

Dacă un sistem de ecuaţii liniare satisface condiţia 2) de convergenţă, sistemul se numeşte diagonal; sistemele diagonale se bucură, deci, de proprietatea că sunt rezolvabile prin metoda iterativă Jacobi.

În ceea ce priveşte criteriul de oprire a căutării soluţiei, dacă x(k) şi x(k+1) sunt două aproximaţii succesive ale vectorului soluţie, generate potrivit relaţiei (5),

, respectiv ,

⎟⎟⎟⎟⎟

⎠

⎞

⎜⎜⎜⎜⎜

⎝

⎛

=

)k(n

)k(2

)k(1

)k(

x

xx

xM

⎟⎟⎟⎟⎟

⎠

⎞

⎜⎜⎜⎜⎜

⎝

⎛

=

+

+

+

+

)1k(n

)1k(2

)1k(1

)1k(

x

xx

xM

atunci vectorul x(k+1) poate fi declarat soluţie a sistemului dacă

( )( ) erxxmaxerr )k(i

1ki

ni1<−= +

≤≤, (8)

unde er reprezintă eroarea maximă admisibilă, stabilită la începutul rezolvării sistemului. 2. Algoritmul metodei program Metoda Jacobi

real er ; eroarea maximă admisibilă real err ; eroarea la pasul curent real s ; suma parţială întreg n ; dimensiunea sistemului întreg i, j, k ; contori tablou real A(n,n) ; matricea sistemului

45 tablou real B(n) ; vectorul-coloană a termenului liber

METODE NUMERICE

tablou real X(n) ; vectorul-soluţie la iteraţia anterioară tablou real Xn(n) ; vectorul-soluţie curentă

; introducere date de intrare şi iniţializare a vectorului-soluţie

citeşte n, er pentru i = 1, n repetă citeşte b(i) pentru j = 1, n repetă citeşte a(i,j) x(i) = 0

; iteraţii

repetă err = 0 pentru i = 1, n repetă s = b(i) pentru j = 1, n repetă s = s - a(i, j)·x(j) s = s + a(i,i)·x(i) xn(i) = s / a(i,i) s = ( ) ( )ixixn − dacă err < s atunci err = s scrie err pentru i = 1, n repetă x(i) = xn(i) scrie x(i) k = k + 1 până când err < er stop 3. Exemplu de aplicare

Să se rezolve, prin metoda iterativă a lui Jacobi, sistemul

⎪⎩

⎪⎨

⎧

−=⋅−⋅+=−⋅+⋅−=⋅+−⋅

.11x6x3x;4xx5x3;11x2xx7

321

321

321

Vom începe prin a constata că matricea sistemului îndeplineşte condiţia (7):

- pentru prima linie, 173

721

aaa

11

1312 <=+

=+

;

46

CAPITOLUL 7 – Rezolvarea sistemelor de ecuaţii liniare prin metoda Jacobi

- pentru a doua linie, 154

513

aaa

22

2321 <=+

=+

;

- pentru a treia linie, 164

631

aaa

33

3231 <=+

=+

şi, prin urmare, sistemul poate fi rezolvat prin metoda lui Jacobi. Pe baza relaţiei (6), se pot scrie formulele de recurenţă necesare pentru

generarea şirului de vectori-aproximaţii succesive ale soluţiei:

( )( )( )

( ) ( ) ( )( )⎪⎪⎪

⎩

⎪⎪⎪

⎨

⎧

⋅−−−⋅−=

+⋅+⋅=

⋅−+⋅=

+

+

+

.x3x1161x

;xx3451x

;x2x1171x

k2

k1

1k3

k3

)k(1

)1k(2

)k(3

)k(2

)1k(1

Utilizându-se, în continuare, un program de calcul bazat pe algoritmul metodei Jacobi, adoptându-se eroarea maximă admisibilă er = 0.001 şi valoarea

iniţială , la rezolvarea sistemului de mai sus se obţin, succesiv,

valorile din tabelul 1.

( )

⎟⎟⎟

⎠

⎞

⎜⎜⎜

⎝

⎛=

000

x 0

Tab.1

Nr. iteraţiei

x1 x2 x3 err

1 1.571 0.800 1.833 1.833 2 1.162 2.109 2.495 1.309 3 1.160 1.996 3.082 0.586 4 0.976 2.112 3.025 0.184 5 1.008 1.991 3.052 0.122 6 0.983 2.016 2.997 0.055 7 1.003 1.990 3.005 0.026 8 0.997 2.003 2.995 0.013 9 1.002 1.997 3.001 0.006

10 0.999 2.001 2.999 0.004 11 1.0004 1.9993 3.0005 0.0018 12 0.9997 2.0004 2.9998 0.0011 13 1.0001 1.9998 3.0001 0.0006

47

METODE NUMERICE

48

Ţinând cont că soluţiile exacte ale sistemului sunt x1 = 1, x2 = 2 şi x3 = 3, se poate observa, din tabel, modul cum aproximaţiile succesive converg către valorile corecte.

4. Chestiuni de verificare

Să se rezolve, prin metoda Jacobi, sistemele de ecuaţii de mai jos, după ce se va fi verificat, în prealabil, îndeplinirea criteriului de convergenţă; să se compare rezultatele cu cele obţinute, la rezolvarea aceloraşi sisteme, la lucrarea 6.

4.1. ⎪⎩

⎪⎨

⎧

−=−+=++−−=+−

.8z6yx2;3z2y4x;5z2y3x7

4.2.

⎪⎪⎩

⎪⎪⎨

⎧

−=−++−=+−+

−=++−=−+−

.1x4xxx;1x2x6x2x

;2xx3x7x;3x2xxx5

4321

4321

4321

4321

CAPITOLUL 8 – Rezolvarea sistemelor de ecuaţii liniare prin metoda Gauss - Seidel

Capitolul 8

REZOLVAREA SISTEMELOR DE ECUAŢII LINIARE PRIN METODA GAUSS - SEIDEL.

1. Principiul metodei

Metoda Gauss-Seidel este, ca şi metoda Jacobi, tot o metodă iterativă

de rezolvare a sistemelor de ecuaţii liniare; mai mult, se poate afirma că metoda Gauss-Seidel este, în fapt, o variantă îmbunătăţită a metodei Jacobi.

Spre deosebire de aceasta, care nu foloseşte în iteraţia curentă vectorul x(k), până când toate componentele sale vor fi fost determinate, metoda Gauss-Seidel foloseşte aceste componente, pe măsură ce ele sunt calculate. Astfel, la pasul (k+1) se folosesc în calculul necunoscutei xi

(k+1) valorile necunoscutelor x1

(k+1), x2(k+1), . . . xi-1

(k+1), deja determinate în iteraţia curentă.; formula

⎟⎟⎟

⎠

⎞

⎜⎜⎜

⎝

⎛⋅−= ∑

≠=

+ n

ij1j

)k(jiji

ii

)1k(i xab

a1x , i = 1, 2, . . . n. (1)

devine, în acest caz

( ) ( )⎟⎟⎠

⎞⎜⎜⎝

⎛⋅−⋅−= ∑∑

+=

−

=

++ n

1ij

kjij

1i

1j

1kjiji

ii

)1k(i xaxab

a1x , i = 1, 2, . . . n. (2)

49

METODE NUMERICE

Condiţia de convergenţă

,1aan

ij1j ii

ij <∑≠=

i = 1, 2, . . . n, (3)

şi criteriul de oprire a căutării soluţiei

( )( ) erxxmaxerr )k(i

1ki

ni1<−= +

≤≤, (4)

sunt aceleaşi ca şi pentru metoda lui Jacobi dar, la acelaşi nivel de precizie, metoda de faţă este sensibil mai rapidă. 2. Algoritmul metodei program Metoda Gauss-Seidel

real er ; eroarea maximă admisibilă real err ; eroarea la pasul curent

real s ; suma parţială întreg n ; dimensiunea sistemului întreg i, j, k ; contori tablou real A(n,n) ; matricea sistemului tablou real B(n) ; vectorul-coloană a termenului liber tablou real X(n) ; vectorul-soluţie

; introducere date de intrare şi iniţializare a vectorului-soluţie

citeşte n, er pentru i = 1, n repetă citeşte b(i) pentru j = 1, n repetă citeşte a(i,j) x(i) = 0 ; iteraţii

repetă err = 0 pentru i = 1, n repetă s = b(i) pentru j = 1, n repetă s = s - a(i, j)·x(j) 50

CAPITOLUL 8 – Rezolvarea sistemelor de ecuaţii liniare prin metoda Gauss - Seidel

s = (s + a(i,i)·x(i)) / a(i,i) dacă err < )i(xs − atunci err = )i(xs − x(i) = s scrie x(i) scrie err k = k+1 până când err < er stop 3. Exemplu de aplicare

Să se rezolve, prin metoda iterativă Gauss - Seidel, sistemul

⎪⎩

⎪⎨

⎧

−=⋅−⋅+=−⋅+⋅−=⋅+−⋅

.11x6x3x;4xx5x3;11x2xx7

321

321

321

S-a considerat acelaşi sistem ca în cazul metodei lui Jacobi, pentru a

putea evidenţia mai bine ceea ce metoda Gauss-Seidel are specific. Astfel, relaţiile utilizate, de această dată, pentru generarea şirului de vectori - aproximaţii succesive ale soluţiei, devin

( )( )( )

( ) ( ) ( )( )⎪⎪⎪

⎩

⎪⎪⎪

⎨

⎧

⋅−−−⋅−=

+⋅+⋅=

⋅−+⋅=

+++

++

+

.x3x1161x

;xx3451x

;x2x1171x

1k2

1k1

1k3

k3

)1k(1

)1k(2

)k(3

)k(2

)1k(1

Utilizând un program de calcul, bazat pe algoritmul metodei Gauss-

Seidel, cu aceleaşi date iniţiale ca în exemplul prezentat la lucrarea 7, se obţin rezultatele din tabelul 1.

Tab. 1 Nr.

iteraţiei x1 x2 x3 err

1 1.571 1.743 2.967 2.967 2 0.973 1.977 2.984 0.599 3 1.001 1.997 2.999 0.028 4 0.9999 1.9997 2.9999 0.0022 5 1.0000 1.9999 2.9999 0.0002

51

METODE NUMERICE

52

Evident, şirul de aproximaţii succesive converge mult mai rapid către soluţia exactă - pentru aceeaşi eroare maximă admisibilă ca în cazul metodei lui Jacobi, sunt necesare doar 5 iteraţii, în loc de 13, pentru găsirea soluţiei aproximative a sistemului. 4. Chestiuni de verificare

Să se rezolve, prin metoda Gauss - Seidel, sistemele de ecuaţii de mai jos; să se compare rezultatele cu cele obţinute, la rezolvarea aceloraşi sisteme, la lucrările 6 şi 7.

4.1. ⎪⎩

⎪⎨

⎧

−=−+=++−−=+−

.8z6yx2;3z2y4x;5z2y3x7

4.2.

⎪⎪⎩

⎪⎪⎨

⎧

−=−++−=+−+

−=++−=−+−

.1x4xxx;1x2x6x2x

;2xx3x7x;3x2xxx5

4321

4321

4321

4321

CAPITOLUL 9 – Polinomul de interpolare al lui Lagrange

Capitolul 9

POLINOMUL DE INTERPOLARE AL LUI LAGRANGE.

1. Noţiuni teoretice

Fie funcţia , dată sub forma [ ] Rb,a:f →

xk x0 x1 ............... xn

f(xk) f(x0) f(x1) ............... f(xn) unde bxxxxa m210 =<<<<= K este o diviziune, nu neapărat echidistantă a intervalului [a, b]. Se aleg ca noduri de interpolare punctele diviziunii de mai sus şi se impune coincidenţa funcţiei f cu un polinom de aproximare algebric de gradul n în nodurile de interpolare, adică

( ) ( ) n,2,1,0k,xPxf knk K== . (1)

Metoda lui Lagrange propune determinarea polinomului de interpolare de forma

( ) ( ) ( ) xyxyxyxL nn1100n ( )ϕ⋅++ϕ⋅+ϕ⋅= K , (2)

unde s-au făcut notaţiile yk = f(xk), k = 0, 1, 2, . . . n iar φ0(x), φ1(x), . . . φn(x) sunt polinoame de gradul n care urmează a fi determinate din impunerea setului de condiţii (1). Rezultă, astfel, sistemul 53

METODE NUMERICE

( ) ( ) ( )( ) ( ) ( )

( ) ( ) ( )⎪⎪⎩

⎪⎪⎨

⎧

=ϕ⋅++ϕ⋅+ϕ⋅

=ϕ⋅++ϕ⋅+ϕ⋅=ϕ⋅++ϕ⋅+ϕ⋅

.yxyxyxy

;yxyxyxy;yxyxyxy

nnnnn11n00

11nn111100

00nn011000

K

M

K

K

(3)

O soluţie a sistemului este, în mod evident, cea pentru care

( ) ( ) ( ) ( )( ) ( ) ( ) ( )

( ) ( ) ( ) ( )

( ) ( ) ( ) ( )⎪⎪⎪

⎩

⎪⎪⎪

⎨

⎧

=ϕ=ϕ=ϕ=ϕ

=ϕ=ϕ=ϕ=ϕ

=ϕ=ϕ=ϕ=ϕ=ϕ=ϕ=ϕ=ϕ

.1x;0x;0x;0x

;0x;1x;0x;0x

;0x;0x;1x;0x;0x;0x;0x;1x

nnnkn1n0

knkkk1k0

1n1k1110

0n0k0100

KK

M

KK

M

KK

KK

(4)

Rezultă că polinomul φk, spre exemplu, se anulează în n puncte şi ia valoarea 1 într-un singur punct - xk. De aici, se poate adopta o formă simplă pentru φk şi anume

( ) ( )( ) ( )( ) ( )n1k1k10kk xxxxxxxxxxax −−−−−⋅=ϕ +− KK . (5)

Constanta ak se determină impunând condiţia ca în punctul xk

polinomul φk să ia valoarea 1; se obţine

( )( ) ( )( ) ( )nk1kk1kk1k0kk xxxxxxxxxx

1a−−−−−

=+− KK

, (6)

de unde, rezultă

( ) ( )( ) ( )( ) ( )( )( ) ( )( ) ( )nk1kk1kk1k0k

n1k1k10k xxxxxxxxxx

xxxxxxxxxxx

−−−−−−−−−−

=ϕ+−

+−

KK

KK. (7)

În fine, se obţine forma polinomului de interpolare Lagrange de gradul n:

( ) ( )( ) ( )( ) ( )( )( ) ( )( ) ( )∑

= +−

+−

−−−−−−−−−−

⋅=n

0k nk1kk1kk1k0k

n1k1k10kn xxxxxxxxxx

xxxxxxxxxxyxL

KK

KK (8)

Dacă se face notaţia

( ) ( )( )∏

≠= −

−=

n

ki0i ik

ik xx

xxxa , k = 0, 1, 2, . . . n, (9)

54

CAPITOLUL 9 – Polinomul de interpolare al lui Lagrange

expresia polinomului Lagrange se poate scrie mai simplu sub forma

( ) ( )∑=

⋅=n

0kkkn xayxL . (10)

Polinomul de interpolare Lagrange coincide cu funcţia f în nodurile de interpolare, dar între nodurile de interpolare poate să se abată destul de mult de la valorile acesteia. Pentru a putea evalua ordinul de mărime al acestor abateri, se consideră funcţia abatere

( ) ( ) ( )xLxfx nn −=ε ; (11)

pentru a determina forma analitică a acestei funcţii, se consideră, mai întâi, o funcţie auxiliară h, de forma

( ) ( ) ( ) ( )( ) ( )n10n xxxxxxKxLxfxh −−−⋅−−= K , (12)

în care K este o constantă care urmează a fi determinată. Se observă că h(xk) = 0 pentru k = 0, 1, 2, . . . n, deci funcţia h are

(n+1) zerouri, oricare ar fi valoarea lui K; această valoare se poate determina impunând ca funcţia h să se anuleze şi într-un punct ( ) kn,0 x,xx ≠α∈α ;se obţine, de aici

( ) ( )

( )( ) ( )n10

n

xxxLf

K−α−α−α

α−α=

K. (13)

Dacă se ţine cont şi de rădăcina α, se observă că funcţia h are (n+2) zerouri pe intervalul [x0, xn] şi, în plus, satisface condiţiile teoremei lui Rolle pe fiecare dintre cele (n+1) subintervale care se formează folosind nodurile interne x1, x2, . . . xn-1, α.

Rezultă că prima derivată a funcţiei h, h’, are cel puţin (n+1) zerouri pe [x0, xn]; raţionând asemănător, se observă că a doua derivată, h”, are cel puţin n zerouri pe [x0, xn] ş.a.m.d. până la derivata de ordin n, care trebuie să aibă cel puţin un punct ( n0 x,x )∈ξ în care să se anuleze, ceea ce conduce la relaţia

( ) ( )( )!1n

fK1n

+ξ

=+

. (14)

Egalând cele două valori obţinute pentru K, conform relaţiilor (13) şi (14), rezultă

( ) ( ) ( )( ) ( )( ) ( )( ) ( .xxx

!1nfLf n10

1n

nn −α−α−α⋅+

ξ=α−α=αε )

+L (15)

Deoarece valoarea lui α este oarecare în intervalul [x0, xn], kx≠α , (dar pentru α = xk, evident, εn(α) = 0) pentru funcţia eroare rezultă expresia: 55

METODE NUMERICE

( )( ) ( )( ) ( )( ) ( n10

1n

n xxxxxx!1n

fx −−−⋅+

ξ=ε

+K ), (16)

unde ξ este un punct oarecare din intervalul (x0, xn). Dacă funcţia f are derivata de ordinul (n+1) mărginită, adică dacă

există M > 0 astfel încât pentru orice ( )n0 x,xx∈ să fie satisfăcută relaţia ( ) ( ) Mxf 1n ≤+ , atunci, conform (16), eroarea de aproximare prin polinomul de

interpolare al lui Lagrange este, la rândul său mărginită:

( ) ( ) n10n xxxxxx!1n

Mx −⋅⋅−⋅−⋅+

≤ε K . (17)

2. Algoritmul metodei program Polinomul Lagrange

întreg n ; gradul polinomului întreg j, k ; contori

tablou real x(n), y(n) ; nodurile de interpolare real x0 ;abscisa punctului de interpolare real y0 ; valoarea interpolată a funcţiei real p ; variabilă auxiliară

; introducerea datelor de intrare

citeşte n citeşte x0 pentru k = 0, n repetă citeşte x(k) citeşte y(k)

; construcţia polinomului Lagrange

y0 = 0 pentru k = 0, n repetă pentru j = 0, n repetă dacă j≠ k atunci

( )( ) ( )jxkx

jxxpp−−

⋅=0

y0 = y0 + y(k)·p scrie y0

stop

56

CAPITOLUL 9 – Polinomul de interpolare al lui Lagrange

3. Exemplu de aplicare

Fie următoarea funcţie, dată tabelar:

xk -2 0 2 4 6 8 10 12 f(xk) 0.25 1 4 16 64 256 1024 4096

Să se calculeze, prin interpolare cu un polinom Lagrange, valoarea funcţiei f în punctul x = 5.

Trebuie observat, mai întâi, că punctele din tabel au fost extrase din mulţimea de valori ale funcţiei f(x) = 2x; astfel, atunci când se va calcula valoarea aproximată prin interpolare, se va putea aprecia şi calitatea aproximării, valoarea exactă fiind cunoscută (25 = 32).

Dacă se utilizează pentru interpolare doar cele două puncte din mijlocul tabelei de valori (x0 = 4 şi x1 = 6), polinomul de interpolare este de gradul întâi (interpolare liniară):

( ) 80x24464x64

646x16xL1 −⋅=

−−

⋅+−−

⋅= ;

rezultă, de aici, L1(5) = 40. Eroarea relativă de aproximare este, în acest

caz ( ) %2510032

324051 =⋅−

=ε - evident, o valoare inacceptabil de mare.

Dacă se utilizează patru noduri de interpolare (x0 = 2, x1 = 4, x2 = 6 şi x3 = 8) polinomul Lagrange de gradul trei rezultă

( ) ( )( )( )

( )( )( )( )( )( )( )( )( )

( )( )( )( )( )( )

( )( )(( )( )( )

) ,6848286x4x2x256

8646268x4x2x64

8464248x6x2x16

8262428x6x4x4xL3

−−−−−−

⋅+−−−−−−

⋅+

+−−−−−−

⋅+−−−−−−

⋅=

respectiv ( ) .80x78x245x

49xL 23

3 −⋅+⋅−⋅= Se obţine, astfel,

L3(5) =28.750, cu o eroare relativă de -10.156 % faţă de valoarea exactă.

În cazul utilizării a şase noduri de interpolare (de la x0 = 0 la x6 = 10), cu ajutorul unui program de calcul special destinat pentru construcţia polinomului de interpolare Lagrange, se obţine L5(5) = 33.496 cu ε5(5) = 4.675 %.

Dacă se utilizează ca noduri de interpolare a tuturor celor opt puncte din tabel, cu acelaşi program se obţine L7(5) = 31.271, respectiv ε7(5) = -2.278 %.

57

METODE NUMERICE

58

Evident, în exemplul dat nodurile de interpolare au fost alese relativ rare, iar punctul în care s-a aproximat cu polinomul Lagrange valoarea funcţiei este amplasat cât mai defavorabil posibil (la mijlocul distanţei dintre două noduri succesive); cu toate acestea, se poate remarca, în cazul utilizării mai multor noduri de interpolare, o precizie relativ buna de aproximare. 4. Chestiuni de verificare 4.1. Fie următoarea funcţie, dată tabelar:

xk -2 0 2 6 8 10 12 f(xk) 0.25 1 4 64 256 1024 4096

Să se calculeze, prin interpolare cu un polinoame Lagrange, de grad 3, 4 şi 5, valoarea funcţiei f în punctele x0 = 1, 2, 3, 4, 5. 4.2. Fie următoarea funcţie, dată tabelar:

xk -2 0 2 4 6 f(xk) 0.11 1 9 81 729

Să se calculeze, prin interpolare cu un polinoame Lagrange, de grad 2, 3, şi 4, valoarea funcţiei f în punctele x0 = -1, 1, 3, 5.

CAPITOLUL 10 – Polinomul de interpolare al lui Newton

Capitolul 10

POLINOMUL DE INTERPOLARE AL LUI NEWTON.

1. Noţiuni teoretice

Fie funcţia f tabelata în puncte echidistante şi fie h pasul tabelei de

valori (xk+1 - xk = h, pentru orice k = 0, 1, 2, . . . n-1). Se numeşte diferenţă finită la dreapta (progresive), de ordinul întâi, pentru funcţia f expresia

( ) ( ) ( )xfhxfxf −+=Δ . (1)

Diferenţele finite de ordin superior se definesc în mod recurent:

( ) ( )( ) ( )( )xfxfxf 1n1nn ΔΔ=ΔΔ=Δ −− (2)

k xk yk Δyk Δ2yk Δ3yk Δ4yk 0 x0 y0 Δy0 1 x1 y1 Δ2y0 Δy1 Δ3y0 2 x2 y2 Δ2y1 Δ4y0 Δy2 Δ3y1 3 x3 y3 Δ2y2 Δy3 4 x4 y4

59

METODE NUMERICE

În aplicaţii, calculul diferenţelor finite se face mai comod, prin simple operaţii de scădere, pornind de la tabela de valori a funcţiei f dispusă vertical; spre exemplu, pentru o tabelă în cinci puncte, reprezentarea este cea de mai sus.

Diferenţele finite, de un anumit ordin, se calculează scăzând diferenţele finite de ordin inferior între care se găsesc, situate la stânga.

În cazul în care se doreşte aproximarea unei funcţii prin interpolare, dacă nodurile de interpolare sunt echidistante se poate utiliza polinomul de interpolare al lui Newton, construit pe baza diferenţelor finite. În cele ce urmează se vor utiliza diferenţele finite la dreapta (regresive).

Se pune deci problema determinării unui polinom algebric de grad mai mic sau egal cu n, care să satisfacă condiţiile de interpolare P(xk) = yk, k = 0, 1, 2, . . . n. Se construieşte un polinom de forma

( ) ( ) ( )( )( )( ) ( .xxxxxxc

xxxxcxxccxN

1n10n

102010dn

−−−−+++−−+−+=

K

K

) (3)

Impunând condiţiile de interpolare polinomului de mai sus, rezultă sistemul

( )( ) ( )( )

( ) ( )( ) ( ) ( )⎪⎪⎪

⎩

⎪⎪⎪

⎨

⎧

=−−+−−+−+

=−−+−+=−+

=

− .yxxxxcxxxxcxxcc

;yxxxxcxxcc;yxxcc

;yc

n1nn0nn1n0n20n10

2120220210

10110

00

K

M

(4) Dacă se ţine seama că nodurile de interpolare sunt echidistante, cu

pasul h, sistemul (4) se poate scrie mai simplu ca

( )⎪⎪⎪

⎩

⎪⎪⎪

⎨

⎧

=⋅⋅⋅⋅⋅++⋅⋅−+⋅⋅+

=⋅⋅⋅+⋅⋅+

=⋅+=

.ychn21ch1nnchnc

;ych21ch2c

;ychc;yc

nnn

22

10

222

10

110

00

KK

M

(5)

Sistemul (5) este un sistem de (n+1) ecuaţii cu (n+1) necunoscute - c0, c1, c2, . . . cn - compatibil determinat. Pentru rezolvare, se înlocuieşte valoarea lui c0 din prima ecuaţie în a doua şi se află c1, apoi se înlocuiesc valorile lui c0 şi c1 în a treia ecuaţie şi se determină c2 ş.a.m.d. Rezultă:

c0 = y0; hy

c 01

Δ= ;

n0

n

n20

2

2h!ny

c;h!2y

c⋅

Δ=

⋅

Δ= K . (6)

60

CAPITOLUL 10 – Polinomul de interpolare al lui Newton

Cu ajutorul acestor valori ale coeficienţilor se obţine formula explicită a polinomului de interpolare Newton de gradul n, construit cu diferenţe finite la dreapta:

( ) ( ) ( ) ( )

( )( ) ( ).xxxxxxh!ny

xxxxh!2y

xxh!1

yyxN

1n10n0

n

1020

2

00

0dn

−−−−⋅

Δ+

++−⋅−⋅⋅

Δ+−

⋅Δ

+=

K

K

(7)

În mod cu totul analog se poate construi si polinomul de interpolare Newton de gradul n, pe baza diferenţelor finite la stânga,

( ) ( ) ( ) ( )

( )( ) ( ).xxxxxxh!ny

xxxxh!2y

xxh!1

yyxN

11nnnn

n

1nn2n

2

nn

nsn

−−−⋅

∇+

++−⋅−⋅⋅

∇+−

⋅∇

+=

−

−

K

K

(8)

Polinomul lui Newton construit pe baza diferenţelor finite la dreapta

se recomandă a fi utilizat când se doreşte estimarea valorilor funcţiei y = f(x) în puncte mai apropiate de x0 decât de xn (situate în prima jumătate a tabelei de valori), în timp ce varianta construită cu diferenţe finite la stânga va fi folosită pentru aproximarea valorilor funcţiei în puncte din vecinătatea lui xn (din a doua jumătate a tabelei de valori).

Pentru a se pune în evidenţă eroarea aproximării prin cele două polinoame de interpolare Newton, se pot face aproximările

( ) ( )1n

01n

1n

hy

f+

++ Δ

≅ξ , respectiv ( ) ( )1n

n1n

1n

hy

f+

++ ∇

≅ξ . (9)

Din cele de mai sus rezultă expresiile funcţiilor eroare pentru cele

două polinoame Newton (cu diferenţe finite la dreapta, respectiv la stânga):

( )( )

( )( ) ( n101n0

1ndn xxxxxx

h!1ny

x −−−⋅+

Δ≅ε

+

+

K ) , (10)

respectiv

( )( )

( )( ) ( n101nn

1nsn xxxxxx

h!1ny

x −−−⋅+

∇≅ε

+

+

K ) . (11)

61

METODE NUMERICE

2. Algoritmul metodei program Polinomul Newton

întreg np ; gradul polinomului întreg n ; ordinul ultimului nod de

interpolare întreg i, k ; contori

tablou real x(n), y(n) ; nodurile de interpolare tablou real D(n,n) ; tabloul diferenţelor finite real x0 ;abscisa punctului de interpolare real y0 ; valoarea interpolată a funcţiei real p ; variabilă auxiliară real h ; pasul reţelei de noduri

; introducerea datelor de intrare

citeşte n; citeşte np; citeşte x0; pentru k = 0, n repetă citeşte x(k) citeşte y(k) D(0,k) = y(k)

; calculul diferenţelor finite

pentru i = 1, n repetă pentru k = 0, n-i repetă D(i,k) = D(i-1,k+1) - D(i-1,k)

; construcţia polinomului Newton

h = x(1) - x(0) y0 = D(0,0) pentru i = 1,np repetă p = 1 pentru k = 0, i-1 repetă

( )( ) hk

kxxpp⋅+

−⋅=

10

y0 = y0 + D(i,0)·p scrie y0 stop

62

CAPITOLUL 10 – Polinomul de interpolare al lui Newton

3. Exemplu de aplicare Fie următoarea funcţie, dată tabelar:

xk -2 0 2 4 6 8 10 12 f(xk) 0.25 1 4 16 64 256 1024 4096

Să se calculeze, prin interpolare cu polinoame Newton, bazate pe diferenţe finite progresive - respectiv regresive - de diferite grade valoarea funcţiei f în punctele x = 1, 5, 9.

Se constată că şi de această dată punctele din tabelă sunt extrase din graficul funcţiei f(x) = 2x, cu scopul declarat de a putea evalua mărimea erorilor de aproximare; se ştie deci, apriori, că f(1) = 21 = 2, f(5) = 25 = 32 iar f(9) = 29 = 512.

Se construieşte tabelul diferenţelor finite, pe baza valorilor funcţiei:

k xk yk Δyk Δ2yk Δ3yk Δ4yk Δ5yk Δ6yk Δ7yk 0 -2 0.25 0.75 1 0 1 2.25 3 6.75 2 2 4 9 20.25 12 27 60.75 3 4 16 36 81 182.25 48 108 243 546.754 6 64 144 324 729 192 432 972 5 8 256 576 1296 768 1728 6 10 1024 2304 3072 7 12 4096

Şirul de valori din căsuţele marcate situate în partea superioară a tabelului reprezintă diferenţele finite la dreaptă (progresive) - Δy0, Δ2y0 . . . Δ7y0 - pe baza cărora se poate construi cu relaţia (4.41) polinomul Newton de grad mai mic sau egal cu 7.

Spre exemplu, polinomul de gradul 2 rezultă de forma

( ) ( ) ( ) 1x9375.0x28125.0x2x2225.22x

275.025.0xN 2

2d2 +⋅+⋅=⋅+

⋅+++= .

63

METODE NUMERICE

Şirul de valori din căsuţele marcate situate în partea inferioară a tabelului reprezintă diferenţele finite la stânga (regresive) - ∇yn, 2yn, . . .

7yn = Δ7y0 - pe baza cărora se poate construi cu relaţia (4.42) polinomul Newton de grad mai mic sau egal cu 7. Astfel polinomul de gradul 3 rezultă de forma

∇∇

( ) ( ) ( )( )

( )( )( ) .14336x11904x792x368x10x12x86

1728

10x12x22

230412x2

30724096xN

23

2s3

−⋅+⋅−⋅=−−−⋅

+

+−−⋅

+−⋅+=

Prin intermediul a două programe de calcul special concepute, se construiesc polinoamele Newton cu diferenţe finite progresive, respectiv regresive, cu gradul n modificabil, succesiv, între 2 şi 7 (maximum posibil) şi se calculează valorile aproximate ale funcţiei în punctele cerute. Rezultă valorile din tabelele 1, respectiv 2.

Tab. 1

n ( )1Ndn ( )1nε

[%] ( )5Nd

n ( )5nε [%]

( )9Ndn ( )9nε

[%] 2 2.219 10.95 12.719 -60.25 32.219 -93.71 3 1.719 -14.05 27.484 -14.11 129.672 -74.67 4 2.271 13.55 33.021 3.19 312.396 -38.98 5 1.560 -22.00 31.360 -2.00 476.849 -6.86 6 2.805 40.25 32.606 1.89 517.962 1.16 7 0.403 -79.85 31.271 -2.28 509.152 -0.56

Tab .2

n ( )1Nsn ( )5Ns

n ( )9Nsn

2 15712 3424 352 3 -9236 -356 460 4 2456 -1.625 490.375 5 -172.859 24.953 501.766 6 -8.407 29.937 506.749 7 0.403 31.271 509.152

După cum se poate observa din rezultatele obţinute, în anumite cazuri eroarea relativă este suficient de mică pentru a putea considera interpolarea rezonabil de precisă; sunt însă şi situaţii în care erorile sunt inacceptabile.

64

CAPITOLUL 10 – Polinomul de interpolare al lui Newton

65

4. Chestiuni de verificare 4.1. Fie următoarea funcţie, dată tabelar:

xk -2 0 2 4 6 f(xk) 0.11 1 9 81 729

Să se calculeze, prin interpolare cu un polinoame Newton, de grad 2, 3, şi 4, valoarea funcţiei f în punctele x0 = -1, 1, 3, 5. 4.2. Fie următoarea funcţie, dată tabelar:

xk 0 2 4 6 8 f(xk) 0 1.41 2 2.45 2.83

Să se calculeze, prin interpolare cu un polinoame Newton, de grad 2, 3, şi 4, valoarea funcţiei f în punctele x0 = 1, 3, 5 şi 7.

CAPITOLUL 11 – Aproximarea cu abatere medie pătratică minimă

Capitolul 11

APROXIMAREA CU ABATERE MEDIE PĂTRATICĂ MINMĂ.

1. Noţiuni teoretice

Dacă funcţia de aproximat este cunoscută tabelar,

xk x0 x1 ............... xm

yk = f(xk) y0 y1 ............... ym

atunci distanţa dintre aceasta şi polinomul de aproximare cu abatere medie pătratică minimă Pn se calculează ca

( ) ( )[ ]∑=

−⋅+

=m

0i

2iinn yxP

1m1f,Pd . (1)

Minimizarea abaterii medii pătratice este echivalentă cu determinarea

valorii minime a sumei de sub radical; rezultă, că se caută minimul funcţiei , RR:S 1n →+

. (2) ( ) ( )[ ]∑=

−=m

0i

2iinn10 yxPa,a,aS K

67

METODE NUMERICE

Determinarea coeficienţilor polinomului Pn se face impunând condiţia ca derivatele parţiale ale funcţiei S să se anuleze. Cum

(∑=

=⋅−⋅++⋅+=∂∂ m

0i

kii

nini10

k0xyxaxaa2

aS

K ) , k = 0, 1, 2, . . . n, (3)

rezultă sistemul de ecuaţii

, k = 0, 1, 2, . . . n. (4) ( ) ∑∑==

⋅=⋅⋅++⋅+m

0i

kii

m

0i

ki

nini10 xyxxaxaa K

Introducând notaţiile

∑=

=m

0i

kik xu , respectiv , (5) ∑

=⋅=

m

0i

kiik xyv

sistemul (4) poate fi scris matriceal sub forma

. (6)

⎟⎟⎟⎟⎟

⎠

⎞

⎜⎜⎜⎜⎜

⎝

⎛

=

⎟⎟⎟⎟⎟

⎠

⎞

⎜⎜⎜⎜⎜

⎝

⎛

⋅

⎟⎟⎟⎟⎟

⎠

⎞

⎜⎜⎜⎜⎜

⎝

⎛

++

+

n

2

1

n

1

0

nn1nn

1n21

n10

v

vv

a

aa

uuu

uuuuuu

MM

K

M

K

K

În cazul particular în care polinomul de aproximarea cu abatere medie

pătratică minimă este de gradul întâi, p1(x) = a0 + a1·x, acesta poartă numele de „dreaptă de regresie”, coeficienţii acesteia calculându-se cu formulele simple

2

120

11020

uuu

vuvua

−⋅

⋅−⋅= , respectiv

2120

01101

uuu

vuvua

−⋅

⋅−⋅= , (7)

obţinute prin rezolvarea sistemului (6). 2. Algoritmul metodei program ampm întreg n ; numărul punctelor din tabel întreg k ; contor

tablou real x(n), y(n) ; punctele din tabelul de valori real u0, u1, u2, u3, u4 ; sume de puteri ale absciselor

68

CAPITOLUL 11 – Aproximarea cu abatere medie pătratică minimă

real v0, v1, v2 ; sume mixte real a0, a1, a2 ; coeficienţii polinomului de aproximare

; introducerea datelor de intrare citeşte n citeşte x0 pentru k = 0, n repetă citeşte x(k) citeşte y(k)

; calculul sumelor de puteri şi determinarea funcţiei de aproximare

u0 = n + 1; u1 = 0; u2 = 0; u3 = 0; u4 = 0; v0 = 0; v1 = 0; v2 = 0 pentru k = 0, n repetă u1 = u1 + x(k) u2 = u2 + x(k)2 u3 = u3 + x(k)3

u4 = u4 + x(k)4 v0 = v0 + y(k) v1 = v1 + x(k)·y(k) v2 = v2 + x(k)2·y(k)

a0, a1, a2 = Metoda Gauss(u0, u1, u2, u3, u4, v0, v1, v2) scrie a0, a1, a2 stop 3. Exemplu de aplicare

Să se găsească polinomul de gradul doi P2(x) = a0 + a1·x + a2·x2 care aproximează cu abatere medie pătratică minimă funcţia dată prin tabelul de mai jos:

xk 0 1 2 3 4 yk = f(xk) 1 2 4 8 16

Se calculează cu (4.99) sumele ; ;

∑=

==4

0i

0i0 5xu

∑=

=4

0

4i 354xu ; 0

∑=

==4

0i

1i1 10xu

∑=

=⋅=4

0ii

0i 31yxv ;∑

===

4

0i

2i2 30xu ;

∑=

=⋅=4

0iii1 yx

∑=

==4

0i

3i3 ;100xu

98v ; ∑=

⋅=4

0ii

2i2 yx

=i

4

= 346v .

69

METODE NUMERICE

70

Rezultă sistemul , care, rezolvat

printr-un program de calcul bazat pe metoda Gauss, conduce la soluţiile: a0 = 1.2857; a1 = -0.9714; a2 = 1.1429. S-a obţinut, astfel,

⎟⎟⎟

⎠

⎞

⎜⎜⎜

⎝

⎛=

⎟⎟⎟

⎠

⎞

⎜⎜⎜

⎝

⎛⋅

⎟⎟⎟

⎠

⎞

⎜⎜⎜

⎝

⎛

3469831

aaa

35410030100301030105

2

1

0

( ) 22 x1429.1x9714.02857.1xP ⋅+⋅−= .

Pentru a putea evalua precizia de aproximare, se calculează valorile lui P2 în abscisele punctelor din tabel şi erorile faţa de funcţia dată în aceste puncte; rezultă:

P2(0) = 1.2857; ε2(0) = 0.2857;

P2(1) = 1.4572; ε2(1) = -0.5428;

P2(2) = 3.9145; ε2(2) = -0.0855;

P2(3) = 8.6576; ε2(3) = -0.3424;

P2(4) = 15.6865; ε2(4) = -0.3135.

4. Chestiuni de verificare

Să se determine funcţiile de gradul al doilea care aproximează cu abatere medie pătratică minimă următoarele funcţii date tabelar:

4.1. xk -2 0 2 4 6

f(xk) 0.11 1 9 81 729 4.2. xk 0 2 4 6 8

f(xk) 0 1.41 2 2.45 2.83

CAPTOLUL 12 – Derivarea cu ajutorul polinomului de interpolare Lagrange

Capitolul 12

DERIVAREA CU AJUTORUL POLINOMULUI DE INTERPOLARE LAGRANGE.

1. Noţiuni teoretice

În cazul unei funcţii dată tabelar prin (n+1) puncte, se poate construi

polinomul de interpolare Lagrange şi se poate scrie

( ) ( ) ( ) ( )( )

( ) ( )[ ]xf!1n

xxxxxayxf 1n

n

0k

n0kk ξ⋅

+−−

+⋅= +

=∑

K. (1)

Prin derivare în raport cu x în relaţia (1), în condiţiile în care

x = xj, oricare ar fi j cuprins între 0 şi n rezultă

( ) ( )( ) ( )[ ]( ) (∏∑

≠=

)+

=−⋅

+

ξ+′⋅=′

n

jk0k

kjj

1nn

0kjkkj xx

!1nxf

xayxf . (2)

Relaţia (2) poartă numele de formulă de aproximare pentru f’(xj) în

(n+1) puncte, deoarece utilizează o combinaţie liniară de (n+1) valori ale funcţiei f.

71

În general, utilizarea unui număr mai mare de puncte în ecuaţia (5.11) va genera o precizie mai mare, dar numărul mare al evaluărilor funcţiei şi creşterea erorilor de rotunjire limitează acest lucru.

METODE NUMERICE

Se pot deduce, în continuare, câteva formule utile de derivare în trei puncte; dacă se calculează derivatele termenilor ak se obţine

( ) ( )( ) ( ) ( )( ) ( )

( )( ) ( ) ( )[ ] ( ),xxxf61xf

xxxxxxx2

xfxxxxxxx2

xfxxxx

xxx2xf

2

ji0i

ijj21202

10j

12101

20j0

2010

21jj

∏≠=

−⋅ξ′′′⋅+⋅⎥⎦

⎤⎢⎣

⎡−−

−−⋅+

+⋅⎥⎦

⎤⎢⎣

⎡−−

−−⋅+⋅⎥

⎦

⎤⎢⎣

⎡−−

−−⋅=′

(3)

pentru j = 0, 1, 2.

În cazul nodurilor egal depărtate, adică x1 = x0 + h şi x2 = x0 + 2·h, (4)

relaţia (3) devine, succesiv, pentru j = 0, 1, respectiv, 2:

( ) ( ) ( ) ( ) ( 0

2

0000 f3

hh2xf21hxf2xf

23

h1xf ξ′′′⋅+⎥⎦

⎤⎢⎣⎡ ⋅+⋅−+⋅+⋅−⋅=′ ) , (5)

( ) ( ) ( ) ( )1

2

001 f6

hh2xf21xf

21

h1xf ξ′′′⋅−⎥⎦

⎤⎢⎣⎡ ⋅+⋅+⋅−⋅=′ , (6)

( ) ( ) ( ) ( ) ( 2

2

0002 f3

hh2xf23hxf2xf

21

h1xf ξ′′′⋅+⎥⎦

⎤⎢⎣⎡ ⋅+⋅++⋅−⋅⋅=′ ). (7)

Dacă în ultimele două relaţii se fac substituţiile x0 + h = x0, respectiv,

x0 + 2·h = x0, după care, pentru generalizare, se înlocuieşte x0 cu xk (inclusiv în (5), rezultă formulele în trei puncte pentru calculul primei derivate a unei funcţii, într-un punct oarecare xk:

( ) ( ) ( ) ( )[ ]h2xfhxf4xf3h2

1xf kkkk ⋅+−+⋅+⋅−⋅⋅

≅′ , (8)

( ) ( ) ([ hxfhxfh2

1xf kkk ++−−⋅⋅

≅′ )], (9)

( ) ( ) ( ) ([ ]kkkk xf3hxf4h2xfh2

1xf ⋅+−⋅−⋅−⋅⋅

≅′ ) . (10)

72

CAPTOLUL 12 – Derivarea cu ajutorul polinomului de interpolare Lagrange

Cu totul asemănător se pot deduce formule pentru calculul primei derivate în mai mult de trei puncte - spre exemplu, mai jos este prezentată una dintre formulele în cinci puncte:

( ) ( ) ( ) ( ) ([ ]h2xfhxf8hxf8h2xfh12

1xf kkkkk ⋅+−+⋅+−⋅−⋅−⋅⋅

≅′ ) . (11)

2. Algoritmul metodei Program Derivare Lagrange

real x0 ; nodul iniţial real h ; pasul reţelei real x ; abscisa punctului curent întreg n ; dimensiunea tablourilor întreg i ; contor tablou real y(n) ; tabloul valorilor funcţiei tablou real dy(n) ; tabloul valorilor derivatei

citeşte h, x0, n pentru i = 0, n repetă citeşte y(i)

dy(0) = ( ) ( ) ( )h

yyy⋅

−⋅+⋅−2

21403 ; formula (8)

dy(n) = ( ) ( ) ( )h

nynyny⋅

⋅+−⋅−−2

3142 ; formula (10)

pentru i = 1, n -1 repetă

dy(i) = ( ) ( )h

iyiy⋅

−−+2

11 ; formula (9)

pentru i = 0, n repetă x = x0 + i·h scrie x, dy(i) stop 3. Exemplu de aplicare Fie următorul tabel, xk 0 1 2 3 4

f(xk) 1 2 4 8 16 conţinând cinci puncte din graficul funcţiei f(x) = 2x. Să se calculeze valoarea primei derivate a funcţiei f în punctul x2 = 2. 73

METODE NUMERICE

74

Evident, pentru rezolvarea problemei se poate utiliza oricare dintre

formulele în trei sau în cinci puncte; pentru a putea analiza oportunitatea utilizării uneia sau a alteia dintre acestea, în continuare se calculează f’(2), succesiv, cu formulele (8) . . . (11); se obţin, astfel, următoarele rezultate:

( ) ( ) 2168443212f =−⋅+⋅−⋅≅′ ;

( ) ( ) 382212f =+−⋅≅′ ;

( ) ( ) 5.243241212f =⋅+⋅−⋅≅′ ;

( ) ( ) 75.21688281212f =−⋅+⋅−⋅≅′ .

Cum valoarea exactă a derivatei pentru x=2 este , se poate remarca faptul că formulele prezentate dau, în

general rezultate satisfăcătoare, iar în cazul utilizării a cinci puncte în loc de trei aduce o îmbunătăţire simţitoare a preciziei de aproximare.

( ) 772.22ln42f ≅⋅=′

4. Chestiuni de verificare

4.1. Fie următoarea funcţie, dată tabelar:

xk -2 0 2 4 6 f(xk) 0.11 1 9 81 729

Să se calculeze, cu ajutorul polinomului de interpolare Lagrange, derivatele funcţiei f în punctele -2, 0, 2, 4, 6. 4.2. Fie următoarea funcţie, dată tabelar:

xk 0 2 4 6 8 f(xk) 0 1.41 2 2.45 2.83

Să se calculeze, cu ajutorul polinomului de interpolare Lagrange, derivatele funcţiei f în punctele 0, 2, 4, 6, 8.

CAPITOLUL 13 – Metode numerice pentru calculul integralei definite

Capitolul 13

METODE NUMERICE PENTRU CALCULUL INTEGRALEI DEFINITE.

1. Noţiuni teoretice

Fie o diviziune a = x0 < x1 < x2 < . . . < xn = b a intervalului [a, b], cu

noduri echidistante. Folosind proprietatea de aditivitate a integralei, se poate scrie

( ) ( )∫ ∑∫−

=

+==b

a

1n

0k

x

x1k

kdxxfdxxfI . (1)

Fie aproximarea funcţiei f pe intervalul [xk, xk+1] printr-un polinom de

interpolare Newton de gradul întâi

( ) ( )kk

k xxh!1

yyxf −⋅

Δ+≅ , (2)

unde h reprezintă pasul diviziunii iar yk = f(xk); înlocuind f din (1) cu valoarea aproximativă dată de (2) rezultă

( ) ( )∫ ∫+ +

⎥⎦⎤

⎢⎣⎡ −

⋅Δ

+≅= 1k

k

1k

k

x

x

x

x kk

kk dxxxh!1

yydxxfI . (3)

75

Prin schimbarea de variabilă thxx k =

− , din (3) se obţine

METODE NUMERICE

( ) ⎟⎠⎞

⎜⎝⎛ Δ

+=Δ⋅+⋅≅ ∫ 2yyhdtytyhI k

k1

0 kkk . (4)

Cum , rezultă că k1kk yyy −=Δ + ( 1kkk yy2hI ++≅ ) şi, în final,

( ) (∫ ∑ )−

=++⋅≅

b

a

1n

0k1kk yy

2hdxxf . (5)

În formă desfăşurată, relaţia (5) se scrie

( ) ( )n1n210b

ayy2y2y2y

2hdxxf +⋅++⋅+⋅+≅ −∫ K . (5’)

Formula astfel dedusă poartă numele de formula trapezelor,

deoarece, aşa după cum se poate observa din figura 1, se aproximează integrala I, reprezentată geometric prin aria suprafeţei cuprinse între curba y = f(x), axa Ox şi ordonatele ridicate în x = a şi x = b, prin suma ariilor suprafeţelor trapezelor de forma XkXk+1Mk+1Mk, construite cu ajutorul diviziunii considerate.

Pentru evaluarea erorii pe care o implică aplicarea formulei trapezelor la calcularea unei integrale, fie F o primitivă a funcţiei integrant f; conform formulei Leibniz - Newton, pentru intervalul [xk, xk+1] se poate scrie

( ) ( ) ( ) ( )[ ] ( ) ( ) ( ) ( )[ ]1kkk1k1kkx

xk xfxf2hxFxFxfxf

2hdxxff 1k

k+++ +−−=+−=ε ∫ +

(6)

Fig.1 - Metoda trapezelor - interpretare geometrică

Dacă se consideră dezvoltările în serie Taylor ale funcţiilor f şi F într-o vecinătate a punctului xk şi se ţine cont de faptul că F’ = f, rezultă: 76

CAPITOLUL 13 – Metode numerice pentru calculul integralei definite

( ) ( ) ( ) ( )k

2

kk1k f!2

hxf!1

hxfxf ξ′′⋅+′⋅+=+ ; (7)

( ) ( ) ( ) ( ) ( )k

3

k

2

kk1k f!3

hxf!2

hxf!1

hxFxF ξ′′⋅+′⋅+⋅+=+ . (8)

În relaţiile de mai sus ξk este un punct oarecare al intervalului

[xk, xk+1]. Prin înlocuirea relaţiilor (7) şi (8) în (6) se obţine

( ) ( )k

3

k f12hf ξ′′⋅−=ε . (9)

Însumând erorile de aproximare de pe fiecare subinterval (dând valori

lui k de la 0 la n-1 în (9)) şi având în vedere că n

abh −= , se poate calcula

eroarea totală generată de aplicarea formulei trapezelor:

( ) ( ) ( ) ( ) ( )∑ ∑−

=

−

=ξ′′⋅

−−=ξ′′⋅−=ε=ε

1n

0k

1n

0k2

3

k

2

kT fn12abf

12hff , (10)

unde punctul ξ a fost ales astfel încât ( ) ( )∑−

=ξ′′⋅=ξ′′

1n

0kkf

n1f .

În cazul în care derivata a doua a funcţiei admite un majorant, adică astfel încât ( ) 0M >∃ ( ) ( ) ( )b,ax,Mxf ∈∀≤′′ , rezultă, în final

( ) ( )2

3

T n12abMf −

≤ε . (11)

Analizând formula (11) se observă că singura cale de reducere a erorii de aproximare este creşterea numărului n (utilizarea unui număr cât mai mare de puncte intermediare).

Formula trapezelor a fost dedusă aproximând funcţia integrant printr-un polinom de interpolare Newton de gradul întâi, pe fiecare subinterval. Dacă gradul polinomului Newton de aproximare se ia doi, procedând analog ca mai sus se obţine, în cazul utilizării unei diviziuni pentru care n este par, formula de cuadratură Cavalieri - Simpson:

( ) [ ]∫ +++++++≅ −−b

a n1n2n3210 yy4y2y4y2y4y3hdxxf K . (12)

Eroarea pe care o implică aplicarea formulei Cavalieri - Simpson este:

( ) ( )4

5

T n180abMf −

≤ε . (13)

77

METODE NUMERICE

Comparând expresiile erorilor pentru cele două formule, (11) şi (13), se poate vedea că formula Cavalieri - Simpson conduce, pentru acelaşi număr de noduri, la o precizie mai bună decât formula trapezelor.

Prin utilizarea unor polinoame de interpolare de grad superior, se pot pune în evidenţă şi alte formule de cuadratură, cu precizie şi mai mare; cu toate acestea, având în vedere rezultatele foarte bune care se pot obţine prin aplicarea formulelor (5) şi (12), considerăm suficientă cunoaşterea acestora. 2. Algoritmii metodelor Program Metoda trapezelor

real a, b ; intervalul de integrare real h ; pasul diviziunii real I ; valoarea integralei calculate întreg n ; numărul de subintervale întreg j ; contor funcţie f(x:real):real f = ................. ; expresia analitică a funcţiei integrant returnează f

citeşte a, b, n

n

abh −=

I = 0 pentru j = 0, n repetă dacă j = 0 sau j = n atunci I = I + f(a + j·h) altfel I = I + 2·f(a + j·h)

I = I2h⋅

scrie I stop Program Metoda Cavalieri_Simpson

real a, b ; intervalul de integrare real h ; pasul diviziunii real I ; valoarea integralei calculate întreg n ; numărul de subintervale întreg j ; contor funcţie f(x:real):real f = ................. ; expresia analitică a funcţiei

78 integrant

CAPITOLUL 13 – Metode numerice pentru calculul integralei definite

returnează f

citeşte a, b, n

n

abh −=

I = 0 pentru j = 0, n repetă dacă j = 0 sau j = n atunci I = I + f(a + j·h) altfel dacă j mod 2 = 1 atunci

I = I + 4·f(a + j·h) altfel I = I + 2·f(a + j·h)

I = I3h⋅

scrie I stop 3. Exemplu de aplicare

Să se calculeze, cu ajutorul formulelor de cuadratură, integralele

şi ∫=4

0x

1 dx2I ∫=10

12 dxx1I .

În primul rând trebuie observat că ambele integrale se pot calcula analitic, cu ajutorul formulei Leibniz - Newton. Astfel,

6404256133.212ln

152ln

2I4

0

x

1 === , respectiv 3025850930.210lnxlnI 1012 === .

Se calculează apoi, cu ajutorul unui program special destinat, valorile aproximative ale integralelor, atât cu formula trapezelor cât şi cu formula Cavalieri - Simpson, pentru diferite numere de subintervale, n. Rezultatele obţinute sunt următoarele:

Tab. 1

Formula n I1 I2

10 21.77887776 2.36521443 1 02 21.64181188 2.30325280 103 21.64043948 2.30259177

rapezelor t

104 21.64042575 2.30258516 10 21.64112961 2.31197800 102 21.64042568 2.30258723

Cavalieri - Simpson

21.6404256133 2.3025850930 103

79

METODE NUMERICE

80

După cum se poate observa din tabelul 1, rezultatele obţinute confirmă precizia (m

. Chestiuni de verificare

ult) mai bună a formulei Cavalieri - Simpson; de asemenea, se poate remarca şi calitatea deosebită a rezultatelor - în cazul formulei amintite, pentru n = 103 diviziuni ale intervalului, primele zece zecimale fiind identice cu cele ale valorii exacte. 4

Să se calculeze, cu ajutorul formulei trapezelor şi apoi cu formula Cavalieri

4.1.

– Simpson, următoarele integrale:

∫5

1dxx

4.2.

.3.

Comparaţi rezultatele obţinute prin cele două metode.

∫5

0x dx3

4 ( )∫ +−7

02 dx5x2x3

CAPITOLUL 14 – Exemple de aplicare a metodelor numerice

Capitolul 14

EXEMPLE DE APLICARE A METODELOR NUMERICE ÎN INGINERIA TEHNOLOGICĂ.

1. Profilarea sculei cilindro-frontală pentru prelucrarea filetelor

pătrate

În cele ce urmează este prezentată un exemplu de aplicare a metodelor numerice de rezolvare a ecuaţiilor, în domeniul profilării sculelor aşchietoare. Astfel, fie cazul generării suprafeţelor elicoidale (în particular, un filet pătrat) cu ajutorul unei scule cilindro-frontală.

Pentru a determina ecuaţiile parametrice ale suprafeţei filetului, conform notaţiilor din figura 1, se porneşte de la dreapta (Δ), care printr-o mişcare elicoidală de pas p descrie suprafaţa filetului:

81

( ),bZ;0Y;uX

===

Δ (1)

în care u este un parametru, variind între umin = Ri şi umax = Re.

Fig.1 - Ecuaţiile parametrice ale

suprafeţei filetului

(Δ)

METODE NUMERICE

Imprimând acestei drepte o mişcare de rotaţie, de parametru unghiular φ, în jurul axei OZ, concomitent cu o translaţie în lungul aceleiaşi axe, coordonatele punctului curent de pe suprafaţa astfel generată satisfac ecuaţia matriceală

ϕ⋅

+⋅ϕϕϕ−ϕ

=p

00

b0u

1000cossin0sincos

ZYX

. (2)

Prin dezvoltarea relaţiei (2) rezultă ecuaţiile parametrice ale suprafeţei

filetului:

( ).pbZ;sinuY;cosuX

ϕ⋅+=ϕ⋅=ϕ⋅=

Σ (3)

Determinarea profilului sculei cilindro-frontală, de axă A , reciproc înfăşurătoare suprafeţei (Σ), se poate face prin impunerea condiţiei Nikolaev,

( ) 0A,r,N =Σ , (4)

în care ΣN reprezintă normala la suprafaţa (2.36), scrisă în punctul curent iar r - vectorul de poziţie al punctului curent.

Dacă în produsul mixt (4) se înlocuiesc

kujcospisinppcosusinu0sincoskji

N ⋅+⋅ϕ⋅−⋅ϕ⋅=ϕ⋅ϕ⋅−

ϕϕ=Σ , (5)

iA = , respectiv ( ) kpbjsinuicosur ⋅ϕ⋅++⋅ϕ⋅+⋅ϕ⋅= ,

atunci condiţia se rezumă la

0pbsinucosu001ucospsinp

=ϕ⋅+ϕ⋅ϕ⋅

ϕ⋅−ϕ⋅ (6)

Dezvoltând determinantul de mai sus după a doua linie, condiţia Nikolaev rezultă, în final, sub forma

. (7) ( ) 0sinucospbp 2 =ϕ⋅+ϕ⋅ϕ⋅+⋅

Ansamblul format de ecuaţiile (3) şi (7) dă profilul căutat al sculei cilindro-frontale. 82

CAPITOLUL 14 – Exemple de aplicare a metodelor numerice

Trebuie remarcat faptul că, dând o valoare oarecare parametrului u, condiţia (7) se transformă într-o ecuaţie transcendentă de necunoscută φ. Pentru rezolvarea acestei ecuaţii se va utiliza, în cele ce urmează, metoda generală prezentată în Lucrarea 2.

Mai concret, pentru determinarea profilului sculei, se dau parametrului u valori într-un număr nu de puncte, pentru fiecare dintre acestea rezolvându-se numeric ecuaţia (7) şi determinând valoarea corespunzătoare a lui φ; fiecare cuplu de valori (u, φ) astfel găsit se înlocuieşte în (3), ceea ce conduce la obţinerea coordonatele unui punct al profilului sculei.

Pentru rezolvarea unei aplicaţii concrete, a fost conceput un program de calcul, potrivit algoritmului de mai jos.

program Filet pătrat

real p, b, Ri, Re ; parametrii geometrici ai filetului real fi, u ; parametrii analitici ai suprafeţei filetului real pasfi, pasu ; incremenţi real X, Y, Z ; coordonatele punctelor profilului sculei real er ; eroarea maximă admisibilă la rezolvarea

ecuaţiei (3) întreg nu, nfi ; numărul punctelor de discretizare pentru

parametrii u şi fi întreg i, j ; contori

citeşte p, b, Ri, Re, er, nu, nfi

1−−

=nu

RiRepasu

11−

=nfi

pasfi

pentru i = 1, nu repetă u = Ri + (i-1)·pasu j = 1 repetă fi = (j-1)·pasfi j = j + 1 până când ( ) ersinucospbp <ϕ⋅+ϕ⋅ϕ⋅+⋅ 2

ϕ⋅+=ϕ⋅=ϕ⋅=

pbZsinuYcosuX

scrie X, Y, Z stop

83

METODE NUMERICE

Cu ajutorul programului menţionat, pentru p = 10 mm, b = 3 mm, Ri = 15 mm, Re = 20 mm, nu = 101 puncte, nfi = 100001 puncte, er = 0.01 s-au obţinut rezultatele din tabelul de mai jos.

Tab. 1

i u X Y Z 1 15.000 14.936 -1.380 2.079 2 15.050 14.987 -1.378 2.083 3 15.100 15.037 -1.376 2.087 4 15.150 15.088 -1.375 2.091 5 15.200 15.138 -1.373 2.096 . . . . . . . . . . . . . . .

99 19.900 19.864 -1.201 2.396 100 19.950 19.914 -1.200 2.398 101 20.000 19.964 -1.198 2.401

2. Profilarea sculei cremalieră pentru generarea canelurilor

triunghiulare

În cele ce urmează este prezentată o aplicaţie pentru metodele numerice de rezolvare a sistemelor de ecuaţii, în domeniul profilării sculelor aşchietoare şi anume profilarea sculei cremalieră utilizată pentru generarea canelurilor triunghiulare prin metoda traiectoriilor [15].

Metoda traiectoriilor se bazează pe o exprimare discretă a spaţiului sculei, introducând noţiunea de „nor de puncte” asociat spaţiului. Se acceptă un mod de reprezentare a spaţiului asociat sculei, ξη, în care se presupun cunoscute (determinabile) un număr suficient de mare (dar finit) de puncte, care formează ceea ce se va numi „nor de puncte”. Acestei reprezentări a spaţiului îi corespunde, în planul ξη, reprezentarea matriceală de forma

⎪⎭

⎪⎬⎫

⎪⎩

⎪⎨⎧

ηηηηξξξξ

=ξnjijj2j1

njijj2j1

KK

KK, j = 1, 2, . . . m (8)

Evident, un caz particular al acestui mod de reprezentare a spaţiului îl

reprezintă ceea ce se va numi „nor ordonat” de puncte, în care distribuţia punctelor corespunzătoare norului se face în baza unei legităţi, punctele formând o matrice cu distribuţii ale punctelor de forma

84

CAPITOLUL 14 – Exemple de aplicare a metodelor numerice

{ }{ } ,n,2,1i,j

,m,2,1j,i

0ij

0ij

K

K

=ηΔ⋅±η=η

=ξΔ⋅±ξ=ξ (9)

în care Δξ şi Δη reprezintă mărimi incrementale, suficient de mici, pentru o acoperire a spaţiului în concordanţă cu precizia dorită a rezolvării problemelor suprafeţelor în înfăşurare.

Această modalitate de reprezentare a spaţiului facilitează operarea cu aceste puncte şi simplifică algoritmii de calcul.

Pentru cazul generării cu scula-cremalieră, se poate imagina un algoritm care să permită determinarea profilului primar al acestei scule, pornind de la cunoaşterea formei suprafeţei de generat şi de la o reprezentare a spaţiului căreia îi aparţine cremaliera, în forma unui nor ordonat de puncte.

Fie figura 2, în care sunt definite următoarele sisteme de referinţă: ξη, care reprezintă un sistem mobil, asociat spaţiului sculei; XY - sistem mobil solidar cu semi-fabricatul; xy - sistem de referinţă fix.

În sistemul de referinţă mobil XY, este definit profilul de generat

Fig.2 - Generarea cu scula

cremalieră

( )( ).uYY

,uXX==

Σ (10)

Un punct oarecare al „norului de puncte” Mij(ξij , ηij) asociat spaţiului

ξη, în mişcarea relativă faţă de sistemul mobil XY ( ) [ ]aX 3 +ξ⋅ϕω= , (11)

în care ( )ϕϕ−ϕϕ

=ϕωcossinsincos

3 iar ϕ⋅−

−=

r

r

RR

a , descrie o traiectorie „T”, care

după înlocuiri şi calcule în (11), rezultă de forma

[ ] [ ][ ] [ ] .cosRsinR

;sinRcosRT

rijrij

rijrij

ϕ⋅ϕ⋅−η+ϕ⋅−ξ−ϕ⋅ϕ⋅−η+ϕ⋅−ξ

(12)

85

METODE NUMERICE

Condiţia ca punctul Mij(ξij , ηij) al norului de puncte să aparţină totodată şi unei traiectorii „T”, (12), care să aibă un singur punct de contact (de tangenţă) cu profilul Σ, conduce la ecuaţiile de identificare:

( ) [ ] [ ]( ) [ ] [ ]( ) [ ]( ) [ ]⎪

⎪

⎩

⎪⎪

⎨

⎧

ϕ⋅ϕ⋅−η−ϕ⋅ξ−=′

ϕ⋅ϕ⋅−η+ϕ⋅ξ−=′

ϕ⋅ϕ⋅−η+ϕ⋅−ξ−=

ϕ⋅ϕ⋅−η+ϕ⋅−ξ=

.sinRcosuY

;cosRsinuX

;cosRsinRuY

;sinRcosRuX

rijij

rijij

rijrij

rijrij

(13)

Acesta este, de fapt, un sistem de patru ecuaţii cu patru necunoscute:

u, φ, ξij, ηij ; prin rezolvarea unui număr n de astfel de sisteme, vor rezulta coordonatele a n puncte aparţinând profilului căutat al sculei-cremalieră.

Profilul unui flanc al canelurilor triunghiulare (a se vedea fig.2) poate fi descris ca

( )( ) .sinuuY

;cosuRuX r

ε⋅=ε⋅+−=

Σ (14)

Derivând în raport cu u relaţiile (14) şi înlocuind în (13), prin

renotarea lui ξij cu x şi a lui ηij cu y rezultă

(15)

[ ] [ ][ ] [ ]

[ ][ ]⎪

⎪⎩

⎪⎪⎨

⎧

=ε−ϕ⋅ϕ⋅−−ϕ⋅−=ε−ϕ⋅ϕ⋅−+ϕ⋅−

=ε⋅−ϕ⋅ϕ⋅−+ϕ⋅−−=ε⋅−+ϕ⋅ϕ⋅−+ϕ⋅−

.0sinsinRycosx;0coscosRysinx

;0sinucosRysinRx;0cosuRsinRycosRx

r

r

rr

rrr

Ecuaţiile din sistemul (15) sunt transcendente; sistemul este, în

principiu, de forma

( )( )( )( )⎪

⎪⎩

⎪⎪⎨

⎧

=ϕ=ϕ=ϕ=ϕ

.0,u,y,xf;0,u,y,xf;0,u,y,xf;0,u,y,xf

4

3

2

1

(16)

Pentru rezolvarea sistemului se poate utiliza o metodă generală, bazată

pe discretizarea variabilelor x, y, u şi φ,

, xixx 0i Δ⋅+= yjyy 0j Δ⋅+= , ukuu 0k Δ⋅+= , ϕΔ⋅+ϕ=ϕ l0l (17) şi declararea ca soluţie a cvartetului (xi, yj, uk, φl) pentru care

86

CAPITOLUL 14 – Exemple de aplicare a metodelor numerice

( )( )( )( )⎪

⎪

⎩

⎪⎪

⎨

⎧

<ϕ

<ϕ

<ϕ

<ϕ

,er,u,y,xf;er,u,y,xf;er,u,y,xf;er,u,y,xf

lkji4

lkji3

lkji2

lkji1

(18)

în care er are o valoare convenabil aleasă (suficient de mică).

Algoritmul programului de calcul care determină profilul sculei-cremalieră pentru generarea canelurilor triunghiulare este prezentat, în pseudo-cod, în cele ce urmează. program Caneluri triunghiulare

real Rr ; raza de rulare a semifabricatului real eps ; unghiul de înclinare al flancurilor canelurii real umax ; lungimea profilului canelurii real fi, u ; parametrii analitici ai suprafeţei generate real dfi, du, dx, dy ; incremenţi real x, y ; coordonatele punctului curent din nor real f1, f2, f3, f4 ; valorile numerice absolute ale condiţiilor

de coincidenţă real er ; eroarea maximă admisibilă la rezolvarea

sistemului (16) întreg nu, nfi, nx,ny ; numărul punctelor de discretizare pentru

parametrii u, fi, x şi y întreg i, j, k, l ; contori

citeşte Rr, eps, umax, nu, nfi, nx, ny, dfi, dx, dy, er

1−=

numaxudu

pentru k = 1, nu repetă u = (k -1)·du pentru l = 1, nfi repetă fi = (l -1)·dfi pentru i = 1, nx repetă x = (i -1)·dx j = 1 repetă y = y + (j -1)·dy calculează f1, f2, f3, f4 conform (15) j = j +1 până când (condiţia 18) sau (j > ny) scrie x, y

87stop

METODE NUMERICE

3. Aproximări ale profilurilor determinate prin puncte

Sunt numeroase situaţiile în care profilurile sculelor aşchietoare sunt determinate numeric, cu ajutorul teoriei suprafeţelor în înfăşurare, sub forma unei mulţimi finite de puncte. Pentru o mai uşoară reprezentare grafică a unui astfel de profil, ca şi pentru rezolvarea diverselor probleme practice, prezintă interes găsirea expresiei analitice a unei funcţii care să aproximeze profilul determinat prin puncte.

În acest sens, în cele ce urmează se va prezenta o aplicaţie în tehnică a noţiunilor teoretice prezentate în capitolul de faţă: aproximarea, prin metoda abaterii medii pătratice minime, a profilului sculei-cremalieră utilizată pentru generarea canelurilor dreptunghiulare.

Fie schema de generare a unei caneluri exterioare dreptunghiulare cu ajutorul unei scule-cremalieră (fig. 3); se consideră următoarele sisteme de referinţă: xy - sistem de referinţă fix; XY - sistem de referinţă mobil, solidar cu suprafaţa generată; ξη - sistem de referinţă mobil, solidar cu scula-cremalieră.

Ecuaţiile suprafeţei generate (în planul z = 0) sunt

.bY;uX

=−=

Σ (19)

88

Mişcarea relativă între sculă şi

semifabricat este descrisă prin relaţia

( ) aXT3 −⋅ϕω=ξ , (20)

în care

( )ϕϕϕ−ϕ

=ϕωcossinsincosT

3 , ϕ⋅−

−=

r

r

RR

a .

Fig.3 - Generarea canelurilor cu scula-cremalieră

(21) Reunind relaţiile (19), (20) şi (21) rezultă, după efectuarea calculelor,

ecuaţiile familiei de suprafeţe:

.Rcosbsinu

;Rsinbcosu

r

r

ϕ⋅+ϕ⋅+ϕ⋅−=η+ϕ⋅−ϕ⋅−=ξ

(22)

Pentru determinarea înfăşurătoarei familiei de suprafeţe (22), adică a

profilului sculei-cremalieră, se pot utiliza mai multe metode.

CAPITOLUL 14 – Exemple de aplicare a metodelor numerice

Conform metodei Gohman, ecuaţiile (22) se asociază cu condiţia de înfăşurare, care în acest caz are forma

ϕ⋅= cosRu r ; (23)

înlocuind (23) în (22) se obţine

.RcosbcossinR

;RsinbcosR

rr

r2

r

ϕ⋅+ϕ⋅+ϕ⋅ϕ⋅−=η+ϕ⋅−ϕ⋅−=ξ (24)

Ecuaţiile (24), în care singurul parametru este φ, reprezintă profilul

sculei-cremalieră; domeniul de variaţie pentru φ se determină cu ajutorul relaţiei (23), ţinând cont că u poate lua valori între

22

imin bRu −= şi 22emax bRu −= , (25)

Ri şi Re fiind raza interioară, respectiv exterioară, pentru arborele canelat considerat.

Pentru determinarea efectivă a punctelor profilului sculei, se discretizează mai întâi u într-un număr n de puncte de la umin la umax; pentru fiecare valoare a lui u se determină, apoi, valoarea corespunzătoare a lui φ, rezolvând ecuaţia (23), după care pentru perechea (u, φ) găsită se calculează punctul de pe profilul sculei, (ξ,η), cu (24).

În acest stadiu, profilul sculei rezultă sub forma unei mulţimi de puncte, ( ) ( ) ({ nn2211 ,,,,, )}ηξηξηξ K . Se poate determina, acum, o funcţie polinomială, P, care să aproximeze cu abatere medie pătratică minimă profilul sculei.

Astfel, dacă se urmăreşte aproximarea printr-un polinom de gradul doi, , se calculează ( ) 01

222 axaxaxP +⋅+⋅=

, k = 0, 1 . . . 4, respectiv , j = 0, 1, 2, (26) ∑=ξ=

n

1i

kiku ∑

=ξ⋅η=

n

1i

jiijv

după care se rezolvă - de exemplu prin metoda Gauss - sistemul

. (27) ⎟⎟⎟

⎠

⎞

⎜⎜⎜

⎝

⎛=

⎟⎟⎟

⎠

⎞

⎜⎜⎜

⎝

⎛⋅

⎟⎟⎟

⎠

⎞

⎜⎜⎜

⎝

⎛

2

1

0

2

1

0

432

321

210

vvv

aaa

uuuuuuuuu

Algoritmul unui program de calcul care, pornind de la profilul

canelurilor de generat, găseşte forma polinomului P2, este prezentat în continuare, redactat în pseudo-cod. 89

METODE NUMERICE

program Caneluri dreptunghiulare

real Rr ; raza de rulare a semifabricatului real Ri, Re, b ; parametrii geometrici ai canelurii de

generat real u, umin, umax ; parametrul u şi limitele sale de variaţie real fi, fimin, fimax ; parametrul fi şi limitele sale de variaţie real dfi, du, dx, dy ; incremenţi real csi, eta ; coordonatele punctului curent de pe

profilul generat real m, S ; parametri auxiliari pentru rezolvarea

sistemului prin metoda Gauss real a0, a1, a2 ; coeficienţii polinomului de aproximare întreg n ; numărul punctelor de discretizare pentru

parametrul u întreg nn ; rangul matricii sistemului (27) întreg i, j, k ; contori tablou real A(nn,nn) ; matricea sistemului tablou real B(nn) ; vectorul-coloană a termenului liber tablou real X(nn) ; vectorul-soluţie

citeşte Rr, b, Ri, Re, n, nn 22 bRiminu −= 22 bRemaxu −=

Rrminuarccosminfi =

Rrmaxuarccosmaxfi =

a(1,2) = 0 a(1,3) = 0 a(2,3) = 0 a(3,3) = 0 c(1) = 0 c(2) = 0 c(3) = 0 pentru i = 1, n repetă

( )1

1−−

⋅−+=n

minfimaxfiiminfifi

( ) ( ) RrfisinbficosRrcsi +⋅−⋅−= 2 ( ) ( ) ( ) fiRrficosbficosfisinRreta ⋅+⋅+⋅⋅−= a(1,2) = a(1,2) + csi a(1,3) = a(1,3) + csi2 a(2,3) = a(2,3) + csi3 a(3,3) = a(3,3) + csi4 c(1) = c(1) + eta

90 c(2) = c(2) + eta·csi

CAPITOLUL 14 – Exemple de aplicare a metodelor numerice

c(3) = c(3) + eta·csi2 a(1,1) = n a(2,1) = a(1,2) a(2,2) = a(1,3) a(3,1) = a(1,3) a(3,2) = a(2,3)

pentru k = 1, nn-1 repetă pentru i = k+1, nn repetă m = a(i,k) / a(k,k) pentru j = k+1, nn repetă a(i,j) = a(i,j) - m·a(k,j) b(i) = b(i) - m·b(k)

x(nn) = b(nn) / a(nn,nn) pentru i = nn-1, 1, -1 repetă s = b(i) pentru j = nn, i+1, -1 repetă s = s - a(i,j)·x(j) x(i) = s / a(i,i) a0 = x(1) a1 = x(2) a2 = x(3) scrie a0, a1, a2 stop

Cu ajutorul unui program de calcul redactat pe baza algoritmului prezentat, s-a obţinut, pentru Ri = 25 mm, Re = 30 mm, Rr = 30 mm, b = 3 mm, n = 1001 puncte, nn = 3,

( ) 98031.2x19017.0x03297.0xP 22 +⋅+⋅= .

Dacă pentru determinarea înfăşurătoarei familiei de suprafeţe (22) se utilizează metoda distanţei minime, condiţia de înfăşurare se scrie sub forma

( ) 0R uru =′η⋅ϕ⋅−η+′ξ⋅ξ , (28)

în care ξ şi η sunt daţi de (22). Înlocuind în (28) ϕ−=′ξ cosu şi ϕ−=′η sinu , condiţia de înfăşurare

rezultă sub forma unei ecuaţii transcendente. Rezolvând numeric, pentru n valori cuprinse între umin şi umax, ecuaţia

(28) - prin metoda generală - se obţin şi în acest caz n cupluri (u, φ), care înlocuite în (22) conduc la n puncte ale profilului sculei. Evident, în acest caz, cele n puncte ale profilului sculei diferă de cele obţinute prin metoda Gohman.

Pentru a putea compara soluţiile problemei considerate, obţinute prin cele două metode, printr-un program asemănător cu cel de mai sus se determină funcţia polinomială de gradul al doilea care aproximează cu abatere medie pătratică minimă profilul dat de punctele rezultate prin metoda distanţei minime. Rezultă, în acest caz

. ( ) 97107.2x19655.0x03219.0xP 22 +⋅+⋅=

91

METODE NUMERICE

Valorile coeficienţilor acestui polinom s-au obţinut pentru o eroare maximă admisibilă, la rezolvarea ecuaţiei (28), er = 0.001. Este evident că formele celor două polinoame sunt aproape identice, aşa după cum era de aşteptat. 4. Determianrea secţiunii aşchiei la prelucrarea arborilor cu profil

hexagonal cu ajutorul sculei-cremalieră

În continuare se va prezenta o aplicaţie a metodelor numerice de integrare în domeniul teoriei generării suprafeţelor în înfăşurare reciprocă: determinarea secţiunii aşchiei la prelucrarea arborilor cu profil hexagonal cu ajutorul sculei-cremalieră.

În acest scop, în conformitate cu notaţiile din figura 4, se consideră următoarele sisteme de referinţă: xy - sistem de referinţă fix; XY - sistem de referinţă mobil, solidar cu suprafaţa generată; ξη - sistem de referinţă mobil, solidar cu scula-cremalieră.

Ecuaţiile suprafeţei generate (în planul z = 0) sunt

.uY

;bX=−=

Σ (29)

Mişcarea relativă între sculă şi semifabricat este descrisă prin relaţia

( ) aXT3 −⋅ϕω=ξ , (30)

92

în care

( )ϕϕϕ−ϕ

=ϕωcossinsincosT

3 , ϕ⋅−

−=

r

r

RR

a

(31)

Reunind relaţiile (29), (30) şi (31) rezultă, după efectuarea calculelor, ecuaţiile familiei de suprafeţe:

Fig.4 - Schema de generare a arborelui hexagonal

cu scula-cremalieră

.Rcosusinb

;Rsinucosb

r

r

ϕ⋅+ϕ⋅+ϕ⋅−=η+ϕ⋅−ϕ⋅−=ξ

(32)

Pentru determinarea înfăşurătoarei familiei de suprafeţe (32), adică a

profilului sculei-cremalieră, se scrie condiţia de înfăşurare conform metodei distanţei minime,

CAPITOLUL 14 – Exemple de aplicare a metodelor numerice

′η

′η=′ξ

′ξ

ϕϕ

uu , (33)

în care ξ şi η sunt conform (32). După derivare şi efectuarea calculelor, rezultă ϕ⋅= sinRu r . (34)

Înlocuind condiţia de înfăşurare, sub forma (34), în relaţia (32), se

obţine profilul sculei-cremalieră:

.RcossinRsinb

;cosRcosb

rr

2r

ϕ⋅+ϕ⋅ϕ⋅+ϕ⋅−=ηϕ⋅+ϕ⋅−=ξ (35)

Dacă se notează cu θ parametrul mişcării de rulare, mişcarea sculei faţă de semifabricat poate fi descrisă de ecuaţia

( )[ ]aX 3 +ξθω= , (36)

93

în care matricele ω3 şi a au semnificaţia dată de relaţia (31).

Înlocuind în (36) pe ω3, a şi pe ξ, conform (35), şi dând parametrului θ valorile corespunzătoare la doi dinţi consecutivi ai sculei- cremalieră, θ1 şi θ2, se obţin poziţiile relative T1 şi T2 ale acesteia faţă de piesa prelucrată, potrivit fig. 5.

Conform notaţiilor alăturate, cele două traiectorii se intersectează cu cercul de secţiune a semifabricatului în C1,

respectiv C2, iar între ele în punctul I.

Fig.5 - Determinarea secţiunii aşchiei

Secţiunea aşchiei este, în aceste condiţii, poligonul curbiliniu IC1C2. Dacă se determină fiecare din arcele IC1, C1C2 şi C2I printr-un număr oarecare de puncte, pe cale numerică, în sistemul de referinţă X0Y0, atunci aria secţiunii aşchiei se poate calcula, în principiu, cu formula

A = A1 - A2 + A3, (37)

în care A1 reprezintă aria cuprinsă între arcul IC1, axa X0 şi dreptele X0 = XI, respectiv X0 = XC1, A2 - aria cuprinsă între arcul IC2, axa X0 şi dreptele X0 = XI,

METODE NUMERICE

respectiv X0 = XC2 iar A3 - aria cuprinsă între arcul C1C2, axa X0 şi dreptele X0 = XC1, respectiv X0 = XC2.

Fiecare dintre cele trei arii de mai sus se poate calcula numeric, prin metoda trapezelor.

Algoritmul unui program de calcul pentru determinarea coordonatelor punctelor I, C1 şi C2 este prezentat, în pseudo-cod, în cele ce urmează.

Program Intersecţii real a, b ; latura şi apotema hexagonului real Rr ; raza de rulare real fi1, t1, u ; parametrii ecuaţiilor (35), (36)

pentru traiectoria T1 real fi2, t2, v ; parametrii ecuaţiilor (35), (36)

pentru traiectoria T2 real X1, Y1 ; coordonatele punctelor de pe T1 real X2, Y2 ; coordonatele punctelor de pe T2 real csi1, csi2, eta1, eta2 ; variabile intermediare real d ; distanţa curentă între două puncte real er ; eroarea maximă admisibilă la

rezolvarea numerică a ecuaţiilor real XI, YI, XC1, YC1,

XC2, YC2 ; coordonatele punctelor de intersecţie întreg n ; numărul de puncte de discretizare a

fiecărei traiectorii întreg i, j ; contori citeşte n, a, Rr, er, t1, t2

23ab ⋅=

pentru i = 1, n repetă

( )1n

a1i2au

−⋅−+−=

Rruarcsinfi =1

calculează csi1, eta1 cu relaţia (35) calculează X1, Y1 cu relaţia (36), pentru teta = t1 j = 1 repetă

( )1n

a1j2av

−⋅−+−=

94

CAPITOLUL 14 – Exemple de aplicare a metodelor numerice

Rrvarcsin2fi =

calculează csi2, eta2 cu relaţia (35) calculează X2, Y2 cu relaţia (36), pentru teta = t2 ( ) ( )22 2Y1Y2X1Xd −+−= j = j + 1

până când (d < er) sau (j > n) dacă d < er atunci XI = X1 YI = Y1 salt la eticheta 1

1 i = 1 repetă

( )1n

a1i2au

−⋅−+−=

Rruarcsin1fi =

calculează csi1, eta1 cu relaţia (35) calculează X1, Y1 cu relaţia (36), pentru teta = t1 222 a1Y1Xd −+= i = i + 1 până când d < er XC1 = X1 YC1 = Y1 i = 1 repetă

( )1n

a1i2au

−⋅−+−=

Rruarcsin2fi =

calculează csi2, eta2 cu relaţia (35) calculează X2, Y2 cu relaţia (36), pentru teta = t2 222 a2Y2Xd −+= i = i + 1 până când d < er

XC2 = X2 YC2 = Y2

stop

Conform celor prezentate mai sus şi considerând cunoscute coordonatele punctelor de intersecţie cunoscute, cu ajutorul unui alt program, bazat pe algoritmul de mai jos, se poate calcula aria secţiunii aşchiei.

95

METODE NUMERICE

Program Secţiunea aşchiei

c

întreg i, j, k

real a, b ; latura şi apotema hexagonului real Rr ; raza de rulare real fi, u ; parametrii ecuaţiei (35)

real t1, t2 ; parametrul θ din ecuaţia (36) pentru traiectoria T , respectiv 1 T2 ; coordonatele punctelor de pe T1 real X1, Y1

real X2, Y2 ; coordonatele punctelor de pe T2 real X, Y ; coordonatele punctelor de pe cer real csi, eta ; variabile intermediare real d1, d2 ; unghiurile la centru pentru C1 şi C2 real d ; unghiul la centru curent

real I1, I2, I3 ; variabile intermediare real XI, YI, XC , Y 1 C1,

XC2, YC2 ; coordonatele punctelor de intersecţie A1, A2, A3 ; ariile parţiale din formula (37) real

real A ; aria secţiunii aşchiei întreg n ; numărul de puncte de discretizare a

fiecărui arc de curbă ; contori

citeşte n, a, Rr, t1, t2, XI, YI, XC1, YC1, XC2, YC2

23ab ⋅=

I1 = XI + XC1 + 2·b I2 = XI + XC2 + 2·b I3 = 0 j = 0 k = 0 pentru i = 1, n repetă

96

( ) ( )1n2a1iu−

⋅−=

Rruarcsinfi =

calculează csi, eta cu relaţia (35) pentru teta = t1 calculează X1, Y1 cu relaţia (36),

dacă (Y1 > YI) şi (Y1 < YC1) atunci I1 = I1 + b1X2 −⋅ j = j + 1 calculează X2, Y2 cu relaţia (36), pentru teta = t2 dacă (Y2 > YI) şi (Y2 < YC2) atunci I2 = I2 + b2X2 −⋅ k = k + 1

CAPITOLUL 14 – Exemple de aplicare a metodelor numerice

97

d1 = 1XC1YCarctg d2 =

2XC2YCarctg

pentru i = 1, n repetă

d = d1+ ( )1n1d2d1i

−−

⋅−

X = a·cosd Y = a·sind dacă (i = 1) sau (i = n) atunci I3 = I3 + (x - b) altfel I3 = I3 + 2·(x - b)

A1 = ( )1j2YI1YC1I

−−

⋅ A2 = ( )1k2YI2YC2I

−−

⋅ A3 = ( )1n21YC2YC3I

−−

⋅

A = A1 - A2 + A3 scrie A stop

Dacă se consideră un exemplu numeric concret, pentru a = 20 mm, Rr = 20 mm, n = 1000 puncte, θ1 = 0˚, θ2 = 10˚, cu primul dintre cele două programe de calcul prezentate rezultă coordonatele punctelor de intersecţie

I(-17.38693; 1.73372), C1(-18.58483; 7.40034) şi C2(-17.98355; 8.74291)

Pe baza datelor de intrare deja menţionate şi a acestor coordonate, cu ajutorul celui de-al doilea program se calculează aria secţiunii aşchiei; rezultă A = 3.15785 mm2.