Curs 9: METODE NUMERICE UTILIZATE N SIMULAREA SISTEMELOR DINAMICE

1 Noiunea de sistem dinamic. Clasificri. Noiunea de simulare

Un sistem dinamic este o entitate care se caracterizeaz printr-un mod specific de transformare a unui semnal de intrare (cauza) ntr-un semnal de ieire (efectul). Matematic, aceste semnale sunt funcii de timp (figura 1).

Fig. 1 Reprezentarea de tip cutie neagr a unui sistem dinamic

Dac mulimea momentelor de timp n care evolueaz un sistem este o submulime a: mulimii numerelor reale, , sistemul este cu timp continuu (sistem continuu); mulimii numerelor ntregi, , atunci sistemul se numete cu timp discret (sistem

discret sau sistem cu eantionare). n modelarea sistemelor dinamice se utilizeaz legi specifice de conservare de mas i

energie, descriind fenomene de acumulare i de dezacumulare, exprimate matematic sub forma ecuaiilor de bilan, care sunt ecuaii difereniale. Aadar, se poate spune c prima imagine asupra comportamentului temporal al unui sistem dinamic este dat de modelul de tip ecuaie diferenial. Ecuaiile difereniale sunt fie ordinare la care variabila de derivare este timpul, t fie cu derivate pariale la care mai apare nc cel puin o variabil de derivare pe lng timp (de exemplu, o coordonat spaial). Fiecare tip de ecuaie diferenial este asociat respectiv unei clase de sisteme:

sistemele cu parametri concentrai sunt descrise de ecuaii difereniale ordinare (exemplul 1);

sistemele cu parametri distribuii sunt descrise de ecuaii difereniale cu derivate pariale (exemplul 2.).

Exemplul 1: Se consider sistemul mecanic din figura 2 a), unde un corp de mas M, legat de un perete printr-un resort de constant de elasticitate k, se deplaseaz cu frecare vscoas, caracterizat de coeficientul de frecare fv, sub aciunea unei fore variabile n timp, f(t).

Modelarea ca sistem dinamic ncepe cu identificarea drept cauz (intrare) a forei f(t) i drept efect (ieire) a deplasrii x(t) (figura 2 b)). Cele dou mrimi sunt legate prin legea de conservare a energiei mecanice, descris de ecuaia diferenial ordinar de ordinul al II-lea:

(1) 2

2d ( ) d ( ) ( ) ( )

dd vx t x tM f k x t f t

tt + + = ,

care este un model al sistemului.

SISTEM DINAMIC (MODEL)

u(t) y(t)

2

a)

b)

Fig. 2 Sistem mecanic cu resort: a) reprezentare fizic; b) simbolizare ca sistem dinamic Sistemul mecanic descris este, deci, un sistem dinamic cu parametri concentrai.

Exemplul 2: Se consider o bar de lungime L, a crei temperatur se modific drept urmare a aplicrii la un moment dat, t0, la capetele ei a unor temperaturi cu variaie cunoscut n timp, 0(t) i L(t) (figura 3). Se cunosc constantele de material: densitatea materialului barei, [kg/m3], cldura lui specific, c [kcalgrd/kg] i conductibilitatea termic, [kcalgrd/m/h], precum i distribuia de temperaturi de-a lungul barei la momentul t0, 0(x).

Fig. 3 Sistem guvernat de transferul unidirecional de cldur

Temperatura n punctele barei este funcie i de timp, dar i de deprtarea de capetele barei, deci de coordonata spaial, notat cu x. Lund ca origine a coordonatei spaiale captul din stnga al barei, se noteaz cu (x,t) temperatura n punctul aflat la abscisa x la momentul de timp t (figura 4.3). (x,t) este efectul datorat cauzelor 0(t) i L(t) (numite i foraje de temperatur). Evoluia temperaturii (x,t) este descris de ecuaia transferului unidirecional de cldur, care este o ecuaie diferenial cu derivate pariale:

(2) 2

2( , ) ( , )x t x t

at x

=

,

unde ( )a c= se numete difuzabilitate termic [m2/s]. Sistemul descris este, deci, un sistem dinamic cu parametri distribuii.

