Gabriela Ciuprina

Algoritmi numerici
- teme laborator -

2009
Cuprins

5 Derivarea numerică. Metoda diferenţelor finite. 1

5.1 Formule de derivare ı̂n cazul unidimensional . . . . . . . . . . . . . . . . . 1

5.2 Metoda diferenţelor finite . . . . . . . . . . . . . . . . . . . . . . . . . . . 5

5.2.1 Studiul regimului tranzitoriu al unui circuit de ordinul 1 . . . . . . 5

Formularea corectă a problemei . . . . . . . . . . . . . . . . . . . . 6

Rezolvarea cu diferenţe finite . . . . . . . . . . . . . . . . . . . . . 8

5.2.2 Studiul regimului electrocinetic staţionar al unui conductor 2D . . . 14

Formularea corectă a problemei . . . . . . . . . . . . . . . . . . . . 15

Rezolvarea cu diferenţe finite . . . . . . . . . . . . . . . . . . . . . 16

5.2.3 Studiul propagării undei scalare . . . . . . . . . . . . . . . . . . . . 26

i
Capitolul 5

Derivarea numerică. Metoda

diferenţelor finite.

Necesitatea evaluării derivatelor funcţiilor este importantă ı̂n problemele de inginerie.

De exemplu, ele apar ı̂n evaluarea senzitivitătilor unor mărimi de interes ı̂n raport cu
anumiţi parametri consideraţi variabili, senzitivităti utile ı̂n proiectarea optimală a unor
dispozitive. Un alt exemplu este cel ı̂n care modelul matematic al problemei conduce la
ecuaţii sau sisteme de ecuaţii diferenţiale, foarte puţine dintre acestea putând fi rezolvate
prin metode analitice.

Derivarea numerică se bazează pe interpolare. Formulele de derivare numerică provin

din derivarea polinomului de interpolare şi ele se reduc la calcule aritmetice ce folosesc
formule relativ simple. Prima parte a acestui capitol ilustrează deducerea acestor formule
şi erorile care apar. Rezolvarea ecuaţiilor diferenţiale folosind formule de derivare numerică
pentru derivatele ce intervin este cunoscută sub numele de metoda diferenţelor finite.
Partea a doua a acestui capitol urmăreşte implementarea şi testarea acestei metode pentru
câteva cazuri simple de test.

5.1 Formule de derivare ı̂n cazul unidimensional

Formulele de derivare numerică se pot deduce fie pornind de la polinomul de inter-

polare, fie prin intermediul seriilor Taylor. Această din urmă abordare are avantajul că
furnizează şi informaţii legate de eroarea de trunchiere.

Vom reaminti câteva rezultate pentru cazul aproximării derivatelor unei funcţii unidi-
mensionale f : [a, b] → IR.

1
2 Capitolul 5. Derivarea numerică. Metoda diferenţelor finite.

Considerând xi punctele diviziunii intervalului [a, b], vom folosi următoarele notaţii

∂f 2 ∂2f
fi = f (xi ), Dfi = (xi ), D fi = (xi ). (5.1)
∂x ∂x2

Cazul unei diviziuni uniforme: xi+1 − xi = h

• Aproximarea derivatelor pornind de la polinomul de interpolare de gradul 1:

fi+1 − fi
(progresivă) Dfi = + O(h) (5.2)
h
fi − fi−1
(regresivă) Dfi = + O(h) (5.3)
h

• Aproximarea derivatelor pornind de la polinomul de interpolare de gradul 2:

fi+1 − fi−1
(centrată) Dfi = + O(h2 ) (5.4)
2h
−3fi + 4fi+1 − fi+2
(progresivă) Dfi = + O(h2 ) (5.5)
2h
fi−2 − 4fi−1 + 3fi
(regresivă) Dfi = + O(h2 ) (5.6)
2h
fi+1 − 2fi + fi−1
(centrată) D2 fi = + O(h) (5.7)
h2

Cazul unei diviziuni neuniforme: xi+1 − xi = αh, xi − xi−1 = βh

• Aproximarea derivatelor de ordinul 1 pornind de la polinomul de interpolare de

gradul 2:

β 2 fi+1 + (α2 − β 2 )fi − α2 fi−1

(centrată) Dfi = + O(h2 ) (5.8)
αβ(α + β)h
−α(α + 2β)fi + (α + β)2 fi+1 − β 2 fi+2
(progresivă) Dfi = + O(h2 ) (5.9)
αβ(α + β)h
α2 fi−2 − (α + β)2 fi−1 + β(2α + β)fi
(regresivă) Dfi = + O(h2 ) (5.10)
αβ(α + β)h

• Aproximarea derivatei de ordinul doi cu o eroare de ordinul1 h:

βfi+1 − (α + β)fi − αfi−1

D2 fi = 2 + O(h) (5.11)
αβ(α + β) h2
1
Derivata de ordinul doi nu poate fi aproximată cu o eroare de ordinul h2 dacă sunt folosite numai
punctele xi−1 , xi şi xi+1 .

 [Link], Draft din 31 decembrie 2009

5.1. Formule de derivare ı̂n cazul unidimensional 3

Exerciţiul 5.1:
Pornind de la polinomul de interpolare de grad doi, ce trece prin punctele (xi−1 , fi−1 ),
(xi , fi ), (xi+1 , fi+1 ) deduceţi una din formulele de derivare de mai sus

Exerciţiul 5.2:
a) Scrieţi ı̂ntr-un fişier următorul script:

clear all;
% calculeaza derivata numerica a functiei sinus in pct x0
x0 = pi/4;
h = pi/10; % pasul de derivare

dp = (sin(x0+h) - sin(x0))/h; % derivata progresiva de ordinul 1

dr = (sin(x0) - sin(x0-h))/h; % derivata regresiva de ordinul 1
dc = (sin(x0+h) - sin(x0-h))/(2*h); % derivata centrata de ordinul 2
dexact = cos(x0);

erp = abs(dp - dexact);

err = abs(dr - dexact);
erc = abs(dc - dexact);
disp(sprintf(’ ==== Der fct sin in x0 = %e, h = %e’,x0,h));
disp(sprintf(’valoarea exacta a der = %e’,dexact));
disp(sprintf(’derivata nr progresiva = %e, eroare = %e’,dp,erp));
disp(sprintf(’ regresiva = %e, eroare = %e’,dr,err));
disp(sprintf(’ centrata = %e, eroare = %e’,dc,erc));

b) Comentaţi valorile obţinute pentru cele trei formule. Schimbaţi pasul de derivare şi
comentaţi rezultatele.

Este foarte probabil ca ı̂n acest moment să fiţi tentaţi să credeţi că, cu cât pasul de
derivare este mai mic, cu atât eroarea e mai mică. Pentru a vedea dacă această afirmaţie
este adevărată, este util să reprezentaţi grafic eroarea ı̂n funcţie de pasul de derivare.

Exerciţiul 5.3:
Scrieţi un script care să permită:
a) Reprezentarea grafică a erorii formulei de derivare progresivă de ordinul 1, ı̂n funcţie de
pasul de derivare. Pentru un pas h cuprins ı̂ntre 1e-10 şi 1e-5 trebuie să obţineţi un grafic
similar celui din partea stângă a fig. 5.1. Acest grafic a fost trasat folosind 1e5 valori.
Repetaţi experienţa pentru derivarea regresivă de ordinul 1. Comentaţi rezultatul.
c) Reprezentarea grafică a erorii formulei de derivare centrată de ordinul 2, ı̂n funcţie de

Document disponibil la [Link] gabriela

4 Capitolul 5. Derivarea numerică. Metoda diferenţelor finite.

