Capitolul 2[1]Metode numerice

53
Capitolul 2. Metode numerice Capitolul 2. Metode numerice 2.1. Metoda înjumătăţirii intervalului Descrierea matematică Fie dată ecuaţia . (1) Vom considera cazul, când funcţia este continuă pe [a, b] şi . Suplimentar vom considera că pe [a, b] semnul derivatei 1 a funcţiei este constant, deci avem doar o singură soluţie. Pentru determinarea soluţiei ecuaţiei (1) vom detecta mijlocul segmentului [a, b], şi vom calcula valoarea funcţiei în acest punct. Dacă , atunci este soluţia ecuaţiei. În caz contrar cercetăm segmentele [a, ] şi [, b]. Pentru aproximarea următoare vom selecta acel segment, pentru care valoarea funcţiei în extremităţi are semne opuse. Dacă Procesul de construire consecutivă a sistemului de segmente încluse, care conţîn soluţia funcţiei pe segmentul [a,b]. La fiecare iteraţie se calculează valoarea

Transcript of Capitolul 2[1]Metode numerice

Page 1: Capitolul 2[1]Metode numerice

Capitolul 2. Metode numerice

Capitolul 2. Metode numerice

2.1. Metoda înjumătăţirii intervalului

Descrierea matematică

Fie dată ecuaţia

. (1)

Vom considera cazul, când funcţia este continuă pe [a, b] şi . Suplimentar vom considera că pe [a, b] semnul derivatei 1 a funcţiei este constant, deci avem doar o singură soluţie.

Pentru determinarea soluţiei ecuaţiei (1) vom detecta mijlocul segmentului [a, b],

şi vom calcula valoarea funcţiei în acest punct. Dacă , atunci este soluţia ecuaţiei. În caz contrar cercetăm segmentele [a, ] şi [, b].

Pentru aproximarea următoare vom selecta acel segment, pentru care valoarea funcţiei în extremităţi are semne opuse. Dacă

, atunci vom continua cercetarea pe segmentul [a1, b1], unde ,

. În caz contrar, extremităţile vor fi , . Noul segment [a1, b1] din nou se divizează, apoi se repetă cercetarea semnelor valorilor funcţiei în extremităţi şi în mijlocul segmentului. Procedura se repetă până când se obţine soluţia exactă sau devierea soluţiei aproximative de la cea exactă nu devine suficient de mică.

În procesul de construcţie a segmentelor succesive obţinem consecutivitatea segmentelor

1

Procesul de construire consecutivă a sistemului de segmente încluse, care conţîn soluţia funcţiei pe segmentul [a,b]. La fiecare iteraţie se calculează valoarea funcţiei la extremităţi şi la mijlocul segmentului.

Page 2: Capitolul 2[1]Metode numerice

Capitolul 2. Metode numerice

[a, b], [a1, b1], [a2, b2], [a3, b3], ... [an, bn]... .

Pentru fiecare din ele are loc relaţia

, . (2)

Deoarece lungimea fiecărui segment următor este egală cu ½ din lungimea celui precedent, putem exprima lungimea oricărui segment prin cea a segmentului iniţial:

. (3)

Din construcţie şi proprietăţile funcţiei , rezultă că şirul extremităţilor stângi a, a1, a2, ..., an , ... este monoton crescător, mărginit în partea superioară, iar şirul extremităţilor drepte b, b1, b2, ..., bn , ... este monoton descrescător, mărginit inferior. De aici rezultă convergenţa ambelor şiruri şi existenţa limitei pentru fiecare din ele.

Trecând la limită în egalitatea (3) obţinem:

.

Trecând la limită în inegalitatea (2) din continuitatea obţinem . Prin urmare, , deci este soluţia ecuaţiei (1).

Deoarece este un punct din intervalul [an, bn], rezultă că

.

În cazul când semnul derivatei întâi alternează pe segmentul [a, b] , adică rădăcinile ecuaţiei nu sunt separate, metoda permite determinarea doar a unei singure soluţii.

Algoritmizarea metodei

Datorită simplităţii sale metoda este uşor de realizat în formă algoritmică:1. Ecuaţia se reduce la forma . Acest lucru simplifică programarea metodei.2. Se introduc valorile a, b – extremităţile segmentului, – şi e – exactitatea cu care

trebuie obţinută soluţia.3. Se calculează f(a), f(b), f(c).4. Dacă , atunci vom considera că a trece în c. În caz contrar b

va primi valoarea lui c.5. Pentru noile valori a şi b repetăm paşii 2 –3 atât timp cât .6. Afişăm în calitate de soluţie mijlocul ultimului segment [a, b].

Exemplu

Elaboraţi un program ce calculează soluţia ecuaţiei

2

Page 3: Capitolul 2[1]Metode numerice

Capitolul 2. Metode numerice

pe segmentul [0, 1] pentru valoarea erorii =0,0001. În program nu veţi utiliza funcţiile standard ale limbajului PASCAL cos(x) şi exp(x). In calcule veţi folosi descompunerea in serie Taylor:

;

.

Valorile funcţiilor vor fi calculate cu aceeaşi exactitate.

Rezolvare. Pentru a calcula funcţiile ex şi cos cu exactitatea cerută, avem nevoie de formula termenului rest în descompunerea funcţiilor respective în serie Taylor. Vom folosi formula Lagrange a termenului rest:

.

Pentru simplitate vom înlocui prin pe segmentul dat:

.

Pentru funcţia ex avem

.

Derivata funcţiei ex de orice ordin este ex şi pe segmentul [0, 1] valoarea ei maxima este 2,71 3. Mai ţinem cont că . Deci, cum atingem un indice n pentru care

,

oprim calculul funcţiei ex. Pentru termenul rest al funcţiei avem:

.

3

Page 4: Capitolul 2[1]Metode numerice

Capitolul 2. Metode numerice

Derivata funcţiei cos(x) de orice ordin impar este sin(x) cu variaţie de semn şi pe segmentul [0, 1] valoarea ei maxima nu depăşeşte 1. Atingând un indice n pentru care

,

oprim calculul funcţiei cos(x).

Program P1;var m1,m2,x,e,a,b,c,a1,b1,c1,c2 : real; i,nb,nc :integer;function expo(x:real) : real;{calculul functiei e la puterea x}

var t,s : real; n : integer; begin s:=0; t:=1; n:=1; while t>e do begin s:=s+t; t:=t*x/n; n:=n+1; end; expo:=s; end;function coso(x:real):real;{calculul functiei cos(x)}var t,s : real; semn,n:integer; begin s:=0;t:=1;n:=2; semn:=-1; while abs(t)>e do begin s:=s+t; t:=semn*t*x*x/(n-1)/n; n:=n+2; end; coso:=s; end;{calculul functiei din ecuatie}function f(x:real):real;begin f:=2*coso(x)*coso(x)-expo(x/2);end;

begin{introducerea datelor: fixam valorile direct} e:=0.0001; a:=0; b:=1;nb:=0; a1:=a; b1:=b; {metoda bisectiei} while (b-a) >e do begin nb:=nb+1; c:=(a+b)/2; if f(c)*f(a)>0 then a:=c else b:=c; end; writeln(' x=',c:8:6,' f(x)=',f(c):10:6,' nb=',nb); {metoda coardelor} a:=a1; b:=b1; m1:=-1;m2:=-5; {m1 -supremul f'(x), m2 - infimul f'(x)}

4

Page 5: Capitolul 2[1]Metode numerice

Capitolul 2. Metode numerice

c1:=a; c2:=b; nc:=0;while abs((m1-m2)/m2)*abs(c1-c2)>e do begin c:=a-f(a)/(f(b)-f(a))*(b-a); c2:=c1; c1:=c; if f(c)*f(a)>0 then a:=c else b:=c; inc(nc); end; writeln(' x=',c:8:6,' f(x)=',f(c):10:6,' nc=',nc); if nb>nc then writeln('Numar minim de iteratii - metoda

coardelor') else writeln('Numar minim de iteratii - metoda bisectiei');readln;end.

Întrebări şi exerciţii

1. Care sînt avantajele şi dezavantajele metodei înjumătăţirii intervalului?2. Daţi formularea matematică a metodei înjumătăţirii intervalului.3. Elaboraţi o procedură care realizează metoda înjumătăţirii intervalului.4. Formulaţi criteriul de convergenţă al metodei înjumătăţirii intervalului.5. Care este condiţia de terminare a iteraţiilor în metoda înjumătăţirii intervalului?6. Reprezentaţi pe un desen procesul de căutare a soluţiei prin metoda înjumătăţirii

intervalului.7. Estimaţi eroarea metodei înjumătăţirii intervalului.8. Elaboraţi un program ce calculează soluţia aproximativă a ecuaţiei

3cos2x2 = ex/2

pe intervalul [0,1; 1] pentru cinci erori ε diferite , egale respectiv cu 0,1; 0,02; 0,003; 0,0004; 0,00005. Pentru fiecare valoare a erorii se va afişa la ecran o linie ce va conţine: eroarea ε, soluţia aproximativă x şi numărul de iteraţii n.

9. Elaboraţi un program ce calculează soluţia aproximativă a ecuaţiei

pe intervalul [0, 1] cu eroarea 0,0001.10. Calculaţi prin metoda înjumătăţirii intervalului soluţiile următoarelor ecuaţii:

a) ex + x – 2 = 0, x[-1, 1];b) x5 - 4x4 + x – 9 = 0, x[0, 4];c) ex – x –cos x – 1 = 0, x[0, 3].

5

Page 6: Capitolul 2[1]Metode numerice

Apropierea succesivă de soluţia ecuaţiei prin metoda coardei

Capitolul 2. Metode numerice

2.2. Metoda coardei

Descrierea matematică

Fie dată ecuaţiaf(x) = 0.

Pentru utilizarea metodei coardei, vom cere îndeplinirea aceloraşi condiţii asupra funcţiei f(x) ca şi pentru metoda precedentă. Ca şi în cazul metodei precedente, în calitate de soluţie vom căuta un punct interior al intervalului, dar de această dată el nu va fi mijlocul segmentului, ci punctul intersecţiei dreptei ce trece prin punctele (a, f(a)) şi (b, f(b)) cu axa 0x.