Dac un sistem rspunde printr-o ieire mrginit unui semnal de intrare mrginit, atunci se spune c sistemul este stabil (intrare-ieire). O clas important de sisteme sunt cele liniare, crora li se asociaz ecuaii difereniale liniare (a se vedea exemplul 1). Sistemele liniare stabile sunt caracterizate prin aceea c las neschimbat forma semnalelor armonice; cu alte cuvinte, dac la intrarea unui sistem liniar stabil se aplic semnal sinusoidal, atunci ieirea lui va fi tot sinusoidal. Majoritatea sistemelor reale nu verific aceast proprietate; se spune c sunt neliniare.

Termenul de simulare este asociat evoluiei temporale a unui sistem dinamic descris printr-un model adecvat, atunci cnd sistemului i se aplic un semnal de intrare cunoscut; n particular, se pune problema determinrii ieirii sistemului, numite i rspuns (la o anume intrare cunoscut). De exemplu, simularea sistemului din exemplul (1) const n a afla cum

0 x L x

0(t) (x,t) L(t)

( )x t

( )f t

vf

k

M

SISTEM DINAMIC

(MODEL) f(t) x(t)

3

variaz deplasarea corpului n timp, x(t), atunci cnd se cunoate f(t), forma de variaie a forei aplicate la un moment dat, t0. Pentru aceasta este necesar integrarea numeric a

ecuaiei difereniale (1) n condiii iniiale cunoscute, ( )0x t i ( ).

0x t . Simularea sistemului din exemplul (2) presupune integrarea numeric a ecuaiei (2) n condiii iniiale cunoscute, (x,t0)=0(x), pentru x [0;L].

A simula un sistem dinamic este echivalent cu a rezolva (integra) ecuaia diferenial care l descrie, n condiii iniiale cunoscute, care reprezint starea iniial a sistemului. Acest capitol este dedicat metodelor de integrare numeric a ecuaiilor difereniale, ca prim metod de simulare numeric a sistemelor dinamice.

O mare parte din problematica determinrii rspunsului sistemelor se reduce la gsirea soluiei ecuaiei difereniale (sistemului de ecuaii difereniale) ca soluie (unic) a problemei cu timp iniial (problemei Cauchy).

Deoarece ecuaiile difereniale care modeleaz procesele din lumea real sunt adesea neliniare i de ordin superior, gsirea unei soluii analitice (exacte) a problemei Cauchy este dificil. Alternativa o constituie gsirea unei soluii aproximative, i anume prin folosirea unui algoritm executabil pe o main de calcul. Astfel, soluia gsit prin aplicarea unei metode numerice este un ir de aproximaii ale valorilor soluiei exacte, calculate la momente discrete de timp, de obicei egal distanate.

2 Simularea sistemelor cu parametri concentrai. Integrarea numeric a ecuaiilor difereniale ordinare i a sistemelor de ecuaii difereniale ordinare

2.1 Fundament teoretic n aceast seciune sunt rezumate principalele rezultate teoretice care stau la baza

metodelor numerice de rezolvare a ecuaiilor difereniale ordinare (EDO) i a sistemelor de ecuaii difereniale ordinare.

n vederea rezolvrii numerice a EDO sau a sistemelor de EDO, indiferent de metoda folosit, mai nti se procedeaz la scrierea acestora sub forma standard:

(3) ( ) ( )0 0 0d , , [ ; ], d fy f t y t t t y t yt

= = ,

unde [ ]1 ... Tny y y= i condiia iniial, 0 00 1 ...T

ny y y = , sunt vectori de

dimensiune n, unde n este ordinul ecuaiei sau al sistemului considerat, [ ]1 ... Tnf f f= este o funcie vectorial de aceeai dimensiune, iar 0[ ; ]ft t este intervalul de integrare.

Scrierea (3) este posibil oricare ar fi ordinul ecuaiei sau al sistemului considerat, n virtutea unui rezultat cunoscut din matematic: orice ecuaie diferenial de ordinul n>1, printr-o schimbare de variabil, poate fi scris ca un sistem de n ecuaii difereniale de ordinul I. Aceast echivalare se face introducnd n noi variabile, al cror mod de alegere nu este unic, dup cum se arat mai jos.

Pentru o ecuaie diferenial de ordin n>1, dat sub forma:

4

(4) ( ) ( ) ( )0 0 0

. .. ( 1)0

. ( 1)0 1 0 2 0

d, , , ,... , ; ,

d

, ,...,

nn fn

nn