pasul de derivare. Pentru un pas h cuprins ı̂ntre 1e-8 şi 1e-3 trebuie să obţineţi un grafic
similar celui din partea dreaptă a fig. 5.1. Acest grafic a fost trasat folosind 1e5 valori.

−5 −6
10 10
Eroarea absoluta a derivatei progresive de ord 1

−7

Eroarea absoluta a derivatei centrate de ord 2

10
−6
10
−8
10

−9
−7 10
10
−10
10
−8
10 −11
10

−12
10
−9
10
−13
10

−10 −14
10 10
−10 −9 −8 −7 −6 −5 −8 −7 −6 −5 −4 −3
10 10 10 10 10 10 10 10 10 10 10 10
Pasul de derivare h Pasul de derivare h

Figura 5.1: Eroarea ı̂n funcţie de pasul de derivare pentru formula de derivare progresivă
de ordinul 1 (stânga) şi formula de derivare centrată de ordinul 2 (dreapta).

Zona de instabilităţi numerice, care apare pentru paşi mai mici decât o anumită valoare,
se datorează erorilor de rotunjire care devin mai mari decât erorile de trunchiere. Pasul
de la care ı̂ncep să predomine erorile de trunchiere se numeşte pas optim de derivare
numerică. Pasul optim nu este pasul pentru care eroarea este minimă. După cum se
poate observa, erori mai mici se pot obţine pentru paşi aflaţi ı̂n zona de instabilitate.
Totuşi lucrul ı̂n această zonă nu se recomandă. Este de preferat ca h să fie ı̂n zona ı̂n care
predomină erorile de trunchiere, zonă ı̂n care erorile de trunchiere se pot estima.

Exerciţiul 5.4:
Inspectând graficele obţinute
a) comparaţi paşii optimi de derivare numerică pentru cele trei formule pe care le-aţi
folosit până acum.
b) verificaţi că eroarea de trunchiere la formulele progresivă şi regresivă de ordinul 1 este
O(h) iar la formula centrată de ordinul 2 este O(h2 ).

Derivata centrată (5.4) oferă o acurateţe mai bună decât formulele progresivă, regresivă
de ordinul 1 (formulele (5.2) şi (5.3)). O comparaţie numerică este sugerată ı̂n exerciţiul
următor.

Exerciţiul 5.5:
Reprezentaţi pe acelaşi grafic eroarea derivatei numerice progresive de ordinul 1 şi centrate
de ordinul 2 pentru paşi cuprinşi ı̂ntre 1e-8 şi 1e-3. Trebuie să obţineţi un grafic similar
celui din figura 5.2.

 [Link], Draft din 31 decembrie 2009

5.2. Metoda diferenţelor finite 5

−2
10

−4 Progresiva de ord. 1
10

Eroarea absoluta a derivatei numerice

Centrata de ord. 2

−6
10

−8
10

−10
10

−12
10

−14
10
−8 −7 −6 −5 −4 −3
10 10 10 10 10 10
Pasul de derivare h

Figura 5.2: Acest grafic trebuie obţinut la exerciţiul 5.5.

5.2 Metoda diferenţelor finite

Aplicarea metodei diferenţelor finite pentru rezolvarea oricărei ecuaţii diferenţiale pre-
supune realizarea următorilor paşi:

• Pasul 1: Se alege o schemă de diferenţe finite pentru aproximarea derivatelor din

ecuaţii şi se rescrie ecuaţia ca ecuaţie cu diferenţe finite.

• Pasul 2: Se stabileşte grila de discretizare şi se scrie ecuaţia discretizată pentru

fiecare nod.

• Pasul 3: Se rezolvă sistemul de ecuaţii pentru determinarea valorilor necunoscute

ı̂n nodurile grilei.

• Pasul 4: Calculul altor valori, ı̂n afara celor din noduri, se face prin interpolare.

Se obişnuieşte să se spună că paşii 1 şi 2 reprezintă etapa de preprocesare, pasul 3
reprezintă rezolvarea, iar pasul 4 postprocesarea.

5.2.1 Studiul regimului tranzitoriu al unui circuit de ordinul 1

Cel mai simplu exemplu pe care ni-l putem imagina este cel care ne conduce matematic
la o ecuaţie diferenţială ordinară de ordinul 1. Vom considera un condensator de capacitate
C = 4µF, iniţial descărcat, care se ı̂ncarcă de la o sursă de tensiune continuă E = 20 mV
printr-un rezistor de rezistenţă R = 10 Ω.

Document disponibil la [Link] gabriela

6 Capitolul 5. Derivarea numerică. Metoda diferenţelor finite.

Formularea corectă a problemei

Dacă notăm cu i curentul prin circuit şi cu u tensiunea la bornele condensatorului,

atunci teorema Kirchhoff II scrisă pe bucla circuitului este

Ri(t) + u(t) = E. (5.12)

Înlocuind ı̂n această relaţie dependenţa dintre tensiunea şi curentul prin condensator
du(t)
i(t) = C , (5.13)
dt
rezultă ecuaţia diferenţială de ordinul unu, satisfăcută de u:
du(t)
RC + u(t) = E. (5.14)
dt
Soluţia acestei ecuaţii este unică deoarece se specifică starea iniţială a condensatorului:

u(0) = u0 = 0. (5.15)

Avantajul acestui exemplu este acela că are o soluţie analitică, uşor de găsit. Vom
reaminti modul de calcul al acesteia.

Conform teoriei cunoscute de la matematică, soluţia acestei ecuaţii se scrie

u(t) = uo (t) + up , (5.16)

unde uo (t) este soluţia ecuaţiei ”omogene” (cum zic matematicienii) adică a ecuaţiei
diferenţiale cu termen liber nul:
duo (t)
RC + uo (t) = 0. (5.17)
dt
Soluţia ecuaţiei omogene este
uo (t) = A exp(xt), (5.18)
unde x este soluţia ecuaţiei caracteristice

RCx + 1 = 0. (5.19)

Dacă notăm produsul dintre rezistenţă şi capacitate cu

τ = RC, (5.20)

mărime care se numeşte constanta de timp a circuitului, rezultă soluţia ecuaţiei omogene

uo (t) = A exp(−t/τ ). (5.21)

 [Link], Draft din 31 decembrie 2009

5.2. Metoda diferenţelor finite 7

Mărimea up este o soluţie particulară, de forma termenului liber al ecuaţiei, deci o con-
stantă. Să notăm această constantă cu B. Soluţia particulară (constantă) trebuie să
satisfacă ecuaţia (5.14), deci
dB
RC + B = E, (5.22)
dt
de unde rezultă
B=E (5.23)
deoarece derivata unei constante este nulă. Soluţia generală a ecuaţiei este ı̂n consecinţă

u(t) = A exp(−t/τ ) + E. (5.24)

Constanta A se determină din impunerea condiţiei iniţiale u(0) = u0 , de unde

A = u0 − E. (5.25)

Expresia finală a soluţiei analitice este

u(t) = (u0 − E) exp(−t/τ ) + E. (5.26)

Exerciţiul 5.6:
a) Care este valorea spre care tinde tensiunea pe condensator?
b) Scrieţi un script Matlab main RC.m cu următorul conţinut

clear all;
% rezolva cu MDF un circuit RC simplu
E = 20e-3;
R = 10;
C = 4e-6;
u_initial = 0;
tau = R*C; % constanta de timp
tmax = 10*tau; % timp maxim de simulare

% solutia analitica - pentru referinta

t_ref = linspace(0,tmax,1e6);
uc_ref = (u_initial - E)*exp(-t_ref/tau) + E;
figure(1);
plot(t_ref,uc_ref,’linewidth’,2);
leg{1} = ’analitic’;

Executaţi-l şi verificaţi că aţi răspuns corect la punctul a.

Document disponibil la [Link] gabriela

8 Capitolul 5. Derivarea numerică. Metoda diferenţelor finite.