Ideea generală a metodei este următoarea: prin punctele (a, f(a)) şi (b, f(b)) se trasează o dreaptă. Se determină punctul c în care ea intersectează axa 0x. Apoi se determină semnul f(c). Extremitatea în care semnul funcţiei coincide cu semnul f(c), trece în c. Procesul se repetă până nu obţinem o apropiere suficientă de soluţia exactă. Această interpretare a metodei o face să fie tot atât de lentă ca şi metoda bisecţiei.

Dacă însă, suplimentar, vom cere prezenţa şi semnului constant al derivatei de ordin 2 pentru funcţia f(x) pe tot segmentul [a, b] ( pe segmentul dat funcţia să fie concavă sau convexă.), putem aplica o variaţie mai puternica a acestei metode. Ecuaţia dreptei ce trece prin 2 puncte date este determinată de relaţia:

.

Prin urmare,

.

În cazul intersecţiei cu axa x avem y = 0, deci putem calcula poziţia punctului c pe axa 0x. Pentru a realiza calculul soluţiei vom crea şirul {xn} conform următoarelor reguli.

Iniţial

.

La paşii următori, valorile elementelor şirului se vor forma, ţinând cont de următoarele situaţii:

6

Page 7: Capitolul 2[1]Metode numerice

Capitolul 2. Metode numerice

a) dacă funcţia este convexă (f ’’(x) > 0 ) (secanta este plasată mai sus de graficul funcţiei), atunci extremitatea segmentului iniţial, în care valoarea funcţiei este pozitivă, rămâne nemişcată (fixă), iar cealaltă extremitate trece de fiecare dată în punctul nou calculat.

b) dacă funcţia este concavă (f ’’(x) < 0 ) (secanta este plasată mai jos de graficul funcţiei), atunci extremitatea segmentului iniţial, în care valoarea funcţiei este negativă, rămâne nemişcată (fixă), iar cealaltă extremitate trece de fiecare dată în punctul nou calculat.

Cu alte cuvinte, metoda secantelor are un punct fix, celălalt modificîndu-se de fiecare dată. Observăm că fixă este extremitatea pentru care are loc relaţia

f (x) f ’’(x) > 0.

Pentru a calcula cealaltă extremitate a noului segment, trebuie să ţinem cont de următoarele cazuri posibile:

1) f ’(x) > 0 şi f ’’(x) > 0;2) f ’(x) > 0 şi f ’’(x) < 0;3) f ’(x) < 0 şi f ’’(x) > 0;4) f ’(x) < 0 şi f ’’(x) < 0.

Vom cerceta cazul 1 în care funcţia este convexă, crescătoare, celelalte cazuri fiind cercetate în mod similar.

Deoarece a doua derivată e mai mare ca 0, graficul funcţiei este inferior coardei. Pornind din extremitatea în care funcţia este pozitivă, vom obţine intersecţia cu axa 0x într-un punct interior al [a, b] (mai întâi se intersectează axa, apoi se ajunge la extremitatea opusă.) Repetând procedura, vom obţine un şir de puncte interioare, care tind spre soluţia exactă. Dacă însă vom fixa extremitatea negativă, coarda generată la pasul 2 va intersecta mai întâi graficul, apoi axa 0x, ceea ce poate genera părăsirea segmentului [a, b].

Deoarece funcţia este crescătoare, vom fixa extremitatea B, luând în calitate de x0

extremitatea A.Conform formulei pentru secantă,

determinăm că a: = x1.Repetând procedura, obţinem:

Vom demonstra acum, că şirul de aproximări a = x0 , x1 , …, xn … converge către soluţia ecuaţiei f(x) = 0

Din formula recurentă

7

Page 8: Capitolul 2[1]Metode numerice

Capitolul 2. Metode numerice

rezultă că xn+1 > xn n (numitorul fracţiei este pozitiv deoarece funcţia este crescătoare, şi negativ, deoarece xn, este punctul de intersecţie a coardei cu 0x, iar graficul funcţiei în acest punct este inferior coardei, deci negativ, (b - xn) – pozitiv, deci tot produsul este negativ, cu “–“ în faţă, trece în “+”).

Aşadar, avem şirul x0 < x1 < … <xn …< 1 < b. Şirul este crescător şi mărginit, deci există limita lui 1. Trecând la limită în formula

recurentă, obţinem

sau, după înlocuire

sau

,

de unde rezultă 0 = f(1) (numitorul fracţiei e diferit de 0, al doilea factor – la fel). Deci 1

este soluţia ecuaţiei iniţiale. Prin urmare, 1 = , ceea ce şi trebuia de demonstrat.Celelalte trei cazuri se demonstrează în mod analogic.

Algoritmizarea metodei

1. Determinăm care va fi extremitatea fixă a segmentului. În acest scop calculăm punctul de intersecţie (x1) a coardei ce trece prin f(a)) şi (b, f(b)) cu axa 0x. Drept extremitate fixă vom lua acel punct pentru care semnul funcţiei diferă de semnul f(x1).

2. Extremitatea “flotantă” primeşte valoarea calculată x1.3. Calculăm următoarea aproximaţie conform formulei:

4. Repetăm pasul 3 până obţinem soluţia cu exactitatea cerută.

Estimarea erorii

Faptul că şirul aproximărilor succesive prin metoda coardelor converge către soluţia exactă implică următoarea concluzie: cu cât mai multe iteraţii ale metodei vor fi realizate, cu atât mai bine va fi aproximată soluţia. Totuşi, am putea determina o formulă, care permite estimarea exactă a erorii de calcul, şi, prin urmare, exactitatea soluţiei obţinute.

Vom considera că prima derivată a funcţiei f(x) este continuă pe segmentul [a, b].Fie atunci m1 şi M1 două numere pozitive, pentru care are loc relaţia

8

Page 9: Capitolul 2[1]Metode numerice

Capitolul 2. Metode numerice

0 < m1 | f ’(x) | M1 < + .

Pentru simplitate, vom cerceta cazul, când aproximările succesive se fac după formula (2) (cazul, când e folosită formula (1) se demonstrează în mod analogic). Aşadar,

Exprimăm xn prin xn-1 ţinând cont de faptul că 0 = f(), pe care îl adăugăm la numărătorul fracţiei, apoi estimăm diferenţa dintre soluţia exactă şi cea aproximativă.

Pentru a deduce formula finală vom folosi următoarea teoremă (se demonstrează în cursul de analiză matematică)

Teorema Lagrange. Fie f:[a, b] R, continuă şi derivabilă pe [a, b] . Atunci c (a, b) astfel încât

f(b) – f(a) = (b -a) f’(c).

Vom aplica formula din teoremă aparte pentru partea stângă şi partea dreaptă a ultimei egalităţi:

Înlocuind, obţinem:

.

Deci, dacă e necesar să obţinem soluţia cu exactitatea dată , vom repeta aproximările până la respectarea inegalităţii

.

Aici xn şi xn+1 sunt două aproximaţii succesive, iar M1 şi m1, respectiv, marginea superioară şi inferioară a f ’(x) pe [a, b].

Exemplu. Elaboraţi un program ce calculează soluţia aproximativă a ecuaţiei

9

Page 10: Capitolul 2[1]Metode numerice

Capitolul 2. Metode numerice

pe intervalul [0,5; 1,5] pentru patru erori Е diferite, egale respectiv cu 0,1, 0,01, 0,001, 0,0001. Pentru fiecare valoare a erorii se va afişa la ecran o linie ce va conţine: eroarea E, soluţia aproximativă x şi numărul de iteraţii n.

Mai întîi verificăm posibilitatea de a aplica metoda secantelor pe intervalul propus:

f ’(x)=-2cos(x)sin(x) – ¼;f ’’(x) = -2 cos(2x);f ’’(x) = 0.

Rezultă că x = /4 + k/2.Unul din zerourile derivatei 2 este situat în intervalul în care se caută soluţia, deci pe

segmentul propus funcţia nu este monotonă, de acea vom împărţi mai întâi segmentul în 2 părţi, pe care le vom cerceta, apoi vom reduce segmentul la fragmentul în care se află soluţia, funcţia fiind monotonă:

[0,5 , /4]; f(0,5)= cos2(0,5) - 1/8 > cos2( /6) - 1/8 = ¾ - 1/8 >0;f( /4)= cos2( /4) - /16 = 1/2 - /16 >0.

Analiza ulterioară nu este necesară, întrucît este evident că soluţia se află pe segmentul [ /4, 1.5]. Pe acest segment funcţia este descrescătoare, iar derivata 2 - pozitivă. Deci extremitatea fixă va fi a = /4, iar formula recurentă

.

Mai rămâne să estimăm supremul şi infimul derivatei pe acest interval. Din expresia pentru derivata 1 rezultă că supremul nu poate depăşi 3/4, iar infimul valoarea 5/4.

Program P2;var e, a,b,c,c1,c2,m1,m2:real; i,n:integer;{calculăm valoarea funcţiei în punctul x}function f(x:real):real; begin f:=cos(x)*cos(x)-1/4; end;

begine:=1;for i:=1 to 4 dobegin{pentru fiecare epsilon repetăm iniţializarea datelor}a:=pi/4;b:=1.5; n:=0;m1:=3/4; m2:=-5/4;e:=e*0.1;c1:=b; c2:=a;while abs((m1-m2)/m2)*abs(c1-c2)>e do {atât timp cât nu a fost atnsă exactitatea} begin c:=c1-f(c1)/(f(c1)-f(a))*(c1-a); {calculăm următoarea aproximare} c2:=c1; c1:=c; inc(n);

10

Page 11: Capitolul 2[1]Metode numerice

Capitolul 2. Metode numerice