y f t y y y y t t tt

y t y y t y y t y

=

= = =

,

se poate opta pentru alegerea noilor variabile astfel:

(5) ( 1) ( ) ( ), 1,i iy t z t i n = = , obinndu-se un sistem de ecuaii difereniale n variabila z:

(6)

( )

.

1 2.

2 3

.

1 2

...

, , ,...,n n

z z

z z

z f t z z z

=

=

=

,

cu condiia iniial ( ) 0 0 0 00 1 1 ...T

n nz t z y z y = = = . Ca urmare a acestei schimbri de variabil, sistemul (6) este de tipul (1):

(7) ( ) ( )0 0d , , d zz f t z z t zt

= = ,

unde:

(8) ( )2 3 1 2... , , ,..., Tz nf z z f t z z z = Observa ie :

Schimbarea de variabil (5) reprezint numai una dintre posibilitile de a transforma o ecuaie diferenial de ordin superior ntr-un sistem de ecuaii difereniale de ordinul I, obinndu-se astfel forma standard (1), care face obiectul metodelor numerice de simulare a rspunsului sistemelor. Se pot utiliza i alte schimbri de variabil care au acelai efect, evident cu expresii diferite pentru funcia vectorial zf . Teorema 1 (unicitatea soluiei problemei (3)):

Fie mulimea convex ( ){ }0, , nfD t y t t t y= i funcia vectorial ( ),f t y continu pe D. Dac exist o constant L>0 astfel nct ( ) ( )1 2 1 2, ,f t y f t y L y y adic f este lipschitzian (notaia desemneaz norma vectorial) atunci soluia problemei (3), ( )y t , este unic pentru 0; ft t t .

Se pune acum problema consecinelor pe care le au asupra soluiei ecuaiei (4.3) eventualele perturbaii datorate erorilor de rotunjire aprute n condiiile iniiale sau/i n funcia f. Se consider problema perturbat:

(9) ( ) ( ) ( ) ( )0 0 0 0d , , ; , d fz f t z t t t t z t y tt

= + = + ,

care admite soluia (exact) z(t). Se demonstreaz urmtorul rezultat. Teorema 2 (soluia problemei perturbate):

5

Dac f este continu i lipschitzian pe 0; ft t , atunci pentru orice 0 > , cu { }0max , ( )t > , exist ( ) 0k > pentru care este ndeplinit relaia:

(10) ( ) ( ) ( ) 0, ; fz t y t k t t t < , unde y(t) este soluia (exact a) problemei neperturbate. Se spune c problema (4.3) este bine formulat.

Metodele numerice de integrare a ecuaiilor difereniale furnizeaz soluii ale problemelor perturbate, aceasta datorit erorilor de rotunjire.

2.2 Elemente necesare pentru integrarea numeric a EDO Aceste elemente trebuie vzute drept date de intrare care sunt furnizate unei metode

numerice de integrare. Ele sunt listate mai jos: 1) EDO de rezolvat trebuie mai nti adus la forma (3) printr-o schimbare de variabil

adecvat (n particular se poate utiliza schimbarea de variabil (5)); rezult astfel funcia vectorial zf , de dimensiune egal cu ordinul ecuaiei, n, care descrie legtura dintre

derivatele celor n noi variabile, .

iz , i valorile lor nederivate, iz ; 2) intervalul de integrare dat prin capetele acestuia, 0t i ft ;

3) condiiile iniiale, sub forma unui vector de dimensiune n, 0 00 1 ...

Tny y y = ;

4) pasul de integrare, notat cu h n toate metodele de integrare se consider o diviziune de ordinul m a intervalului 0; ft t (o mprire a acestui interval n intervale de aceeai lungime, h), n ale crei puncte, { }0 , 0,...,kt t k h k n= + , se calculeaz valorile soluiei numerice a ecuaiei, notate cu ky , care sunt aproximri ale valorilor exacte, ( )ky t ; este intuitiv evident c soluia numeric se apropie cu att mai mult de soluia exact cu ct diviziunea este mai fin (pasul de integrare, h, este mai mic);

5) un subprogram (funcie) care calculeaz valoarea funciei zf pentru orice moment de timp, t, i care va fi apelat la fiecare pas de integrare, kt , de programul unde este implementat metoda de integrare numeric;

6) alte date care se constituie ca parametri constani n funcia zf .