Rezolvarea cu diferenţe finite

Vom rezolva acum ecuaţia diferenţială cu diferite scheme de derivare numerică. Ecuaţia
(5.14) o rescriem ca
du(t) 1 1
+ u(t) = E. (5.27)
dt τ τ
Vom urmări calculul numeric ı̂n intervalul de timp [0, tmax ] unde tmax = 10τ ı̂ntr-o
reţea echidistantă de N puncte tk , unde pasul de discretizare este

tk+1 − tk = h, pentru k = 1, . . . , N − 1. (5.28)

Vom nota valorile discrete obţinute prin rezolvare numerică cu uk . Ele vor fi aproximaţii
ale mărimii reale u.
uk ≈ u(tk ). (5.29)

Varianta I - utilizarea formulei de diferenţe finite progresive de ordinul 1

Vom discretiza ecuaţia (5.27) prin scrierea ei ı̂n momentul de timp tk şi folosind pentru
derivată o formulă de diferenţe finite progresive de ordinul 1:
uk+1 − uk 1 1
+ uk = E, (5.30)
h τ τ
de unde
h h
uk+1 − uk + uk = E. (5.31)
τ τ
Se observă că mărimea uk+1 poate fi calculată explicit cu formula
 
h h
uk+1 = uk 1 − + E. (5.32)
τ τ
Această formulă este extrem de uşor de implementat.

Exerciţiul 5.7:
a) Completaţi scriptul de mai sus cu

%%%%%%%%%%%%%%%% preprocesare %%%%%%%%%%%%%%%%%%%%%

h = tau/2; % pasul gridului de discretizare pe axa timpului
NI = floor(tmax/h); % nr de intervale de timp
N = NI + 1; % numarul de puncte
t = linspace(0,NI*h,N); % vectorul timp

%%%%%%%%%%%%%%% rezolvare %%%%%%%%%%%%%%%%%%%%%%%%

up = mdf_circuitRC(u_initial,N,h,tau,E,’progresiv1’);

 [Link], Draft din 31 decembrie 2009

5.2. Metoda diferenţelor finite 9

b) Scrieţi o funcţie mdf circuitRC.m de următorul tip

function u = mdf_circuitRC(u_initial,N,h,tau,E,metoda)
% rezolva ecuatia diferentiala du/dt + 1/tau u = 1/tau E cu
% metoda diferentelof finite

u = zeros(N,1);
u(1) = u_initial;
switch metoda
case ’progresiv1’
% calcul explicit
for k = 1:N-1
................. completati...........
end
otherwise
error(’Metoda de discretizare necunoscuta’);
end

c) Completati scriptul main RC.m cu

%%%%%%%%%%%%%%%% postprocesare %%%%%%%%%%%%%%%%%%%%

figure(1);
hold on;
plot(t,up,’r-*’,’linewidth’,2);
leg{2} = ’prog ord 1’;
title(sprintf(’h/tau = %f’,h/tau));

Experimentaţi comportamentul metodei pentru h = 2τ, τ, τ /10. Observaţi că metoda

este instabilă pentru h < τ . Această metodă este cunoscută şi sub numele de metoda
Euler explicită pentru rezolvarea ecuaţiilor diferenţiale ordinare.
Varianta a II-a - utilizarea formulei de diferenţe finite regresive de ordinul
1
Vom discretiza ecuaţia (5.27) prin scrierea ei ı̂n momentul de timp tk şi folosind pentru
derivată o formulă de diferenţe finite regresive de ordinul 1:
uk − uk−1 1 1
+ uk = E, (5.33)
h τ τ
de unde
h h
uk − uk−1 + uk = E. (5.34)
τ τ

Document disponibil la [Link] gabriela

10 Capitolul 5. Derivarea numerică. Metoda diferenţelor finite.

Se observă că mărimea uk poate fi calculată explicit cu formula

   
h h
uk = uk−1 + E / 1 + . (5.35)
τ τ

Şi această formulă este extrem de uşor de implementat.

Exerciţiul 5.8:
a) Completaţi funcţia mdf circuitRC.m cu cazul ”regresiv1”:

case ’regresiv1’
% calcul explicit
for k = ..............completati
u(k) = ...........completati
end

b) Completaţi scriptul cu apelul funcţiei ı̂n cazul regresiv

ur = mdf_circuitRC(u_initial,N,h,tau,E,’regresiv1’);

şi partea de postprocesare cu reprezentarea grafică a soluţiei obţinute, pe acelaşi grafic ca

soluţia analitică şi soluţia obţinută ı̂n prima variantă.

Experimentaţi comportamentul metodei pentru h = 2τ, τ, τ /10. Observaţi că metoda

nu mai este instabilă pentru h < τ . Această metodă este cunoscută şi sub numele de
metoda Euler implicită2 pentru rezolvarea ecuaţiilor diferenţiale ordinare.

Varianta a III-a - utilizarea formulei de diferenţe finite centrate de ordinul

2

Vom discretiza ecuaţia (5.27) prin scrierea ei ı̂n momentul de timp tk şi folosind pentru
derivată o formulă de diferenţe finite centrate de ordinul 2 dată de relaţia (5.4):
uk+1 − uk−1 1 1
+ uk = E, (5.36)
2h τ τ
de unde
2h 2h
−uk−1 + uk + uk+1 = E. (5.37)
τ τ
Se observă că mărimea uk nu mai poate fi calculată explicit. Relaţia (5.37) poate fi scrisă
pentru k = 2, . . . , N − 1 şi reprezintă N − 2 ecuaţii cu N necunoscute. Pentru a obţine
2
Ecuaţia generală este dy/dt = f (t, y) unde f este o funcţie ı̂n general neliniară. Folosirea derivatei
regresive pentru f conduce la o relaţie implicită pentru calculul lui uk , relaţie ce reprezintă o ecuaţie
neliniară. În cazul particular studiat, f este liniară şi de aceea soluţia se poate calcula explicit.

 [Link], Draft din 31 decembrie 2009

5.2. Metoda diferenţelor finite 11

un sistem bine formulat, trebuie să adăugăm ı̂ncă două relaţii. Acestea sunt relaţiile la
capete. La t = 0 vom impune condiţia iniţială:

u(1) = u0 , (5.38)

iar pentru ultimul punct k = N vom scrie ecuaţia discretizată dar ı̂n care vom folosi
pentru derivată o formulă de diferenţe regresive de ordinul 2 dată de relaţia (5.6):
uN −2 − 4uN −1 + 3uN 1 1
+ uN = E, (5.39)
2h τ τ
de unde  
2h 2h
uN −2 − 4uN −1 + 3 + uN = E. (5.40)
τ τ

Exerciţiul 5.9:
a) Completaţi funcţia mdf circuitRC.m cu cazul ”centrat2”:

case ’centrat2’
%% asamblare sistem
A = sparse(N,N);
tl = zeros(N,1);
A(1,1) = 1; % conditia initiala
tl(1) = u_initial;
for k = 2:N-1 % noduri interioare
A(k,k-1) = -1;
.................completati...........
end
A(N,N-2) = 1; % conditia la tmax
....................completati.............
tl(N) = 2*h*E/tau;
%% rezolvare
u = A\tl;

b) Completaţi scriptul cu apelul funcţiei ı̂n cazul centrat şi partea de postprocesare cu
reprezentarea grafică a soluţiei obţinute, pe acelaşi grafic ca soluţia analitică şi soluţia
obţinută ı̂n prima variantă.

Experimentaţi comportamentul metodei pentru h = 2τ, τ, τ /10. Observaţi că metoda

nu este instabilă şi este mai precisă decât implementările anterioare. Acest lucru era de
aşteptat deoarece formulele de derivare numerică de ordinul 2 sunt mai precise decât cele
de ordinul 1. Creşterea acurateţii se face ı̂nsă pe seama creşterii complexităţii algoritmului.