end; writeln('e=',e:6:4,' x =',c:8:6,' functia=',f(c):10:6,'

iteratii:', n);end; end.

Întrebări şi exerciţii

1. Care sînt avantajele şi dezavantajele metodei coardei?2. Daţi formularea matematică a metodei coardei.3. Elaboraţi o procedură care realizează metoda coardei.4. Formulaţi criteriul de convergenţă pentru metoda coardei.5. Care este condiţia de terminare a iteraţiilor în metoda coardei?6. Reprezentaţi pe un desen procesul de căutare a soluţiei prin metoda coardei.7. Estimaţi eroarea metodei coardei.8. Elaboraţi un program ce calculează prin metoda coardei soluţia aproximativă a

ecuaţiei

3 cos2 x2 = ex/2

pe intervalul [0,1; 1] pentru cinci erori ε diferite , egale respectiv cu 0,1; 0,02; 0,003; 0,0004; 0,00005. Pentru fiecare valoare a erorii se va afişa la ecran o linie ce va conţine: eroarea ε, soluţia aproximativă x şi numărul de iteraţii n.

9. Elaboraţi un program ce calculează prin metoda coardei soluţia aproximativă a ecuaţiei

pe intervalul [0, 1] cu eroarea 0,0001.10. Calculaţi prin metoda coardei soluţiile următoarelor ecuaţii:

a) ex + x – 2 = 0, x[-1, 1];b) x5 - 4x4 + x – 9 = 0, x[0, 4];c) ex – x –cos x – 1 = 0, x[0, 3].

11. Se cunoaşte faptul că funcţia y = x sin(x/2) are pe intervalul [3, 5] un maximum local unic. Elaboraţi un program cu care, fiind dată eroarea ε = 0,0001, se află valoarea xmax

pentru care funcţia ia valoarea maximă, utilizînd metoda coardelor pentru aflarea zerourilor derivatei. La ecran se va afişa xmax, ymax, şi numărul de iteraţii n.

12. Compuneţi un program ce calculează soluţia aproximativă a ecuaţiei

sin2 x = x/2

11

Page 12: Capitolul 2[1]Metode numerice

Capitolul 2. Metode numerice

pe intervalul [0,5, 1] pentru trei erori ε diferite, egale respectiv cu 0,1, 0,01, 0,001. În program nu veţi utiliza funcţiile trigonometrice standard ale limbajului Pascal. În calcule veţi folosi descompunerea în serie Taylor

sin x = x – x3/3! + x5/5! – x7/7! + ….

Valoarea funcţiei se va calcula cu aceeaşi exactitate. Pentru fiecare valoare a erorii se va afişa la ecran o linie ce va conţine eroarea ε, soluţia aproximativă x şi numărul de iteraţii n.

13. Compuneţi un program ce calculează soluţia aproximativă a ecuaţiei

x2 = sin 5x

pe intervalul [0,5, 0,6] pentru nouă erori ε diferite, egale respectiv cu 0,1, 0,02, 0,003,…, 0,000000009. În program nu veţi utiliza funcţiile trigonometrice standard ale limbajului Pascal. În calcule veţi folosi descompunerea în serie Taylor

sin x = x – x3/3! + x5/5! – x7/7! + ….

Valoarea funcţiei se va calcula cu aceeaşi exactitate. Pentru fiecare valoare a erorii se va afişa la ecran o linie ce va conţine eroarea ε, soluţia aproximativă x şi numărul de iteraţii n.

14. Elaboraţi un program ce va calcula soluţia aproximativă a ecuaţiei

x/4 = cos2 x

pe intervalul [0,5, 1,5] prin două metode: metoda coardelor şi metoda înjumătăţirii intervalului pentru valoarea erorii ε egală cu 0,000001. În program nu veţi utiliza funcţiile trigonometrice standard ale limbajului Pascal. În calcule veţi folosi descompunerea în serie Taylor

cos x = 1 - x2/2! + x4/4! – x6/6! +… .

Valoarea funcţiei se calculează cu aceeaşi exactitate. Programul va afişa la ecran eroarea ε, soluţia aproximativă x şi denumirea metodei prin care soluţia cu exactitatea cerută se calculează printr-un număr mai mic de iteraţii.

2.3. Metoda lui Newton

Descrierea matematică

Ideea generală a metodei este următoarea: prin punctul (b, f(b)) se duce o dreaptă tangentă la graficul functiei. Se determină punctul c în care ea intersectează axa 0x. Acest punct se considera noua extremitate, prin care se duce tangenta. Procesul se repetă, până obţinem o apropiere suficientă de soluţia exactă. Pentru această metodă vom cere suplimentar existenţa semnului constant al derivatei de ordin 2 pentru funcţia f(x) pe tot segmentul [a, b] sau, cu alte cuvinte, pe segmentul dat funcţia să fie sau concavă sau convexă.

12

Page 13: Capitolul 2[1]Metode numerice

Apropierea succesivă de soluţia ecuaţiei prin metoda tangentelor

Capitolul 2. Metode numerice

Pentru a determina extremitatea din care pornim, trebuie să ţinem cont de următoarele cazuri posibile:

1) f ’(x) > 0, f ’’(x) > 0; 2) f ’(x) > 0, f ’’(x) < 0;

3) f ’(x) < 0, f ’’(x) > 0; 4) f ’(x) < 0, f ’’(x) < 0.

Este adevărată afirmaţia: selectarea în calitate de punct iniţial al metodei tangentelor extremităţii pentru care se îndeplineşte relaţia f(x) · f ’’(x) > 0 permite determinarea soluţiei cu orice grad de exactitate e.

Vom cerceta cazul 1 – funcţia este convexă, crescătoare. (Celelalte cazuri se cercetează la fel). Deoarece a doua derivată e mai mare ca 0, graficul funcţiei este superior coardei. Pornind din extremitatea în care funcţia este pozitivă, vom obţine intersecţia cu axa 0x într-un punct interior al [a, b]. Vom considera acest punct noua extremitate a intervalului. Repetând procedura, vom obţine un şir de puncte interioare, care tind spre soluţia exactă. Dacă însă vom fixa extremitatea negativă, tangenta generată poate părăsi segmentul [a, b]. Deoarece funcţia este crescătoare, vom fixa extremitatea B.

Teoremă. Dacă f(a) · f(b )<0, f ’(x) şi f ’’(x) 0 şi îşi păstrează semnul pe [a, b], atunci, pornind de la o aproximare iniţială x0 [a, b] ( f (x0) · f ’’(x0) > 0), soluţia unică a ecuaţiei pe [a, b] poate fi calculată cu orice grad de exactitate .

Demonstraţie. Fie f(a) < 0, f(b) >0, f’(x) > 0 f ’’(x) > 0. Considerăm x0 = b, f(x0)>0. Pentru demonstrare vom folosi metoda inducţiei matematice, verificând ipoteza cum că toate valorile şirului xn sunt > , n.

Fie, pentru un oarecare n ipoteza este adevărată, xn > . Considerăm = xn + ( - xn), apoi aplicăm formula Tailor (pentru descompunerea funcţiei în sumă de polinoame) şi obţinem:

0 = f () = f(xn) + f’(xn) ( - xn) +1/2 f’’(cn) ( - xn)2 unde <cn < xn.

Deoarece f’’(x) > 0, rezultăf(xn) + f’(xn) ( - xn) < 0 sau f’(xn) ( - xn) < - f(xn) ( - xn) < - f(xn)/ f’(xn) < xn -

f(xn)/ f’(xn) = xn+1.

Din formula recurentă rezultă direct că x0 > x1 > … >xn …> 1 - limita xn.

Trecând la limită, obţinem

1 = 1 – f(1)/ f’(1),

de unde f(1) = 0, deci 1 = , ceea ce şi trebuia de demonstrat.

Algoritmizarea metodei

13

Page 14: Capitolul 2[1]Metode numerice

Capitolul 2. Metode numerice

1. Calculăm semnul derivatei 2 pe segmentul [a, b]. 2. Fixăm punctul iniţial x0 conform formulei: f(x0) · f’’(x0) > 0.3. Calculăm următoarea aproximaţie conform formulei:

4. Repetăm pasul 2 până obţinem soluţia cu exactitatea cerută.

Estimarea erorii

Faptul că şirul aproximărilor succesive prin metoda tangentelor converge către soluţia exactă implică următoarea concluzie: cu cât mai multe iteraţii ale metodei vor fi realizate, cu atât mai bine va fi aproximată soluţia. Totuşi, am putea determina o formulă, care ar permite estimarea exactă a erorii de calcul şi, prin urmare, exactitatea soluţiei obţinute.

Vom considera că prima derivată a funcţiei f(x) este continuă pe segmentul [a, b]. Fie atunci m1 şi M1 două numere pozitive, pentru care are loc relaţia

0 < m1 | f’(x) | M1 < + ;

.Exprimăm xn prin xn-1, ţinând cont de faptul că 0 = f(), şi îl adăugăm în partea stângă a

egalităţii, apoi estimăm diferenţa între soluţia exactă şi cea aproximativă. Pentru a deduce formula finală vom folosi teorema ce urmează.

Teorema Lagrange. Fie f : [a, b] R, continuă şi derivabilă pe [a, b] . Atunci c(a, b) astfel încât:

f(b) –f(a) = (b -a) f’(c).Prin urmare,

Înlocuind, obţinem:

Deci dacă e necesar să obţinem soluţia cu exactitatea dată , vom repeta aproximările până la respectarea inegalităţii

.

Aici xn şi xn+1 sunt două aproximaţii succesive, iar M1 şi m1 respectiv marginea superioară şi inferioară a derivatei f’(x) pe intervalul [a, b].

14

Page 15: Capitolul 2[1]Metode numerice

Capitolul 2. Metode numerice

Exemplul 1.

Elaboraţi un program ce calculează soluţia aproximativă a ecuaţiei

pe segmentul [0,5, 1] pentru 4 erori diferite, egale cu: 0,1, 0,01, 0,001, 0,0001. Pentru fiecare valoare a erorii pe ecran se va afişa o linie ce va conţine: eroarea e, soluţia aproximativă x, numărul de iteraţii n.

Rezolvare. Vom verifica mai întîi dacă metoda lui Newton poate fi aplicată pe intervalul dat. Pentru aceasta vom determina semnul derivatei 2 pe [0,5, 1]:

f(x) = sin2(x) – x/2;f’(x) = 2sin(x)cos(x) – 1/2 = sin(2x) – 1/2;f’’(x) = 2cos(2x);sin(2x) = 1/2; 2x = /6 sau 2x = 5/6 => x = /12 sau x = 5/12.

Zerouri pentru derivata 1 pe interval nu avem. Intervalul [0.5 , 1] se conţine în [/12, 5/12]. Pentru x=/4 f ’(x)=sin(/2) – ½ = ½ > 0. Deci funcţia este crescătoare pe acest interval.

