Post on 06-Mar-2016
description
a.ch. Noiembrie 2008
1
CURS 6 - Rezumat
METODE NUMERICE PENTRU ECUAII DIFERENIALE ORDINARE
6-I METODE NUMERICE PENTRU ECUAII DIFERENIALE DE
ORDINUL NTI
n aceast seciune se vor prezenta metode numerice pentru ecuaii i sisteme de
ecuaii difereniale de ordinul nti problema cu valori iniiale. O ecuaie sau
sistem de ordin mai mare dect unu se pot reduce la un sistem echivalent de
ordinul unu, prin adugarea de funcii necunoscute.
Exemplu: Fie sistemul de ordinul doi,
),,,,(
),,,,(
yxyxtgy
yxyxtfx
Punnd yvxu , , sistemul devine
),,,,(
),,,,(
vuyxtgv
vuyxtfu
vy
ux
Pentru sistemele de ordinul doi, se vor da metode speciale n Seciunea 7-II.
1 Problema cu valori iniiale (consideraii generale)
Fie ecuaia
),( xtfdt
dx (1)
cu condiia iniial
00 )( xtx (1')
Ecuaia (1) cu condiia iniial (1') constituie o problem cu valori iniiale (sau o
problem Cauchy). Dac funcia f ndeplinete urmtoarele condiii pe domeniul
D = I, unde I este definit de att || 0 , iar de bxx || 0 :
a.ch. Noiembrie 2008
2
1) f este definit i continu pe D;
2) f este lipschitzian n raport cu x, adic: exist o constant pozitiv A, astfel c
pentru orice It i orice *, xx avem
|||),(),(| * xxAxtfxtf ,
atunci: notnd cu M marginea superioar a funciei |f | pe D, problema are o soluie
unic x(t) definit pe intervalul || 0tt , unde )/,min( Mba .
n particular, condiia 2 este ndeplinit dac f are derivat parial n raport cu x,
mrginit n D (sau, mai mult, continu pe D).
Pentru un sistem de m ecuaii difereniale cu m funcii necunoscute, fie
T
mxx ][ 1x , T
m tftf )],(),([ 1 xxf , T
mxx ][)0()0(
1
)0( x , i sistemul
),,,( 1 mxxtdt
df
x (2)
cu condiia iniial
)0(
0 )( xx t (2')
Cu domeniul D = I, unde unde I este definit de att || 0 , iar de
mibxx iii ,1,||)0( , condiiile 1 i 2 devin:
1') f este definit i continu pe domeniul D;
2') f este lipschitzian n raport cu argumentele mxx ,,1 , adic: exist
constantele pozitive mjAj ,1, , astfel c pentru orice It i orice *, xx
avem:
mixxAtftfm
j
jjjii ,1,|||),(),(|1
**
xx .
Notm cu iM marginea superioar a funciei || if pe D, i cu imi
MM,1
max
. Dac
1' i 2' sunt ndeplinite, atunci exist o soluie unic x(t) definit pe intervalul
a.ch. Noiembrie 2008
3
|| 0tt , unde )/,,/,min( 1 MbMba m . n particular, condiia 2' este
ndeplinit dac f are derivate pariale n raport cu mjx j ,1, , continue pe I.
Not: Condiia Lipschitz pentru funcia f se poate considera i sub forma:
||||||),(),(|| ** xxxfxf Att ,
iar marginea M este dat de Mt ||),(|| xf , pentru It ),( x . Norma
considerat este norma-
n ceea ce urmeaz vom considera probleme cu valori iniiale (1, 1') au (2, 2'),
pentru care vom presupune ndeplinite condiiile de existen i unicitate ale
soluiei. Considerm calculul soluiei pentru un interval de integrare ],[ 0 TTt ,
inclus n intervalul de existen a soluiei. Metodele numerice vor fi prezentate
pentru o singur ecuaie diferenial (1), i vor fi generalizate la sisteme (2).
2 Operatori de integrare numeric (intr-un singur pas, n mai muli
pai, explicii, implicii)
Gsirea soluiei ecuaiei (1) printr-o metod numeric se va numi integrare
numeric sau integrare pas cu pas. Metoda const n urmtoarele:
a) Intervalul de integrare ],[ 0 TTt se divizeaz prin punctele niti ,0, , unde
TTtn .
b) Ecuaia (1) se cere s fie satisfcut n punctele it , iar ntre aceste puncte,
variaia funciei x(t) se estimeaz.
Vom nota n ceea ce urmeaz:
)( itx = soluia exact;
ix = soluia calculat n it ;
iii extx )( ,
unde ie este eroarea de trunchiere global a metodei, pe pasul i.
a.ch. Noiembrie 2008
4
Un operator de integrare numeric este reprezentat de o formul care d soluia la
momentul 1it n funcie de soluia calculat la k momente anterioare
11 ,,, kiii ttt , i anume:
),,,( 111 kiiii xxxgx (3)
- Dac n membrul doi din (3), g este funcie numai de ix , i eventual 1ix ,
operatorul se zice ntr-un pas, altfel se zice n mai muli pai (i anume, n k
pai). Adic: ),( 11 iii xxgx ; sau )(1 ii xgx .
- Dac n membrul doi din (3) apare i 1ix , operatorul se zice implicit, n caz
contrar se zice explicit.
Integrarea prin operatori implicii conduce la rezolvarea ecuaiei (3) n
necunoscuta 1ix , printr-o metod pentru ecuaii neliniare. O comparaie ntre
operatori ntr-un singur pas i n mai muli pai se va face n 4.8.
Distana dintre dou puncte succesive de diviziune a intervalului de integrare se
zice pas de integrare:
iii tth 11 .
Cazul comun este acela n care pasul este constant: hhi 1 . Avem
htt ii 1 (h = constant).
Exist ns, algoritmi care utilizeaz pai variabili.
3 Operatori ntr-un singur pas (Taylor, Euler, Runge-Kutta)
3.1 Serii Taylor, eroare de trunchiere, ordin al metodei
Se desvolt x(t) n serie Taylor n jurul lui t pn la termenul de ordinul p. De
exemplu pentru p = 3, avem:
)(!3
)(!2
)()()( )3(32
txh
txh
txhtxhtx (4)
Ecuaia (1) este
),( xtfx
a.ch. Noiembrie 2008
5
i prin derivare succesiv obinem:
fxxffx xt ;
xxfxffx xxttt ;)3(
Eroarea n desvoltarea (4) este dat de restul seriei Taylor
),();(!4
1 )4(44 httxhT
Eroarea 4T se numete eroarea de trunchiere local. Derivata )4(x n se poate
aproxima prin derivata n t, i aceasta din urm prin diferena divizat, obinnd
estimarea )]()([!4
1 )3()3(34 txhtxhT .
n general, considernd desvoltarea pn la ordinul p 1, eroarea de trunchiere
local este
),();()!1(
1 )1(1 httxhp
T ppp
sau
)( 1 pp hOT
Eroarea de trunchiere global pe este eroarea produs de eroarea local n calculul
lui )( ntx , adic eroarea dup n pai unde httn n /)( 0 i ea va fi de ordinul
pnT , adic de ordinul ph . Avem urmtoarea
Definiie: Ordin
(1) Dac eroarea de trunchiere global este de ordinul ph , metoda (sau
operatorul) se zice de ordinul p
Definiii echivalente ale ordinului sunt urmtoarele:
(2) Metoda este de ordinul p dac formula metodei coincide cu seria Taylor
trunchiat pn la termenul de ordinul p inclusiv
(3) Metoda este de ordinul p dac formula metodei este exact pentru un
a.ch. Noiembrie 2008
6
polinom de gradul p (i nu mai este exact, pentru un polinom de gradul
1p ).
Formula (3) a metodei se zice exact pentru o funcie )(tx dac, din
ipoteza c n membrul doi avem 1,),( kiijtxx jj , rezult ca avem
i )( 11 ii txx
n cazul de fa, formula metodei este chiar seria Taylor (4) trunchiat, scris
pentru 1, ii thttt , i anume:
)()3(32
1!!3!2
p
i
p
iiiii xp
hx
hx
hhfxx , i 0 (5)
n care ),( iii xtff , iar ,,)3(
ii xx reprezint derivatele calculate n it .
Avantaje i dezavantaje ale metodei seriei Taylor
c) Avantajele sunt simplitatea metodei i precizia mare care poate fi atins.
Precizia crete cu ordinul p, dar calculul cere evaluarea a mai multor derivate.
d) Dezavantajul principal const n calculul derivatelor de ordin superior. Mai
mult, trebuie ca funcia f s aib derivate pn la ordinul p, ceea ce, n general,
nu este cerut pentru existena soluiei. Totui, pentru multe din problemele
practice, aceast condiie este realizat.
3.2 Metoda Euler
Metoda Euler corespunde cazului n care p = 1. Formula metodei este, cf. (4),
),(1 iiii xthfxx (6)
Metoda are avantajul c nu cere dect calculul lui f. Ordinul ei este p = 1, i pentru
a atinge o precizie convenabil, pasul h trebuie luat foarte mic. Metoda are mai
degrab o importan teoretic. Ea servete la demonstrarea teoremelor de
existen, i la exemplificarea noiunilor de convergen i stabilitate pe exemplul
unei metode simple.
a.ch. Noiembrie 2008
7
3.3 Metode Runge-Kutta
3.3.1 Construcia metodelor Runge-Kutta
Metodele Runge-Kutta (abreviat RK) utilizeaz desvoltarea n serie Taylor, dar
nlocuiesc calculul derivatelor de ordin superior, cu calculul funciei f n puncte
de forma ),( hxht , unde i sunt definii de coeficienii metodei.
Relund desvoltarea Taylor cu rest, avem:
)(!!2
)()( 1)1(2
1
pp
i
p
iiii hOfp
hf
hhftxtx (7)
n care s-a inut cont de ))(,()( txtftx , iar itti
ff | i itt
nnn
i dtdff )/( )()( .
Reamintim c notm prin ix soluia calculat n it , prin )( itx soluia exact, i c
punem condiia )( ii txx (pn la termenul de ordinul p n h).
O caracteristic a metodei este numrul de evaluri al membrului doi al ecuaiei
(1) sau sistemului (2), pe un pas. Acest numr este numit numrul de evaluri de
funcii. O metod RK care face q evaluri de funcii va fi numit cu q-trepte
(q-stage). Pentru a obine o metod cu q-trepte, punem:
),,(1 hxthxx iiii (8)
n care
q
m
mmii khxt1
),,( (9)
unde m sunt coeficieni ai metodei, iar ),,( hxtkk iimm . Se obine
q
m
mmii khxx1
1 (10)
n (9, 10) funciile mk se definesc astfel:
a) Pentru o metod explicit:
1
1
),(m
j
jmjimim khxhtfk (11)
i 01 , astfel c avem:
a.ch. Noiembrie 2008
8
),(1 ii xtfk ,
),( 12122 khxhtfk ii ,
etc.
b) Pentru o metod implicit:
q
j
jmjimim khxhtfk1
),( (12)
Coeficienii m se mai zic noduri, iar m se mai zic ponderi.
Se obinuiete ca coeficienii mjm , i m , s se dea n tabloul Butcher:
B
n care: Tq ][ 21 , ][ mjB , i ][ 21 q .
Pentru o metod explicit ( 01 , i 0mj pentru 1 mj ) tabloul Butcher
este:
qqqqq
121
1,21
32313
212
0
(13)
Condiii pentru coeficienii metodei:
- Coeficienii m ndeplinesc condiia de consisten:
11
q
m
m (14)
Aceasta asigur convergena metodei v. 3.3.3.
- Coeficienii mjm , sunt supui la condiiile:
m
m
j
mj
1
1
, qm ,2 (14')
a.ch. Noiembrie 2008
9
adic: 221 , 33231 , , qqqqq 1,21 .
Aceste condiii simplific deducerea coeficienilor pentru metodele de ordin mai
mare ca 2. Pentru justificri ale condiiilor (14'), v. Ralston & Rabinowitz (1978),
i Isaacson & Keller (1966).
Ordin:
Eroarea de trunchiere local 1iT , pe pasul 1i , se definete ca eroarea formulei
(8) a metodei, cnd nlocuim aproximaiile jx cu soluia exact )( jtx . Adic,
definim 1iT prin:
11 )),(,()()( iiiii Thtxthtxtx , (15)
Dac
)( 11
p
i hOT (15')
metoda se zice de ordinul p. (Mai precis, p este cel mai mare ntreg pentru care
avem (15')). Aceasta revine la condiia ca ca formula (8) s coincid cu seria
Taylor a lui )()( 1 htxtx ii , trunchiat pn la termenii de ordinul p n h
inclusiv. Pentru a obine o metod de ordin p, coeficienii mjm , i m se
determin din condiia de mai sus, cu respectarea condiiilor (14, 14').
Eroarea de trunchiere global pe pasul 1i , este eroarea aproximaiei 1ix , adic
111 )( iii xtxe .
n 3.3.6 se va arta c )( 11
p
i hOT )(1p
i hOe . Astfel, o metod RK de
ordinul p are o eroare global de ordinul ph .
n ceea ce urmeaz vom analiza numai metodele RK explicite. Pentru metodele
implicite trimitem la Hairer & Wanner (1991).
Exemplu:
Metoda RK explicit, cu 2-trepte i de ordinul 2 ( 2q i 2p ). Se obine o
familie cu un parametru, de metode explicite RK cu 2-trepte, de ordinul doi,
definite de formulele:
a.ch. Noiembrie 2008
10
])1[( 22121 kkhxx ii
),(1 ii xtfk
)2
,2
( 122
2 kh
xh
tfk ii
Metode cunoscute se obin cu 1,,43
21
2 . De exemplu, pentru 12 metoda se
zice metoda Runge de ordinul 2, iar tabloul Butcher (13) este:
10
0
21
21
3.3.2 Ordin i numr de trepte (evaluri de funcii / pas )
Se arat c, n general, pentru ca o metod explicit s aib ordinul p, ea trebuie s
aib q p trepte, i anume: pentru p = 1, 2 ,3 4, avem pq min ; pentru p > 4,
pq min . Mai precis, avem urmtoarele rezultate datorate lui Butcher (Hairer,
Nrsett, & Wanner (1987)):
c) Pentru 5p nu exist metode explicite de ordin p, cu pq trepte.
d) Pentru 7p nu exist metode explicite de ordin p, cu 1 pq trepte.
e) Pentru 8p nu exist metode explicite explicite de ordin p, cu 2 pq
trepte.
Aceste rezultate sunt numite barierele Butcher. Pentru p = 9, 10 se cunosc
numai margini pentru minq , iar pentru 10p nu se cunosc evaluri pentru qmin.
Rezultatele anterioare se pot sintetiza n tabloul urmtor (Cartwright & Piro
(1992)):
p 1 2 3 4 5 6 7 8 9 10
)(min pq 1 2 3 4 6 7 9 11 12 17 13 17
Ordinul maxim pentru care avem q = p este p = 4. Din acest motiv, metoda RK de
ordinul 4 este cea mai frecvent utilizat. (Pentru 4p trebuie adugate cel puin
a.ch. Noiembrie 2008
11
dou trepte, ceea ce mrete timpul de calcul i introduce erori de rotunjire
suplimentare.). n ceea ce privete metodele RK implicite, pentru orice numr de
trepte q, exist metode de ordinul p = 2q. V. Hairer, Nrsett, & Wanner (1987).
3.3.3 Convergen i consisten
Metoda RK se zice convergent dac, pentru 0h , soluia calculat tinde la
soluia exact (pe fiecare it ). Considernd intervalul de integrare ],[ 0 itt , i notnd
0ttc i , numrul de pai de integrare va fi hci / , sau ih = c. Astfel, condiia
se exprim prin limita:
)(lim0
ii
cihh
txx
Metoda RK, definit de (8) se zice consistent (cu problema cu valori iniiale)
dac avem
))(,()0),(,( iiii txtftxt (20)
Cu (20) i expresiile (11, 12) ale funciilor mk , avem
q
m
miihm
q
m
mii txtfktxt1
0
1
))(,(|)()0),(,(
i condiia de consisten (20) este echivalent cu condiia
11
q
m
m (21)
Se demonstreaz c, consistena este o condiie necesar i suficient pentru
convergen (Cartwright & Piro (1992)).
3.3.4 Metode RK de ordinul 4
O metod RK explicit, de ordinul 4 (abreviat RK4), este definit de (10) cu q = 4:
)( 443322111 kkkkhxx ii
n care, conform (11), avem:
),(1 ii xtfk ,
),( 12122 khxhtfk ii ,
))(,( 23213133 kkhxhtfk ii
a.ch. Noiembrie 2008
12
))(,( 34324214144 kkkhxhtfk ii
Deducerea coeficienilor metodei conduce la o familie cu doi parametri v.
Hairer, Nrsett, & Wanner (1987), Ralston & Rabinowitz (1978). Cele mai uzuale
metode RK4 sunt definite de urmtoarele tablouri Butcher:
Metoda RK4 Regula 3/8
61
62
62
61
21
21
21
21
1001
0
0
81
83
83
81
21
31
32
31
31
1111
0
Se verific condiia de consisten (21) i condiiile (14'). Prima metod este cea
mai uzual, fiind denumit Metoda RK de ordinul 4. A doua este ceva mai
precis dect prima (Hairer et al. (1987)).
Explicit, Metoda RK4 este dat de formulele:
)22(6
43211 kkkkh
xx ii (22)
n care:
),(
),(
),(
),(
34
221
21
3
121
21
2
1
hkxhtfk
hkxhtfk
hkxhtfk
xtfk
ii
ii
ii
ii
(23)
Pentru un sistem de ecuaii difereniale (de ordinul nti), formulele metodei RK4
sunt similare cu (19, 20), variabilele scalare x, f, k, nlocuindu-se cu vectorii x, f,
k:
)22(6
43211 kkkkxx h
ii (22a)
),(
),(
),(
),(
34
221
21
3
121
21
2
1
kxfk
kxfk
kxfk
xfk
hht
hht
hht
t
ii
ii
ii
ii
(23a)
a.ch. Noiembrie 2008
13
n programarea formulelor (22a, 23a) vectorii se reprezint prin tablouri:
x(0:n), f(1:m), k(1:m). Reamintim c m desemneaz numrul de ecuaii,
iar n numrul pailor de integrare.
Metode RK de ordin mai nalt
Cel mai nalt ordin pentru care s-au construit metode RK explicite este 10: Curtis
(18 trepte, 1975) i Hairer (17 trepte, 1978).
3.3.5 Metode RK mbricate
Fie o metod RK de ordin p, cu q trepte, care calculeaz soluia
q
m
mmii khxx1
1 (24)
Funciile mk sunt definite de (11) i revin la calculul funciei f n puncte de forma
1
1
),(m
j
jmjimi khxht .
Ideea metodei mbricate este de a o combina metoda (24) cu o metod RK de
ordin p' (uzual 1 pp sau 1 pp ), cu acelai numr de trepte q, i care s
calculeze funcia f pe aceleai puncte ca (24) adic avnd aceeai coeficieni
mjm , . Fie cea de-a doua metod, care calculeaz soluia
q
m
mmii khxx1
1 . (25)
n (24) i (25), q este numrul de trepte din metoda de ordin mai mare. Pentru
fixarea ideilor s presupunem c pp : atunci )(min pqq , iar metoda de ordin
p va avea numrul de trepte )(min pqq . Astfel, metoda de ordin mai mic are
trepte sau grade de libertate suplimentare. Coeficienii metodei imbricate se
determin astfel ca ei s minimizeze coeficienii care definesc eroarea n una din
cele dou metodele. Soluia 1 ix se utilizeaz pentru estimarea erorii de trunchiere
prin (citete egal prin estimare):
11
iip xxT , (26)
O astfel de metod se va nota RK )( pp exemplu RK 4(5).
Explicit, eroarea de trunchiere local se estimez prin
a.ch. Noiembrie 2008
14
q
m
mmmp khT1
)( (26')
Metode Runge-Kutta-Fehlberg:
Fehlberg a construit astfel de metode de ordine )1( pp , care s minimizeze
coeficienii erorii n metoda de ordin mai mic. Ele sunt numite metode Runge-
Kutta-Fehlberg (RKF). Cele mai cunoscute sunt metodele RKF 4(5) i 7(8). Cea
mai utilizat dintre acestea este metoda de ordinul 4 cu 6 trepte (p = 4, p' = 5,
6)( pq ), definit de urmtorul tablou Butcher n ultima linie sunt dai
coeficienii m :
Metoda RKF 4(5)
552
509
5643028561
128256656
13516
51
41042197
25651408
21625
4011
41041859
25653544
278
21
4104845
5133680
216439
21977296
21971932
1312
329
323
83
41
41
0
00
2
81
0
Metoda RKF 4(5) se gsete implementat n multe pachete de programe pentru
integrarea numeric a ecuaiilor difereniale. Implementarea conduce la o metod
RKF cu pas variabil: estimarea (26) se utilizeaz pentru a controla eroarea
metodei (24) i a modifica pasul dac eroarea depete o toleran impus v.
mai jos.
Metode Dormand-Prince (DOPRI):
Dormand & Prince au construit metode mai precise, de ordine )(1 pp , n care se
minimizeaz coeficienii erorii n metoda de ordin mai mare. Soluia calculat este
dat de metoda cu ordinul 1p , iar metoda de ordinul p se utilizeaz numai
pentru controlul pasului. Acestea sunt metodele DOPRI 5(4) ordin 5, cu 7 trepte,
i DOPRI 8(7) ordin 8, cu 12 trepte. Coeficienii metodelor, ca i coduri Fortran,
sunt date n tratatul Hairer, Nrsett, & Wanner (1987). Codurile Fortran gsesc i
a.ch. Noiembrie 2008
15
la adresa: http://www.unige.ch/math/folks/hairer/software.html.
Metode DOPRI sunt prezentate n Dormand (1996). Coduri Fortran sunt date la
adresa: ftp://ftp.tees.ac.uk/pub/j.r.dormand/.
Metodele DOPRI sunt cele mai precise metode explicite pentru integrarea
numeric a ecuaiilor difereniale de ordinul nti, existente n momentul de fa.
Alegerea pasului
Considerm o metod imbricat cu pp , i un pas curent, notnd pentru
simplificare ixx 0 i 11 ixx . Eroarea de trunchiere local estimat conform
(26) este:
11 xxTp
Avem:
)()()()( 11'1001 ppp hOhOxhtxhtxxT
sau, cu pp , avem eroarea absolut:
1|| pp ChTerr
Pasul optim este cel pentru care eroarea este aproximativ egal cu tolerana tol
specificat de utilizator, adic:
1 poptChtol ,
Eliminnd C ntre ultimele dou realii, rezult:
1
p
opt
h
h
err
tol
din care,
herr
tolh
p
opt
1
1
Pentru siguran, n program se pune:
herr
tolh
p
opt
1
1
9.0
a.ch. Noiembrie 2008
16
n formula anterioar, |||| 11 pTxxerr unde pT este estimarea (26') a erorii.
Pentru un sistem, modulul se nlocuiete cu norma: |||| 11 xx err .
Observaii
1) Dac se cere specificare toleranei tolrel la eroarea relativ n modul rel, atunci
avem || 1 errx
errrel
i rezult
|| 1
1
errx
Chrel
p
, || 1
1
errx
Chtolrel
p
opt
, de unde
1
p
opt
h
h
rel
tolrel,
hrel
tolrelh
p
opt
1
1
n formula anterioar, rel este dat de expresia de mai sus n care err este estimarea
(26') n modul. Pentru un sistem, avem ||
max1 jj
j
j errx
errrel
.
2) Eroarea err se mai estimeaz i prin aa numita extrapolare Richardson,
calculnd n paralel, soluia 2x cu doi pai de mrime h i soluia 2X cu un pas
dublu 2h, i estimnd eroarea prin diferena celor dou soluii. Pentru o metod de
ordinul p se obine (Hairer et al. (1987)):
)(12
)2( 22220
p
phO
Xxxhtx ,
12
22
pp
XxT .
Soluia
pTxx 22
este o aproximaie a lui )2( 0 htx , cu o eroare de ordinul p+1. Pentru controlul
erorii avem:
a.ch. Noiembrie 2008
17
12
|| 22
p
Xxerr ,
|| 2x
errrel .
Pentru sistem, n err modulul se nlocuiete cu norma, iar ||
max2 j
j
j x
errrel , unde
jerr este estimarea err pentru coordonata j a soluiei. Estimrile err, rel se pot
folosi n formulele anterioare pentru opth .
3) Codurile care implementeaz metode cu pas variabil utilizeaz fie tol, fie
tolrel, fie ambele, mpreun cu alte mecanisme de control al pasului care previn
creterea sau scderea excesiv a pasului. De exemplu, n unul din cele mai noi
coduri v. RKSUITE n Brankin and Gladwell (1994), utilizatorul specific
tolerana TOL a erorii relative, iar testul de eroare cere ca pe fiecare pas i:
))(),(max(|)(| jpragjmagTOLjeroare
unde )( jmag este o mrime medie a coordonatei j a soluiei ix pe pasul
considerat, iar prag este un tablou specificat de utilizator. Astfel, dac prag(j) >
mag(j) rezult un test de eroare absolut cu tolerana tol = TOLprag(j), iar pentru
prag(j) < mag(j) rezult un test de eroare relativ cu tolrel = TOL.
3.3.6 Estimarea erorii de trunchiere globale
Notm acum cu 1ie i 1iT , modulul erorii de trunchiere local, i respectiv
global, pe pasul i+1. Dup definiiile din 3.3.1, avem:
|)(| 111 iii xtxe , (27)
i definim 00 e , i
|)),(,()()(| 11 htxthtxtxT iiiii (28)
Se arat c: eroarea de trunchiere global este )( pi hOe .
Eroarea de trunchiere local (n modul) se poate scrie sub forma
)())(( 211
pp
ii hOhtxT (30)
Primul termen se zice eroarea de trunchiere local principal.
a.ch. Noiembrie 2008
18
Pentru un sistem, n expresiile anterioare modulul se nlocuiete cu norma.
3.3.7 Stabilitatea metodelor RK (stabilitatea absolut liniar)
Ne vom limita la stabilitatea liniar a metodei. Aceasta se studiaz prin
liniarizarea ecuaiei (1) n jurul unei soluii a acesteia. Fie ecuaia diferenial
),()( xtftx (31)
i )(tx o soluie neted a acesteia, adic
))(,()( ttft . (32)
Considerm o perturbaie )(tx a soluiei (provenind dintr-o perturbare a condiiei
iniiale), unde |)(| tx :
)()()(),()()( txttxttxtx
i scznd relaia (32) din (31) rezult
))(,())()(,()( ttftxttftxdt
d (33)
Desvoltm membrul doi n (33) n jurul lui )(t pn la termenul de ordinul nti
n )(tx . Obinem:
)()()())(,()( txtJtxtt
x
ftx
dt
d (34)
n care ))(,(|)/()( ttxftJ . Ecuaia (34) liniarizat se obine neglijnd termenii
nescrii, i anume:
)()()( txtJtxdt
d
n fine, n prim aproximaie, considerm JtJ )( = constant (i anume
)( *tJJ , unde ),(* httt ) i avem
)()( txJtxdt
d (35)
Ecuaia (35) poate fi scalat, astfel c perturbaia s fie de mrime arbitrar.
Punem
)()( txCty , unde C este o constant, i ecuaia (35) devine
a.ch. Noiembrie 2008
19
Jyy . (35')
n fine, notnd x n loc de y , ecuaia (35') devine o ecuaie de tipul
xx (36)
n care, n general, va fi considerat complex (v. mai jos, cazul unui sistem).
Ecuaiei (36) i atam o condiie iniial arbitrar:
)0(
0 )( xtx (36')
Problema (36, 36') constituie testul pentru stabilitatea liniar a metodei numit i
testul Dalquist. Soluia exact a problemei este:
)()0( 0)(tt
extx (37)
Dac 0)Re( , atunci avem 0)( txt . Se zice c problema are un
punct fix stabil, n x = 0.
Observaie
Punctele fixe ale unei ecuaii (31) sunt valorile x pentru care avem ,0),( xtf
0tt . Punctele fixe ale unei metode numerice explicite )(1 ii xgx , sunt date
de ii xx 1 (adic de soluiile ecuaiei )(xgx ). Punctele fixe ale unei metode
RK definit de (8) sunt date de
0),,(1
q
m
mmii khxt
n care, pentru o metod explicit:
1
1
),(m
j
jmjimim khxhtfk
Dac X este un punct fix al ecuaiei ),( xtfx , atunci avem 0,0),( ttXtf ,
i rezult: 0),(1 Xtfk , 0),(2 Xtfk , , sau qmkm ,1,0 . Pentru o
metod implicit avem acelai rezultat. Urmeaz c punctele fixe ale ecuaiei sunt
i puncte fixe ale metodei RK
Definiie
O metod numeric este stabil liniar dac, aplicat ecuaiei (36) avem:
a.ch. Noiembrie 2008
20
0 ixi ,
adic metoda pstreaz stabilitatea punctului fix x = 0
Pentru un sistem de m ecuaii difereniale (3)
),( xfx t
avem, analog cu cazul unei singure ecuaii: fie soluia )(t
))(,()( ttt f ,
punem
)()()( ttt xx
i rezult
)()(),(),()( tttttdt
dxJfxfx
n care ))(,(|]/[)( ttki xft J este jacobianul funciei f n raport cu x. Aproximm
J(t) = A = constant. Cu aceasta, schimbnd notaia xx , modelul liniar este
)0(
0 )( xx
Axx
t (38)
n care A este o matrice constant mm. Presupunem, pentru simplificare, c A c
are valori proprii mjj ,1, distincte (i, n general, complexe). Presupunem c
valorile proprii au partea real negativ: atunci avem un punct fix stabil n x = 0.
ntruct valorile proprii sunt distincte, exist o baz ortogonal format din
vectorii proprii n care matricea A se diagonalizeaz, iar ecuaiile (38) se
decupleaz, sistemul reducndu-se la m ecuaii independente de forma (31). ntr-
adevr, dac vectorii proprii sunt },1,{ mjj v , definii de jjj vAv , punem
j
jjy vx i avem j
jjy vx )( , j
jjj y vAx . nlocuind n (35) rezult
mjyy jjj ,1,
la care adugm condiii iniiale arbitrare, de exemplu
mjyty jj ,1,)()0(
0
a.ch. Noiembrie 2008
21
Astfel, n ipotezele fcute, analiza stabilitii pentru sistemul (38) se poate face pe
o singur ecuaie de forma (36)
Revenind la ecuaia (36) s-i aplicm metoda explicit RK de ordinul 2,
considerat n 3.3.1:
])1[( 22121 kkhxx ii .
Cu xxtf ),( , rezult ixk 1 , )2
(2
2 ii xh
xk
, i
)2
()]2
()1[( 2
2
221 iiiiiiii xh
xhxxh
xxhxx
Avem:
ii xhRx )(1 ,
unde
2
)(1)(
2
hhhR
S observm c, cu 10 x , avem )(1 hRx .
Definiie
R(h) se numete funcia de stabilitate a metodei. Ea poate fi considerat ca
soluia numeric dup un pas, a problemei liniare de test (36, 36'), cu 1)0( x
Regiunea de stabilitate absolut pentru o metod, este mulimea valorilor h i
(h = real i nenegativ, = complex), pentru care avem 0ix pentru i ,
adic punctul fix x = 0 (originea) este stabil. Pentru aceasta este necesar i
suficient ca s avem 1|| R . Punnd hz , regiunea de stabilitate este mulimea
}1|)(|;{ zRzS C .
Uneori, regiunea de stabilitate este definit mpreun cu frontiera sa, prin condiia
1|| R . Pentru metoda explicit RK de ordinul 2, regiunea de stabilitate va fi dat
de:
1|2/1| 2 zz
a.ch. Noiembrie 2008
22
Pentru un sistem, va fi valoarea proprie de modul maxim a matricii jacobian A.
S considerm acum, cazul general al unei metode explicite cu p trepte, de ordinul
p (adic, p 4). Considerm desvoltarea lui x(t) n serie Taylor, pn la ordinul p.
Cu x' = x, rezult xxx 2 , i n general, prxx rr ,)( , astfel c
avem: )()(!
)(!2
)()()( 122
1
p
i
pp
iiii hOtxp
htx
htxhtxtx
n fine, cu )( ii txx , i omind restul )(1phO , avem:
i
pp
i xp
hhhx )
!!21( 2
2
1
care arat c funcia de stabilitate este
!
)(
!2
)(1
2
p
hhhR
p ; p 4. (39)
Pentru o metod explicit de ordin p, cu pq trepte, funcia de stabilitate va fi
q
pj
j
j
p
hp
hhhR
1
2
)(!
)(
!2
)(1
, (40)
unde j sunt definii de coeficienii metodei. De exemplu, pentru metoda DOPRI
5(4) cu 6 trepte (treapta 7 se utilizeaz numai pentru estimarea erorii) termenul
adiional n (40) este (h)6/600. (Hairer & Wanner (1991).
Din aceasta rezult c funcia de stabilitate a unei metode explicite cu q trepte este
un polinom de gradul q n h. Condiia | R | < 1 conduce la o regiune de stabilitate
mrginit. (Dac aceasta ar fi nemrginit nu putem avea | R | < 1, ntruct
|||| Rh ). Metodele RK implicite pot avea regiuni de stabilitate
nemrginite. Aceste metode se aplic pentru ecuaii difereniale rigide v. 5, la
care metodele explicite nu mai convin.
Pentru reprezentarea regiunilor de stabilitate liniar n cazul = complex, pentru
metodele RKp, p =1, 2, 3, 4 v. Hairer & Wanner (1991), Cartwright & Piro
(1992). Interseciile regiunilor cu axa real dau intervalele de stabilitate pentru
cazul = real. Aceste intervale se gsesc din condiia | R | < 1, unde R este definit
de (37) cu = real i negativ (conform ipotezei Re() < 0). Rezult:
a.ch. Noiembrie 2008
23
p = 1, 2: -2 < h p = 3: -2.512745 < h p = 4: -2.785296 < h
3.3.8 Stabilitatea absolut neliniar
Stabilitatea neliniar este o problem mult mai complex. Ea are conexiune cu
dinamica haotic. n cazul unei probleme neliniare, regiunea de stabilitate a unei
metode RK poate fi diferit de regiunea ei de stabilitate liniar. Cea mai
important diferen const n aceea c, pentru o problem neliniar, metodele RK
pot conine pe lng punctele fixe ale problemei v. Observaia din 3.3.7 i
puncte fixe adiionale. Excepie face metoda Euler care are numai punctele fixe
ale problemei. Punctele fixe adiionale sunt numite puncte fixe fantom. Recent
(1991) s-a artat c, n unele cazuri, puncte fixe fantom pot exista la orice
lungime a pasului (diferit de zero), adic la pai pentru care h este n regiunea
de stabilitate liniar absolut. Dac un asemenea punct fix este stabil la pai orict
de mici, atunci o traiectorie (calculat) poate converge la un punct fix care nu
exist n dinamica problemei originale. Diferena ntre problemele liniare i
neliniare const n aceea c, pentru probleme neliniare bazinul de atracie este
mrginit, n timp ce pentru o problem liniar acesta este nemrginit. Astfel,
pentru o problem liniar exist convergen pentru orice condiii iniiale, cu
condiia ca h s fie n interiorul regiunii de stabilitate, n timp ce pentru o
problem neliniar este necesar, n plus, ca condiiile iniiale s fie coninute n
bazinul de atracie. Pentru desvoltri, trimitem la Cartwright & Piro (1992).
n practica de calcul s-a constatat c pentru un rspuns haotic, unde calculaia se
face pe un mare numr de pai (sute de mii sau milioane), codul Runge-Kutta de
ordinul 4 formulele (22a, 23a) este foarte sensibil la mici schimbri ca:
utilizarea de variabile locale, asocierea n operaiile aritmetice, vectorizarea
ciclurilor DO n subrutina de integrare a sistemului dat, opiunile de build (ca
optimizarea codului), etc. V. raportul Chisli A. & al. (1998).
3.3.9 Exemplu de test Problema celor dou corpuri
Urmtoarea problem, constituit de problema celor dou corpuri n cazul micrii
eliptice, este luat ca test pentru metodele de integrare numeric a problemei cu
valori iniiale v. Dormand and Prince (1978), Brankin and Gladwell (1994).
Problema consider micarea relativ a dou puncte materiale care interacioneaz
a.ch. Noiembrie 2008
24
prin legea atraciei universale, i este descris, n coordonate carteziene, de
sistemul de ecuaii difereniale:
33 /,/ ryyrxx ,
n care 2/122 )( yxr . Se consider conditiile iniiale pentru cazul micrii
eliptice )1/()1()0(,0)0(,0)0(,1)0( eeyyxex ,
n care e < 1. Soluia analitic este dat de:
ue
uey
ue
uxueyeux
cos1
cos1,
cos1
sin,sin1,cos
22
n care u se determin din ecuaia lui Kepler: tueu sin . Soluia este periodic
cu perioada minim T = 2, iar orbita este o elips cu excentricitatea e i semi-axa
mare egal cu 1. Problema reprezint un test sever, datorit periodicitii soluiei.
Pentru rezolvarea numeric, sistemul dat se transform ntr-un sistem echivalent
de 4 ecuaii de ordinul nti:
2/12233 )(;/,/,, yxrrywrxvwyvx
cu condiiile iniiale: )1/()1()0(,0)0(,0)0(,1)0( eewvyex .
Calculm soluia pe intervalul [0, 20], adic peste trei perioade, pentru valorile
e = 0.1 i e = 0.9 ale excentricitii. Calculul este fcut n dubl precizie, cu
metodele:
(a) RK4 (pas constant) v. codul n ANA_EcDif.
(b) Runge-Kutta-Verner 5(6), cu subrutina DIVPRK din IMSL (pas variabil) v.
IMSL Libraries Reference (1998) cu argumentele: tol = 1D-7 2.23 D-
16, param(10) = 1 (se utilizeaz norma- a erorii). Subrutina se bazeaz pe
codul scris de Hull, Enright i Jackson (1976), care utilizeaz formulele lui
Verner de ordinul 5 i 6 v. DVERK, in site-ul:
http://www.cs.toronto.edu/NA/index.html. Rutina poate utiliza
pai n plaja 2.22D-15 2.0 (valori implicite). Intervalul de integrare s-a
mprit n 20, i respectiv n 200, sub-intervale. Rezultatele mai precise se
obin pentru mprirea n 200 sub-intervale n acest caz pasul maxim posibil
este 0.1 (egal cu lungimea sub-intervalului).
a.ch. Noiembrie 2008
25
(c) RK 8(7), cu subrutina DIVMRK din IMSL (pas variabil). S-a utilizat apelul
subrutinei DI2MRK, cu specificarea argumentelor. Tolerana tol s-a luat n
plaja 1D-7 2.23 D-15. Subrutina implementeaz codul din RKSUITE
metodele RK de ordine 3(2), 5(4), i metoda Dormand i Prince de ordin 8(7)
v. Brankin and Gladwell (1994). Ordinul metodelor este 3, 5, i 8, respectiv.
Intervalul de integrare s-a mprit n 20, i respectiv n 200 sub-intervale.
Rezultatele mai precise se obin pentru mprirea n 20 sub-intervale n
acest caz pasul maxim posibil este 1.0 (lungimea sub-intervalului)..
Pentru soluia exact, ecuaia lui Kepler se rezolv prin metoda punctului fix, cu
tolerana eps =1D-13. n tabelele urmtoare este dat eroarea absolut maxim i
minim a soluiei calculate, la timpul cel mai apropiat de 3 perioade. n paranteze
se indic funcia dintre x, y, x , y pentru care are loc extremul erorii absolute.
e = 0.1: Erori absolute extreme la t = 18.84 (RK4); 18.60 (RKV); 18.0 (RK 8(7)).
Extrem eroare
absolut
Metoda
RK4
h = 0.01
RKV 5(6)
tol = 1D-10
RK 8(7)
tol = 1D-10
Maxim
Minim
7.71 D-9 ( x )
4.38 D-11 ( x )
2.68 D-9 ( y )
4.49 D-10 ( x )
1.57 D-10 ( y )
9.63 D-11 ( y )
Numr pai
Nr. apeluri FCN
2000
8000
439
3512
138
2090
e = 0.9, Metoda RK4: Erori absolute extreme la t = 18.84 (h = 0.01; 0.005) i t =
18.849 (h = 0.001; 0.0005)
Pasul
h
Eroarea absolut
Maxim ( x ) Minim ( x )
Numr
de pai
0.01 3.52 D0 3.54 D-1 2000
0.005 7.12 D-1 4.26 D-3 4000
0.001 6.02 D-4 3.33 D-7 20000
0.0005 3.28 D-5 1.82 D-8 40000
Eroarea maxim are loc n y
a.ch. Noiembrie 2008
26
e = 0.9, Metoda RKV 5(6) 200 sub-intervale:
Erori absolute extreme la t = 18.60
Tolerana
tol
Eroarea absolut
Maxim ( x ) Minim ( y )
Numr
de pai
Numr
apeluri FCN
1 D-7 3.98 D-4 6.97 D-6 315 2779
1 D-10 1.28 D-6 2.23 D-7 622 5165
1 D-13 1.92 D-9 3.35 D-10 1800 14435
2.23 D-15 2.88 D-14 6.11 D-16 2244 17959
e = 0.9, Metoda RK 8(7) 20 sub-intervale:
Erori absolute extreme la t = 18.00
Tolerana
tol
Eroarea absolut
Maxim ( x ) Minim ( y )
Numr
de pai
Numr
apeluri FCN
1 D-7 2.16 D-6 1.31 D-7 152 2542
1 D-10 1.29 D-9 8.55 D-11 356 4984
1 D-13 9.00 D-13 5.96 D-14 847 11223
2.23 D-15 2.21 D-13 1.38 D-14 1380 18464
Tolerana minim admis = 2.22 D-15.
Urmtorul grafic d o comparaie a eficienei metodelor de mai sus, pentru cazul
e = 0.9. n ordonat este reprezentat numrul r = log10(Eroarea absolut minim).
El indic cea mai bun precizie atins de metod (eroarea minim este de ordinul
r10 ) i este reprezentat n funcie de numrul de evaluri de funcii. Metoda este
cu att mai eficient, cu ct realizeaz o precizie dat cu un numr mai mic de
evaluri de funcii.
a.ch. Noiembrie 2008
27
Problema celor dou corpuri, e = 0.9: Eficiena metodelor RKV 5(6), RK 8(7)
Observaii
- Testul cel mai sever este cazul e = 0.9. Cu acelai pas (n RK4), sau aceeai
toleran (n RKV 5(6), RK 8(7)), metodele dau rezultate cu o precizie
inferioar cazului e = 0.1. Din acest motiv, s-au efectuat integrri cu pai,
respectiv tolerane, mai mici. S remarcm c pasul h = 0.01 reprezint
aproximativ T/628, unde T este perioada micrii. n cazul e = 0.1, pentru
metodele RKV i RK 8(7), s-a ales tolerana 1D-10 pentru a avea erori
comparabile cu cele din metoda RK4.
- n subrutina DIPVRK (metoda RKV 5(6)), argumentul tol servete pentru
controlul normei erorii locale, n scopul de a se ncerca meninerea erorii
globale aproximativ proporional cu valorea tol (v. referinele citate mai sus).
- n subrutina DI2MRK (metoda RK 8(7)), argumentul tol servete la controlul
erorii relative, i la alegerea ordinului metodei astfel: 1D-4 < tol 1D-2, 1D-6
< tol 1D-4, i tol < 1D-6, produc alegerea metodei de ordinele 3(2), 5(4) i
8(7), respectiv.
- Coloana Numr apeluri FCN d numrul de apeluri ale subrutinei FCN care
calculeaz membrii doi ai sistemului de ecuaii. Acest numr este referit ca
a.ch. Noiembrie 2008
28
numrul de evaluri de funcii al metodei. Metoda RK4 face 4 apeluri ale
subrutinei FCN pe un pas (n codul din Anexa, 4.1, FCN este DERIVS.).
- Se remarc creterea preciziei odat cu micorarea pasului (RK4), sau a
argumentului tol (RKV 5(6), RK 8(7)), dar cu preul mririi numrului de pai
sau a numrului total de evaluri de funcii. Se remarc creterea preciziei cu
cretera ordinului metodei. Din comparaia eficienei celor trei metode rezult
c, pentru problema considerat, metoda RK 8(7) ofer cel mai bun raport
precizie/numr de evaluri de funcii cu excepia cazului unei tolerane
apropiat de cea minim admis (2.22D-15), cnd metoda RKV 5(6) este
superioar
a.ch. Noiembrie 2008
29
4 Operatori n mai muli pai
4.1 Definiii. Operatori liniari
Considerm din nou, problema cu valori iniiale (1,2)
),( xtfx ; )0(0 )( xtx (41)
pentru care se cere soluia pe intervalul ],[ 0 TTt , i nodurile Njt j ,0, (unde
TTtN ). Notm, ca nainte, cu )( jtx soluia exact i cu jx soluia calculat n
1, jt j , iar )0(
0 xx . Un operator n k-pai este o formul de forma
),,,,( 1111 kiiiii xxxxgx , (42)
care calculeaz soluia 1ix n funcie de k valori 1,,1,, kiiijx j
calculate anterior. La primul pas al metodei, aceste valori trebuie determinate
printr-o procedur special de start. Dac 1ix apare n membrul doi operatorul
este implicit, altfel este explicit. n ceea ce urmeaz vom considera numai cazul
operatorilor multi-pas liniari i cu pas constant h, adic:
- 1,0 jjhtt j ;
- Funcia g este liniar n jx i n ),( jjj xtff , 1,,,1 kiiij .
Pentru convenien, ecuaia (42) se scrie sub forma
)),(),(),(( 110111
101211
kikiiikiik
kiikiki
xtfbxtfbxtfbh
xaxaxax
(43)
n (43) vom presupune c cel puin unul dintre 0a , 0b este diferit de zero
(operatorul are k pai); n rest, oricare alt coeficient poate fi zero. Dac 0kb
operatorul este implicit, iar dac 0kb el este explicit. Condensat, (43) se scrie:
1,1
0 0
111
kifbhxaxk
l
k
l
lkillkili (43a)
a.ch. Noiembrie 2008
30
Pentru simplificarea notaiei indexate, s notm valorile curente cum urmeaz:
valoarea care se calculeaz cu kx ( 1,1 kiik ), i cele k valori anterioare cu
021 ,,, xxx kk . (Aceasta nu restrnge generalitatea, coeficienii ll ba ,
nedepinznd de i). Astfel, relaia (44) se scrie:
1
0 0
k
l
k
l
llllk fbhxax (44)
Forma general a unei metode multi-pas (44) este
k
l
k
l
llll fbhxa0 0
(45)
n care s-a pus 1,0, klaa ll .
n (45) se pun condiiile:
1ka (mai general, 0ka , pentru explicitare n raport cu kx )
000 ba (metoda are k pai).
Explicit:
)( 111100111100 kkkkkkkk xbfbfbfbhxaxaxaxa
Exemple:
1) Metoda mijlocului este definit de formula
1),,(211 ixthfxx iiii
i este un operator explicit n 2 pai.
2) Metoda trapezului este definit de
0)),,(),((2
111 ixtfxtfh
xx iiiiii
i este un operator implicit ntr-un singur pas
a.ch. Noiembrie 2008
31
4.2 Ordin
Formula (44) sau (45) se zice exact pentru o funcie )(tx dac, din ipoteza c
n membrul nti avem 1,),( kiiltxx ll , rezult ca avem i )( 11 ii txx
n limita erorilor de rotunjire.
Definiie
Dac formula (43) sau (44) este exact pentru polinoamele de grad p, zicem c
operatorul are ordinul p (sau, formula are ordinul de precizie p)
Lucrm pe forma (45)
k
l
k
l
llll fbhxa0 0
S presupunem c cerem ca (45) s aib ordinul p. n acest caz )())(,( txtxtf
este un polinom de grad 1p . Condiia pus revine la condiia c formula s fie
exact pentru polinoamele pqttx q ,0,)( (acestea alctuiesc o baz pentru
polinoamele de grad p). Avem 1)( qqttx . Putem presupune 00 t
(coeficienii nu pot depinde de 0t ), avem lhtl , i rezult:
0,)( qhltxx qqll ;
1,)( 11 qhqltxf qqll , i 0lf , pentru 0q .
innd cont de faptul c termenul al doilea n (44) conine ql qlhhf , urmeaz c
qh se va simplifica. Putem pune atunci, h = 1, i avem:
0, qlx ql ;
1,1 qqlf ql ; 0lf , pentru 0q .
nlocuind n (45), rezult, pentru q = 0, 1 i q = 2, , p:
1) q = 0 ( 0,1 ll fx ):
00
k
l
la (46a)
2) q = 1 ( 1, ll flx ):
a.ch. Noiembrie 2008
32
Rezult, innd cont de (45a):
000
k
l
l
k
l
l bl (46b)
3) 1;,,2 ppq ( ql lx , 1 ql ljf ):
00 0
1
k
l
k
l
l
q
l
q blqal , pq ,,2 (46c)
0)1(0 0
1
k
l
k
l
l
p
l
p blpal )1( pq (46d)
Rezumnd, avem condiiile:
)1(0)1(
_______________________________________
,,2,0
)1(0
)0(0
0 0
1
1
0 0
1
00
1
0
0
pqblpald
pqblqald
qblad
qad
k
l
k
l
l
p
l
p
p
k
l
k
l
l
q
l
q
q
k
l
l
k
l
l
k
l
l
(46)
n (46) avem 0,10 ll .
n particular, formulele explicite pentru q = 0, 1, 2, 3 sunt cele de mai jos, n care
sumele se opresc la termenii de indice k, inclusiv:
432100 aaaaad
)(432 4321043211 bbbbbaaaad
)432(21694 432143212 bbbbaaaad
)1694(364278 432143213 bbbbaaaad
Observaie
Cu coeficienii din (43), ll aa , 1,1 kl , i 1 kk aa condiiile (46) sunt:
a.ch. Noiembrie 2008
33
)1(0)1(
_______________________________________
,,2,0
)1(0
)0(01
1
1 0
11
1
0 0
1
0
1
0
1
1
0
0
pqblpalkd
pqblqalkd
qbalkd
qad
k
l
k
l
l
p
l
pp
p
k
l
k
l
l
q
l
q
k
l
l
k
l
l
k
l
l
Cu cele de mai sus avem urmtoarea propoziie:
Propozitie
Ordinul unei metode (45) este numrul natural p pentru care avem
0,,0,0 10 pddd i 01 pd .
Coeficienii qd sunt definii de (46)
Polinoamele generatoare ale metodei
Polinoamele obinute prin nlocuirile lll rfx , n cei doi membri din (45) se zic
polinoamele generatoare ale metodei, i anume:
k
l
l
l
k
l
l
l
rr
rr
0
0
)(
)(
(47)
Observai c:
)1(0 d ; )1()1(1 d
Consisten
O metod pentru care coeficienii verific primele dou relaii (46) adic,
condiiile pentru q = 0, 1:
0,1 10 dd ,
a.ch. Noiembrie 2008
34
se zic consistent. Aceasta echivaleaz cu condiia ca metoda s fie exact pentru
polinoame de gradul unu. Condiia de consisten se poate exprima i sub forma
urmtoare, utiliznd polinoamele generatoare:
)1()1(,0)1( . (48)
Exemple-2
1) Relum metoda mijlocului (Exemple-1), n care k = 2:
1,211 ihfxx iii
Avem: 2;1,0,1 0210 baaa . Rezult:
02100 aaad ; 0222 0211 baad
04)0(24 0212 baad
Astfel, ordinul metodei este p = 1.
2) Metoda trapezului (k = 1):
0),(21
121
1 iffhxx iiii
Avem: 21
121
010 ,;1,1 bbaa , i
00 d ; 0)(1 21
21
1 d
0)(21)(221
112 bad
031)(321
113 bad
Rezult p = 2
4.3 Construcia operatorilor n mai muli pai
Coeficienii n (44), (45) se determin astfel ca formula s fie exact pentru
polinoamele de un grad dat. Dac, din condiiiile puse, rmn coeficieni liberi
(parametri), acetia se vor determina astfel ca s avem ndeplinite una sau mai
multe din condiiile:
- Eroarea de trunchiere s fie ct mai mic;
- Propagarea erorilor s fie ct mai mic;
- Formula s fie ct mai simpl, de exemplu unii coeficieni s fie zero.
a.ch. Noiembrie 2008
35
n afar de aceste condiii, vom cere ca metoda s fie stabil, v. mai jos.
Determinarea coeficienilor n (45) se face prin una din urmtoarele metode:
- metoda coeficienilor nedeterminai
- prin integrare numeric
- prin derivare numeric
Acestea se expun n continuare.
4.3.1 Metoda coeficienilor nedeterminai
Relaiile (46) reprezint un sistem de 1p ecuaii n cel mult 12 k coeficieni
ll ba , (conform 1ka ). Dac numrul coeficienilor este egal cu 1p atunci,
sistemul poate fi rezolvat n ll ba , . Dac acest numr este mai mare dect 1p ,
unii coeficieni rmn ca parametri.
Exemple-3
1) S determinm metodele 1-pas, de ordin doi (k = 1, i p = 2):
)( 01101 iiii fbfbhxax .
n forma (45), metoda se scrie:
)( 01101 iiii fbfbhxax
Cu 11 a , ecuaiile (46) sunt: 010 a , 0)(1 10 bb , 021 1 b . Din
acestea rezult 10 a , 2/101 bb , astfel c metoda cutat este metoda
trapezului:
)(2
11 iiii ffh
xx .
2) Metode 2-pas, explicite, de ordin 1 (k = 2, p = 1; 01 b )
iiii fhbxaxax 01011 .
Sau, n forma (45): iiii fhbxaxax 01011 . Condiiile (46) sunt
110 01 aa , 01 0 b , astfel c metoda este
iiii hfxaxax 1001 )1( ,
a.ch. Noiembrie 2008
36
n care 00 a rmne un parametru
4.3.2 Metode bazate pe integrare numeric (metodele Adams i Milne)
Integrnd ecuaia (41) pe intervalul ],[ 1ii tt , avem:
dttxtftxtxi
i
t
t
ii
1
))(,()()( 1
Cu k aproximaii cunoscute 11 ,,, kiii xxx , pe nodurile 11 ,,, kiii ttt , gsim
aproximaiile funciei f pe aceste noduri:
1,),,( kiijxtff jjj .
Metode Adams explicite
n membrul doi, nlocuim funcia necunoscut ))(,( txtf cu polinomul de
interpolare Newton pe nodurile 11 ,,, kiii ttt . (k noduri, polinom de grad 1k ).
Se obine metodele Adams explicite, referite i ca metode Adams-Bashforth:
)( 3322110
1
0
1
iiiiik
l
lilii fcfcfcfchxfchxx (49)
Coeficienii lc din (49) sunt dai n tabelul urmtor, pentru k = 1, 2, 3, 4.
Coeficienii pentru lc pentru metodele Adams explicite ecuaia (49)
k if 1if 2if 3if 4if
1 1
2 3/2 -1/2
3 23/12 -16/12 5/12
4 55/24 -59/24 37/24 -9/24
5 1901/720 -2774/720 2616/720 -1274/720 251/720
a.ch. Noiembrie 2008
37
De exemplu, metoda pentru k = 4 este:
)9375955(24
3211 iiiiii ffffh
xx . Formulele (49) sunt de tipul (44),
cu 1,0,10 laa l , i 01 b , ll cb .
Ordin:
Metodele Adams explicite au ordinul p = k .
Metode Adams implicite
Analog, cu aproximaiile ),( jjj xtff pe nodurile 11 ,,, kiii ttt , utilizm
polinomul de interpolare pentru f .( 1k noduri, polinom de grad k).
Se obin metodele Adams implicite, referite i ca metode Adams-Moulton:
)( 23121100
11
iiiiik
l
lilii fcfcfcfchxfchxx (50)
Pentru k = 0, 1, 2, 3, 4, coeficienii din (50) se dau n tabelul urmtor:
Coeficienii lc pentru metodele Adams implicite ecuaia (50)
k 1if if 1if 2if 3if
0 1
1 1/2 1/2
2 5/12 8/12 -1/12
3 9/24 19/24 -5/24 1/24
4 251/720 646/720 -264/720 106/720 -19/720
Cazul special k = 0 produce metoda Euler implicit: 11 iii hfxx . Formulele
obinute sunt de tipul (44) cu 1,0,10 laa l , i 1 ll cb . Metodele
implicite au o precizie mai mare dect cele explicite. Determinarea lui 1ix din
formulele de mai sus, se face cu o metod pentru ecuaii neliniare (de exemplu,
pentru h = mic, metoda punctului fix).
Ordin:
a.ch. Noiembrie 2008
38
Ordinul metodei este 1 kp .
Metode Milne-Simpson (implicite)
Ecuaia (41) se integreaz pe intervalul ],[ 11 ii tt , obinnd
dttxtftxtxi
i
t
t
ii
1
1
))(,()()( 11
Sub integral, nlocuim f cu polinomul kp utilizat n metodele Adams implicite
(polinomul de interpolare pe nodurile 11 ,, kii tt ).
Se obin metodele:
k
l
lilii fchxx0
111 (51)
Coeficienii lc n (51), pentru k = 0, 1, 2, 3, 4, sunt dai mai jos.
Coeficienii lc pentru metodele Milne-Simpson ecuaia (51)
K 1if if 1if 2if 3if
0 2
1 0 2
2, 3 1/3 4/3 1/3
4 29/90 124/90 24/90 4/90 -1/90
Metoda pentru k = 2 este numit metoda Milne, i este dat de
)4(3
1111 iiiii fffh
xx .
Ordin:
Metodele Milne-Simpson au ordinul 1 kp
a.ch. Noiembrie 2008
39
4.3.3 Metode bazate pe derivare numeric (BDF)
n metodele anterioare s-a integrat ecuaia (41) i s-a utilizat polinomul de
interpolare pentru funcia f(t, x(t)). Considerm acum ecuaia (41), i polinomul de
interpolare pentru funcia x(t), pe nodurile 11 ,,, kiii ttt i valorile 11 ,, kii xx
(polinom de grad k):
Avem:
1
1
1
1
ik
j
i
j hfxj
(54)
Formulele (54) se numesc formule de derivare napoi (backward differentiation
formulae) i metodele cu aceste formule se zic metode BDF. Ele se utilizeaz n
integrarea numeric a ecuaiilor difereniale rigide v. 5.
Formulele (54) explicitate n termenii lix 1 sunt:
k
l
ilil hfxc0
11 (54')
Pentru 6,1k , coeficienii lc sunt dai mai jos. Pentru k > 6, metodele BDF sunt
instabile (Hairer et al. (1987)).
Coeficienii lc pentru metodele BDF ecuaia (54')
k 1ix ix 1ix 2ix 3ix 4ix 5ix
1 1 -1
2 3/2 -2 1/2
3 11/6 -3 3/2 -1/3
4 25/12 -4 3 -4/3 1/4
5 137/60 -5 5 -10/3 5/4 -1/5
6 147/60 -6 15/2 -20/3 15/4 -6/5 1/6
Ordin: Metodele BDF au ordinul p = k
a.ch. Noiembrie 2008
40
4.4 Stabilitatea metodelor n mai muli pai
Stabilitatea metodei analizeaz comportarea soluiei pentru n i 0h , cu
condiia nh = constant ( 0tTTnh = lungimea intervalului de integrare). Pentru
0h , din (43') se obine:
0111101 kikiii xaxaxax (55)
1ix din aceast ecuaie, poate fi interpretat ca soluia dat de metod, pentru
ecuaia diferenial 0x .
Pentru a rezolva ecuaia liniar (i omogen) cu diferene (55), cutm soluii de
forma jj rx . nlocuind, obinem
0111
10
1 ki
k
iii rararar ,
i mprind cu 1kir rezult
0)( 12
1
1
0
k
kkk arararr . (56)
Remarcm c )(r este primul polinom generator definit n (47). Fie kr ,1,
rdcinile polinomului )(r . Dac rdcinile sunt simple, un sistem de soluii
fundamentale este },,{ 1j
k
jrr , i soluia general este o combinaie liniar de
acestea:
1,,1,1
kiijrcxk
j
jj
.
Dac exist rdcini multiple r , cu ordinul de multiplicitate 1m , un sistem de
soluii fundamentale corespunznd rdcinii r sunt },,,{1
rjjrr
m . Soluia
general are forma
k
j
jj rjpx1
)(
,
n care )( jp j sunt polinoame de grad 1m n j. n ambele cazuri, pentru ca
soluia s rmn mrginit pentru n , trebuie ca rdcinile ecuaiei (56) s se
situeze n discul unitate ( 1|| r ), iar rdcinile de modul 1 s fie simple.
a.ch. Noiembrie 2008
41
Remarcm c, pentru metode consistente polinomul )(r are ntotdeauna o
rdcin 1r , conform (48).
Definiie
Metoda multi-pas (43) se zice stabil, dac polinomul generator )(r
satisface condiia de rdcini, i anume:
a) Rdcinile sunt situate n discul unitate: 1|| r ;
b) Rdcinile de modul 1 sunt simple: dac 1|| r , atunci 0)( r
Exemple
1) Stabilitatea metodelor Adams ( 4.2.2):
Polinomul generator pentru metodele Adams (explicite i implicite) este
1)( kk rrr . Radcinile sunt r = 1 simpl, i r = 0 multipl de ordinul
(k-1). Polinomul satisface condiia de rdcini i, n consecin, metodele Adams
sunt stabile.
2) Stabilitatea metodelor Milne-Simpson ( 4.2.2):
Polinomul generator 2)( kk rrr satisface condiia de rdcini, i deci,
metodele sunt stabile. Existena rdcinii r = -1 duce ns la fenomenul de
instabilitate slab v. Hairer & Wanner (1991).
3) Stabilitatea metodelor BDF ( 4.2.3):
Se arat c metodele BDF sunt stabile pentru k 6, i instabile pentru k 7.
Ordinul maxim al unei metode stabile
Ordinul maxim pentru care metoda este stabil este dat de urmtorul rezultat
datorat lui Dalquist, i numit prima barier Dalquist.
Propoziie
Ordinul p al unei metode liniare n k-pai stabil, satisface:
a) p k+2 pentru k = par
b) p k+1 pentru k = impar
c) p k pentru 01 b (n particular, pentru o metod explicit)
a.ch. Noiembrie 2008
42
4.5 Convergena metodelor n mai muli pai
Considerm din nou, problema cu valori iniiale (41), n care presupunem c
),( xtf satisface condiiile din seciunea 6-I, 1, i n consecin problema are
soluie unic. Considerm integrarea numeric pe intervalul ],[ 0 TTt , inclus n
intervalul de existen a soluiei. Reamintim c notm prin )( jtx i jx , respectiv
soluia exact i calculat pe punctul jhtt j 0 . Pentru ceea ce urmeaza vom
presupune c vrem s calculm soluia pe punctul fixat ],[ 0 TTtt , cu pai h din
ce n ce mai mici. Punem atunci: nhtt 0 constant, i avem httn /)( 0 ,
astfel c nh 0 . Notm cu )(hxt , soluia calculat pe punctul fixat t,
cu pasul h. Convergena va cere ca, sub anumite ipoteze, s avem limita:
)()(lim
fixat0
txhxtth
Practic, pentru a calcula )(hxt , utilizm metoda (43) n care punem 1 itt ,
1 in , i 01)1( tthi i constant. n acest caz vom scrie )(1 hxi n loc de
)(1
hxit
.
Definiie
1) Metoda liniar multi-pas (43) se zice convergent dac, pentru orice
problem cu valori iniiale (41), urmtoarea condiie este ndeplinit:
Dac valorile de start )(hx j (pe punctele jt , 10 kj ), satisfac
0)()( hxtx jj pentru 0h , 1,0 kj
atunci avem i
],[ 0 TTtt , t = fixat, 0)()( hxtx t pentru 0h .
2) Metoda se zice convergent de ordinul p dac, pentru orice problem (41)
cu f suficient derivabil, exist constantele pozitive CCh ,, 00 , astfel ca s
avem:
Dac valorile de start satisfac
a.ch. Noiembrie 2008
43
p
jj hChxtx 0||)()(|| pentru 0hh , 1,0 kj
atunci
],[ 0 TTtt , t = fixat, p
t Chhxtx ||)()(|| pentru 0hh .
Rezultatul principal din urmtoarele teoreme este c proprietile de consisten +
stabilitate ale unei metode, sunt condiii necesare i suficiente pentru convergena
acesteia: Convergen Consisten + Stabilitate.
Teorema 1
Dac metoda (43) este convergent, atunci ea este stabil i consistent
Teorema 2
1) Dac metoda (43) este stabil i de ordinul p = 1 (adic, consistent),
atunci ea este convergent.
2) Dac metoda (43) este stabil i de ordinul p, atunci ea este convergent de
ordinul p
4.6 Stabilitate relativ i stabilitate slab
Considerm ecuaia de test
1)0(, xxx
care are soluia tetx )( .
Definiie
Fie o metod (43) consistent, i 0r rdcina principal a polinomului
caracteristic (57). Metoda se zice:
- Relativ stabil pe intervalul ],[ 0, dac pentru orice h n acest
interval, rdcinile polinomului caracteristic satisfac condiiile:
1,1|,)(||)(| 0 khrhr (61a)
i, dac |||| 0rr , atunci r este rdcin simpl. (61b)
- Absolut stabil pe intervalul ],[ dac, pentru orice h n acest interval:
a.ch. Noiembrie 2008
44
1,0,1|)(| khr (62)
- Metoda satisface condiia tare de rdcini, dac:
1,1,1|)0(| kr (63)
O metod stabil dar care nu este relativ stabil, se zice slab stabil.
Observaii
- Stabilitatea absolut poate avea loc numai n cazul 0 (mai general, are
partea real negativ) v. relaia (60). n acest caz, stabilitatea absolut
echivaleaz cu condiia ca 0 jj xt .
- Condiia tare de rdcini implic stabilitatea relativ: cu 1)0(0 r , condiia
(63) se scrie )0(|)0(| 0rr , iar din continuitatea rdcinilor ca funcii de h
rezult c aceasta are loc pe o vecintate a lui 0. Reciproca nu este, n general
adevrat.
- ntruct definiiile anterioare se aplic i la sisteme, se consider, n general,
complex. Mulimea valorilor h pentru care au loc (60), respectiv (61), se zic
regiunea de stabilitate relativ, respectiv absolut, ale metodei. Determinarea
regiunii de stabilitate absolut este mai simpl dect cea de stabilitate relativ,
i ea se poate face pe criterii algebrice (Hurwitz-Routh, Schur v. Ralston &
Rabinowitz (1978)). Regiunile de stabilitate absolut (n planul complex),
pentru metodele Adams-Bashforth i Adams-Moulton sunt reproduse n
Atkinson (1978). Din aceste diagrame rezult c, cu ct ordinul este mai mare,
regiunea de stabilitate este mai mic v. i Exemple-3. Totui, chiar pentru
ordine mari, | h | rmne relativ mare, i nu introduce restricii majore asupra
lui h, cu excepia cazului n care | | este mare
Exemple
1. Metoda mijlocului este dat de 1),,(211 ixthfxx iiii , i cu
xxtf ),( devine iii xhxx 211 . Polinomul caracteristic este:
01)(22 rhr . Punem K , unde K > 0, i avem
a.ch. Noiembrie 2008
45
1)()( 21,0 hKhKhKr . Metoda este stabil (| r | < 1) pentru hK < 1, dar
avem |||| 01 rr , deci metoda nu este relativ stabil.
2. Metodele Adams: Verificm condiia (62). Pentru = 0, polinomul
caracteristic este 1)( kk rrr , cu rdcinile 10 r i 1,1,0 kr .
Condiia tare (63) este satisfcut i deci, metoda este relativ stabil.
3. S determinm intervalul (presupunem real, i negativ) de stabilitate
absolut pentru metodele Adams-Moulton k = 2, k = 3 (de ordinele 3, 4) v.
Tabelul din 4.2.2. jf se nlocuiete cu jx . Pentru metoda k = 2 avem:
12/)85( 111 iiiii xxxhxx . Punem h = z (z < 0) i
,,1 12 rxx ii , rezult zrzrzzp )812()512()(2 . Condiiile
pentru | r | < 1 sunt: 0 ; p(-1) > 0; p(1) > 0; -1 < (6+4z)/(12-5z) < 1, care
conduc la 6 < z < 0.
Pentru metoda k = 3 avem 24/)5199( 2111 iiiiii xxxxhxx . Cu
notaiile anterioare avem zzrrzrzzp 5)1924()924()( 23 . Condiii
necesare pentru | r | < 1 sunt (avem z < 0): p(-1) < 0; p(1) > 0, care conduc la
3 < z < 0. Se arat c acesta este rezultatul final. (Exerciiu: verificai c
0)( rp , 0z .)
4.6 Eroarea de trunchiere
Considerm metoda n k-pai (44), n care punem acum ))(,()( txtftx
1,1
0
1
1
1
kixbhxaxk
l
k
l
lillili (64)
n (64), jx i jx sunt aproximaii pentru )( jtx i )( jtx . nlocuind jx , jx cu
valorile exacte, formula va avea o eroare pe care o notm 1iT i care reprezint
eroarea de trunchiere local pe pasul i+1:
1
1
0
1
1
1 )()()(
ik
l
k
l
lillili Ttxbhtxatx (65)
Presupunem c avem 1iT = 0, pentru cazul cnd x este un polinom de grad p
(ordinul metodei este p).
a.ch. Noiembrie 2008
46
Se arat c eroarea 1iT are forma
),(),( 11)1(1
1
ikipp
pi ttxhdT (69)
Expresia coeficientului pd se d mai jos.
2
1
2
0
11 )1()1()1()!1(k
l
p
l
k
l
p
l
p
p lkbplkakdp (74)
Explicit:
]1)2()1()[1(
1)2()1()!1(
2101
2
1
1
1
0
1
k
ppp
k
ppp
p
bkbkbkbp
akakakdp
(74')
Exemple
1. Metode Adams-Bashforth (p = k):
Nr. pai
k
Ordin
p
Coeficient
pd
1 1 1/2
2 2 5/12
3 3 3/8
4 4 251/720
2. Metode Adams-Moulton (p = k+1):
Nr. pai
k
Ordin
p
Coeficient
pd
1 2 -1/12
2 3 -1/24
3 4 -19/720
4 5 -3/160
a.ch. Noiembrie 2008
47
4.7 Metode predictor-corector
4.7.1 Predictori i corectori
Eroarea unei metode de ordin p este dat de (69). Pentru o ecuaie dat, la acelai
pas i acelai ordin, mrimea erorii este dat de || pd , unde pd este dat de (74). n
general, || pd este mai mic pentru o metod implicit dect pentru una explicit.
Un exemplu l constituie metodele Adams compar valorile din exemplele 1 i 2
de mai sus. Astfel, este preferabil a se determina soluia cu o metod implicit.
Aceasta are forma general ),,( 111 kiii xxgx i determin soluia pe un pas
cu o metod iterativ pentru rezolvarea ecuaiei n .1ix Astfel, se cere o estimare
a aproximaiei iniiale (la fiecare pas). Cea mai bun cale este de a calcula aceast
aproximaie cu o metod explicit care va fi numit predictor. Apoi se va
corecta valoarea prin iteraie n metoda implicit aceasta va fi numit
corector. Metoda obinut prin cuplarea unui predictor i a unui corector, se va
numi o metod predictor-corector. Criterii pentru alegerea predictorului i
corectorului se vor discuta n 4.7.4.
4.7.2 Convergena iteraiei de punct fix
Explicitnd n membrul doi termenul n 1ix , metoda (64) se scrie:
)(),()( 1111
1
0
1
iiik
l
lillili xgxtfhbxhbxax (75)
S rezolvm (75) prin metoda punctului fix. Fie la pasul 1i (fixat) estimarea
)0(
1ix a lui 1ix , cu aceasta calculm ),()0(
11 ii xtf i nlocuim n (74) obinnd )1(
1ix ,
etc. n general, la pasul j al iteraiei avem:
0),,()( )( 111
1
0
)1(
1
jxtfhbxhbxax jiik
l
lillil
j
i (76)
i remarcm c de la pasul j la pasul j+1, suma din membrul doi nu se modific.
Condiia de convergen n (76) este 1|)(| xg , pe o vecintate a lui )0( 1ix .
Avem
a.ch. Noiembrie 2008
48
x
xtfhbxg
),()( 1 ,
Presupunnd c derivata xf / este mrginit ntr-o vecintate I a lui ),( )0( 11 ii xt
care conine punctele ),( )( 11j
ii xt :
Ixtx
xtf
),(,
,( (77)
rezult c trebuie s avem:
11 hb (78)
Mai mult, rata convergenei este 1hb , astfel c pentru o iteraie rapid vom cere
s avem 11 hb . n (78) s-a presupus 01 b , ceea ce are loc pentru toate
metodele considerate anterior. Pentru convergena pe orice pas (i+1) corespunznd
lui ],[ 01 TTtti , n (77) se va lua ],[ 0 TTtI .
Observaie
Pentru ecuaii difereniale rigide (v. 5), unde este mare, nu putem avea (77)
dect pentru un pas h excesiv de mic. n acest caz se va utiliza, n loc de iteraia de
punct fix, metoda Newton
4.7.3 Estimarea de tip Milne a erorii. Modificarea pasului
Diferena ntre valoarea corectat )( 1c
ix i valoarea prezis )0(
1ix , constituie o
estimare a erorii pe pasul i+1, numit estimare de tip Milne:
)0(
1
)(
11 ic
ii xx
Dac n raport cu o toleran tol impus, avem:
- toli || 1 : )(
1
c
ix este acceptat (ca aproximaie pentru 1ix ) i calculul
continu. Dac toli || 1 , pasul h poate fi mrit pentru calculul valorii
urmtoare 2ix .
- toli || 1 : )(
1
c
ix nu este acceptat i se recalculeaz cu un pas h mai mic;
Utilizarea estimrii de mai sus este ns supus condiiei ca predictorul i
corectorul s aib acelai ordin.
a.ch. Noiembrie 2008
49
4.7.4 Eroarea de trunchiere a metodei predictor-corector
Fie predictorul i corectorul definii de (44), explicit i respectiv implicit:
1
0
1
0
)0(
1
k
l
k
l
lillili fhxx - predictor (79a)
1
0
)0(
111
1
0
1 ),(k
l
ii
k
l
lillili xtfhbfbhxax - corector (79b)
i presupunem c facem o singur iteraie n corector adic, aplicm (79b) cu
valoarea )0( 1ix furnizat de (79a), obinnd )1(
11 ii xx . Eroarea de trunchiere a
metodei (79a, b) se obine cum urmeaz.
Eroarea de trunchiere a metodei predictor-corector este de ordinul erorii de
trunchiere a corectorului. Pentru aceasta trebuie s avem:
1 CP pp
unde Pp i Cp sunt ordinele predictorului, respectiv corectorului.
4.7.5 Alegerea predictorului i corectorului. Exemple
Alegerea predictorului este legat de o ct mai bun estimare a lui )0( 1ix . ntre
predictori de acelai ordin, criteriul este coeficientul pd al erorii (ct mai mic).
Alegerea predictorului este mai puin critic dect cea a corectorului. Alegerea
corectorului este dictat n principal de proprietile lui de stabilitate. ntre
corectori de acelai ordin, criteriile de alegere sunt: coeficientul erorii (ct mai
mic) i regiunea de stabilitate absolut (ct mai mare). Dac predictorul este
suficient de exact v. 4.7.4 , ordinul metodei predictor-corector este ordinul
corectorului (dac ar fi utilizat singur).
Vom da cteva exemple de metode predictor-corector, dintre cele mai utilizate. n
predictor, notm:
),( )( 11)(
1
j
ii
j
i xtff
Primul exemplu este cel al unei metode de cel mai mic ordin (corector 1-pas,
2p ).
0. Metod de ordin 2:
a.ch. Noiembrie 2008
50
Predictor (Euler, ordin 1):
iii hfxx )0(
1
Corector (Metoda trapezului, ordin 2):
)(2
)(
1
)1(
1
j
iii
j
i ffh
xx
1. Metodele Milne i Hamming
Predictor (ordin 4):
)2(3
4213
)0(
1 iiiii fffh
xx
Corector (Milne, ordin 4):
)4(3
1
)(
11
)1(
1
iij
ii
j
i fffh
xx , 0j
Coeficienii erorii sunt: 14/45 (predictor) i -1/90 (corector).
Corectorul nu este ns relativ stabil, pentru valori > 0 (ecuaia de test este
xx v. 4.4).
O modificare a corectorului conduce la metoda Hamming care este stabil pentru
h 0.69. Corectorul devine:
)(~ )0(121112)0(
1
)0(
1 iiii xxxx - modificare
)2~
(8
3)9(
8
11
)(
12
)1(
1
iij
iii
j
i fffh
xxx - corector
unde f~
este calculat n valoarea modificat x~ .
2. Metode Adams ordin 4:
Predictor (Adams-Bashforth, ordin 4):
)9375955(25
321
)0(
1 iiiiii ffffh
xx
Corector (Adams-Moulton, ordin 4):
)5199(24
21
)(
1
)1(
1
iiij
ii
j
i ffffh
xx
a.ch. Noiembrie 2008
51
3. Metode Adams ordin 5:
Predictor (Adams-Bashforth, ordin 5):
)2511274261627741901(720
4321
)0(
1 iiiiiii fffffh
xx
Corector (Adams-Moulton, ordin 5):
)19106264646251(720
321
)(
1
)1(
1
iiiij
ii
j
i fffffh
xx
4.7.6 Implementarea metodelor predictor-corector
Implementarea clasic urmeaz, n principal, dou scheme:
a) PECE: Predicie Evaluare Corectare (1 iteraie) Evaluare
Aceasta nseamn, pe un pas 1i :
1. )0( 1ix este estimat de predictor; (P)
2. Se calculeaz ),( )0( 11 ii xtf ; (E)
3. Se face o singur iteratie n corector, rezult )1( 1ix ; (C)
4. Se calculeaz ),( 111 iii xtff (dac pasul este acceptat). (E)
Schema conduce la 2 evaluri de funcii pe un pas sub-paii 2, 4. Dac
predictorul este de ordin (ordinul corectorului 1), eroarea de trunchiere a
schemei este de ordinul erorii corectorului.
b) P(EC)M:
1) )0( 1ix este estimat de predictor; (P)
2.0) Testare tolxx jij
i
||)(
1
)1(
1 . Dac este satisfcut, se trece la pasul i+2.
Dac nu, se face:
2.1) Evaluare: ),( )( 11)(
1
j
ii
j
i xtff ; (E)
2.2) Iterare n corector: rezult )1( 1
j
ix , i se revine la pasul 2.0. (C)
M este numrul de iteraii pe pasul 1i , i poate varia de la un pas la altul.
Schema (b) conduce la M evaluri de funcii pe un pas, dar numrul total de
a.ch. Noiembrie 2008
52
evaluri de funcii (pentru ntregul interval de integrare) poate fi mai mic dect n
schema (a), i n acest caz, schema (b) este mai eficient.
Observaie
O alt schem (c), const n a impune un numr fixat de iteraii M n schema (b).
n acest caz, la pasul 2.0 se testeaz dac numrul de iteraii = M. n cazul
satisfacerii testului 2.0, pasul este urmat de evaluarea 2.1. n acest caz, schema se
simbolizeaz prin P(EC)ME. Cazul M = 1 revine la schema (a)
Stabilitate:
n schema (a) PECE, utilizarea unei singure iteraii modific caracteristicile de
stabilitate ale metodei n raport cu cele ale schemei (b), i stabilitatea este
influenat de predictorul utilizat. Concret, schema (a) micoreaz regiunea de
stabilitate pe care o are corectorul utilizat singur. n schema (b) P(EC)M, iterarea
pn la convergen nu modific stabilitatea metodei i aceasta nu este influenat
de predictorul utilizat. V. o discuie mai ampl n Ralston & Rabinowitz (1978), i
reprezentri ale regiunii de stabilitate pentru schema (a) n Hairer et al. (1991).
Observaie
Codurile care implementeaz aceste scheme utilizeaz procedeul pas variabil
ordin variabil (abreviat VSVO), adic n funcie de un test pe un pas 1i acceptat,
la pasul urmtor se poate varia att ordinul metodei ct i mrimea pasului
4.7.7 Determinarea valorilor de start
La primul pas ( 0i ), metoda n k-pai cere k valori de start 110 ,,, kxxx .
)0(
0 xx este dat de condiiile iniiale, dar pentru celelalte valori trebuie utilizat
o procedur de determinare. Aceasta poate fi:
a) Seria Taylor
b) Metode Runge-Kutta
c) Metode multi-pas de ordin mai mic
(a) Se utilizeaz desvoltarea lui )(tx n serie Taylor, n jurul lui 0t :
0
2
000!2
)()( x
jhxjhxjhtx 1,,1 kj
a.ch. Noiembrie 2008
53
Derivata 0x se calculeaz, cu condiiile iniiale, din ecuaia dat ),( xtfx ,
iar derivatele )( 0)()(
0 txxnn prin derivarea ecuaiei. Desvoltatea n serie se
face pn la termenul de ordinul p inclusiv, unde p este ordinul metodei.
(b) Metodele Runge-Kutta sunt auto-start i pot astfel furniza valorile de start. Se
va utiliza o metod al crei ordin este egal cu ordinul metodei multi-pas.
(c) Procedurile a, b, au fost utilizate nainte de apariia implementrii metodelor
VSVO (pas variabil - ordin variabil). n acestea din urm, se ncepe integrarea
cu o metod 1-pas (v. Exemplul 0 n 4.7.5) care calculeaz 1x , i se continu
cu metode n 2, 3, , (k-1) pai. Cu cele k valori calculate, se continu cu
metoda n k-pai. Astfel, aceste metode devin metode auto-start.
4.8 Comparaia metodelor predictor-corector (PC) cu metodele Runge Kutta
(RK)
Comparaia se face lund n considerare urmtoarele criterii:
a) Necesitatea valorilor de start;
b) Precizie;
c) Numr de evaluri de funcii / pas;
d) Numrul de evaluri suplimentare de funcii, necesare pentru controlul erorii
locale de trunchiere.
(a) Metodele R-K sunt auto-start, pe cnd metodele PC cer o procedur pentru
valorile de start. Totui, n implementarea predictor-corector, ordin variabil-
pas variabil, metodele PC devin auto-start.
(b) La ordine egale, precizia metodelor RK este uor superioar. Compar
rezultatele din 4.9 (v. mai jos) cu cele din 3.3.9.
(c) Pentru metodele RK, numrul de evaluri/pas este ordinul metodei. Pentru
metodele PC, n implementarea PECE acest numr este 2, dar n
implementarea P(EC)M
, respectiv P(EC)ME, acesta depinde de numrul M de
iteraii (fiind egal cu M, respectiv M+1).
(d) Controlul erorii locale de trunchiere: cere evaluri suplimentare n metodele
RK, n timp ce la metodele PC nu cere evaluri suplimentare.
a.ch. Noiembrie 2008
54
Codurile actuale se orienteaz mai mult spre metodele RK, sau metode de tip
similar pentru ecuaii difereniale de ordinul doi (metodele Nystrm) cu
excepia ecuaiilor difereniale rigide, pentru care se utilizeaz metode BDF i
metode RK implicite.
4.9 Exemple numerice
Vom relua problema celor dou corpuri din 3.3.9, pentru e = 0.9, calculnd soluia
pe [0, 20] n dubl precizie, cu urmtoarele metode predictor-corector:
- Adams de ordinele 4 i 5 (4.7 Exemplele 2 i 3), cu pas constant i cu 1
iteraie n corector utiliznd codul din Anexa, 4.2.
- Adams, ordin variabil-pas variabil, ordinul 12, lucrnd cu subrutina
DIVPAG din IMSL. S-au considerat dou cazuri: ordinul 5, i 12.
ntervalul de integrare s-a mprit n 20, respectiv 200, sub-intervale.
Rezultate mai precise se obin pentru 20 sub-intervale (pasul maxim este astfel
1.).
Rezultatele se dau n tabelele urmtoare.
e = 0.9, Metoda Adams ordin 4, pas constant: Erori absolute extreme la t = 18.84
Pasul
h
Eroarea absolut
Maxim ( x ) Minim ( x )
Nr. apeluri
DERIVS
0.01 3.65 D0 3.30 D-1
4010
0.001 2.70 D-2 1.14 D-5 40010
0.0005 2.09 D-3 1.14 D-6 80010
Eroarea maxim are loc n y i cea minim n y.
e = 0.9, Metoda Adams ordin 5, pas constant: Erori absolute extreme la t = 18.84
Pasul
h
Eroarea absolut
Maxim ( x ) Minim ( x )
Nr. apeluri
DERIVS
0.01 4.29 D0 2.97 D-1
4013
0.001 6.64 D-4 3.69 D-7 40013
0.0005 3.33 D-5 1.86 D-8 80013
Eroarea maxim are loc n y i cea minim n x .
a.ch. Noiembrie 2008
55
Observaie
Exemplele din tabelele anterioare s-au rulat i iternd n corector pn la
convergen, cu tolerana tol = 1D-10 (pentru iteraia de punct fix). Ordinul erorii
a fost aproximativ acelai, iar numrul de iteraii a crescut nesemnificativ: de
exemplu, pentru Adams ordinul 5, h = 0.001, valorile din tabelul de mai sus sunt:
1.71 D-4, 9.85 D-8, 40942
e = 0.9, Metoda Adams, ordin maxim 5, ordin variabil-pas variabil (DIVPAG):
Erori absolute extreme la t = 18.0
Tolerana
TOL
Eroarea absolut
Maxim ( x ) Minim ( y )
Numr
de pai
Numr
apeluri FCN
1 D-10 1.11 D-6 1.05 D-8 2437 2577
1 D-15 6.84 D-11 4.44 D-13 16276 16424
e = 0.9, Metoda Adams, ordin maxim 12, ordin variabil-pas variabil (DIVPAG):
Erori absolute extreme la t = 18.0
Tolerana
TOL
Eroarea absolut
Maxim ( x ) Minim ( y )
Numr
de pai
Numr
apeluri FCN
1 D-10 2.07 D-07 1.44 D-08 1611 1760
1 D-15 6.28 D-12 3.25 D-13 4228 4417
Observaii
- Pasul maxim care a putut fi utilizat de DIVPAG a fost 1, adic lungimea sub-
intervalului (nu s-au impus nici pasul maxim, nici pasul minim, care pot fi
specificai n parametrii de intrare hmax, hmin). Aceasta explic numrul
redus de pai cu care lucreaz rutina chiar pentru TOL foarte mic.
- TOL este un parametru care servete la controlul normei erorii locale, astfel ca
eroarea global s fie proporional cu TOL. Norma aleas n exemplele de
mai sus a fost norma-. TOL i tipul de norm sunt parametri de intrare ai
rutinei DIVPAG (n tabloul param). A nu se confunda parametrul TOL cu
tolerana tol pentru iteraia de punct fix n 4.7.6.
a.ch. Noiembrie 2008
56
- Se observ superioritatea implementrii n forma ordin variabil-pas variabil,
cu controlul erorii, fa de implementarea cu pas fix.
- Comparnd rezultatele cu cele obinute prin metoda Runge-Kutta 8(7) v.
3.3.9, se constat o precizie mai mare a acesteia din urm, cu preul unui
numr mai mare de evaluri de funcii: de exemplu, pentru tolerana 2.22D-15,
ordinul erorii este n plaja D-13 D-14, la un numr de 1380 pai i cca.
18000 evaluri
5 Ecuaii difereniale rigide
Vom face n continuare o scurt introducere n problematica ecuaiilor difereniale
rigide. Pentru o tratare extensiv trimitem la tratatul Hairer et al. (1991), dedicat
integrrii ecuaiilor difereniale rigide.
n aplicarea unei metode numerice, dimensiunea pasului se alege dintr-o condiie
de precizie a soluiei, impunnd o toleran pentru eroarea de trunchiere local. De
obicei, pasul rezultat este n regiunea de stabilitate a metodei. Dac ns
dimensiunea pasului este dictat mai degrab de condiia de stabilitate dect de
condiia de precizie, zicem c avem o problem rigid.
Exemplu 1
S considerm urmtoarea problem:
0)0(,1)0(;0100101 xxxxx
Punnd ecuaia sub forma unui sistem de ordinul nti, avem
uxuux 101100, ; 0)0(,1)0( ux ,
sau matriceal, Axx , unde Tux ][x , iar
101100
10A .
Valorile proprii ale lui A sunt 100,1 21 . Soluia exact este de forma
tt eCeCtx 10021)( (84)
i cu condiiile iniiale date rezult
tt eetx 100991
99100)( . (85)
Termenul n te 100 tinde rapid spre 0 (se amortizeaz): de exemplu, pentru t = 0.1
avem 510 1054.4 e , astfel c pentru t > 0.1 avem:
a.ch. Noiembrie 2008
57
tetx 99
100)( (86)
Soluia dat de al doilea termen din (85) zice tranzitorie. Soluia (86) reprezint
soluia staionar. Pentru t suficient de mare (de exemplu, t > 0.1), soluia
tranzitorie nu mai contribuie la soluia (85), dar determin stabilitatea metodei.
Chiar i n cazul n care termenul al doilea dispare din soluia (84) condiiile
iniiale sunt astfel c rezult 02 C valoarea proprie 2 determin intervalul de
stabilitate al metodei. ntr-adevr, s presupunem c integrm sistemul cu metoda
RK4. Intervalul de stabilitate absolut a metodei este (v. 3.3.9): -2.785296 < h <
0, sau 0 < h(- < 2.785296, ceea ce conduce (pentru 2 ) la h < 0.02785. ntruct
ordinul erorii globale a metodei RK4 este 4h , dac am vrea s calculm soluia cu
o eroare de ordinul 10-4
, ar fi suficient un pas de ordinul h = 0.1 care este ns n
afara intervalului de stabilitate. ntr-adevr, calculnd cu h = 0.025 (pas constant)
se obine x(10) = 4.5858514950 E-5 care are o eroare de 6.21E-13, n timp ce cu
h = 0.28 i 0.3, rezult respectiv x(9.996) = -2.74E+1, i x(9.990) = -
1.1459E+44 care probeaz instabilitatea
Fie un sistem liniar Axx , de m ecuaii, i i valorile proprii ale matricii A.
Faptul c sistemul este rigid se definete prin urmtoarele condiii:
mii ,1,0)Re(
|)Re(|min|)Re(|max,1,1
imi
imi
Cu alte cuvinte, sistemul este rigid dac are un punct fix stabil, i A are valori
proprii de mrimi foarte diferite. Definiia de mai sus are mai multe
inconveniente, i anume:
- Nu mai convine n cazul n care matricea A are o valoare proprie egal cu
zero;
- Nu se poate aplica pentru un sistem neliniar, sau pentru o singur ecuaie.
Se d atunci urmtoarea definiie mai puin precis, dar aplicabil att sistemelor
liniare ct i celor neliniare, ct i pentru o singur ecuaie (Lambert, v. Cartwright
and Piro (1992)):
a.ch. Noiembrie 2008
58
Dac o metod numeric este forat s utilizeze, ntr-un anumit interval, un pas
excesiv de mic pentru o problem a crei soluie exact este neted n acel interval,
atunci problema se zice rigid n acel interval
(O funcie este neted n [a, b], dac are derivat continu n [a, b].)
O problem poate fi rigid pe unele sub-intervale ale soluiei i non-rigid pe
altele. Definiia de mai sus permite codurilor pentru integrarea numeric, s
recunoasc rigiditatea problemei pe intervalul pe care se calculeaz soluia (sau
pe sub-intervale ale acestuia), prin faptul c rutina este forat s micoreze
excesiv pasul, pentru a satisface tolerana impus erorii de trunchiere. Utilizarea
unui pas foarte mic poate crea probleme datorate acumulrii erorilor de rotunjire
sau creterii timpului de calcul.
Pentru o problem rigid controlat de un parametru, s-ar cere ca metoda de
integrare s fie stabil pentru orice dimensiune a pasului, la orice valoare a
parametrului pentru care problema este stabil. De exemplu, problema de test
pentru stabilitatea absolut, xx , este stabil pentru 0)Re( , iar metoda ar
trebui s fie stabil pentru orice h , oricare ar fi cu 0)Re( . Regiunea de
stabilitate absolut este atunci tot