Document disponibil la [Link] gabriela

12 Capitolul 5. Derivarea numerică. Metoda diferenţelor finite.

Pentru verificare, comparaţi graficele pe care le-aţi obţinut cu cele din figurile 5.3 ÷
5.5.

Exerciţiul 5.10:
Completaţi scriptul main RC.m cu comenzi astfel ı̂ncât, ı̂ntr-o altă fereastră grafică să
obţineţi graficele erorilor absolute ı̂n funcţie de pasul de timp.

h/tau = 0.500000 −3 h/tau = 0.500000

x 10
2.5
prog ord 1
0.02 regr ord 1
centr ord 2
0.018 2
0.016

Eroare absoluta [V]

0.014
1.5
0.012
u [V]

0.01
1
0.008

0.006

0.004 analitic 0.5

prog ord 1
0.002 regr ord 1
centr ord 2
0 0
0 1 2 3 4 0 1 2 3 4
t [s] x 10
−4 t [s] x 10
−4

Figura 5.3: Cazul h = τ /2.

h/tau = 1.000000 −3 h/tau = 1.000000

x 10
8
prog ord 1
0.02 regr ord 1
7
centr ord 2
0.018
6
0.016
Eroare absoluta [V]

0.014 5

0.012
u [V]

4
0.01

0.008 3

0.006
2
0.004 analitic
prog ord 1
1
0.002 regr ord 1
centr ord 2
0 0
0 1 2 3 4 0 1 2 3 4
t [s] x 10
−4 t [s] x 10
−4

Figura 5.4: Cazul h = τ .

În capitolul 3, aţi văzut că asamblarea matricelor ı̂n Matlab se poate face mult mai
rapid dacă se construieşte formatul pe coordonate al matricei, după ce, ı̂n prealabil, a fost
alocat un spaţiu de memorie exact cât este necesar pentru stocarea matricei.

Exerciţiul 5.11:
a) Completaţi funcţia mdf circuitRC.m cu cazul ”centrat2b”:

case ’centrat2b’
%% asamblare sistem

 [Link], Draft din 31 decembrie 2009

5.2. Metoda diferenţelor finite 13

h/tau = 2.000000 h/tau = 2.000000

0.04 0.025
analitic
prog ord 1
0.035
regr ord 1
centr ord 2 0.02
0.03

Eroare absoluta [V]

0.025 prog ord 1
0.015 regr ord 1
centr ord 2
u [V]

0.02

0.01
0.015

0.01
0.005
0.005

0 0
0 1 2 3 4 0 1 2 3 4
t [s] x 10
−4 t [s] x 10
−4

Figura 5.5: Cazul h = 2τ .

tic;
no_nnz = 3*(N-1) + 1;
r_idx = zeros(no_nnz,1);
c_idx = zeros(no_nnz,1);
val = zeros(no_nnz,1);
tl = zeros(N,1);
m = 1;
% A(1,1) = 1 % conditia initiala
r_idx(m) = 1;
c_idx(m) = 1;
val(m) = 1;
tl(1) = u_initial;
for k = 2:N-1 % noduri interioare
%A(k,k-1) = -1;
m = m + 1;
r_idx(m) = k;
c_idx(m) = k-1;
val(m) = -1;
.................completati...........
end
....................completati.............
tl(N) = 2*h*E/tau;
A = sparse(r_idx,c_idx,val,N,N);
t2 = toc;
disp(sprintf(’centrat2b - timp asamblare matrice = %e’,t2));
%% rezolvare

Document disponibil la [Link] gabriela

14 Capitolul 5. Derivarea numerică. Metoda diferenţelor finite.

u = A\tl;

b) Completaţi scriptul cu apelul funcţiei ı̂n acest caz. Pentru a verifica faptul că nu
s-au strecurat greşeli, reprezentaţi grafic ı̂ntr-o altă fereastră soluţia obţinută ı̂n cazul
”centrat2b” cu soluţia obţinută ı̂n cazul ”centrat”.
c) Introduceţi instrucţiuni pentru contorizarea timpului de asamblare al matricei ı̂n ambele
variante de implementari. Comparaţi aceşti timpi pentru cazul h = τ /1000.

5.2.2 Studiul regimului electrocinetic staţionar al unui conduc-

tor bidimensional

Exemplul din acest paragraf va conduce la o ecuaţie diferenţială cu derivate parţiale

de ordinul 2, de tip eliptic.

y
σ=0
1V Q V=1 R dV/dn=0
b P
dV/dn=0
b/2 dV/dn=0 ∆V = 0 S
σ
0V V=0

O a/2 a x O dV/dn=0 N

Figura 5.6: Problema 2D de regim electrocinetic: domeniul de calcul (stânga), condiţii

de frontieră (dreapta).

Vom considera un conductor omogen, de conductivitate σ, situat ı̂ntr-un mediu per-

fect izolant. Conductorul are o dimensiune după axa Oz mult mai lungă decât celelalte
două. Figura 5.6 reprezintă o secţiune perpendiculară pe direcţia dimensiunii foarte mari.
Această secţiune este un dreptunghi, de dimensiuni a, b. Conductorul are două borne
supraconductoare, una aflată la potenţial V0 = 1V, iar cealaltă aflată la potenţial nul.
Bornele sunt plasate ca ı̂n figură. Dorim să reprezentăm liniile echipotenţiale şi spectrul
liniilor de curent. Ar fi interesant să calculăm şi mărimi globale cum ar fi energia câmpului
acumulată ı̂n domeniu precum şi rezistenţa pe unitatea de lungime a acestui rezistor. La
aceste aspecte ne vom ı̂ntoarce ı̂nsă după ce vom discuta problema integrării numerice.

 [Link], Draft din 31 decembrie 2009

5.2. Metoda diferenţelor finite 15

Formularea corectă a problemei

Înainte de orice, o problemă de câmp electromagnetic trebuie corect formulată, ı̂n

conformitate cu teorema de unicitate demonstrată pentru regimul studiat.

Problema fiind de regim electrocinetic staţionar, formularea se va face ı̂n V – potenţial

electrocinetic, unde E = −grad V , E fiind intensitatea câmpului electric.

Ecuaţia de ordinul doi satisfăcută de potenţialul electrocinetic este:

−div (σ grad V ) = 0 (5.41)

unde V : D → IR este potenţialul definit pe domeniul bidimensional D = ONPQ. Ecuaţia

(5.41) este o ecuaţie eliptică, de tip Laplace generalizată.

Deoarece domeniul este omogen (conductivitatea σ are acceaşi valoare ı̂n orice punct
al domeniului), ecuaţia se simplifică la o ecuaţie de tip Laplace:

∆V = 0, (5.42)

unde ∆ = div grad este operatorul Laplace, care ı̂n coordonate carteziene (cazul 2D) are
expresia
∂ 2V ∂2V
∆V = + . (5.43)
∂x2 ∂y 2

În consecinţă, ecuaţia diferenţială ce trebuie rezolvată este o ecuaţie cu derivate parţiale:

∂2V ∂2V
+ = 0, (5.44)
∂x2 ∂y 2
unde V = V (x, y) : D → IR.

Pentru buna formulare a problemei, trebuie impuse pentru potenţial condiţiile de fron-
tieră. Frontierele supraconductoare (SN şi QR) au potenţial constant, şi ı̂n consecinţă pe
ele trebuie impuse codiţii Dirichlet ı̂n concordanţă cu valorile potenţialului. Pe frontierele
de lângă izolant (NOQ şi RPS), trebuie impuse condiţii Neumann nule3 .

Condiţiile de frontieră impuse potenţialului sunt:

V = V0 pe QR (Dirichlet) (5.45)
V = 0 pe SN (Dirichlet) (5.46)
∂V
= 0 pe NOQ ∪ RPS (Neumann) (5.47)
∂n
3
În izolant nu există curent de conducţie, ı̂n consecinţă J · n = 0 la frontiera cu izolantul, deci
−σgrad V · n = 0, echivalent cu −σ ∂V
∂n = 0 pe acele frontiere.

Document disponibil la [Link] gabriela

16 Capitolul 5. Derivarea numerică. Metoda diferenţelor finite.

Şi pentru o astfel de problemă se pot găsi rezolvări analitice bazate de exemplu pe
metoda separării variabilelor sau pe aproximarea liniilor de câmp. Deşi nu este dificilă,
prezentarea acestor posibile soluţii nu este atât de simplă ca cea din exemplul anterior.
Scopul nostru este ı̂nsă de a prezenta modul ı̂n care se poate rezolva această problemă cu
metoda diferenţelor finite aşa ı̂ncât ı̂n cele ce urmează ne vom concentra exclusiv asupra
acestui aspect.

Rezolvarea cu diferenţe finite

Faţă de cazul unidimensional discutat ı̂n exemplul anterior, aici gridul de discretizare
este unul bidimensional, iar ı̂n loc de aproximarea unei derivate totale de ordinul 1, acum
avem de aproximat derivate parţiale de ordinul 2.
Gridul de discretizare bidimensional poate fi privit ca fiind obţinut din produsul
cartezian dintre un grid de discretizare pe axa Ox şi un grid de discretizare pe axa Oy:
0 = x1 < x2 < . . . < xnx = a
0 = y1 < y2 < . . . < yny = b.
astfel ı̂ncât, un nod al gridului este identificat de poziţia proiecţiilor sale pe cele două axe:
(xi , yj ) i = 1, . . . , nx , j = 1, . . . , ny .

Metoda diferenţelor finite va determina valorile potenţialului ı̂n nodurile acestui grid:
not
V (xi , yj ) = Vi,j i = 1, nx , j = 1, ny . (5.48)

Forma discretizată (aproximativă) a ecuaţiei potenţialului este obţinută aproximând

prin diferenţe finite ecuaţia (5.44). Forma specifică depinde de tipul nodului: interior
sau pe frontiera Neumann. Nodurile de pe frontiera Dirichlet nu ridică probleme, ecuaţia
asociată unui astfel de nod constând ı̂n simpla atribuire a valorii potenţialului nodului
respectiv.
Ecuaţia asociată unui nod interior:
Ecuaţia aproximativă pentru un nod interior se deduce ı̂nlocuind derivatele parţiale
după x şi după y cu formule de tipul (5.7).
Pentru a simplifica scrierea formulelor, vom presupune că cele două griduri pe Ox şi
Oy sunt uniforme, cu paşi hx şi respectiv hy . Derivatele parţiale se vor ı̂nlocui cu:
∂ 2V Vi+1,j − 2Vi,j + Vi−1,j
2
(xi , yj ) = , (5.49)
∂x h2x
∂2V Vi,j+1 − 2Vi,j + Vi,j−1
2
(xi , yj ) = . (5.50)
∂y h2y

 [Link], Draft din 31 decembrie 2009

5.2. Metoda diferenţelor finite 17

Ecuaţia discretizată asociată unui nod interior devine

Vi+1,j − 2Vi,j + Vi−1,j Vi,j+1 − 2Vi,j + Vi,j−1

+ = 0, (5.51)
h2x h2y
sau  
1 1 1 1 1 1
2 2
+ 2 Vi,j − 2
Vi+1,j − 2 Vi−1,j − 2 Vi,j+1 − 2 Vi,j−1 = 0. (5.52)
hx hy hx hx hy hy
Relaţia (5.52) arată că potenţialul ı̂ntr-un nod interior este o combinaţie liniară a potenţialelor
nodurilor ı̂nvecinate. Dacă paşii de discretizare pe cele două axe ar fi egali (hx = hy ),
atunci potenţialul ı̂ntr-un nod interior este media aritmetică a potenţialelor celor patru
noduri vecine.

y
n
C
C
D1
D O B O B

A A
x
x
Figura 5.7: Nodurile marcate intervin ı̂n scrierea ecuaţiilor: stânga - cazul unui nod
interior; dreapta - cazul unui nod pe frontieră Neumann dreaptă.

Exerciţiul 5.12:
În cazul unor griduri neuniforme, folosind notaţiile din fig. 5.7, arătaţi că ecuaţia dis-
cretizată asociată unui nod interior este
 
1 1 1 1
VO + − VA − VB −
hB hD hA hC hA (hA + hC ) hB (hB + hD )
1 1
−VC − VD = 0, (5.53)
hC (hA + hC ) hD (hB + hD )

unde hA = kOAk, hB = kOBk, hC = kOCk, hD = kODk.

Ecuaţia asociată unui nod pe frontieră dreaptă:

În cazul unui punct pe frontieră O (fig. 5.7- dreapta), să notăm cu g = −σ ∂V
∂n
condiţia
de frontieră Neumann şi cu gO valoarea ei medie ı̂n punctul O:

gOA hA + gOC hC
gO = (5.54)
hA + hC

Document disponibil la [Link] gabriela

18 Capitolul 5. Derivarea numerică. Metoda diferenţelor finite.

unde gOA este condiţia de frontieră asociată segmenului OA iar gOC este condiţia de
frontieră asociată segmentului OC.

Pentru deducerea ecuaţiei aproximative pentru punctul O, vom considera un nod “fan-
tomă” D1 , situat ı̂n exterior, la o distanţă hD = kOD1 k de punctul O. Pentru simplitate
putem considera hD = hB .

Din condiţia de frontieră Neumann, deducem valorea potenţialului fantomă ı̂n funcţie
de valoarea acestei condiţii. Pentru cazul din figură (normala exterioară ı̂n sens contrar
axei x) rezultă:
∂V g
= . (5.55)
∂x σ
Scriind relaţia aproximativă (5.4) pentru condiţia de frontieră Neumann, rezultă că:
gO
VD1 = VB − 2hB . (5.56)
σ

Pentru nodul O se scrie acum o ecuaţie de tipul (5.53), ca pentru un nod interior şi
ı̂nlocuind expresia (5.56) rezultă ecuaţia finală (5.57):
 
1 1 1 1 1 gO
VO 2
+ − VA − VB 2 − VC =− . (5.57)
hB hA hC hA (hA + hC ) hB hC (hA + hC ) σhB

În cazul unei condiţii Neumann nule, relaţia devine:

 
1 1 1 1 1
VO 2
+ − VA − VB 2 − VC = 0. (5.58)
hB hA hC hA (hA + hC ) hB hC (hA + hC )

Ecuaţia asociată unui nod pe frontieră - colţ exterior:

Nodurile situate ı̂n colţuri reprezintă o problemă pentru metoda diferenţelor finite.
Aceasta deoarece pentru un colţ nu poate fi definită direcţia normalei4 . Dacă notăm
g = −σ ∂V ∂n
, şi notăm cu gOA valoarea lui g pe segmentul OA (fără capătul O) şi cu gOB
valoarea lui g pe segmentul OB (fără capătul O), atunci vom presupune că ı̂n nodul O (fig.
5.8) este cunoscută derivata după o direcţie care rezultă din prelungirea ı̂n O a valorilor
funcţiei g.

Pentru deducerea ecuaţiei vom considera două noduri “fantomă” C1 şi D1 , situate ı̂n
exterior, la distanţele hD = kOD1 k şi hC = kOC1 k de punctul O. Pentru simplitate putem
considera hD = hB şi hC = hA .
4
În realitate, colţul reprezintă o modelare brutală a realităţii. Trecerea de la o suprafaţă la alta se
face prin racordări.

 [Link], Draft din 31 decembrie 2009