cos(2x)=0.2x = /2 + k sau x = /4+ k/2.

Pe intervalul [0,5, 1] avem zeroul derivatei 2 în x = /4. Urmează să stabilim pe care din segmentele [0,5, /4] sau [/4, 1] se află soluţia ecuaţiei în studiu.

f(/4) = sin2(/4) – /8 = 1/2 – /8 >0.

Deoarece f(/4) >0 şi funcţia este crescătoare, soluţia trebuie căutată în intervalul [0,5, /4]. În calitate de punct iniţial vom selecta extremitatea /4, deoarece pentru ea se îndeplinesc condiţiile teoremei: f’’(/4) ·f(/4) >0.

Valorile M şi m pentru derivata 1 le vom stabili din expresia derivatei 1: M =1/2, m = -3/2, deoarece sin(x) variază între 1 şi –1.

Program P3;var x,e,a,b,c,c1,c2,m1,m2:real; i,n:integer;{x – valoarea curentă, e – eroarea, a,b – extremităţile segmentului, c1,c2 - 2 aproximări succesivem1, m2 – supremul şi infimul derivatei pe segment}function f(x:real):real; begin f:=sin(x)*sin(x) - x/2; end;

function fd1(x:real):real; begin fd1:=sin(2*x) - 1/2; end;

begine:=1;for i:=1 to 4 do {calculul se repetă pentru 4 aproximări succesive}beginb:=pi/4; a:=0.5; n:=0; {iniţializăm datele}m1:=1/2; m2:=-3/2;

15

Page 16: Capitolul 2[1]Metode numerice

Capitolul 2. Metode numerice

e:=e*0.1;c1:=b; c2:=a; c:=b;while abs((m1-m2)/m2)*abs(c1-c2)>e do {Repetăm calculul până la obţinerea exactităţii cerute} begin c:=c1-f(c1)/fd1(c1); c2:=c1; c1:=c; inc(n); end; writeln('e=',e:6:4,' x =',c:8:6,' functia=',f(c):10:8,' iteratii:', n);end;end.

Exemplul 2.

Elaboraţi un program ce determină dacă metodă lui Newton poate fi aplicată la rezolvarea ecuaţiei

x3 - 2x2 + x - 3 = 0

în punctul iniţial x0=2,2 şi care va rezolva, în caz afirmativ, ecuaţia dată pentru 9 erori diferite, egale cu: 0,1, 0,02, 0,003,... 0,0000000009. Pentru fiecare valoare a erorii se va afişa pe ecran o linie ce va conţine: eroarea e, soluţia aproximativă x, numărul de iteraţii n.

Rezolvare. Evident,

f(x) = x3 – 2x2 + x – 3;f’(x) = 3x2 – 4x + 1; x1=1/3, x2= 1 –zerourile derivatei 1f’’(x) = 6x – 4; x = 2/3 –zeroul derivatei 2.

În problema dată nu se indică segmentul pe care se caută soluţia. Vom începe prin o analiză matematică a problemei şi separarea soluţiilor:

f(1/3) < 1/3 + 1/3 + 1/3 +3 < 0;f(1) = 1 – 2 + 1 – 3 <0;f(2,2) = 2,22(2,2 -2) – 0,8 = 4,84 x 0,2 – 0,8 = 0,968 – 0,8 > 0.

Deci soluţia se află în intervalul [1,2, 2] si poate fi calculată prin metoda Newton, deoarece derivata 2 îşi păstrează semnul constant anume pe acest interval. Mai mult decât atât, observăm că şi f(2) = 1 < 0, deci restrîngem intervalul la [2, 2,2]. Verificarea prin program a posibilităţii începutului calculului prin metoda Newton constă în programarea condiţiilor teoremei respective.

Program r4;var x,e,e1, a,b,c,c1,c2,m1,m2:real; i,n:integer;{x – valoarea curentă, e – eroarea, a,b – extremităţile segmentului, c1,c2 - 2 aproximări succesivem1, m2 – supremul şi infimul derivatei pe segment}function f(x:real):real; begin f:=x*x*x - 2*x*x +x -3; end;function fd1(x:real):real; begin fd1:=3*x*x - 4*x +1; end;function fd2(x:real):real; begin fd2:=6*x - 4; end;

16

Page 17: Capitolul 2[1]Metode numerice

Capitolul 2. Metode numerice

beginx:=2.2;{verificarea posibilităţii de aplicare a metodei}if f(x)*fd2(x)>0 then begine:=1;for i:=1 to 9 do {Repetarea calculelor pentru următoarea eroare}beginb:=2.2; a:=2; n:=0; {iniţializarea datelor}m1:=fd1(2.2); m2:=fd1(2);e:=e*0.1; e1:=e*i;c1:=b; c2:=a; c:=b;while abs((m1-m2)/m2)*abs(c1-c2)>e1 do {cât nu a fost atinsă exactitatea cerută } begin c:=c1-f(c1)/fd1(c1); c2:=c1; c1:=c; inc(n); end; writeln('e=',e1:10:9,' x =',c:8:6,' functia=',f(c):10:6,' iteratii:', n);end;endelse writeln('Conditii initiale incompatibile');end.

Întrebări şi exerciţii

1. Care sînt avantajele şi dezavantajele metodei lui Newton?2. Daţi formularea matematică a metodei lui Newton.3. Elaboraţi o procedură care realizează metoda lui Newton.4. Formulaţi criteriul de convergenţă pentru metoda lui Newton.5. Care este condiţia de terminare a iteraţiilor în metoda lui Newton?6. Reprezentaţi pe un desen procesul de căutare a soluţiei prin metoda lui Newton.7. Estimaţi eroarea metodei lui Newton.8. Elaboraţi un program ce calculează soluţia aproximativă a ecuaţiei

3cos2x2 = ex/2

pe intervalul [0,1; 1] pentru cinci erori ε diferite , egale respectiv cu 0,1; 0,02; 0,003; 0,0004; 0,00005. Pentru fiecare valoare a erorii se va afişa la ecran o linie ce va conţine: eroarea ε, soluţia aproximativă x şi numărul de iteraţii n.

9. Elaboraţi un program un program ce calculează soluţia aproximativă a ecuaţiei

pe intervalul [0, 1] cu eroarea 0,0001.10. Calculaţi soluţiile următoarelor ecuaţii:

a) ex + x – 2 = 0, x[-1, 1];

17

Page 18: Capitolul 2[1]Metode numerice

Capitolul 2. Metode numerice

b) x5 - 4x4 + x – 9 = 0, x[0, 4];c) ex – x –cos x – 1 = 0, x[0, 3].

11. Se cunoaşte faptul că funcţia y = x sin(x/2) are pe intervalul [3, 5] un maximum local unic. Elaboraţi un program care, avînd eroarea ε = 0,0001, află valoarea xmax pentru care funcţia ia valoarea maximă, utilizînd metoda coardelor pentru aflarea zerourilor derivatei. La ecran se va afişa xmax, ymax, şi numărul de iteraţii n.

12. Compuneţi un program ce calculează soluţia aproximativă a ecuaţiei

sin2 x = x/2

pe intervalul [0,5, 1] pentru trei erori ε diferite, egale respectiv cu 0,1, 0,01, 0,001. În program nu veţi utiliza funcţiile trigonometrice standard ale limbajului Pascal. În calcule veţi folosi descompunerea în serie Taylor

sin x = x – x3/3! + x5/5! – x7/7! + ….

Valoarea funcţiei se va calcula cu aceeaşi exactitate. Pentru fiecare valoare a erorii se va afişa la ecran o linie ce va conţine eroarea ε, soluţia aproximativă x şi numărul de iteraţii n.

13. Compuneţi un program ce calculează soluţia aproximativă a ecuaţiei

x2 = sin 5x

pe intervalul [0,5, 0,6] pentru nouă erori ε diferite, egale respectiv cu 0,1, 0,02, 0,003,…, 0,000000009. În program nu veţi utiliza funcţiile trigonometrice standard ale limbajului Pascal. În calcule veţi folosi descompunerea în serie Taylor

sin x = x – x3/3! + x5/5! – x7/7! + ….

Valoarea funcţiei se va calcula cu aceeaşi exactitate. Pentru fiecare valoare a erorii se va afişa la ecran o linie ce va conţine eroarea ε, soluţia aproximativă x şi numărul de iteraţii n.

14. Elaboraţi un program care va determina dacă poate fi aplicată metoda lui Newton pentru rezolvarea ecuaţiei

x3 - 2x2 + x –3 = 0

pornind din punctul iniţial x0 = 2,2 şi care va rezolva, în caz afirmativ, ecuaţia dată pentru nouă erori ε diferite, egale respectiv cu 0,1, 0,02, 0,003,…, 0,000000009. Valoarea funcţiei se va calcula cu aceeaşi exactitate. Pentru fiecare valoare a erorii la ecran se va afişa o linie ce va conţine eroarea ε, soluţia aproximativă x şi numărul de iteraţii n.

2.4. Noţiunea de determinant. Proprietăţile de bază ale determinanţilor

Vom cerceta o matrice pătratică arbitrară de rang n:

18

Page 19: Capitolul 2[1]Metode numerice

Capitolul 2. Metode numerice

.

Fiecărei din matricele de acest tip îi vom asocia o valoare numerică, numită determinant. Vom defini această mărime în mod inductiv.

Fie că rangul matricei A este 1. Matricea este formată dintr-un singur element a11. Determinantul matricei A va fi chiar valoarea elementului a11.

.

Pentru o matrice de rangul 2 determinantul va fi egal cu valoarea expresiei

a11a22 – a12a21,

adică cu diferenţa produselor elementelor de pe diagonala principală şi cea secundară.Pentru a defini noţiunea de determinant al unei matrice de rang n vom avea nevoie de

următoarele noţiuni.Vom numi minor de ordin n-1 al elementului aij al matricei de rang n A determinantul

de ordin n-1, care se obţine din matricea A la excluderea din ea a rândului i şi a coloanei j. Vom nota minorul elementului aij prin , unde i indică rândul, iar j coloana care corespund elementului aij.