2.3 Metode unipas (directe) Metodele unipas de integrare numeric a EDO sunt metodele Taylor, care nu se

folosesc, fiind complicate, i metodele Runge-Kutta (RK). Aceste metode au ca baz de pornire expresia dezvoltrii polinomiale n serie Taylor a funciei y(t), presupus derivabil de ordinul (p+1) i cu derivatele continue:

(11) ( ) ( ) ( ) ( ) ( ) ( ) ( )2 1

. .. ( ) ( 1)1 2! ! 1 !

p pp p

k k k k k kh h hy t y t h y t y t y t y

p p

++

+ = + + + + + + ,

unde ( )1,k k kt t + . nlocuind n (4.11) expresia derivatei de ordinul i a funciei ( )y t : (12) ( ) ( )( )( ) ( 1) , , 1, ( 1)i iy t f t y t i p= = +

6

n punctul kt , se obine:

(13) ( ) ( ) ( )( ) ( )( ) ( ) ( )( )2 1

. ( )1 , , ... ,2! 1 !

pp

k k k k k k k kh hy t y t hf t y t f t y t f y

p

+

+ = + + + + +

Relaia (13) este o ecuaie n diferene ce permite calculul soluiei n punctul 1kt + pe baza soluiei n punctul kt . Transpunerea numeric a acestei relaii conduce la calculul soluiei la un anume pas de integrare, 1ky + , n funcie de soluia la pasul anterior, ky . Aceast dependen, care ia n considerare istoria integrrii cu un singur pas n urm justific numele acestei clase de metode (unipas).

n general, scopul final al oricrei metode numerice este s determine cu suficient acuratee o soluie aproximativ a problemei cu minimum de resurse de memorie i numr de operaii. Cea mai important msur a acurateii o constituie eroarea local de trunchiere. Aceasta reprezint la un anumit pas msura n care soluia exact a ecuaiei difereniale satisface ecuaia n diferene utilizat pentru calcul.

Se demonstreaz c la metodele Taylor, care sunt bazate pe relaia (13), eroarea local de trunchiere este mrginit superior de un polinom de grad p n pasul de integrare, ceea ce nseamn c este de ordinul ph ( ( )pO h ). Aceasta arat c, ntr-adevr, acurateea metodei crete odat cu micorarea lui h. Pe de alt parte, micorarea lui h antreneaz creterea numrului de operaii efectuate i, deci, a erorilor de rotunjire. Se demonstreaz c exist o valoare a lui h pentru care se atinge minimul erorii de trunchiere n prezena erorilor de rotunjire. Concluzia este c pasul de integrare, h, rezult dintr-un compromis ntre acurateea i precizia metodei.

Pentru p=0 se obine metoda Taylor de ordinul 1, cunoscut i ca metoda Euler sau metoda Runge-Kutta de ordinul 1:

(14) ( )1 ,k k k ky y h f t y+ = + , la care minimizarea erorii de trunchiere se obine pentru:

2h M= , unde este eroarea de rotunjire la fiecare pas de integrare, iar M este o constant care mrginete superior valorile derivatei a doua ale soluiei exacte a ecuaiei pe intervalul de

integrare: ( ).. 0, ; fy t M t t t . Metodele Taylor de ordinul p>1 nu sunt folosite ca atare, deoarece necesit calculul

expresiei analitice a derivatelor funciei din membrul drept din (4.3). Metodele RK elimin acest dezavantaj, pstrnd avantajul erorii locale de trunchiere mici.

Metoda Runge-Kutta de ordinul 2 sau metoda punctului median care are o eroare de trunchiere de ( )2O h , are la baz expresia:

(15) ( )1 , ,2 2k k k k k kh hy y h f t y f t y+ = + + +

Alte dou metode cunoscute sub numele general de metode Runge-Kutta de ordinul 2 sunt metoda Euler modificat, care are relaia n diferene:

7

(16) ( ) ( )( )1 1, , ,2k k k k k k k khy y f t y f t y h f t y+ + = + + +

i metoda Heun, cu relaia n diferene:

(17) ( ) ( )1 1 2 2, 3 , ,4 3 3k k k k k k k khy y f t y f t h y h f t y+ +

= + + + +

Metoda RK de ordinul 3 este puin utilizat. Una dintre cele mai utilizate metode unipas este metoda Runge-Kutta de ordinul 4,