5.2. Metoda diferenţelor finite 19

C1
n
D1
O B

A
n
x

Figura 5.8: Nodurile marcate intervin ı̂n scrierea ecuaţiei unui nod pe frontieră colţ
exterior.

Din condiţia de frontieră Neumann pe segmentul OB, scrisă ı̂n nodul O, putem deduce
valoarea potenţialului fantomă C1 :
gOB
VC1 = VA − 2hA . (5.59)
σ

Din condiţia de frontieră Neumann pe segmenul OA, scrisă ı̂n nodul O, putem deduce
valoarea potenţialului fantomă D1 :
gOA
VD1 = VB − 2hB . (5.60)
σ

Pentru nodul O se scrie acum o ecuaţie de tipul (5.53), ca pentru un nod interior, şi
ı̂nlocuind expresiile (5.59) şi (5.60) rezultă ecuaţia finală (5.61):
 
1 1 1 1 gOB gOA
VO 2
+ 2 − VA 2 − VB 2 = − − (5.61)
hA hB hA hB σhA σhB

În cazul unor condiţii Neumann nule, relaţia devine:

 
1 1 1 1
VO 2
+ 2 − VA 2 − VB 2 = 0 (5.62)
hA hB hA hB

În concluzie, discretizarea problemei conduce la un sistem de ecuaţii algebric liniar,

prin a cărui rezolvare vom obţine valorile potenţialelor ı̂n nodurilei gridului de discretizare.
Pentru asamblarea acestui sistem cu ajutorul unui algoritm, este util ca nodurile să fie
numerotate astfel ı̂ncât necunoscutele să reprezinte un vector, nu o matrice cum sugerează
notaţiile de până acum. Vom numerota nodurile de jos ı̂n sus şi de la stânga la dreapta,
ca ı̂n figura 5.9. Se poate verifica uşor că relaţia ı̂ntre numărul nodului k şi poziţiile
proiecţiilor sale pe cele două axe i şi j este

k = (i − 1)ny + j, i = 1, . . . , nx , j = 1, . . . , ny . (5.63)

Document disponibil la [Link] gabriela

20 Capitolul 5. Derivarea numerică. Metoda diferenţelor finite.

ny 2 ny nx ny

j k

1 ny+1 (i−1)ny+1 (nx−1)ny+1

Figura 5.9: Numerotarea nodurilor.

unde nx şi ny reprezintă numărul de puncte de discretizare pe axele x, respectiv y.

Exerciţiul 5.13:
Scrieţi ecuaţiile discretizate pentru cazul nx = ny = 3, hx = hy = 1. Completaţi matricea
coeficienţilor A şi vectorul termenilor liberi t pentru sistemul de ecuaţii de rezolvat:

A= t=

În acest moment putem trece la conceperea algoritmului. Vom scrie funcţii speciale
pentru: citirea datelor problemei, generarea gridului, asamblarea sistemului de ecuaţii de
rezolvat şi pentru postprocesare.
Să ı̂ncepem cu citirea datelor problemei. În principiu, am putea citi datele de la
tastatură sau din fişiere, dar cum scopul nostru este de a construi şi testa partea cea mai
complicata a algoritmului (preprocesarea), vom prefera ca datele problemei să le ”citim”
prin instanţiere, ı̂ntr-o funcţie numită citire date elcin. Datele problemei sunt grupate
ı̂ntr-o structură numită date, ce are mai multe câmpuri.
Exerciţiul 5.14:
a) Scrieţi un script mdf elcin.m cu următorul conţinut:

clear all;

 [Link], Draft din 31 decembrie 2009

5.2. Metoda diferenţelor finite 21

% rezolva ecuatia Laplace V = 0

% intr-un domeniu dreptunghiular de dimensiuni a, b
% prin metoda diferentelor finite
% Conditii de frontiera
% * pe latura de sus, in stanga exista un electrod la potential V1
% * pe latura din dreapta, jos exista un electrod la potential 0
% * in rest cond Neumann nule

%% citirea datelor problemei

date = citire_date_elcin();

b) Scrieţi o funcţie citire date elcin.m cu următorul conţinut

function date = citire_date_elcin()

date.a = 2; % dim pe OX
date.b = 2; % dim pe Oy
date.T1 = date.a/2; % dim terminalului de sus
date.T2 = date.b/2; % dim terminalului din dreapta
date.V1 = 1; % potentialul lui T1
date.V2 = 0; % potentialul lui T2

c) Rulaţi codul scris până acum pentru a observa dacă s-au strecurat greşeli.

Acum ı̂ncepe partea de preprocesare. Mai ı̂ntâi trebuie pregătite informaţiile legate
de gridul de discretizare. Vom presupune că acesta este uniform pe cele două axe.
Exerciţiul 5.15:
a) Completaţi scriptul cu comenzile

%% ====== PREPROCESARE =======

%% alegerea gridului
my_grid = generare_grid_elcin(date);

b) Scrieţi o funcţie generare grid elcin.m cu următorul conţinut

function grid = generare_grid_elcin(date)

a = date.a;
b = date.b;
nx = 3; % nr pct de discretizare pe Ox
ny = 3; % nr pct de discretizare pe Oy

Document disponibil la [Link] gabriela

22 Capitolul 5. Derivarea numerică. Metoda diferenţelor finite.

% pasi de discretizare pe Ox si Oy
hx = ..................... completati............
hy = ..................... completati............
% numarul total de noduri
N = ..................... completati.............
[Link] = nx;
[Link] = ny;
[Link] = hx;
[Link] = hy;
grid.x = linspace(0,a,nx);
grid.y = linspace(0,b,ny);
grid.N = N;

c) Ce reprezintă grid.x şi grid.y?

Urmează acum partea cea mai complicată, asamblarea sistemului de ecuaţii.

Exerciţiul 5.16:
a) Completaţi scriptul cu comenzile

%% asamblare sistem
[A,tl] = asamblare_mdf_elcin(date,my_grid);

b) Scrieţi o funcţie asamblare mdf elcin.m cu următorul conţinut:

function [A,tl] = asamblare_mdf_elcin(date,my_grid)

N = my_grid.N;
nx = my_grid.nx;
ny = my_grid.ny;
hx = my_grid.hx;
hy = my_grid.hy;
a = date.a;
b = date.b;
T1 = date.T1;
T2 = date.T2;
V1 = date.V1;
V2 = date.V2;

%% asamblarea sistemului de ecuatii

 [Link], Draft din 31 decembrie 2009

5.2. Metoda diferenţelor finite 23

A = sparse(N,N);
tl = sparse(N,1);
%% nodurile interioare
for i = 2:nx-1
for j = 2:ny-1
k = (i-1)*ny + j;
ksus = k+1;
kjos = k-1;
kst = k - ny;
kdr = k + ny;
A(k,k) = 2*(1/hx^2+1/hy^2);
A(k,ksus) = -1/hy^2;
................. completati .................
end
end
%% cond. de frontiera Neumann - segm vertical stanga
i = 1;
for j = 2:ny-1
k = (i-1)*ny + j;
ksus = k+1;
kjos = k-1;
kdr = k + ny;
A(k,k) = 2*(1/hx^2+1/hy^2);
A(k,ksus) = -1/hy^2;
A(k,kjos) = -1/hy^2;
A(k,kdr) = -2/hx^2;
end
%% cond. de frontiera Neumann - segm orizontal jos
j = 1;
for i = 2:nx-1
.................... completati ......................
end
%% cond. de frontiera Neumann - segm orizontal sus (partial)
iT1dr = floor(T1/a*nx)+1; % idx pe i atins in dreapta de T1
j = ny;
for i = iT1dr+1:nx-1
.................... completati ......................