Vom numi determinant al matricei A de rang n valoarea expresiei

care se va marca prin simbolul

.

Următoarele proprietăţi ale determinanţilor permit calculul lor fără a recurge la calcularea minorilor algebrici.

Proprietatea 1: Dacă de o parte a diagonalei principale a determinantului sunt numai zerouri, atunci valoarea determinantului este produsul elementelor de pe diagonala principală.

Proprietatea 2: Dacă schimbăm locurile a două linii (coloane) ale determinantului, semnul acestuia se schimbă în opus.

Proprietatea 3: Dacă adăugăm la o linie a determinantului o altă linie înmulţită cu un număr, valoarea determinantului nu se schimbă.

Aceste trei proprietăţi permit elaborarea unui algoritm simplu pentru calculul determinanţilor.

19

Page 20: Capitolul 2[1]Metode numerice

Capitolul 2. Metode numerice

Scopul final este de a aduce determinantul la forma triunghiulară (adică de o parte a diagonalei principale să fie doar zerouri). Pentru a realiza aceasta, vom folosi iterativ proprietatea 3.

Fie pe linia i elementul aii 0. Dacă înmulţim linia i cu – aki/aii şi o adăugăm la linia k, valoarea determinantului nu se va schimba. În acelaşi timp, elementul aki va deveni 0. Într-adevăr, pentru el vom avea aki + aii (– aki/aii ) = aki + (- aki ) = 0. De aici rezultă, în cazul aplicării consecutive a acestei proceduri pentru toţi k > i, apariţia zerourilor în coloana i a determinantului mai jos de diagonala principală. De aici rezultă şi algoritmul integral:

1. Numărul liniei (i) – 0.2. Numărul liniei se măreşte cu 1 (i:=i+1).3. Dacă elementul aii 0, atunci atribuim j:=i+1 şi pentru fiecare din liniile (j) de la

i+1 la n repetăm procedura (P):

(P) linia i o înmulţim cu – aji/aii şi o adăugăm la linia j. Dacă însă aii = 0, atunci mai întîi cercetăm elementul aji pe liniile situate sub linia i (coloana i). Dacă detectăm un element nenul (fie în linia k), schimbăm locurile liniilor i şi k, apoi realizăm procedura (P). Dacă toate elementele în coloana i situate mai jos de rîndul i sînt 0, atunci valoarea determinantului este 0, corespunzător se încheie lucrul algoritmului.

4. Dacă i< n revenim la pasul 1.5. Calculăm produsul elementelor de pe diagonala principală.

Întrebări şi exerciţii

1. Explicaţi noţiunile de minor şi determinant.2. Formulaţi proprietăţile de bază ale determinanţilor.3. Elaboraţi un algoritm şi procedura respectivă pentru calculul determinanţilor.4. Cîte operaţii aritmetice se efectuează pentru calculul unui determinant de ordin n?5. Calculaţi următorii determinanţi:

a) b) c)

d) e) f)

2.5. Metoda lui Cramer de rezolvare a sistemelor de ecuaţii liniare

Definiţie. Matricea inversă matricei A – matricea, care fiind înmulţită la A (fie la stânga, fie la dreapta), dă în rezultat matricea unitară E. Matricea inversă matricei A va fi notată A-1.

20

Page 21: Capitolul 2[1]Metode numerice

Capitolul 2. Metode numerice

.

Definiţie. Matricea pătratică A se numeşte nesingulară, dacă determinantul ei este diferit de 0.

Teoremă. Pentru orice matrice nesingulară există o matrice inversă.

Metoda Cramer de rezolvare a sistemelor de ecuaţii liniare

Fie dat sistemul din n ecuaţii liniare cu n necunoscute:

.

Vom nota prin matricea coeficienţilor liberi ai matricei

sistemului,

prin - vectorul termenilor liberi ai sistemului, iar prin - vectorul

soluţiilor.

În formă vectorială sistemul va avea forma Ax = b. În cazul când determinantul A este diferit de 0, există matricea inversă A-1. Înmulţind ambele părţi ale egalităţii precedente la A-1, obţinem:

A-1Ax = A-1 b sau x = A-1 b,

ceea ce permite calculul soluţiilor sistemului iniţial. În cele ce urmează vom detalia formula obţinută, determinând formulele de calcul pentru componentele vectorului x.

Vom folosi matricea adjunctă matricei A:

în care Aij este minorul cu semn al elementului aij din matricea A. Observăm că minorii elementelor din linia i vor forma coloana cu acelaşi indice, adică are loc operaţia de transpunere a minorilor elementelor matricei A.

Este cunoscut că

,

prin urmare

21

Page 22: Capitolul 2[1]Metode numerice

Capitolul 2. Metode numerice

sau, la nivel de componente,

,

unde

(*).

Într-adevăr, produsul unei matrice cu un vector coloană este un vector, componenta i a căruia se calculează ca suma produselor elementelor liniei i a matricei cu componentele respective ale vectorului coloană. Deoarece elementele liniei i a matricei sunt minorii matricei A, descompuse după coloana i, suma produselor lor şi a elementelor vectorului b vor corespunde descompunerii după coloana i a matricei (*). De aici urmează direct formulele de calcul a soluţiilor sistemului iniţial:

.

Algoritmizarea metodei.

1. Verificăm dacă determinantul al matricei A a coeficienţilor sistemului este diferit de 0. Dacă da, trecem pa pasul 2, în caz contrar – afişăm mesajul de imposibilitate de calcul, apoi SFÂRŞIT.

2. i:=1.3. Pentru coloana activă (cu numărul i) facem o copie (Acop) a matricei coeficienţilor

sistemului în care înlocuim coloana i cu vectorul termenilor liberi.4. Calculăm determinantul matricei Acop, memorăm valoarea lui în componenta i a

vectorului x.5. Dacă i < n, incrementăm i şi revenim la pasul 3, în caz contrar trecem la pasul 66. Componentele vectorului x conţin valorile i . Calculăm soluţia împărţind

componentele x la .7. Afişăm soluţia. SFÂRŞIT.

Structuri de date

1. Doua tablouri bidimensionale cu n linii şi n coloane: primul va păstra matricea coeficienţilor A, al doilea – copia ei cu coloana i înlocuită prin vectorul termenilor liberi.

2. Două tablouri unidimensionale din n elemente reale – pentru păstrarea vectorului termenilor liberi şi vectorului soluţiei obţinute.

Întrebări şi exerciţii

1. Explicaţi termenii matrice inversă, matrice nesingulară şi matrice singulară.2. Elaboraţi un algoritm şi procedura respectivă pentru rezolvarea sistemelor de ecuaţii

liniare prin metoda lui Cramer.

22

Page 23: Capitolul 2[1]Metode numerice

Capitolul 2. Metode numerice

4. Cîte operaţii aritmetice se efectuează pentru rezolvarea prin metoda lui Cramer a unui sistem de ecuaţii liniare de ordin n?

5. Se consideră că este dat modulul SYSEQLIN care conţine funcţia Determinant. Această funcţie returnează valoarea determinantului ca rezultat de tip real. Parametrii funcţiei sînt: matricea pătratică A de tip real şi variabila n de tip integer – dimensiunea matricei. Elaboraţi un program pentru rezolvarea sistemului de n ecuaţii liniare cu n necunoscute folosind regula lui Cramer. Coeficienţii sistemului de ecuaţii se introduc de la tastatură. Soluţia sistemului se va afişa la ecran.

Notă. Nu se cere elaborarea algoritmul propriu-zis de calculare a determinanţilor.6. Scrieţi o funcţie cu numele Det care are drept argument o matrice din 3 linii şi 3

coloane cu elemente de tip real. Funcţia returnează ca rezultat de tip real determinantul matricei.

7. Rezolvaţi prin metoda lui Cramer următoarele sisteme de ecuaţii liniare:

a) b) c)

2.6. Metoda lui Gauss de rezolvare a sistemelor de ecuaţii liniare

Fie dat sistemul din n ecuaţii liniare cu n necunoscute:

.

Sistemul are o soluţie unică doar în cazul în care determinantul matricei coeficienţilor A este diferit de 0. Evident, determinarea acestui fapt necesită calculul determinantului matricei A. Metoda Gauss constă în eliminarea succesivă a necunoscutelor din sistem, până când în ultima ecuaţie nu va rămâne doar o necunoscută (pasul direct). Din această ecuaţie se află necunoscuta, care apoi se înlocuieşte în ecuaţie cu numărul n-1. Procesul se repetă până nu se calculează ultima componentă necunoscută din vectorul soluţiilor (pasul invers).

Pasul directFie a11 0 (dacă a11=0, se alege o altă ecuaţie (i), în care ai1 0). Vom înmulţi

coeficienţii primei ecuaţii la -a21/a11 şi vom aduna ecuaţia obţinută la ecuaţia a doua din sistem (spre deosebire de determinanţi, aici se va transforma şi termenul liber al ecuaţiei 2). Apoi vom înmulţi coeficienţii primei ecuaţii la –a31/a11 şi vom aduna ecuaţia obţinută la ecuaţia cu indicele (3) ... Procesul se repetă până la ultima ecuaţie din sistem, pentru care ecuaţia (1) se înmulţeşte cu -an1/a11.

În rezultat obţinem un sistem nou.

.În acest sistem începem eliminarea necunoscutei x2 din toate ecuaţiile, începând cu a

treia. Pentru aceasta vom repeta algoritmul pasului precedent, dar vom porni de la elementul

23

Page 24: Capitolul 2[1]Metode numerice

Capitolul 2. Metode numerice

a(1)22 Coeficienţii ecuaţiei 2 vor fi pe rând înmulţiţi la - a(1)

k2 / a(1)22, iar ecuaţia obţinută se va

aduna la ecuaţia cu numărul de ordine k. Vom obţine sistemul

.

În caz general algoritmul este următorul. Fie pe linia i elementul aii 0. Dacă înmulţim linia i cu – aki/aii şi o adăugăm la linia k, soluţia sistemului nu se va schimba. În acelaşi timp, elementul aki va deveni 0. Într-adevăr, pentru el vom avea

aki + aii (– aki/aii ) = aki + (- aki ) = 0.