creia i corespund urmtoarele formule de calcul:

(18) ( )1 1 2 3 41 2 26k ky y k k k k+ = + + + + , unde:

(19) ( )

( )1 2 1

3 2 4 3

1, ,

2 21

, ,

2 2

k k k k

k k k k

hk h f t y k h f t y k

hk h f t y k k h f t h y k

= = + +

= + + = + +

Aceast metod are eroarea local de trunchiere ( )4O h , iar funcia y(t) trebuie s fie de clas 5C (continu i cu derivatele pn la ordinul 5 continue).

Exemplu: Sa se scrie un program Matlab care sa realizeze integrarea numerica a sistemului de ecuatii diferentiale:

( )1 2

2 2 11 25 5 252

x x

x x x

=

=

Cerinte: - Se va utiliza metoda de integrare numerica Runge-Kutta de ordinal 4. - Se vor afisa grafic, in aceeasi figura, rezultatele obtinute prin integrare

Rezolvare: % program principal pentru integrare numerica clear all;close all; x=[7 10]; h=0.1; t0=0; tf=10; t(1)=0; i=1; for m=0.1:h:tf, k1=func(t,x(i,:)); k2=func(t+h/2,x(i,:)+h*k1/2); k3=func(t+h/2,x(i,:)+h*k2/2); k4=func(t+h,x(i,:)+h*k3); x(i+1,:)=x(i,:)+h/6*(k1+2*k2+2*k3+k4); t(i+1)=m;

8

i=i+1; end subplot(211); hold on; plot(t,x(:,1),'k'); subplot(212); hold on; plot(t,x(:,2),'k');

% implementare sistem de ecuatii diferentiale function xd=func(t,x) xd(1)=x(2); xd(2)=1/2*(25-5*x(2)-25*x(1));

2.4 Metode directe cu pas variabil n general nu se poate determina eroarea global a unei metode directe de integrare, dar

se tie c exist o strns legtur ntre aceasta i eroarea local de trunchiere. Prin utilizarea adecvat a metodelor directe de diferite ordine se poate predicta eroarea local de trunchiere i, mpreun cu variaia pasului, se poate controla eroarea global.

Metodele directe cu pas variabil se mai numesc i cu pas adaptiv pentru c adapteaz numrul i poziia punctelor n care se face calculul astfel nct eroarea local de trunchiere s se menin inferioar unei limite impuse. Acest paragraf este dedicat manierei n care eroarea local de trunchiere poate fi controlat prin variaia pasului de integrare.

Eroarea local de trunchiere produs prin aplicarea unei metode directe de ordin p poate fi controlat prin folosirea unui pas de integrare modificat. Presupunnd c aceast modificare este de forma nou vechih h= , se demonstreaz c este suficient s se aleag:

(20) 1

1 12p

k ky y+ +

,

unde cu s-a notat limita impus a erorii locale de trunchiere (tolerana impus), iar 1ky + reprezint soluia la pasul k+1 dat de o alt metod de integrare dect prima folosit.

O metod cunoscut care utilizeaz aceast tehnic, de control al erorii globale prin variaia pasului de integrare, este metoda Runge-Kutta-Fehlberg. Aceast metod const n utilizarea unei metode Runge-Kutta cu eroarea local de trunchiere de ordinul cinci:

(21) 1 1 3 4 5 616 6656 28561 9 2135 12825 56430 50 55k k

y y k k k k k+ = + + + + ,

pentru a estima eroarea local de trunchiere a unei metode RK de ordinul 4 dat prin relaia:

(22) 1 1 3 4 525 1408 2197 1

216 2565 4104 5k ky y k k k k+ = + + + ,

unde:

9

(23)

( )1 2 13 1 2

4 1 2 3

5 1 2 3 4

6 1

1, ,

4 43 3 9

,

8 32 3212 1932 7200 7296

,

13 2197 2197 2197439 3680 845

, 8216 513 4104

8, 2

2 27

k k k k

k k

k k

k k

k k

hk h f t y k h f t y khk h f t y k k

hk h f t y k k k

k h f t h y k k k khk h f t y k k

= = + +

= + + +

= + + +

= + + +

= + + 2 3 4 53544 1859 112565 4104 40

k k k

+