Document disponibil la [Link] gabriela

24 Capitolul 5. Derivarea numerică. Metoda diferenţelor finite.

end
%% cond. de frontiera Neumann - segm vertical dreapta (partial)
jT2sus = floor(T2/b*ny)+1; % idx pe j atins sus de T2
i = nx;
for j = jT2sus+1:ny-1
.................... completati ......................
end
%% cond de frontiera Neumann - colt stanga jos
k = 1;
kdr = ny+1;
ksus = 2;
A(k,k) = 2*(1/hx^2+1/hy^2);
A(k,ksus) = -2/hy^2;
A(k,kdr) = -2/hx^2;
%% cond de frontiera Neumann - colt dreapta sus
.................... completati ......................
%% cond de frontiera Dirichlet - T1
j = ny;
for i = 1:iT1dr
k = (i-1)*ny + j;
A(k,k) = 1;
tl(k) = V1;
end
%% cond de frontiera Dirichlet - T2
.................... completati ......................

c) Verificaţi codul scris cu comanda mlint. Executaţi programul scris până acum şi
verificaţi că matricea coeficienţilor şi vectorul termenilor liberi au valorile pe care le-aţi
obţinut la exerciţiul 5.13

În Matlab, rezolvarea acestui sistem este extrem de uşoară din punct de vedere al
utilizatorului:

Exerciţiul 5.17:
Completaţi scriptul cu comenzile

%% ====== REZOLVARE =======

V = A\tl;

 [Link], Draft din 31 decembrie 2009

5.2. Metoda diferenţelor finite 25

Urmează etapa de postprocesare. Vom reprezenta liniile echipotenţiale şi vom desena
gridul de discretizare.

Exerciţiul 5.18:
a) Completaţi scriptul cu comenzile

%% ====== POSTPROCESARE =======

echipotentiale_elcin(my_grid,V);
draw_my_grid(my_grid);

b) Scrieţi o funcţie echipotentiale elcin.m cu următorul conţinut

function echipotentiale_elcin(my_grid,V)
nx = my_grid.nx;
ny = my_grid.ny;
x = my_grid.x;
y = my_grid.y;
v = zeros(nx,ny);
for i = 1:nx
for j = 1:ny
k = (i-1)*ny + j;
v(i,j) = V(k);
end
end
..................completati..........(ceva similar ati facut la capitolul 1)

c) Scrieţi o funcţie draw my grid.m cu următorul conţinut

function draw_my_grid(my_grid)
nx = my_grid.nx;
ny = my_grid.ny;
x = my_grid.x;
y = my_grid.y;
for i = 1:nx
xcrt = x(i);
hold on;
plot([xcrt xcrt],[y(1),y(ny)],’k:’);
end
..................completati..........

Document disponibil la [Link] gabriela

26 Capitolul 5. Derivarea numerică. Metoda diferenţelor finite.

2 1 1 21 1
8 6
9 0. 7 0.
0. 0.

0.7
5 0.9
1.8 0.9 0. 1.8 0.9

0.8 5
1.6 0.8 0.4 1.6 0.
6
0.8 0.
3 0.8 0.
0.
6 0.4
1.4 0.7 1.4 0.7
5
0.7 0.
1.2 2 1.2 0.7 0.3
0.6 0. 1 5
0.
4 0. 0.6 0. 4 0.2
1 1 0.

0 .3

0.1
6 0.6
0. 5
0.8 0. 0.8

0 .3
0.6 0.6 5
0.
0.2

0.2
0.4
0.4 0.4
0.1
0.4

5
0.
0.3

0.2 0.2

0.1
5
0.
0 0
0 0.5 1 1.5 2 0 0.5 1 1.5 2

Figura 5.10: Linii echipotenţiale pentru un grid cu nx = ny = 3 (stânga) şi un grid cu

nx = ny = 21 (dreapta).

d) Verificaţi că aţi lucrat corect, comparând rezultatele cu cele din fig. 5.10.

Exerciţiul 5.19:
a) Observaţi cu atenţie fig. 5.10 - stânga şi descrieţi ideea care stă la baza algoritmului
de desenare al curbelor de echivalori.
b) Propuneţi ı̂mbunătăţiri posibile pentru funcţia de asamblare.
c) Ştiind că E = −grad V , scrieţi formulele cu care se poate calcula câmpul electric ı̂n
nodurile gridului de discretizare.

5.2.3 Studiul propagării undei scalare

În acest paragraf vom rezolva cu metoda diferenţelor finite cel mai simplu exemplu de
ecuaţie hiperbolică, şi anume ecuaţia undei scalare:

∂ 2u 2
2∂ u
= v , (5.64)
∂t2 ∂x2
unde mărimea scalară necunoscută u(x, t) : [0, a] × [0, T ] → IR depinde de o singură
coordonată spaţială şi de timp, iar v este o constantă cunoscută. Soluţia acestei ecuaţii
este o undă care se propagă cu viteză v.

Buna formulare a acestei probleme necesită impunerea

• condiţiei iniţiale u(x, 0) = h0 (x),

• condiţiilor la capete - relaţii ı̂n care intervin mărimile u(0, t) = h1 (t) şi u(a, t) =
h2 (t).

 [Link], Draft din 31 decembrie 2009

5.2. Metoda diferenţelor finite 27

Se poate demonstra că soluţia ecuaţiei undei scalare (5.64) este

u(x, t) = ud(x, t) + ui(x, t) + U0 , (5.65)

unde ud(x, t) se numeşte undă directă şi este o funcţie care se poate scrie sub forma

ud(x, t) = f (x − vt), (5.66)

iar ui(x, t) se numeşte undă inversă şi este o funcţie care se poate scrie sub forma

ud(x, t) = g(x + vt). (5.67)

Funcţiile f şi g rezultă ı̂n mod univoc, din impunerea condiţiilor iniţiale şi condiţiilor
la capete. Mărimea f (x − vt) este o undă direcţă deoarece o anumită valoare a acestei
funcţii, ı̂ntr-un anumit punct x şi ı̂ntr-un anumit moment de timp t se va regăsi după un
interval de timp ∆t ı̂n punctul x + ∆x, unde ∆x = v∆t. Unda directă se propagă deci ı̂n
sensul pozitiv al axei Ox. Un raţionament similar se poate face pentru unda inversă.

Următorul exerciţiul vă permite să vizualizaţi conceptul de undă. El nu reprezintă

rezolvarea numerică a vreunei ecuaţii diferenţiale.

Exerciţiul 5.20:
a) Scrieţi o funcţie cu următorul conţinut

function demo_unde()
a = 1; % domeniul spatial este [0,a]
T = 1; % domeniul temporal este [0, T]
N = 500; % nr pct pentru discretizarea spatiala
M = 500; % nr pct pentru discretizarea temporala

v = 1; % viteza de propagare
dt = T/(M-1); % pas de discretizare temporal - uniform

x = linspace(0,a,N);
t = 0;
u = ud(x,t,v);
plot(x,u);
drawnow;
for j = 1:M
t = t + dt;
u = ud(x,t,v);

Document disponibil la [Link] gabriela

28 Capitolul 5. Derivarea numerică. Metoda diferenţelor finite.

plot(x,u);
drawnow;
end

function rez = ud(x,t,v)

rez = sin(2*pi*t - 2*pi*x/v);

Observaţi rezultatul executării ei.

b) Modificăti funcţia astfel ı̂ncât să vizualizaţi propagarea unei unde inverse.
c) Modificaţi funcţia astfel ı̂ncât să vizualizaţi propagarea sumei dintre unda directă şi
unda inversă. Veţi avea nevoie să fixaţi axele reprezentării grafice. Pentru aceasta, ı̂nainte
de prima comanda drawnow adăugaţi