De aici rezultă, în cazul aplicării consecutive a acestei proceduri pentru toţi k > i, apariţia zerourilor în coloana i a ecuaţiilor sistemului situate mai jos de ecuaţia cu indicele i. Tot în felul acesta rezultă şi algoritmul integral:

1. Numărul liniei (i) – 0.2. Numărul liniei se măreşte cu 1 (i := i+1).3. Dacă elementul aii 0, atunci considerăm j := i+1 şi pentru fiecare din liniile (j) de la i+1

la n repetăm procedura (P):

(P) linia i o înmulţim cu – aji/aii şi o adăugăm la linia j (inclusiv termenii liberi). Dacă însă aii = 0, atunci mai întîi cercetăm elementul aji pe liniile situate sub linia i (coloana i). Dacă detectăm un element nenul (fie în linia k), schimbăm locurile liniilor i şi k, apoi realizăm procedura (P). Dacă toate elementele în coloana i situate mai jos de rîndul i sînt 0, atunci valoarea determinantului care corespunde acestui sistem de ecuaţii este 0 şi respectiv soluţia unică nu există.

4. Dacă i< n revenim la pasul 1.

În final, vom obţine sistemul:

.

Aici se încheie pasul direct al metodei Gauss.

Pasul înapoi (Etapa indirectă)

24

Page 25: Capitolul 2[1]Metode numerice

Capitolul 2. Metode numerice

Din ultima ecuaţie determinăm xn:

xn = b(n-1)n / a(n-1)

nn.

Înlocuind în ecuaţia n-1, obţinem:

În caz general

Conform formulei obţinute toate componentele vectorului de soluţii x se calculează recurent, începând cu xn

Pentru realizarea practică a metodei Gauss vor fi necesare următoarele resurse de memorie:

1. Tabloul bidimensional cu n linii şi n+1 coloane : în fiecare linie coloana n+1 conţine termenii liberi ai ecuaţiilor, primele n – coeficienţii de pe lângă vectorul x în ecuaţia respectivă. Deoarece în procesul transformărilor se realizează împărţirea, tipul datelor va fi real, chiar dacă iniţial coeficienţii sunt întregi.

2. Tablou unidimensional din n elemente reale – pentru păstrarea soluţiei obţinute.

Întrebări şi exerciţii

1. Explicaţi termenii pasul direct şi pasul înapoi.2. Elaboraţi un algoritm şi procedura respectivă pentru rezolvarea sistemelor de ecuaţii

liniare prin metoda lui Gauss.4. Cîte operaţii aritmetice se efectuează pentru rezolvarea prin metoda lui Gauss a unui

sistem de ecuaţii liniare de ordin n?5. Elaboraţi un program de rezolvare a unui sistem de n ecuaţii liniare cu n necunoscute

pentru cazul, când coeficienţii ce se află sub diagonala principală sunt egali cu zero, iar cei situaţi pe diagonala principală sunt diferiţi de zero. Coeficienţii sistemului se introduc de la tastatură. Programul va afişa soluţia sistemului la ecran.

6. Rezolvaţi prin metoda lui Gauss următoarele sisteme de ecuaţii liniare:

a) b) c)

25

Page 26: Capitolul 2[1]Metode numerice

Capitolul 2. Metode numerice

2.7. Formula dreptunghiurilor pentru calculul aproximativ al integralei

Una dintre cele mai des aplicate implementări ale metodelor numerice este calcularea integralei definite prin metode aproximative. Metodele directe nu întotdeauna permit calculul analitic al integralei şi, de multe ori, formula care defineşte funcţia ce trebuie integrată nici nu e cunoscută. De obicei sunt date doar o serie de puncte în care este cunoscută valoarea funcţiei. În aceste cazuri integrala poate fi calculată doar prin metode aproximative (în condiţia în care funcţia de sub integrală este continuă pe segmentul pe care se operează integrarea).

Din cursul de analiză matematică se ştie că sensul geometric al integralei definite este aria trapezului curbiliniu determinat de axa 0X, dreptele y = a şi y = b, şi graficul funcţiei f(x) pe segmentul [a, b].

Aparatul matematic. Fie f o funcţie derivabilă pe [a, b] şi

.

Considerăm o diviziune a segmentului [a, b] în forma

a = x0 < x1 < x2 < … < xn-1 < xn = b.

Pe fiecare din segmentele [xi, xi+1] vom aproxima suprafaţa trapezului curbiliniu cu aria dreptunghiului cu baza h = |xi, - xi+1| şi înălţimea egală cu valoarea funcţiei în unul din punctele:

a) x = xi – atunci metoda se numeşte a dreptunghiurilor de stânga;b) x = xi+1 – atunci metoda se numeşte a dreptunghiurilor de dreapta;c) x=(xi,+ xi+1)/2 – atunci metoda se numeşte a dreptunghiurilor medii.

Pentru simplitate, vom considera segmentele [xi, xi+1] egale. Atunci lungimea segmentului va fi de

h = |b - a| / n,

iar punctele xi vor putea fi determinate conform formulei xi =a + i h.Dreapta

y = (xi,+ xi+1)/2

26

Page 27: Capitolul 2[1]Metode numerice

Capitolul 2. Metode numerice

(corespunzător y = xi pentru metoda dreptunghiurilor de stânga şi y = xi+1 pentru metoda dreptunghiurilor de dreapta) este graficul unei funcţii constante (a polinomului de interpolare Lagrange de grad 0, care aproximează funcţia pe segmentul elementar [xi, xi+1] prin mijlocul lui (corespunzător prin extremitatea stângă – pentru metoda dreptunghiurilor de stânga şi prin extremitatea dreaptă - pentru metoda dreptunghiurilor de dreapta).

În acest caz valoarea integralei definite pe segmentul [a, b] pentru funcţia f se va aproxima conform formulei

Pentru calculul numeric transformarea integralei într-o expresie aritmetică (sumă) care depinde doar de numărul de diviziuni ale segmentului (acest număr sau e dat apriori sau e stabilit în procesul de calcul) şi de valoarea funcţiei în nodurile de interpolare (valoare care este de asemenea dată sau poate fi calculată) este foarte comodă – procedura de calcul devine elementară.

Determinarea erorii de calcul (marginea superioară). Eroarea de interpolare pe un interval elementar se estimează ca modulul diferenţei dintre valoarea exactă a funcţiei pe interval şi funcţia de interpolare. Corespunzător, eroarea de la integrare prin funcţia de interpolare va fi integrala erorii de interpolare. Pe un segment elementar eroarea de interpolare este:

| f(z) - g(z) | M * (z –(xi+1 + xi )/2) ),

unde(xi+1 + xi ) / 2 – punctul de interpolare, iar g(z) este o funcţie constantă:g(z) = f((xi+1 + xi )/2) – pentru metoda mediilor;

27

Page 28: Capitolul 2[1]Metode numerice

Capitolul 2. Metode numerice

g(z)= f(xi) – pentru metoda dreptunghiurilor de stânga;g(z)= f(xi+1) – pentru metoda dreptunghiurilor de dreapta.

Estimarea erorii o facem prin integrarea diferenţei | f(z) - g(z) | utilizând înlocuirea variabilei x prin t:

x = x(t) = (xi+1 + xi )/2 + t*(xi+1 - xi ) / 2.

Prin urmare, eroarea totală va fi

.

Formula dată ne permite să stabilim apriori care este numărul necesar de iteraţii pentru obţinerea unei erori mai mici decât cea prestabilită:

O metodă mai simplă, care ne fereşte de riscul obţinerii unor valori din afara diapazonului admisibil, este cea bazată pe divizarea binară a intervalului de iteraţii. Vom selecta o valoare maximal posibilă pentru n2 – numărul iniţial de iteraţii şi 0 pentru n1 – numărul minim posibil de iteraţii. Deoarece eroarea este o expresie ce depinde de numărul de iteraţii, calculăm eroarea pentru

28

Page 29: Capitolul 2[1]Metode numerice

Capitolul 2. Metode numerice

n=(n1 + n2) div 2,

cea ce reprezintă un număr mediu de iteraţii. Dacă eroarea obţinută e este mai mică decât cea cerută, considerăm n2 <- n (micşorăm numărul maxim de iteraţii), în caz contrar n1 <- n (mărim numărul minim de iteraţii).

Acesta e un proces de divizare binară şi va continua atât timp cât n1 < n2. Va fi nevoie de log2 (maxlongint div 2) verificări pentru a stabili exact numărul minim necesar de iteraţii pentru obţinerea exactităţii date.

Algoritmizarea problemei – varianta mediilor.

1. Se introduc limitele de integrare a, b.2. Se stabileşte numărul necesar de divizări n.3. Se calculează pasul de deplasare h.4. Pornind de la a calculăm mijlocurile segmentelor elementare zi şi ariile dreptunghiurilor

elementare.5. Însumăm ariile elementare.6. Afişăm rezultatul.

Algoritmizarea problemei – varianta dreptunghiurilor de stânga

1. Se introduc limitele de integrare a, b.2. Se stabileşte numărul necesar de divizări n.3. Se calculează pasul de deplasare h.4. Pornind de la a calculăm extremităţile stângi ale segmentelor elementare zi şi ariile

dreptunghiurilor elementare.5. Însumăm ariile elementare.6. Afişăm rezultatul.

Algoritmizarea problemei – varianta dreptunghiurilor de dreapta1. Se introduc limitele de integrare a, b.2. Se stabileşte numărul necesar de divizări n.3. Se calculează pasul de deplasare h.4. Pornind de la a calculăm extremităţile drepte ale segmentelor elementare zi şi ariile

dreptunghiurilor elementare.5. Însumăm ariile elementare.6. Afişăm rezultatul.

Exemplul 1. Elaboraţi un program ce calculează integrala

aplicând formula dreptunghiurilor (metoda mediilor) pentru diferite numere de diviziuni ale intervalului de integrare n=10, n=100, n=1000. Programul va afişa la ecran valoarea integralei I pentru fiecare număr de diviziuni N.

29

Page 30: Capitolul 2[1]Metode numerice

Capitolul 2. Metode numerice

Rezolvare.