Un avantaj al acestei metode este c sunt necesare numai ase evaluri ale funciei f. Metodele RK de ordinul 4 i 5, utilizate mpreun necesit cel puin patru evaluri pentru cea de ordinul 4 i ase pentru cea de ordinul 5, deci n total cel puin zece evaluri. O valoare iniial pentru h la pasul k este utilizat pentru a gsi primele valori ale lui 1ky + i 1ky + care, la rndul lor, sunt utilizate n calculul lui . n mod practic, pentru alegerea unei valori iniiale pentru h se folosete o regul mai mult sau mai puin empiric. Pentru metoda Runge-Kutta-Fehlberg de ordinul 4 aceast valoare se poate obine particulariznd p=4 n relaia (4.20):

1 14 4

1 1 1 10.84

2 k k k ky y y y+ + + +

=

Algoritmul Runge-Kutta-Fehlberg (de calcul al soluiei ecuaiei difereniale (3) cu eroare local de trunchiere avnd o toleran dat, ) Date de intrare: t0, tf, y0, hmin, hmax, Pas 1. 0t t= , 0y y= , maxh h= , 1flag = Pas 2. Ct timp 1flag = repet paii 3-11. Pas 3. ( )1 ,k kk h f t y ( )2 1/ 4, / 4k kk h f t h y k + + ( )3 1 23 /8, 3 / 32 9 / 32k kk h f t h y k k + + + ( )4 1 2 312 /13, 1932 / 2197 7200 / 2197 7296 / 2197k kk h f t h y k k k + + + ( )5 1 2 3 4, 439 / 216 8 3680 / 513 845 / 4104k kk h f t h y k k k k + + + ( )6 1 2 3 4 5/ 2, 8 / 27 2 3544 / 2565 1859 / 4104 11 / 40k kk h f t h y k k k k k + + + Pas 4. 1 3 4 5 6/ 30 128 / 4275 2197 / 75240 / 50 11 / 40R k k k k k h + Pas 5. Dac R , atunci t t h= + 1 3 4 525 / 216 1408 / 2565 2197 / 4104 / 5y y k k k k= + + + Pas 7. Scrie t, y.

Pas 8. 140.84

R

Pas 9. Dac 1.0 , atunci 0.1h h= altfel, dac 4 , atunci 4h h= altfel h h=

10

Pas 10. Dac maxh h> , atunci maxh h= Pas 11. Dac ft t , atunci 0flag = altfel, dac ft h t+ > , atunci fh t t= altfel, dac minh h> , atunci 0flag = altfel scrie Pas minim depit!

2.5 Metode multipas (indirecte) 2.5.1 Introducere

Metodele unipas, discutate n capitolele anterioare, folosesc numai informaia de la momentul curent, kt , pentru a calcula soluia EDO la momentul urmtor, 1kt + . Dei aceste metode folosesc i valori calculate n interiorul intervalului de integrare [ ]1,k kt t + , totui acestea nu sunt reinute pentru utilizare ulterioar i, prin urmare, nu sunt disponibile utilizatorului. Deoarece informaia calculat la momentele 0 1, ,..., kt t t este disponibil i eroarea local de trunchiere, ( )k ky y t , crete odat cu k, atunci apare motivat dezvoltarea de metode care s foloseac informaia de la momentele de integrare anterioare. Aceste metode care, pentru a calcula soluia 1ky + , folosesc istoria integrrii cu cel puin doi pai n urm, se numesc metode multipas sau indirecte. Def in i ia 4.1: O metod multipas de ordinul m folosit pentru a determina soluia

aproximativ a problemei Cauchy (3) presupune calculul lui 1ky + cu relaia:

(24) ( ) ( ) ( )( )1 1 2 1 11 1 1 1 1...

, , ,

k m k m k o k mm k k m k k o k m k m

y a y a y a yh b f t y b f t y b f t y

+ +

+ + + +

= + + ++ + + + ,

unde , , 0,..., 1i ia b i m= sunt constante. Cnd 0mb = , metoda se numete explicit sau deschis, n caz contrar ea se numete

implicit sau nchis. Condiia iniial, 0y , este cunoscut. Folosirea relaiei (24) ntr-o structur de ciclare

