.

1

I.2.10. Controlul pasului

În unele cazuri este necesară pentru obţinerea unei erori date folosirea unui pas variabil în rezolvarea numerică. Metodele numerice care folosesc un pas variabil se numesc metode adaptive. Pentru controlul pasului vom folosi metodele imbricate descrise anterior. Fie o metodă numerică unipas de ordin ( )nhO de forma:

( ) 1...,,2,1,,,1 −=+=+ �ihythyy iiii φ (I.93)

ataşată unei probleme cu valori iniţiale:

=

≤≤=

ayay

btaytfty

)(~),~,()('~

.

2

Se pune problema găsirii unui număr minim de noduri � şi a unei grile

variabile { } 1,1 −= �iih astfel încât pentru o toleranţă dată 0>ε eroarea globală

( ) ii yty −~ să nu depăşească această eroare.

Pentru estimarea erorii considerăm o nouă metodă similară cu prima (de

exemplu ambele metode sunt de tip Runge-Kutta) dar de ordin de

exactitate ( )1+nhO :

( ) 1...,,2,1,,, ***1

* −=+=+ �ihythyy iiii φ (I.94)

O estimare a erorii este dată de (I.91)

( ) ( ) ( )[ ]hythyth

hytrp

,,,,1

,, **φφ −= (I.91)

.

3

iar legătura cu eroarea locală de trunchiere este dată de

( ) ( ) ( ) ( ) ( )11 ,,,,,, ++ +≈+= ppii

ppiiii hOhhytrhOhhythytT τ (I.95)

Dar (I.95) este echivalentă cu

( ) pii hKhytT ≈,, (I.96)

Atunci eroarea locală de trunchiere pentru un nou pas qh poate fi

estimată în modul următor:

( ) ( ) ( ) ( ) ( )

( ) ( )[ ]hythytq

hhytrqhytTqKhqqhKqhytT

iiiip

pii

pii

ppppii

,,,,

,,,,,,

**

(I.91)

(I.95)(I.96)

φφ −=

=≈≈=≈

(I.97)

Dacă impunem ca eroarea lui ( )qhytT ii ,, să fie mărginită de ε atunci din

.

4

(I.97) găsim:

( ) ( )p

ii

p

iiii yy

h

hythytq

1

*

1

** ,,,,

−=

−≤

ε

φφ

ε (I.98)

unde am folosit (I.93) şi (I.94).

De obicei în formula (I.98) se introduce un factor θ (de exemplu 8.0≈θ ),

numit factor de siguranţă, astfel încât:

p

ii yy

hq

1

*

−≤

εθ (I.99)

Pentru a calcula diferenţa ( ) ( )hythyt iiii ,,,, **φφ − se folosesc de obicei

metode imbricate pentru a reduce numărul de evaluări ale funcţiei f . Una

.

5

dintre cele mai cunoscute tehnici este metoda Runge-Kutta-Fehleberg

care combină o metodă de ordin patru cu o metodă de ordin cinci şi care

are tabela Butcher:

0 0

4

1 4

1 0

8

3 32

3 32

9 0

13

12 2197

1932 2197

7200−

2197

7296 0 1 216

439 - 8 513

3680 4104

845− 0

2

1 27

8− 2 2565

3544−

4104

1859 40

11− 0

216

25 0 2565

1408 4104

2197 5

1−

135

16 0 12825

6656 56430

28561 50

9−

55

2

iar pasul variabil este dat de:

.

6

4

1

*84.0

−=

ii yy

hq

ε (I.100)

În general asupra paşilor variabili obţinuţi cu (I.100) se mai pune o

condiţie de mărginire de forma maxmin hhh i ≤≤ pentru a evita consumul

prea mare de timp cu calculul în regiunile în care derivata lui f prezintă

neregularităţi şi pentru a evita paşii prea mari care nu ar putea analiza

variaţii locale mai mici .

O altă tehnică foarte cunoscută este metoda Bogacki-Shampine ce

combină o metodă de ordinul doi cu o metodă de ordin 3 şi stă la baza

.

7

funcţiei ode23 din Matlab. Ea are are tabela Butcher:

0 0

2

1 2

1 0

4

3 0 4

3 0 1 9

2 9

3 9

4 0 9

2 9

3 9

4 24

7 4

1 3

1 8

1

Vom da în continuare un exemplu în care vom folosi o metodă cu pas

variabil calculat pe baza metodei imbricate Euler-Heun cu tabela

.

8

Butcher:

0 0 0 1 1 1 1 0 1/2 1/2

Precizăm că metoda folosită în exemplul următor are un scop pur

didactic, în practcă fiind necesară o tehnică Bogacki-Shampine sau

Runge-Kutta-Fehleberg.

Exemplul 1: Rezolvaţi problema iniţială

.

9

1)0(,40),sin( =≤≤= yttydt

dyπ

folosind un pas variabil cuprins între 0.01 şi 0.1cu o eroare estimată 31 −< eε . Soluţia exactă este ( ) )cos(1 tety −= .

Aplicăm în continuare metoda imbricată Euler-Heun.

11

1 ),(),(

hKyy

ytfytK

ii

iiii

+=

=

+

şi ( )( )21

*1

*

1**

2

**1

2

1

,),,(

),(),(

KKhyy

hKyhtfhytK

ytfytK

ii

iiii

iiii

++=

++=

=

+

deci putem evalua în ecuaţia (I.99) ( )21* 5.0 KKhyy ii −=− .

.

10

În urma aplicării algoritmului s-a obţinut pasul minim, 0.01min =h ,

pasul maxim, 0.1max =h , iar numărul de noduri necesare integrării

este 567=� .

Programul Matlab este:

%exemplu pas variabil Euler-Heun

clear all

format short

a=0;b=2*pi;%capetele intervalului

hmax=0.1;hmin=0.01;

tol=1e-3;

y(1)=1; %conditia initiala

T(1)=a;

i=1;

t=a;

while t<b

.

11

K1=fPasVar(t,y(i));

K2=fPasVar(t+h,y(i));

R=0.5*h*abs(K1-K2);

if R<=tol

T(i+1)=T(i)+h;

y(i+1)=y(i)+0.5*h*(K1+K2);

i=i+1;

t=t+h;

h=hmax;

else

q=0.8*(tol*h/R)^0.25;

h=q*h;

if h<hmin

h=hmin;

end

if h>hmax

h=hmax;

end

end

if(t+h>b)

.

12

h=b-t;

end

end %while

plot(T,y,'.k') %reprezentam grafic solutia numerica

hold on

plot(T,exp(1-cos(T)),'r');

Am rezolvat problema şi pentru pasul cuprins între 0.001 şi 0.1 cu

o eroare estimată 41 −< eε soluţia numerică fiind mult

îmbunătăţită. S-a obţinut 0.0019min =h , 0.0437max =h şi 1736=� .

.

13

Figura 11. Variaţia pasului pentru 1.001.0 ≤≤ h şi 31 −< eε .

.

14

Figura 12. Soluţia numerică şi soluţia analitică pentru

1.001.0 ≤≤ h şi 31 −< eε .

.

15

Figura 13. Soluţia numerică şi soluţia analitică pentru 1.0001.0 ≤≤ h şi 41 −< eε .

.

16

I.2.10. Probleme stiff

.

17

Fie o metodă cu un pas ataşată problemei (I.101) de forma:

(I.101)

(I.101)

(I.101)

(I.101)

(I.103)

.

18

(I.104)

:

(I.105)

(I.106)

.

19

(I.101)

(I.107)

(I.108)

(I.103)

(I.108)

(I.109) (I.102)

.

20

Aproximare Padé

(I.104)

(I.110)

(I.106) (I.110)

:

.

21

(I.111)

(I.112)

(I.112)

:

.

22

:

(I.113)

(I.114)

(I.115)

10: (I.104)

9:

.

23

Metoda implicită a lui Euler

Metoda implicită a lui Euler a fost tratată şi în cursul 1.

(I.116)

(I.117)

(I.101)

(vezi (I.106)) 9

.

24

Metoda implicită a trapezului

(I.118)

(I.119)

(I.101)

.

25

Exemplul 1: Rezolvaţi problema iniţială 1)0(,4/30,2 =≤≤= ytydt

dy .

Soluţia exactă este ( ) )1/(1 tty −= .

Programele Matlab sunt:

%metoda implicita a trpezului

a=0;b=0.75; %capetele intervalului

h=0.05; %pasul retelei

N=round((b-a)/h)+1; %numarul de noduri

ye=zeros(N,1);%initializam vectorul solutie pentru metoda

lui Euler

yt=zeros(N,1);%initializam vectorul solutie pentru

ye(1)=1; yt(1)=1; yex(1)=1; %conditia initiala

t=a:h:b;

.

26

for i=2:N

ye(i)=ye(i-1)+h*ye(i-1)^2; %metoda Euler explicita

yt(i)=fNewt(yt(i-1),yt(i-1),h,1e-5);%metoda trapezului

yex(i)=-1/(t(i)-1); %solutia exacta

end plot(t,y,'.k') %reprezentam grafic solutia numerica(Euler)

hold on

plot(t,yt,'.r')%reprezentam grafic solutia numerica(trapez)

plot(t,yex,'b') %reprezentam grafic solutia exacta

function y1=fNewt(y0,yi,h,tol)

err=1;

y1=y0;

while(err>tol)

y1=y0-ftrap(y0,yi,h)/fptrap(y0,h);

err=abs(y1-y0);

y0=y1;

end

function yp=ftrap(y,yi,h)

yp=y-yi-0.5*h*(yi^2+y^2);

function yp=fptrap(y,h)

yp=1-y*h;

.

27

Figura 14. Soluţii numerice şi soluţia analitică

.

28

Formule Runge-Kutta implicite

O formulă Runge-Kutta implicită cu r stadii are forma:

( ) ( )∑=

=r

siissii hfytKchfyt

1,,,,,,φ (I.120)

rsKahyhtfytKr

jjsjisiiis ,...,2,1,,),(

1=

++= ∑

=µ

Se poate arăta că metoda este de ordinul p, rpr 2≤≤ dacă

[ ]( )dp baCf R×∈ , şi

.

29

Formulele (I.121) şi (I.122) pot fi privite ca nişte formule de cuadratură

(I.121)

(I.122)

(I.121)

(I.121’)

(I.122’)

.

30

Se poate arăta că o formulă Runge-Kutta cu ordinul maxim r2 este A –

stabilă (vezi Trambiţaş şi ceilalţi, 2010).

(I.122’)