Program P4;var x1,S,a,b,h,M :real; j,i,n : integer;function f(x:real):real; begin f:=x*x*sqrt(1-x*x*x); end; begin b:=1; a:=0; n:=10; h:=(b-a)/n; for j:=1 to 3 do {contor utilizat pentru schimabarea numarului de divizari} begin s:=0; for i:=0 to n-1 do {calculul ariilor elementare si sumarea lor} begin x1:=a+i*h+h/2; s:=s+h*f(x1); end; {afisarea rezultatelor pentru numarul dat de divizari} writeln ('n=',n,' s=',s:8:4); n:=n*10; {schimabarea numarului divizarilor} h:=(b-a)/n; {recalcularea pasului} end;end.

Exemplul 2. Elaboraţi un program, ce calculează aria figurii mărginite de liniile ce nu se intersectează definite de funcţiile y=sin(x) + 2, y=x2 + 1 si de dreptele x = 0 si x = 1. Calculul se va efectua utilizând formula dreptunghiurilor de stânga cu exactitatea e = 0,01. La ecran se va afişa valoarea S a ariei întregii figuri şi numărul de diviziuni n a intervalului de integrare.

Rezolvare. In problema data se cere sa se calculeze 2 integrale, apoi sa se calculeze diferenta lor. Aria căutată va fi

.

Menţionăm că în condiţiile problemei este prezent un moment de nedeterminare: exactitatea cerută se poate obţine în un număr diferit de iteraţii pentru fiecare din funcţiile date. Din acest motiv vom afişa numărul de iteraţii pentru fiecare din integralele calculate. Mai întîi vom determina M pentru derivata funcţiilor date pe segmentul propus:

y=sin(x)+2 y=cos(x) 1 pe segmentul [0, 1];y=x2+1 y=2x 2 pe segmentul [0, 1].

Vom considera M1=1 si M2=2.

Program P5;

30

Page 31: Capitolul 2[1]Metode numerice

Capitolul 2. Metode numerice

var x1,S1,s2,a,b,h,M1,m2,e :real; j,i,n,n2 : integer;

function f1(x:real):real; {prima functie integrabila} begin f1:=sin(x)+2; end;

function f2(x:real):real; {a doua functie integrabila} begin f2:=1+x*x; end;

begin {initializarea procesului de calcul. Datele le fixam direct in program, deoarece nu sunt puse restrictii asupra introducerii lor} b:=1; a:=0; e:=0.01; M1:=1; n:=2; h:=(b-a)/n;

{calculul primei integrale} while (b-a)*M1*h >e do {atit timp cat nu e atinsa exactitatea recalculam integrala} begin s1:=0; for i:=0 to n-1 do {calculul ariilor elementare si sumarea lor} begin x1:=a+i*h; s1:=s1+h*f1(x1); end; n:=n+1; {schimabarea numarului divizarilor} h:=(b-a)/n; {recalcularea pasului} end;

{initializarea datelor pentru calculul integralei 2:} M2:=2; n2:=2; h:=(b-a)/n2; while (b-a)*M2*h >e do {atit timp cat nu e atinsa exactitatea recalculam integrala} begin s2:=0; for i:=0 to n2-1 do {calculul ariilor elementare si sumarea lor} begin x1:=a+i*h; s2:=s2+h*f2(x1); end; n2:=n2+1; {schimabarea numarului divizarilor} h:=(b-a)/n2; {recalcularea pasului} end; writeln ('numarul de divizări pentru prima componenta ',n); writeln ('numarul de divizări pentru componenta doi ',n2); writeln ('aria figurii ',abs(s1-s2));end.

Întrebări şi exerciţii

1. Daţi formularea matematică a metodei dreptunghiurilor (dreptunghiuri de stînga, dreptunghiuri medii, dreptunghiuri de dreapta).

2. Elaboraţi o procedură care realizează fiecare din variantele metodei dreptunghiurilor: de stînga, medii, de dreapta.

3. Care este eroarea de calcul a metodei dreptunghiurilor?

31

Page 32: Capitolul 2[1]Metode numerice

Capitolul 2. Metode numerice

4. Se consideră figura mărginită de graficul funcţiei y = |x| cos x + 1, axa 0x şi de dreptele x = -2 şi x = 2. Elaboraţi un program care va determina poziţia dreptei verticale, care împarte figura examinată în două părţi, ariile cărora, respectiv din stânga şi din dreapta, se raportă ca 3:1. Calculele se vor efectua utilizând formula dreptunghiurilor de stânga cu exactitatea ε = 0,001. La ecran se va afişa valoarea S a ariei întregii figuri şi coordonata punctului b de pe axa 0x prin care trece dreapta verticală căutată.

5. Compuneţi un program ce calculează integrala

utilizînd formula dreptunghiurilor (metoda mediilor) cu exactitatea ε = 0,01. Programul va afişa la ecran valoarea integralei I şi numărul de diviziuni n pentru care se realizează exactitatea dată.

5. Elaboraţi un program ce calculează integrala

utilizînd formula dreptunghiurilor de dreapta cu exactitatea ε = 0,001. Programul va afişa la ecran valoarea integralei I şi numărul de diviziuni n pentru care se realizează exactitatea dată.

6. Elaboraţi un program ce calculează integrala

utilizînd formula dreptunghiurilor de stânga pentru diferite numere de diviziuni ale intervalului de integrare n =10, 100, 1000. Programul va afişa la ecran valoarea integralei I pentru fiecare număr de diviziuni n.

7. Elaboraţi un program ce calculează aria figurii mărginite de liniile ce nu se intersectează definite de funcţiile y = sin x + 2, y = x2 + 1 şi de dreptele x = 0 şi x = 1. Calculul se va efectua utilizând formula dreptunghiurilor de stânga cu exactitatea ε = 0,01. La ecran se va afişa valoarea S a ariei figurii şi numărul de diviziuni n ale intervalului de integrare.

8. Elaboraţi un program ce calculează aria figurii mărginite de graficul funcţiei y = e|x|, axa 0x şi de dreptele x = -2 şi x = 2. În program nu se va utiliza funcţia standard a limbajului Pascal exp(x). În calcule se va folosi descompunerea în serie Taylor

ex = 1 + x/1! + x2/2! + x3/3 + x4/4! + … .

Calculele se vor efectua utilizând formula dreptunghiurilor de stânga cu exactitatea ε = 0,001. Cu aceeaşi exactitate se va calcula şi valoarea funcţiei. La ecran se vor afişa valoarea S a ariei figurii şi numărul de diviziuni n ale intervalului de integrare.

9. Elaboraţi un program ce calculează aria figurii mărginite de graficul funcţiei y = |x| cos x + 1, axa 0x şi de dreptele x = -2 şi x = 2. În program nu veţi utiliza funcţiile trigonometrice standard ale limbajului Pascal. În calcule veţi folosi descompunerea în serie Taylor

cos x = 1- x2/2! + x4/4! - x6/6! +… .

32

Page 33: Capitolul 2[1]Metode numerice

Capitolul 2. Metode numerice

Calculele se vor efectua utilizând formula dreptunghiurilor de stânga cu exactitatea ε = 0,001. Cu aceeaşi exactitate se va calcula şi valoarea funcţiei. La ecran se va afişa valoarea S a ariei figurii şi numărul de diviziuni n ale intervalului de integrare.

10. Se consideră figura mărginită de graficul funcţiei y = sin3x + 1, axa 0x şi de dreptele x = -2 şi x = 1. Elaboraţi un program care va determina poziţia dreptei verticale, care împarte figura examinată în două părţi de arie egală. Calculele se vor efectua utilizând formula dreptunghiurilor de dreapta cu exactitatea ε = 0,001. La ecran se va afişa valoarea S a ariei întregii figuri şi coordonata punctului b de pe axa 0x prin care trece dreapta verticală căutată.

11. Compuneţi un program ce calculează integrala

utilizînd formula dreptunghiurilor de stânga cu exactitatea ε = 0,001. Programul va afişa la ecran valoarea integralei I şi numărul de diviziuni n pentru care se realizează exactitatea dată. În program nu se vor utiliza funcţiile trigonometrice standard ale limbajului Pascal. În calcule se vor folosi descompunerile în serie Taylor

sin x = x – x3/3! + x5/5! – x7/7! +… ;cos x = 1 - x2/2! + x4/4! – x6/6! +… .

Valorile funcţiilor se vor calcula cu aceeaşi exactitate.

2.8. Formula trapezelor pentru calculul integralei

În paragraful precedent am folosit pentru aproximarea funcţiei integrale funcţii constante, care coincid cu valoarea funcţiei într-un singur punct (începutul, mijlocul sau sfârşitul segmentului elementar, pe care se operează aproximarea). Este logică dezvoltarea metodei prin interpolarea mai exactă a funcţiei de sub integrală, de exemplu cu polinoame de interpolare Lagrange de grad 1 (funcţia de interpolare va coincide cu funcţia integrală în 2 puncte). În calitate de noduri de interpolare vom folosi extremităţile segmentului elementar [xi, xi+1] pe care realizăm interpolarea.

Grafic procesul de aproximare a funcţiei îl putem urmări în figura alăturată. Menţionăm că funcţia de aproximare va fi o funcţie liniară ce va trece prin punctele (xi, f(xi), xi+1, f(xi+1)).

33

Page 34: Capitolul 2[1]Metode numerice

Capitolul 2. Metode numerice

Formularea matematică a metodei

Pentru a avea posibilitatea de a aplica metoda trapezelor, vom cere suplimentar ca funcţia integrală să fie de două ori derivabilă.

Fie f o funcţie de două ori derivabilă pe [a, b] şi .

Considerăm o diviziune a intervalului [a, b] în forma a = x0 < x1 < x2 < … < xn-1 < xn

= b.Pe fiecare din segmentele [xi, xi+1] vom aproxima suprafaţa trapezului curbiliniu cu aria

trapezului cu baza h = |xi, - xi+1| , laturile laterale cu lungimile f(xi,) şi f(xi+1), iar latura superioară dată de segmentul cu extremităţile (xi, f(xi)) şi (xi+1, f(xi+1)).

Ca şi în cazul formulei dreptunghiurilor vom considera segmentele [xi, xi+1] egale. Lungimea segmentului elementar va fi de h = |b - a|/n, iar punctele xi vor putea fi determinate de formula xi =a + i * h.