necesit ca urmtoarele m-1 valori care urmeaz condiiei iniiale, 1 2 1, ,..., my y y , s fie calculate separat, nainte de a se intra n ciclare. De obicei, aceste valori se calculeaz printr-o metod direct, Runge-Kuta. Aceast operaiune poart denumirea de amorsare. Dup amorsare, valorile 1 1, ,..., k k k my y y + se calculeaz n ciclu utiliznd formula (24) particularizat pentru metoda multipas respectiv (aceast faz se numete integrarea propriu-zis).

Metodele implicite sunt folosite n practic mai ales pentru a mbunti valoarea calculat la un pas de integrare printr-o metod explicit. Combinaia dintre metodele explicite i cele implicite d natere la aa numitele metode predictor-corector sau cu predicie i corecie. La aceste metode faza de integrare propriu-zis se mparte la fiecare pas n dou etape: calculul (predicia) valorii lui 1ky + prin formula explicit i, respectiv, corecia acestei valori cu formula implicit, pn la o anume precizie impus. Metodele cu predicie i corecie sunt singurele metode care permit controlul direct al preciziei de calcul.

11

n cele ce urmeaz, clasa metodelor predictor-corector este tratat separat de restul metodelor multipas, pe care le vom numi simple.

2.5.2 Metode multipas simple Paii majori ai integrrii numerice prin metode multipas simple sunt dai mai jos (notaia

[] desemneaz partea ntreag). Algoritmul general de integrare numeric a EDO prin metode multipas simple Date de intrare: t0, tf, y0, h Pas 1. 0t t= , 0y y= Pas 2. Amorsarea: calculul valorilor 1 2 1, ,..., my y y (printr-o metod direct, Runge-Kutta).

Pas 3. Pentru k de la m-1 la 0 1ft t

h

( )1 1 1AB , ,..., k m k k k my y y y+ + = repet

Expresiile de calcul pentru formulele de integrare multipas se deduc pornind de la

integrala funciei ( ).y t pe intervalul [ ]1,k kt t + :

(25) ( ) ( ) ( ) ( )( )1 1.1 d , dk kk k

t t

k kt t

y t y t y t t f t y t t+ ++ = =

Deoarece nu se poate integra ( )( ),f t y t fr cunoaterea lui ( )y t , se va utiliza un polinom de interpolare, P(t), care trece prin punctele ( ) ( )1 1, , , ,...,k k k kt y t y ( )1 1,k m k mt y+ + . Dac se presupune i ( )k ky t y , atunci relaia (25) devine:

(26) ( ) ( )11 dkk

t

k kt

y t y P t t+

+ +

n principiu, orice metod de interpolare este utilizabil pentru deteminarea lui P(t), dar cea mai des folosit este interpolarea Newton de tip diferen napoi.

Eroarea local de trunchiere corespunztoare metodelor multipas se calculeaz similar metodelor unipas.

Fr a detalia calculele, n tabelul 1 se prezint formulele de integrare i erorile de trunchiere pentru metoda multipas explicit de ordinul m, numit i Adams-Bashforth (AB), n cazurile m=2, 3, 4.

m Formula de integrare Eroarea local de trunchiere la pasul k+1

2 (27) ( ) ( )( )1 1 13 , ,2k k k k k khy y f t y f t y+ = + ( ) ( ) ( )

2(3)

1 2 15

, ,

12k k k k khh y t t+ =

3 (28) ( ) ( )( ( ))1 1 1

2 2

23 , 16 ,12

5 ,k k k k k k

k k

hy y f t y f t yf t y

+

= +

+ ( ) ( ) ( )

3(4)

1 2 13

, ,

8k k k k khh y t t+ =

12

4 (29) ( ) ( )( ( ) ( ))1 1 1

2 2 3 3

55 , 59 ,24

37 , 9 ,k k k k k k

k k k k

hy y f t y f t yf t y f t y

+

= + +

+ ( ) ( ) ( )

4(5)

1 2 1251

, ,

720k k k k khh y t t+ =

Tabel 1 Caracteristicile metodei multipas Adams-Bashforth pentru diferite ordine

Metodele multipas implicite, numite i Adams-Moulton (AM), utilizeaz perechea ( )( )1 1 1, ,k k kt f t y+ + + ca punct adiional al polinomului de interpolare n aproximarea

integralei ( )( )1 , dkk

t

tf t y t t+ .

n tabelul 2 se prezint formulele de integrare i erorile de trunchiere pentru metodele AM de diferite ordine.

m Formula de integrare Eroarea local de trunchiere la pasul k+1