xmin = 0; xmax = a; ymin = -2.2; ymax = 2.2;

axis([xmin, xmax,ymin,ymax]);

şi ı̂nainte de a două comandă drawnow adăugaţi doar comanda axis ca mai sus. Ceea ce
obţineţi sunt unde staţionare.

Vom concepe acum un algoritm bazat pe metoda diferenţelor finite, care rezolvă ecuaţia
(5.64). Ca exemplu numeric vom considera a = 1 m, v = 1 m/s.

Avem nevoie de o discretizare spaţială şi de una temporală. Pentru a simplifica

prezentarea, vom presupune că ambele discretizări sunt uniforme şi vom nota cu ∆z pasul
discretizării spaţiale şi cu ∆t pasul discretizării temporale. Vom nota cu N numărul de
puncte de discretizare spaţiale şi cu M numărul de paşi de timp simulaţi. În consecintă
timpul maxim de simulare este T = M ∆t.

Într-un prim exemplu, vom considera condiţii iniţiale nule şi condiţii Dirichlet la ambele
capete:

• condiţia iniţială nulă u(x, 0) = h0 (x) = 0,

• condiţia la capătul din stânga - excitaţia cu un impuls Gauss: u(0, t) = h1 (t) =

exp(−(t − T /10).2 /(2 ∗ (T /50)2 ))

• condiţia la capătul din dreapta u(a, t) = h2 (t) = 0.

Mărimea u(x, t) este discretizată atât ı̂n spaţiu cât şi ı̂n timp. Prin rezolvarea cu
metoda diferenţelor finite, vom obţine aproximaţii pentru valorile reale u(xk , tj ).
(j)
u(xk , tj ) ≈ uk , k = 1, . . . , N, j = 1, . . . , M. (5.68)

 [Link], Draft din 31 decembrie 2009

5.2. Metoda diferenţelor finite 29

Alegând pentru derivatele ce intervin ı̂n (5.64) formule de derivare ce provin din poli-
nomul de interpolare de ordin doi, rezultă următoarea relaţie discretizată:
(j−1) (j) (j+1) (j) (j) (j)
uk − 2uk + uk 2 uk−1 − 2uk + uk+1
=v , (5.69)
(∆t)2 (∆x)2

de unde rezultă că mărimea u, ı̂ntr-un anumit punct k şi la un anumit moment de timp
j + 1 depinde explicit de valoarea ı̂n acel punct şi ı̂n cele ı̂nvecinate la momentul de timp
anterior j şi de valoarea ı̂n acel punct la momentul j − 1:
 2 
(j+1) v∆t (j) (j) (j)

(j)
uk = uk−1 − 2uk + uk+1 +2uk −ukj−1 , k = 2, . . . , N −1, j = 0, . . . , M −1.
∆x
(5.70)
Momentul −1 ı̂l vom considera identic cu momentul iniţial 0.

Exerciţiul 5.21:
a) Scrieţi un script cu următorul conţinut

clear all;
a = 1; % domeniul spatial este [0,a]
T = 1; % intervalul de timp necesar propagarii pe distanta a
N = 500; % nr pct pentru discretizarea spatiala
M = 500; % nr pct pentru discretizarea temporala a lui T

v = 1; % viteza de propagare
dx = a/(N-1); % pas de discretizare spatial - uniform
dt = T/(M-1); % pas de discretizare temporal - uniform

r = v*dt/dx;
r2 = r^2;

% conditia initiala
x = linspace(0,1,N);
solv = zeros(1,N);
solvv = solv;
plot(x,solv);
xmin = 0; xmax = a; ymin = -2.2; ymax = 2.2;
axis([xmin, xmax,ymin,ymax]);

MM = 2*M;

Document disponibil la [Link] gabriela

30 Capitolul 5. Derivarea numerică. Metoda diferenţelor finite.

for idx_t = 1:MM

t = idx_t*dt;
% conditie de frontiera la stanga
sol(1) = exp(-(t - T/10).^2 / (2 * (T/50)^2));
% conditie de frontiera la dreapta
sol(N) = 0;
for idx_x = 2:N-1
............... completati .................
end
solvv = solv;
solv = sol;
plot(x,solv);
title(sprintf(’tmax = %4.2e, no_t = %d, idx_t = %d’,T,MM,idx_t))
axis([xmin, xmax,ymin,ymax]);
drawnow
end

b) Observaţi ce se ı̂ntâmplă după reflexie (fig. 5.11). În momentul reflexiei, apare o undă
inversă deoarece condiţia Dirichlet la capătul din dreapta impune ca mărimea să fie zero.
În consecinţă, unda inversă la capătul din dreapta este exact opusul undei directe, şi ea
se propagă ı̂n sens contrar axei Ox.
c) Cazul implementat mai sus corespunde situaţiei ∆x = v∆t. Observaţi ce se ı̂ntâmplă
dacă ∆x > v∆t. Alegeţi de exemplu N = 100, M = 500.
c) Observaţi ce se ı̂ntâmplă dacă ∆x < v∆t. Alegeţi de exemplu N = 100, M = 500.

tmax = 1.00e+000, no = 1000, idx = 400 tmax = 1.00e+000, no = 1000, idx = 600
t t t t

2 2

1.5 1.5

1 1

0.5 0.5

0 0

−0.5 −0.5

−1 −1

−1.5 −1.5

−2 −2

0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1

Figura 5.11: Propagarea unui impuls Gauss (stânga) şi rezultatul reflexiei după atingerea
unei frontiere pe care s-a impus condiţie Dirichlet nulă (dreapta).

 [Link], Draft din 31 decembrie 2009

5.2. Metoda diferenţelor finite 31

În concluzie, atunci când o problemă necesită atât o discretizare spaţială cât şi una
temporală, pentru ca soluţia numerică să fie stabilă, este necesar să fie ı̂ndeplinită condiţia
lui Courant
|v|∆t ≤ ∆x. (5.71)

Exerciţiul 5.22:
a) Reveniţi la cazul ∆x = v∆t şi schimbaţi condiţia la capătul din stânga cu u(0, t) =
h1 (t) = sin(10πt).
b) Observaţi apariţia undelor staţionare, de data aceasta obţinute prin rezolvarea ecuaţiei
undei scalare (fig. 5.12).
tmax = 1.00e+000, no = 1000, idx = 400 tmax = 1.00e+000, no = 1000, idx = 600
t t t t

2 2

1.5 1.5

1 1

0.5 0.5

0 0

−0.5 −0.5

−1 −1

−1.5 −1.5

−2 −2

0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1

Figura 5.12: Propagarea unui unde sinusoidale (stânga) şi rezultatul reflexiei după atin-
gerea unei frontiere pe care s-a impus condiţie Dirichlet nulă (dreapta).

În problemele ı̂n care apare propagare, este util uneori ca frontiera domeniului spaţial
să fie modelată ca o frontieră absorbantă, ”invizibilă” din punct de vedere al propagării
mărimilor ı̂n domeniul spaţial modelat. O condiţie de frontieră absorbantă pentru capătul
din dreapta al problemei studiate este
∂u ∂u
+v = 0. (5.72)
∂t ∂z

Exerciţiul 5.23:
a) Implementaţi condiţia de frontieră absorbantă (5.72) folosind o discretizare cu diferenţe
progresive de ordinul 1 pentru derivata temporală şi diferenţe regresive de ordinul 1 pentru
derivata spaţială:
(j) (j)
u − uN −1 u(j+1) − uN
(j)
v N + N = 0. (5.73)
∆z ∆t

b) Verificaţi codul implementat pentru cazul ı̂n care la capătul din stânga este impusă
condiţia Dirichlet folosită la exerciţiul 5.21

Document disponibil la [Link] gabriela