Atunci valoarea integralei definite pe un segment elementar [xi, xi+1] pentru funcţia f(x) se va aproxima conform formulei

,

iar pentru întreg segmentul [a, b] vom obţine

Determinarea erorii de calcul (marginea superioară)

Eroarea de interpolare pe un interval elementar se estimează ca modulul diferenţei între valoarea exactă a funcţiei pe interval şi funcţia de interpolare. Corespunzător, eroarea de la integrare, exprimată prin funcţia de interpolare, va fi integrala erorii de interpolare. Pe un segment elementar eroarea de interpolare (conform formulei pentru estimarea erorii de interpolare Lagrange) va fi:

34

Page 35: Capitolul 2[1]Metode numerice

Capitolul 2. Metode numerice

| f(z) - g(z, xi , xi+1 ) | M * |(z –xi+1 )(z - xi )| / 2!

undexi , xi+1 – punctele de interpolare,M – supremul derivatei 2 pentru funcţia f(z) pe segmentul elementar [xi , xi+1].

Atunci eroarea de integrare pe segmentul elementar e de

Eroarea de integrare pe tot segmentul poate fi interpretată drept suma erorilor de integrare pe segmentele elementare, prin urmare poate fi estimată prin formula:

unde M este supremul derivatei de ordin 2 al funcţiei integrate pe segmentul [a, b];h – lungimea unui segment elementarFormula dată ne permite să stabilim apriori care este numărul necesar de iteraţii pentru

obţinerea unei erori mai mici decât cea prestabilită:

.

Pentru metoda trapezelor această formulă este mai puţin utilă din cauza convergenţei rapide a metodei – de obicei numărul de divizări necesare este destul de mic. În anumite condiţii, totuşi ea poate optimiza procesul de calcul şi poate fi utilizată direct pentru stabilirea numărului de iteraţii necesare.

Algoritmizarea problemei

Varianta A (pentru numărul cunoscut de divizări).1. Se introduc limitele de integrare a, b.2. Se indică numărul de divizări n.

35

Page 36: Capitolul 2[1]Metode numerice

Capitolul 2. Metode numerice

3. Se calculează pasul de deplasare h (pentru cazul când e cunoscut numărul de divizări).

4. Pornind de la a calculăm valorile funcţiei în extremităţile segmentelor elementare şi ariile trapezelor.

5. Însumăm ariile elementare.6. Afişăm rezultatul.

Varianta B (pentru exactitatea prestabilită).1. Se introduc limitele de integrare a, b şi sup(f ’’(x)) pe [a, b].2. Se indică exactitatea cerută .3. Se fixează sistemul de noduri pentru divizarea segmentului [a, b] (iniţial poate fi

considerat din 2 puncte – extremităţile segmentului [a, b] ).4. Pornind de la a calculăm valorile funcţiei în extremităţile segmentelor elementare şi

ariile trapezelor.5. Însumăm ariile elementare.6. Estimăm eroarea conform formulei. Dacă eroarea este mai mare decât exactitatea

cerută, mai adăugăm un nod la divizarea curentă, apoi repetăm paşii 4-5, în caz contrar afişăm rezultatul. Sfîrşit.

Funcţii frecvent utilizate

function itera(a,b,M,e:real):longint;{funcţia determină numărul minim de divizări necesar pentru obţinerea exactităţii cerute - metoda dreptunghiurilor, toate variaţiile}

var n1,n2,n : longint; e1:real;begin{fixăm numărul minim posibil şi maxim acceptabil de divizări} n1:=0;n2:=maxlongint div 2; while(n2-n1)>1 do {divizăm intervalul atât timp, cât lungimea lui e mai mare ca 2} begin n:=trunc((n1+n2) / 2); e1:=(b-a)*M*(b-a)/4/n; if e1<e then n2:=n else n1:=n; {verificăm necesitatea micşorării sau măririi numărului de divizări} end; itera:=n;end;

function expo(x:real):real;{calculul functiei e la puterea x}var t,s:real; n:integer; begin s:=0;t:=1;n:=1; while t>e do begin s:=s+t; t:=t*x/n; n:=n+1; end;

36

Page 37: Capitolul 2[1]Metode numerice

Capitolul 2. Metode numerice

expo:=s; end;

function mysin(x:real):real;{calculul functiei sin x}var t,s:real; semn,n2 : integer; begin s:=0;t:=x;n2:=2; semn:=-1; while abs(t)>e do begin s:=s+t; t:=semn*t*x*x/(n2+1)/n2; n2:=n2+2; end; mysin:=s; end;

function coso(x:real):real;{calculul functiei cos(x)}var t,s:real; semn,n:integer; begin s:=0;t:=1;n:=2; semn:=-1; while abs(t)>e do begin s:=s+t; t:=semn*t*x*x/(n-1)/n; n:=n+2; end; coso:=s; end;

ExempluSe consideră figura mărginită de graficul funcţiei y=e|x/2| , axa 0x şi de dreptele x=-2 şi

x=2. Elaboraţi un program, care va determina poziţia dreptei verticale, care împarte figura dată în două părţi, ariile cărora, respectiv din stânga şi din dreapta, se raportă ca 1:2. Calculul se va efectua utilizând formula trapezelor cu exactitatea e=0,0001.

Rezolvare. Vom determina mai întâi M pentru funcţia dată pe intervalul propus:

y=e|x/2| => y’’=1/4 e|x/2|.

Deoarece funcţia este pară e suficient să determinăm M pe intervalul (0, 2). Pe acest interval derivata 2 e crescătoare şi valoarea maximă se obţine în x=2. De aici

M=e/4 0,7.

Pentru o exactitate mai mare vom considera M=1.

Program P6;var s1,x1,x2,e,S,a,b,h,M:real; {x1,x2 - extremităţile segmentului elementar; i,n: integer; {i – contor, n – numărul de divizări} S - aria totală; S1 aria figurii separate în stânga de linia verticală în raportul cerut ab – extremităţile intervalului; h – pasul;

37

Page 38: Capitolul 2[1]Metode numerice

Capitolul 2. Metode numerice

M – marginea superioară a derivatei 2} function f(x:real):real; begin f:=exp(abs(x/2)); end;

begin {deoarece nu se impun restricţii la modul de introducere a datelor iniţiale, le fixăm valorile direct în program} b:=2; a:=-2; M:=0.8; e:=0.001; n:=2; {numărul iniţial de divizări} h:=(b-a)/n; while (b-a)/12*h >e do {atât timp, cât nu e atinsă exactitatea cerută repetăm procesul} begin s:=0; for i:=0 to n-1 do begin x1:=a+i*h; x2:=a+(i+1)*h; s:=s+h*(f(x1)+f(x2))/2; end; n:=n+1; {Mărim numărul de divizări. Creşterea numărului de divizări poate fi făcută şi mai mare, ex n:=n+10} h:=(b-a)/n; {Recalculăm lungimea segmentului elementar} end; s1:=0; i:=0; n:=0; while s1<s/3 do {Pornind de la extremitatea stângă a figurii, adăugăm trapezele elementare, până nu ajungem la proporţia cerută} begin inc(n); x1:=a+i*h; x2:=a+(i+1)*h; s1:=s1+h*(f(x1)+f(x2))/2; i:=i+1; end;{afişarea rezultatelor} writeln('s=', s:10:6); writeln('x= ', (x1+x2)/2:10:6); {în calitate de punct de trecere luăm mijlocul laturii de jos a ultimului trapez}end.

Întrebări şi exerciţii

1. Daţi formularea matematică a metodei trapezilor.2. Elaboraţi o procedură care realizează metoda trapezilor.3. Care este eroarea de calcul a metodei trapezilor?4. Se consideră figura mărginită de graficul funcţiei y = e|x/2|, axa 0x şi de dreptele x = -

2 şi x = 2. Elaboraţi un program care va determina poziţia dreptei verticale care împarte figura examinată în două părţi, ariile cărora, respectiv din stînga şi din dreapta, se raportă ca 1:2. Calculele se vor efectua utilizând formula trapezelor cu exactitatea ε = 0,001. La ecran se

38

Page 39: Capitolul 2[1]Metode numerice

Capitolul 2. Metode numerice

va afişa valoarea S a ariei întregii figuri şi coordonata punctului b de pe axa 0x prin care trece dreapta verticală căutată.

5. Compuneţi un program ce calculează integrala

utilizînd formula trapezelor cu exactitatea ε =0,001. Programul va afişa la ecran valoarea integralei I şi numărul de diviziuni n pentru care se realizează exactitatea dată.

6. Compuneţi un program ce calculează integrala

aplicând formula trapezelor pentru diferite numere de diviziuni ale intervalului de integrare n = 10, 100, 1000. Programul va afişa la ecran valoarea integralei I pentru fiecare număr de diviziuni n.

7. Elaboraţi un program ce calculează aria figurii mărginite de liniile ce nu se intersectează definite de funcţiile y = x3+3, y = etg 2x şi de dreptele x = -0,3 şi x = 0,1. Calculele se vor efectua utilizând formula trapezelor cu exactitatea ε =0,001. La ecran se va afişa valoarea S a ariei figurii şi numărul de diviziuni n a intervalului de integrare.

8. Elaboraţi un program ce calculează integrala

cu exactitatea ε = 0,001 prin două metode: utilizînd formula trapezelor şi formula dreptunghiurilor (metoda mediilor). La ecran se va afişa valoarea integralei I şi numărul de diviziuni n1 şi n2 pentru care a fost calculată integrala prin fiecare metodă.

9. Elaboraţi un program ce calculează integrala

cu exactitatea ε = 0,01 prin două metode: aplicând formula trapezelor şi formula dreptunghiurilor de stânga. La ecran se va afişa valoarea integralei I şi numărul de diviziuni n1 şi n2 pentru care a fost calculată integrala prin fiecare metodă.

10. Elaboraţi un program ce calculează integrala

39

Page 40: Capitolul 2[1]Metode numerice

Capitolul 2. Metode numerice

cu exactitatea ε = 0,0001 prin două metode: aplicând formula trapezelor şi formula dreptunghiurilor (metoda mediilor). La ecran se va afişa valoarea integralei I şi numărul de diviziuni n1 şi n2 pentru care a fost calculată integrala prin fiecare metodă.

40