0 (30) ( ) ( )( )1 1 1, ,2k k k k k khy y f t y f t y+ + += + + ( ) ( ) ( )

3(3)

1 1, ,12k k k k khh y t t+ + =

1 (31) ( ) ( )( ( ))1 1 1

1 1

5 , 8 ,12

,

k k k k k k

k k

hy y f t y f t yf t y

+ + +

= + +

( ) ( ) ( )3

(4)1 1 1, ,24k k k k k

hh y t t+ + =

2 (32) ( ) ( )( ( ) ( ))1 1 1

1 1 2 2

9 , 19 ,24

5 , ,k k k k k k

k k k k

hy y f t y f t yf t y f t y

+ + +

= + +

+ ( ) ( ) ( )

4(5)

1 2 119

, ,

720k k k k khh y t t+ + =

Tabel 2 Caracteristicile metodei multipas Adams-Moulton pentru diferite ordine

2.5.3 Metode multipas corective (cu predicie i corecie) Pentru a schia algoritmul de integrare numeric folosind metode corective, se introduce

notaia jky , care semnific soluia la pasul k de predicie i pasul j de corecie. Rezult c valoarea predictat la fiecare pas poate fi considerat drept o corecie la pasul 0 i se noteaz, deci, cu 0ky . ky desemneaz valoarea final (dup predicie i corecie) a soluiei la pasul k.

Presupunem c predicia se realizeaz cu o metod AB de ordinul m, iar corecia cu o metod AM de ordinul p (metoda corectiv rezultat se mai numete i Adams-Bashforth-Moulton). Precizia de corecie se noteaz cu , iar notaia desemneaz norma vectorial. Algoritmul general de integrare numeric a EDO prin metode PC Date de intrare: t0, tf, y0, h, Pas 1. 0t t= , 0y y= Pas 2. Amorsarea: calculul valorilor 1 2 1, ,..., my y y (printr-o metod direct, Runge-Kutta).

Pas 3. Pentru k de la m-1 la 0 1ft t

h

Predicia la pasul de integrare k+1: ( )0 1 1 1AB , ,..., k m k k k my y y y+ + =

Corecia la pasul 1 al pasului de integrare k+1:

13

( )1 01 1 1 1AM , , ,..., k p k k k k py y y y y+ + + = j 1 Ct timp 11 1

j jk ky y

+ + >

Corecia la pasul j al pasului de integrare k+1: ( )1 1 11 1AM , , ,..., j jp k k k pk ky y y y y+ + + +=

j j+1 repet

11 1

jk ky y

++ +

repet Acurateea metodelor multipas corective (PC) depinde de efectul pe care l are asupra

erorilor cuplarea unei anumite formule de integrare explicite n faza de predicie cu o formul implicit n faza de corecie. Analiza numeric a acestui efect a condus la perechile predictor-corector optime sintetizate n tabelul 3.

Formulele de integrare Comentarii

PC1 (33) ( ) ( )( )

( ) ( )( )0

1 1 1

11 1 1

( ) : 3 , ,2

( ) : , ,2

k k k k k k

j jk k k k k k

hP y y f t y f t yhC y y f t y f t y

+

++ + +

= +

= + +

predicia cu AB pentru m=2 i corecia cu AM pentru m=0

PC2 (34) ( ) ( )(( ))

( ) ( )(( ))

01 1 1

2 2

11 1 1

1 1

( ) : 23 , 16 ,12

5 ,

( ) : 5 , 8 ,12

,

k k k k k k

k k

j jk k k k k k

k k

hP y y f t y f t yf t y

hC y y f t y f t yf t y

+

++ + +

= +

+ = + +

predicia cu AB pentru m=3 i corecia cu AM pentru m=1

PC3 (35) ( ) ( )(

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

01 1 1

2 2 3 3

11 1 1

1 1 2 2

( ) : 55 , 59 ,24

37 , 9 ,

( ) : 9 , 19 ,24

5 , ,

k k k k k k

k k k k

j jk k k k k k

k k k k

hP y y f t y f t yf t y f t y

hC y y f t y f t yf t y f t y

+

++ + +

= + +

+

= + +

+

predicia cu AB pentru m=4 i corecia cu AM pentru m=2

Tabel 3 Formulele de integrare pentru cele mai utilizate metode predictor-corector