Post on 26-Jun-2018
transcript
Formulare Formule de der. Algoritmi Der.ord.sup Der.part. MDF-1 MDF-2 MDF-3
Derivarea numerica. Introducere în metodadiferentelor finite.
Conf.dr.ing. Gabriela Ciuprina
Universitatea "Politehnica" Bucuresti, Facultatea de Inginerie Electrica,Departamentul de Electrotehnica
Suport didactic pentru disciplina Algoritmi Numerici, 2012
Formulare Formule de der. Algoritmi Der.ord.sup Der.part. MDF-1 MDF-2 MDF-3
Cuprins
Formularea problemei derivarii numerice
Formule de derivare numerica
Algoritmi numerici pentru derivarea functiilor
Derivate de ordin superior
Derivate partiale
Metoda diferentelor finite - exemplul 1 (ODE)
Metoda diferentelor finite - exemplul 2 (PDE, eliptic)
Metoda diferentelor finite - exemplul 3 (PDE, hiperbolic)
Formulare Formule de der. Algoritmi Der.ord.sup Der.part. MDF-1 MDF-2 MDF-3
Importanta evaluarii derivatelor
• relatii utile pentru evaluarea unor marimi
i(t) = −dqDΣ
dt
u(t) = −dφSΓ
dt• evaluarea senzitivitatilor
∂F
∂p
Ideea: derivarea numerica se bazeaza pe interpolare.
• ecuatii sau sisteme de ecuatii diferentiale
⇒ Metoda diferentelor finite
Formulare Formule de der. Algoritmi Der.ord.sup Der.part. MDF-1 MDF-2 MDF-3
Formularea problemei derivarii numericeSe da functia f : [a, b] → IR (cunoscuta prin date sau prin cod)Se cere evaluarea numerica a derivatei f ′(x0), unde x0 ∈ [a, b].Matematic:
f ′(x0) = limx→x0
f (x)− f (x0)
x − x0. (1)
Nu se poate aplica aceasta formula:
• nu se poate trece la limita ⇒ erori de trunchiere;
• nedeterminare 0/0 ⇒ probleme numerice când x este alesprea aproape de x0.
Ideea: g(x) interpolarea sau aproximarea functiei si
g′(x) =dgdx
≈ f ′(x) =dfdx
. (2)
Precizia cu care se face interpolarea sau aproximarea vadetermina precizia de calcul a derivatei numerice.
Formulare Formule de der. Algoritmi Der.ord.sup Der.part. MDF-1 MDF-2 MDF-3
Formule de derivare numerica - cu eroare de ordinul 1
g - interpolare liniara pe portiuni.
0 2 4 6 8 105
10
15
20
25
30
35
x
fg
Functia reala f si interpolarea ei pe portiuni g.
0 2 4 6 8 10−15
−10
−5
0
5
10
15
20
x
f′
g′
Aproximarea derivatei f ′ cu valoarea g′ .
Formulare Formule de der. Algoritmi Der.ord.sup Der.part. MDF-1 MDF-2 MDF-3
Formule de derivare numerica - cu eroare de ordinul 1
Derivare progresiv a de ordinul 1 :
y ′
k =yk+1 − yk
xk+1 − xk, (3)
Derivare regresiv a de ordinul 1 :
y ′
k =yk − yk−1
xk − xk−1. (4)
Ordinul erorii de trunchiere specifice acestei aproximari sepoate estima cu ajutorul dezvoltarii în serie Taylor.
Formulare Formule de der. Algoritmi Der.ord.sup Der.part. MDF-1 MDF-2 MDF-3
Formule de derivare numerica - cu eroare de ordinul 1Derivata progresiva y ′
k = (yk+1 − yk )/(xk+1 − xk )
f (x) = f (xk ) + f ′(xk )(x − xk ) +f ′′(ξ)
2!(x − xk )
2 (5)
x = xk+1 ⇒
yk+1 = yk + f ′(xk )(xk+1 − xk ) +f ′′(ξ)
2!(xk+1 − xk )
2, (6)
Eroarea de trunchiereyk+1 − yk
xk+1 − xk− f ′(xk ) =
f ′′(ξ)2!
(xk+1 − xk ), (7)
Pp. |f ′′(x)| ≤ M2
Notam h = xk+1 − xk
|et | ≤M2
2h, (8)
Pe scurt: |et | =O(h), eroare de ordin 1(1 = puterea lui h, sau gradul polinomului de interpolare folosit)
Formulare Formule de der. Algoritmi Der.ord.sup Der.part. MDF-1 MDF-2 MDF-3
Formule de derivare numerica - cu eroare de ordinul 1
Derivata regresiva y ′
k = (yk − yk−1)/(xk − xk−1)
f (x) = f (xk ) + f ′(xk )(x − xk ) + O((x − xk )2) (9)
x = xk−1 si notând xk − xk−1 = h
yk−1 = yk − f ′(xk )h + O(h2), (10)
yk − yk−1
h− f ′(xk ) = O(h). (11)
|et | =O(h), eroare de ordin 1
Formulare Formule de der. Algoritmi Der.ord.sup Der.part. MDF-1 MDF-2 MDF-3
Formule de derivare numerica - cu eroare de ordinul 2
Dezavantajul g - lpp: nu este derivabila în noduri.Remediu: cresterea gradului polinomului de interpolare.Polinom de interpolare de grad 2:
g(x) =(x − xk )(x − xk+1)
(xk−1 − xk )(xk−1 − xk+1)yk−1 +
(x − xk−1)(x − xk+1)
(xk − xk−1)(xk − xk+1)yk +
+(x − xk−1)(x − xk )
(xk+1 − xk−1)(xk+1 − xk )yk+1, (12)
g′(x) =(x − xk ) + (x − xk+1)
(xk−1 − xk )(xk−1 − xk+1)yk−1 +
(x − xk−1) + (x − xk+1)
(xk − xk−1)(xk − xk+1)yk +
+(x − xk−1) + (x − xk )
(xk+1 − xk−1)(xk+1 − xk )yk+1. (13)
Formulare Formule de der. Algoritmi Der.ord.sup Der.part. MDF-1 MDF-2 MDF-3
Formule de derivare numerica - cu eroare de ordinul 2Notam: h1 = xk − xk−1 si h2 = xk+1 − xk .Evaluând polinomul derivat (34) în xk se obtine formula dederivare centrata de ordinul 2
y ′
k = −h2
h1(h1 + h2)yk−1 −
h1 − h2
h1h2yk +
h1
h2(h1 + h2)yk+1. (14)
Evaluând polinomul derivat (34) în punctele xk−1 si xk+1 seobtine formula de derivare progresiva de ordinul 2
y ′
k−1 = −2h1 + h2
h1(h1 + h2)yk−1+
h1 + h2
h1h2yk −
h1
h2(h1 + h2)yk+1. (15)
si, respectiv, de derivare regresiva de ordinul 2
y ′
k+1 = −h2
h1(h1 + h2)yk−1−
h1 + h2
h1h2yk +
h1 + 2h2
h2(h1 + h2)yk+1. (16)
Formulare Formule de der. Algoritmi Der.ord.sup Der.part. MDF-1 MDF-2 MDF-3
Formule de derivare numerica - cu eroare de ordinul 2
În cazul unui pas echidistant h1 = h2 = hDerivare centrat a
y ′
k =yk+1 − yk−1
2h(17)
Derivare progresiv a
y ′
k−1 =−3yk−1 + 4yk − yk+1
2h(18)
Derivare regresiv a
y ′
k+1 =−yk−1 − 4yk + 3yk+1
2h(19)
Formulare Formule de der. Algoritmi Der.ord.sup Der.part. MDF-1 MDF-2 MDF-3
Formule de derivare numerica - cu eroare de ordinul 2
f (x) = f (xk )+f ′(xk )(x−xk )+f ′′(xk )
2!(x−xk )
2+O((x−xk )3). (20)
Evaluam în punctele xk−1 si xk+1 (pp. pas echidistant)xk+1 − xk = xk − xk−1 = h, rezulta:
f (xk+1) = f (xk ) + f ′(xk )h +f ′′(xk )
2!h2 + O(h3), (21)
f (xk−1) = f (xk )− f ′(xk )h +f ′′(xk )
2!h2 − O(h3). (22)
Scadem relatiile
f (xk+1)− f (xk−1) = 2f ′(xk )h + O(h3), (23)
de undeyk+1 − yk−1
2h− f ′(xk ) = O(h2). (24)
Formula de derivare centrata are eroarea de trunchiere deordinul 2
Formulare Formule de der. Algoritmi Der.ord.sup Der.part. MDF-1 MDF-2 MDF-3
Formule de derivare numericaFormulele de ordinul 2 sunt mai precise decât formulele deordinul 1.
f
panta acestei drepte = f′(xk)
panta acestei drepte = derivata regresiva de ordin 1panta acestei drepte = derivata progresiva de ordin 1panta acestei drepte = derivata centrata de ordin 2
(xk−1
,yk−1
)
(xk+1
,yk+1
)
(xk,y
k)
Semnificatia geometrica a celor mai simple formule de derivare numerica (pas echidistant).
Formulare Formule de der. Algoritmi Der.ord.sup Der.part. MDF-1 MDF-2 MDF-3
Algoritmi numerici pentru derivarea functiilor
Este important modul în care este cunoscuta functia.• Functia data tabelar (prin date)
1. se evalueaza derivatele în punctele din tabel2. se evalueaza derivata într-un punct oarecare, facând o
interpolare liniara pe portiuni a acestor valori.
• Functia data printr-un cod1. se alege pasul optim de derivare pentru o anumita formula
de derivare numerica;2. se evalueaza derivata cu formula de derivare numerica
într-un punct oarecare.
Formulare Formule de der. Algoritmi Der.ord.sup Der.part. MDF-1 MDF-2 MDF-3
Algoritmi numerici - functie data tabelar1. Se evalueaza derivatele în punctele din tabel2. Se evalueaza derivata într-un punct oarecare, facând o
interpolare liniara pe portiuni a acestor valori.procedur a pregatire_derivate(n, x, y, dy ); calculeaza derivatele numerice în noduri; declaratiiîntreg n ; numarul de puncte din tabel minus 1tablou real x [n], y [n]; tabelul de valori, indici de la 0tablou real dy [n] ; valorile derivatelor; alte declaratii· · ·; la capatul din stânga diferente progresive de ordinul 2h1 = x1 − x0h2 = x2 − x1dy0 = −(2h1 + h2)/(h1(h1 + h2))dy0 + (h1 + h2)/(h1h2)dy1 − h1/(h2(h1 + h2))dy2; la capatul din dreapta diferente regresive de ordinul 2h1 = xn−1 − xn−2h2 = xn − xn−1dyn = −h2/(h1(h1 + h2))dyn−2 − (h1 + h2)/(h1h2)dyn−1 + (h1 + 2h2)/(h2(h1 + h2))dyn; în nodurile interioare diferente centrate de ordinul 2pentru k = 1, n − 1
h1 = xk − xk−1h2 = xk+1 − xkdyk = −h2/(h1(h1 + h2))dyk−1 − (h1 − h2)/(h1h2)dyk + h1/(h2(h1 + h2))dyk+1
•
rez = interpolare_lpp(n, x, dy, xcrt)
Formulare Formule de der. Algoritmi Der.ord.sup Der.part. MDF-1 MDF-2 MDF-3
Algoritmi numerici - functie data prin cod• Pasul de discretizare trebuie ales de utilizator.
• Eroarea de trunchiere este cu atât mai mica cu cât pasul estemai mic. Dar, exista si o eroare de rotunjire, care este cu atâtmai mare cu cât pasul este mai mic.
10−10
10−9
10−8
10−7
10−6
10−5
10−10
10−9
10−8
10−7
10−6
10−5
Pasul de derivare h
Ero
area
abs
olut
a a
deriv
atei
pro
gres
ive
de o
rd 1
10−8
10−7
10−6
10−5
10−4
10−3
10−14
10−13
10−12
10−11
10−10
10−9
10−8
10−7
10−6
Pasul de derivare h
Ero
area
abs
olut
a a
deriv
atei
cen
trat
e de
ord
2
Teste numerice pentru sin′(π/4): Eroarea în functie de pasul de derivare pentru formula de derivare progresiva de
ordinul 1 (stânga) si formula de derivare centrata de ordinul 2 (dreapta).
Pas optim de derivare numeric a = pasul de la care încep sapredomine erorile de trunchiere.Atentie: Nu este pasul pentru care eroarea este minima!.
Formulare Formule de der. Algoritmi Der.ord.sup Der.part. MDF-1 MDF-2 MDF-3
Algoritmi numerici - functie data prin codEstimarea hoptim în cazul formulei de derivare progresiva deordinul 1.
y ′
k =yk+1 − yk
xk+1 − xk=
yk+1 − yk
h(25)
et = |y ′
k − f ′(xk )| ≤M2
2h, (26)
unde |f ′′(x)| ≤ M2. Notam marginea erorii de trunchiere:
at =M2
2h (27)
Pp. |eyk/yk | < eps, si eh = 0
|er | ≤|eyk+1 |+ |eyk |
h≤
eps(|yk+1|+ |yk |)
h(28)
Daca |f (x)| ≤ M0 atunci
|er | ≤2M0eps
hnot= ar (29)
Formulare Formule de der. Algoritmi Der.ord.sup Der.part. MDF-1 MDF-2 MDF-3
Algoritmi numerici - functie data prin codEroarea totala e = et + er va fi majorata de
|e| = |et + er | ≤ |et |+ |er | ≤ at + ar = m(h)
m(h) =M2
2h +
2M0epsh2 , (30)
min(m(h)) = m(hoptim). (31)
h
Mar
gine
a er
orii
Marginea erorii de trunchiere a
t(h)
Marginea erorii de rotunjire ar(h)
Marginea erorii totale a(h)
hoptim
Pasul optim de derivare numeric a este pasul de la care încep s a predomine erorile de trunchiere.
Formulare Formule de der. Algoritmi Der.ord.sup Der.part. MDF-1 MDF-2 MDF-3
Algoritmi numerici - functie data prin codConditia de minim
m′(h) = 0
rezulta valoarea optima a pasului de derivare
hoptim = 2
√
M0epsM2
. (32)
Pasul optim de derivare depinde de:
• eroarea de rotunjire
• marginea functiei
• marginea derivatei a doua.
El este pasul pentru care marginea erorii totale este minima sinu trebuie confundat cu pasul pentru care eroarea totala esteminima, pas care s-ar putea sa fie mai mic decât pasul optim,dar care nu poate fi estimat apriori.
Formulare Formule de der. Algoritmi Der.ord.sup Der.part. MDF-1 MDF-2 MDF-3
Algoritmi numerici - functie data prin codAlgoritm pentru derivata progresiva de ordinul 1 a unei functii dateprin cod, folosind un pas de derivare optim. M2 ≈ f (x)−2f (x+h)+f (x+2h)
h2 .
func¸tie derivata_f(x, h0, maxit); calculeaza derivata numerica a functiei f folosind diferente progresive de ordinul 1 si pas optim de derivarereal x ; punctul în care se va evalua derivatareal h0 ; pasul initialîntreg maxit ; numarul maxim de iteratii pentru calculul pasului; alte declaratii· · ·eps = zeroul_masinii ()f 0 = f (x) ; valoarea functiei în punctul de derivareh = h0k = 0repet a
k = k + 1f 1 = f (x + h)f 2 = f (x + 2h)M2 = |f 0 − 2f 1 + f 2|/h2 ; estimarea marginii derivatei a douadac a M2 < eps ; testeaza daca derivata a doua este zero
hoptim = h ; functia e liniara, eroarea de trunchiere e zeroaltfel
M0 = max(|f 0|, |f 1|, |f 2|) ; estimarea marginii functieihoptim = 2
√
M0eps/M2r = h/hoptim ; rata de modificare a pasuluih = hoptim
pân a când (k > maxit) sau (r > 0.5si r < 2)întoarce (f (x + hoptim) − f 0)/hoptim
Formulare Formule de der. Algoritmi Der.ord.sup Der.part. MDF-1 MDF-2 MDF-3
Derivate de ordin superior
Pot fi privite ca aplicatii recursive ale derivatei de ordinul unu.
pregatire_derivate(n, x , dy , d2y )
rez = interpolare_lpp(n, x , d2y , xcrt)
Experimentele numerice arata ca procedând astfel, acuratetearezultatului se pierde pe masura ce ordinul derivatei
Formulare Formule de der. Algoritmi Der.ord.sup Der.part. MDF-1 MDF-2 MDF-3
Derivate de ordin superiorO alta posibilitate: deducerea formulelor de derivare pornind dela polinomul de interpolare initial.
g(x) =(x − xk )(x − xk+1)
(xk−1 − xk )(xk−1 − xk+1)yk−1 +
(x − xk−1)(x − xk+1)
(xk − xk−1)(xk − xk+1)yk +
+(x − xk−1)(x − xk )
(xk+1 − xk−1)(xk+1 − xk )yk+1, (33)
g′(x) =
(x − xk ) + (x − xk+1)
(xk−1 − xk )(xk−1 − xk+1)yk−1 +
(x − xk−1) + (x − xk+1)
(xk − xk−1)(xk − xk+1)yk +
+(x − xk−1) + (x − xk )
(xk+1 − xk−1)(xk+1 − xk )yk+1. (34)
g′′(x) =
2
(xk−1 − xk )(xk−1 − xk+1)yk−1 +
2
(xk − xk−1)(xk − xk+1)yk +
+2
(xk+1 − xk−1)(xk+1 − xk )yk+1, (35)
În cazul particular xk+1 − xk = xk − xk−1 = h
g′′(x) =yk−1 − 2yk + yk+1
h2 . (36)
Formulare Formule de der. Algoritmi Der.ord.sup Der.part. MDF-1 MDF-2 MDF-3
Derivate de ordin superiorEvaluarea erorii - folosind seria Taylor cu un numar suplimentarde termeni.
f (xk+1) = f (xk ) + f ′(xk )h +f ′′(xk )
2!h2 +
f ′′′(xk )
3!h3 + O(h4), (37)
f (xk−1) = f (xk )− f ′(xk )h +f ′′(xk )
2!h2 −
f ′′′(xk )
3!h3 − O(h4). (38)
f (xk+1) + f (xk−1) = 2f (xk ) + f ′′(xk )h2 + O(h4). (39)
yk+1 − 2yk + yk−1
h2 − f ′′(xk ) = O(h2). (40)
• Polinomul de interpolare initial trebuie sa aiba un gradsuficient de mare pentru ca derivata de ordin superior safie nenula.
• Ar putea apare fenomenului Runge care va afecta siacuratetea derivatelor numerice. De aceea, calcululnumeric al derivatelor de ordin superior trebuie evitat.
Formulare Formule de der. Algoritmi Der.ord.sup Der.part. MDF-1 MDF-2 MDF-3
Derivate partiale
f (x , y) : [a, b]× [c, d ] → IR.
a = x0, x1, . . . , xi , . . . , xn = b
c = y0, y1, . . . , yj , . . . , ym = d
x
y
b=xn
a=x0
c=y0
d=ym
xi
yj (x
i,y
j)
Discretizarea domeniului de definitie.
Valorile din"tabel":
f (xi , yj) = zi,j
i = 0, n j = 0,m
Formulare Formule de der. Algoritmi Der.ord.sup Der.part. MDF-1 MDF-2 MDF-3
Derivate partialeEvaluarea numerica a derivatelor partiale se reduce laevaluarea numerica a derivatelor totale.
∂f∂x
∣
∣
∣
∣
x=x0,y=y0
= limx→x0
f (x , y0)− f (x0, y0)
x − x0= lim
x→x0
f1(x)− f1(x0)
x − x0=
df1dx
,
(41)unde f1(x) = f (x , y0) este o functie care depinde de o singura
variabila reala. Similar,
∂f∂y
∣
∣
∣
∣
x=x0,y=y0
= limy→y0
f (x0, y)− f (x0, y0)
y − y0= lim
y→y0
f2(y)− f2(y0)
y − y0=
df2dy
,
(42)unde f2(y) = f (x0, y) este o functie care depinde de o singura
variabila reala.Formulele de derivare progresiva de ordinul 1
∂g∂x
(xi , yj) =zi+1,j − zi,j
xi+1 − xi,
∂g∂y
(xi , yj) =zi,j+1 − zi,j
yj+1 − yj(43)
etc.
Formulare Formule de der. Algoritmi Der.ord.sup Der.part. MDF-1 MDF-2 MDF-3
Derivate partiale
Exemplu: evaluarea numerica a operatorului Laplace aplicatunei functii de doua varibile.
∆f =∂2f∂x2 +
∂2f∂y2 . (44)
Pp. xi+1 − xi = yj+1 − yj = h, pentru orice i = 0, n − 1,j = 0,m − 1
∆g =zi−1,j − 2zi,j + zi+1,j
h2 +zi,j−1 − 2zi,j + zi,j+1
h2 =
=zi−1,j + zi+1,j + zi,j−1 + zi,j+1 − 4zi,j
h2 , (45)
Formulare Formule de der. Algoritmi Der.ord.sup Der.part. MDF-1 MDF-2 MDF-3
Derivate partiale
yj
xi
zi+1,j
zi−1,j
zi,j−1
zi,j+1
zi,j
În discretizarea operatorului Laplace în punctul (xi , yj ) apar numai valorile
nodurilor marcate.
∆g =zi−1,j + zi+1,j + zi,j−1 + zi,j+1 − 4zi,j
h2(46)
Daca astfel de rationamente se folosesc pentru rezolvarea numericaa problemelor de inginerie formulate cu ajutorul unor ecuatii cuderivate totale (ODE) sau partiale (PDE) atunci se spune ca pentrurezolvarea problemei se aplica metoda diferentelor finite.
Formulare Formule de der. Algoritmi Der.ord.sup Der.part. MDF-1 MDF-2 MDF-3
Metoda diferentelor finite
• Pasul 1: Se alege o schema de diferente finite pentruaproximarea derivatelor din ecuatii si se rescrie ecuatia caecuatie cu diferente finite.
• Pasul 2: Se stabileste grila de discretizare si se scrieecuatia discretizata pentru fiecare nod.
• Pasul 3: Se rezolva sistemul de ecuatii pentrudeterminarea valorilor necunoscute în nodurile grilei.
• Pasul 4: Calculul altor valori, în afara celor din noduri, seface prin interpolare.
1 si 2 = preprocesare3 = rezolvare4 = postprocesare
Formulare Formule de der. Algoritmi Der.ord.sup Der.part. MDF-1 MDF-2 MDF-3
Exemplul 1 - rezolvarea unei ecuatii diferentialeordinare
Vom considera un condensator de capacitate C = 4µF, initialdescarcat, care se încarca de la o sursa de tensiune continuaE = 20 mV printr-un rezistor de rezistenta R = 10 Ω.Fie i curentul prin circuit si u tensiunea la bornelecondensatorului.
Ri(t) + u(t) = E . (47)
i(t) = Cdu(t)
dt, (48)
RCdu(t)
dt+ u(t) = E . (49)
u(0) = u0 = 0. (50)
Solutie analitica
u(t) = (u0 − E) exp(−t/τ) + E . (51)
unde τ = RC este constanta de timp a circuitului.
Formulare Formule de der. Algoritmi Der.ord.sup Der.part. MDF-1 MDF-2 MDF-3
Exemplul 1 - rezolvarea unei ecuatii diferentialeordinare
Ecuatia (49) o rescriem ca
du(t)dt
+1τ
u(t) =1τ
E . (52)
Vom urmari calculul numeric în intervalul de timp [0, tmax ] undetmax = 10τ într-o retea echidistanta de N puncte tk , unde pasulde discretizare h este
tk+1 − tk = h, pentru k = 1, . . . ,N − 1. (53)
Vom nota valorile discrete obtinute prin rezolvare numerica cuuk . Ele vor fi aproximatii ale marimii reale u.
uk ≈ u(tk ). (54)
Formulare Formule de der. Algoritmi Der.ord.sup Der.part. MDF-1 MDF-2 MDF-3
Exemplul 1 - rezolvarea unei ecuatii diferentialeordinare
du(t)
dt+
1
τu(t) =
1
τE. (55)
Varianta I - diferente finite progresive de ordinul 1
uk+1 − uk
h+
1τ
uk =1τ
E , (56)
uk+1 − uk +hτ
uk =hτ
E . (57)
⇒ uk+1 poate fi calculata explicit cu formula
uk+1 = uk
(
1 −hτ
)
+hτ
E . (58)
Formulare Formule de der. Algoritmi Der.ord.sup Der.part. MDF-1 MDF-2 MDF-3
Exemplul 1 - rezolvarea unei ecuatii diferentialeordinare
uk+1 = uk
(
1 −hτ
)
+hτ
E . (59)
procedur a mdf_ode_p1(u0,E,tau,h,N,u)
; rezolva ecuatia du(t)dt + 1
τu(t) = 1
τE cu metoda diferentelor finite
; d/dt se discretizeaza folosind diferente progresive de ordinul 1real u0 ; conditia initiala - datareal E ; coeficient în ecuatie - datareal tau ; constanta de timp - datareal h ; pas de discretizare al intervalului de timp - datîntreg N ; numar de valori de timp - dattablou real u[N] ; solutia discreta - rezultatu(1) = u0pentru k = 1,N-1
u(k+1) = u(k)*(1-h/tau) + h*E/tau•retur
Aceasta metoda, cunoscuta si sub numele de Euler explicita,este instabila pentru h > τ .
Formulare Formule de der. Algoritmi Der.ord.sup Der.part. MDF-1 MDF-2 MDF-3
Exemplul 1 - rezolvarea unei ecuatii diferentialeordinare
du(t)
dt+
1
τu(t) =
1
τE. (60)
Varianta a II-a - diferente finite regresive de ordinul 1
uk − uk−1
h+
1τ
uk =1τ
E , (61)
uk − uk−1 +hτ
uk =hτ
E . (62)
⇒ uk poate fi calculata explicit:
uk =
(
uk−1 +hτ
E)
/
(
1 +hτ
)
. (63)
Formulare Formule de der. Algoritmi Der.ord.sup Der.part. MDF-1 MDF-2 MDF-3
Exemplul 1 - rezolvarea unei ecuatii diferentialeordinare
uk =
(
uk−1 +hτ
E)
/
(
1 +hτ
)
. (64)
procedur a mdf_ode_r1(u0,E,tau,h,N,u)
; rezolva ecuatia du(t)dt + 1
τu(t) = 1
τE cu metoda diferentelor finite
; d/dt se discretizeaza folosind diferente regresive de ordinul 1real u0 ; conditia initiala - datareal E ; coeficient în ecuatie - datareal tau ; constanta de timp - datareal h ; pas de discretizare al intervalului de timp - datîntreg N ; numar de valori de timp - dattablou real u[N] ; solutia discreta - rezultatu(1) = u0pentru k = 2,N
u(k) = (h*E/tau + u(k-1))/(1 + h/tau)• retur
Metoda (Euler implicita pentru rezolvarea ODE) nu mai este instabilapentru h > τ .Ecuatia generala este dy/dt = f (t, y) unde f este o functie în general neliniara. Folosirea derivatei regresive
pentru f conduce la o relatie implicita pentru calculul lui uk , relatie ce reprezinta o ecuatie neliniara. În cazul
particular studiat, f este liniara si de aceea solutia se poate calcula explicit.
Formulare Formule de der. Algoritmi Der.ord.sup Der.part. MDF-1 MDF-2 MDF-3
Exemplul 1 - rezolvarea unei ecuatii diferentialeordinare
du(t)
dt+
1
τu(t) =
1
τE. (65)
Varianta a III-a - diferente finite centrate de ordinul 2
uk+1 − uk−1
2h+
1τ
uk =1τ
E , (66)
− uk−1 +2hτ
uk + uk+1 =2hτ
E . (67)
k = 2, . . . ,N − 1 ⇒ N − 2 ecuatii cu N necunoscute.Trebuie adaugate relatiile la capete. La t = 0: u(1) = u0
Pentru k = N diferente regresive de ordinul 2:
uN−2 − 4uN−1 + 3uN
2h+
1τ
uN =1τ
E , (68)
uN−2 − 4uN−1 +
(
3 +2hτ
)
uN =2hτ
E . (69)
Formulare Formule de der. Algoritmi Der.ord.sup Der.part. MDF-1 MDF-2 MDF-3
Exemplul 1 - rezolvarea unei ecuatii diferentialeordinare
procedur a mdf_ode_c2(u0,E,tau,h,N,u)
; rezolva ecuatia du(t)dt + 1
τu(t) = 1
τE cu metoda diferentelor finite
; d/dt se discretizeaza folosind diferente centrate de ordinul 2real u0 ; conditia initiala - datareal E ; coeficient în ecuatie - datareal tau ; constanta de timp - datareal h ; pas de discretizare al intervalului de timp - datîntreg N ; numar de valori de timp - dattablou real u[N] ; solutia discreta - rezultattablou real A[N,N] ; matricea coeficientilor sistemului asamblat - stocata rartablou real tl[N] ; vectorul termenilor liberi ai sistemului asamblatA = 0 ; pseudocod simplificattl = 0A(1,1) = 1 ; conditia initiala tl(1) = u0pentru k = 2,N-1 ; noduri interioare
A(k,k-1) = -1A(k,k) = 2*h/tauA(k,k+1) = 1tl(k) = 2*h*E/tau
•A(N,N-2) = 1 ; conditia la tmax = N*hA(N,N-1) = -4A(N,N) = (3 + 2*h/tau)tl(N) = 2*h*E/tauu = ”A−1tl” ; rezolva sistemul algebric liniar asamblatretur
Formulare Formule de der. Algoritmi Der.ord.sup Der.part. MDF-1 MDF-2 MDF-3
Exemplul 1 - rezolvarea unei ecuatii diferentialeordinare
0 1 2 3 4
x 10−4
0
0.002
0.004
0.006
0.008
0.01
0.012
0.014
0.016
0.018
0.02
t [s]
u [V
]
h/tau = 0.500000
analiticprog ord 1regr ord 1centr ord 2
0 1 2 3 4
x 10−4
0
0.5
1
1.5
2
2.5x 10
−3
t [s]E
roar
e ab
solu
ta [V
]
h/tau = 0.500000
prog ord 1regr ord 1centr ord 2
Cazul h = τ/2.
Formulare Formule de der. Algoritmi Der.ord.sup Der.part. MDF-1 MDF-2 MDF-3
Exemplul 1 - rezolvarea unei ecuatii diferentialeordinare
0 1 2 3 4
x 10−4
0
0.002
0.004
0.006
0.008
0.01
0.012
0.014
0.016
0.018
0.02
t [s]
u [V
]
h/tau = 1.000000
analiticprog ord 1regr ord 1centr ord 2
0 1 2 3 4
x 10−4
0
1
2
3
4
5
6
7
8x 10
−3
t [s]E
roar
e ab
solu
ta [V
]
h/tau = 1.000000
prog ord 1regr ord 1centr ord 2
Cazul h = τ .
Formulare Formule de der. Algoritmi Der.ord.sup Der.part. MDF-1 MDF-2 MDF-3
Exemplul 1 - rezolvarea unei ecuatii diferentialeordinare
0 1 2 3 4
x 10−4
0
0.005
0.01
0.015
0.02
0.025
0.03
0.035
0.04
t [s]
u [V
]
h/tau = 2.000000
analiticprog ord 1regr ord 1centr ord 2
0 1 2 3 4
x 10−4
0
0.005
0.01
0.015
0.02
0.025
t [s]
Ero
are
abso
luta
[V]
h/tau = 2.000000
prog ord 1regr ord 1centr ord 2
Cazul h = 2τ .
În cazul utilizarii diferentelor centrate de ordinul 2, metoda nu esteinstabila si este mai precisa decât implementarile anterioare. Acestlucru era de asteptat deoarece formulele de derivare numerica deordinul 2 sunt mai precise decât cele de ordinul 1. Crestereaacuratetii se face însa pe seama cresterii complexitatii algoritmului.
Formulare Formule de der. Algoritmi Der.ord.sup Der.part. MDF-1 MDF-2 MDF-3
Exemplul 2 - rezolvarea unei ecuatii derivate partialede tip eliptic
O
y
xa/2
b/2
b
a
σ=0
σ0 V
1 V
∆V = 0
O
RQ
N
S
P
dV/dn=0
V=1
dV/dn=0
dV/dn=0
V=0
dV/dn=0
Problema 2D de regim electrocinetic: domeniul de calcul (stânga), conditii de frontiera (dreapta).
Se da un conductor omogen, de conductivitate σ, situat într-un mediu perfect
izolant. Conductorul are o dimensiune dupa axa Oz mult mai lunga decât
celelalte doua. Figura reprezinta o sectiune perpendiculara pe directia
dimensiunii foarte mari. Aceasta sectiune este un dreptunghi, de dimensiuni
a, b. Conductorul are doua borne supraconductoare, una aflata la potential
V0 = 1V, iar cealalta aflata la potential nul. Se cere câmpul în acest domeniu.
Formulare Formule de der. Algoritmi Der.ord.sup Der.part. MDF-1 MDF-2 MDF-3
Exemplul 2 - rezolvarea unei ecuatii derivate partialede tip eliptic
Formularea în V - potential electrocinetic, unde E = −gradV , Eeste intensitatea câmpului electric.Ecuatia de ordinul doi satisfacuta de V :
− div (σ gradV ) = 0 (70)
V : D → IR; D = ONPQ; (70) - ecuatie eliptica, de tip Laplacegeneralizata. Domeniul omogen ⇒ ecuatie de tip Laplace:
∆V = 0, (71)
∆ = div grad este operatorul Laplace; în 2D, xOy:
∆V =∂2V∂x2 +
∂2V∂y2 . (72)
⇒ Ecuatia de rezolvat, unde V = V (x , y) : D → IR:
∂2V∂x2 +
∂2V∂y2 = 0, (73)
Formulare Formule de der. Algoritmi Der.ord.sup Der.part. MDF-1 MDF-2 MDF-3
Exemplul 2 - rezolvarea unei ecuatii derivate partialede tip eliptic
Conditii de frontiera:
V = V0 pe QR (Dirichlet) (74)
V = 0 pe SN (Dirichlet) (75)∂V∂n
= 0 pe NOQ∪ RPS (Neumann) (76)
Gridul de discretizare bidimensional:
0 = x1 < x2 < . . . < xnx = a
0 = y1 < y2 < . . . < yny = b.
Un nod al gridului:
(xi , yj) i = 1, . . . , nx , j = 1, . . . , ny .
MDF determina valorile potentialului în nodurile acestui grid:
V (xi , yj)not= Vi,j i = 1, nx , j = 1, ny . (77)
Formulare Formule de der. Algoritmi Der.ord.sup Der.part. MDF-1 MDF-2 MDF-3
Exemplul 2 - rezolvarea unei ecuatii derivate partialede tip eliptic
Ecuatia asociata unui nod interior: Pp hx , hy pasii dediscretizare:
∂2V∂x2 (xi , yj) =
Vi+1,j − 2Vi,j + Vi−1,j
h2x
, (78)
∂2V∂y2 (xi , yj) =
Vi,j+1 − 2Vi,j + Vi,j−1
h2y
. (79)
Vi+1,j − 2Vi,j + Vi−1,j
h2x
+Vi,j+1 − 2Vi,j + Vi,j−1
h2y
= 0, (80)
2
(
1h2
x+
1h2
y
)
Vi,j −1h2
xVi+1,j −
1h2
xVi−1,j −
1h2
yVi,j+1−
1h2
yVi,j−1 = 0.
(81)Potentialul într-un nod interior este o combinatie liniara a
potentialelor nodurilor învecinate.
Formulare Formule de der. Algoritmi Der.ord.sup Der.part. MDF-1 MDF-2 MDF-3
Exemplul 2 - rezolvarea unei ecuatii derivate partialede tip eliptic
Ecuatia asociata unui nod interior - retea neuniforma:
O
A
B
C
D
y
x
Nodurile marcate intervin în scrierea ecuatiei unui
nod interior.
VO
(
1hBhD
+1
hAhC
)
−
−VA1
hA(hA + hC)−
−VB1
hB(hB + hD)−
−VC1
hC(hA + hC)−
−VD1
hD(hB + hD)= 0, (82)
hA = ‖OA‖, hB = ‖OB‖, hC =
‖OC‖, hD = ‖OD‖.
Formulare Formule de der. Algoritmi Der.ord.sup Der.part. MDF-1 MDF-2 MDF-3
Exemplul 2 - rezolvarea unei ecuatii derivate partialede tip eliptic
Ecuatia asociata unui nod pe frontiera dreapta:
O
A
B
C
x
n
D1
Nodurile marcate intervin în scrierea ecuatiilor unui
nod pe frontiera Neumann dreapta.
Fie g = −σ ∂V∂n si gO valoarea
medie în O:
gO =gOAhA + gOChC
hA + hC(83)
unde gOA este CF asociata luiOA, iar gOC este CF asociata luiOC.Nod “fantoma" D1, hD = ‖OD1‖.Pentru simplitate hD = hB.
∂V∂x
=gσ. (84)
⇒VD1 = VB − 2hB
gO
σ. (85)
Formulare Formule de der. Algoritmi Der.ord.sup Der.part. MDF-1 MDF-2 MDF-3
Exemplul 2 - rezolvarea unei ecuatii derivate partialede tip eliptic
Ecuatia asociata unui nod pe frontiera dreapta:
O
A
B
C
x
n
D1
Nodurile marcate intervin în scrierea ecuatiilor unui
nod pe frontiera Neumann dreapta.
Pentru O se scrie o ecuatieca pentru un nod interior si în-locuind expresia (85) ⇒:
VO
(
1h2
B
+1
hAhC
)
−
−VA1
hA(hA + hC)−
−VB1h2
B
−
−VC1
hC(hA + hC)= −
gO
σhB.(86)
Formulare Formule de der. Algoritmi Der.ord.sup Der.part. MDF-1 MDF-2 MDF-3
Exemplul 2 - rezolvarea unei ecuatii derivate partialede tip eliptic
• Discretizarea problemei conduce la un sistem de ecuatii algebricliniar, prin a carui rezolvare se obtin valorile potentialelor înnodurile gridului.
• Pentru asamblarea sistemului este util ca nodurile sa fienumerotate a.î. necunoscutele sa reprezinte un vector.
2
ny 2 ny nx ny
k
(i−1)ny+1
j
(nx−1)ny+11 ny+1
Exemplu de numerotare a nodurilor. Relatia între numarul nodului k si pozitiile proiectiilor sale pe cele doua axe i si
j este k = (i − 1)ny + j, i = 1, . . . , nx , j = 1, . . . , ny
Formulare Formule de der. Algoritmi Der.ord.sup Der.part. MDF-1 MDF-2 MDF-3
Exemplul 2 - rezolvarea unei ecuatii derivate partialede tip eliptic
0.1
0.1
0.2
0.2
0.3
0.3
0.3
0.4
0.4
0.4
0.5
0.5
0.5
0.5
0.6
0.6
0.6
0.6
0.7
0.7
0.7
0.80.8
0.8
0.90.9
1 1
0 0.5 1 1.5 20
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2
0.1
0.1
0.2
0.2
0.3
0.3
0.4
0.4
0.4
0.5
0.5
0.5
0.5
0.6
0.6
0.6
0.7
0.7
0.7
0.80.8
0.90.9
1 1
0 0.5 1 1.5 20
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2
Linii echipotentiale pentru un grid cu nx = ny = 3 (stânga) si un grid cu nx = ny = 21 (dreapta).
Formulare Formule de der. Algoritmi Der.ord.sup Der.part. MDF-1 MDF-2 MDF-3
Exemplul 3 - rezolvarea unei ecuatii cu derivatepartiale de tip hiperbolic
Ecuatia undei scalare:
∂2u∂t2 = v2∂
2u∂x2 , (87)
u(x , t) : [0, a]× [0,T ] → IR, v este o constanta cunoscuta.Solutia acestei ecuatii este o unda care se propaga cu viteza v .Buna formulare a acestei probleme necesita impunerea
• conditiei initiale u(x , 0) = h0(x),
• conditiilor la capete - relatii în care intervin marimileu(0, t) = h1(t) si u(a, t) = h2(t).
Formulare Formule de der. Algoritmi Der.ord.sup Der.part. MDF-1 MDF-2 MDF-3
Exemplul 3 - rezolvarea unei ecuatii cu derivatepartiale de tip hiperbolic
Se poate demonstra ca
u(x , t) = ud(x , t) + ui(x , t) + U0, (88)
unde ud(x , t) se numeste unda directa si se poate scrie subforma
ud(x , t) = f (x − vt), (89)
iar ui(x , t) se numeste unda inversa si se poate scrie sub forma
ud(x , t) = g(x + vt). (90)
f si g rezulta în mod univoc, din impunerea conditiilor initiale siconditiilor la capete.f (x − vt) - o anumita valoare a acestei functii, într-un anumit punct x si într-un anumit moment de timp t se va regasi
dupa un interval de timp ∆t în punctul x + ∆x , unde ∆x = v∆t (se propaga deci în sensul pozitiv al axei Ox).
Formulare Formule de der. Algoritmi Der.ord.sup Der.part. MDF-1 MDF-2 MDF-3
Exemplul 3 - rezolvarea unei ecuatii cu derivatepartiale de tip hiperbolic
Avem nevoie de o discretizare spatiala si de una temporala.Pp. ca ambele discretizari sunt uniforme si vom nota cu
• ∆z pasul discretizarii spatiale• ∆t pasul discretizarii temporale• N numarul de puncte de discretizare spatiale• M numarul de pasi de timp simulati ⇒ timpul maxim de
simulare este T = M∆tVom considera urmatoarele conditii de unicitate:
• conditia initiala nula u(x , 0) = h0(x) = 0,• conditia la capatul din stânga - excitatia cu un impuls
Gauss: u(0, t) = h1(t) = exp(−(t − T/10).2/(2 ∗ (T/50)2))
• conditia la capatul din dreapta u(a, t) = h2(t) = 0.În problemele în care apare propagare, este util uneori ca frontiera domeniului spatial sa fie modelata ca o frontiera
absorbanta, "invizibila" din punct de vedere al propagarii marimilor în domeniul spatial modelat. Conditia de frontiera
absorbanta pentru capatul din dreapta al problemei studiate este ∂u∂t + v ∂u
∂z = 0.
Formulare Formule de der. Algoritmi Der.ord.sup Der.part. MDF-1 MDF-2 MDF-3
Exemplul 3 - rezolvarea unei ecuatii cu derivatepartiale de tip hiperbolic
Marimea u(x , t) este discretizata atât în spatiu cât si în timp.Prin rezolvarea cu MDF, vom obtine aproximatii:
u(xk , tj) ≈ u(j)k , k = 1, . . . ,N, j = 1, . . . ,M. (91)
Formule de derivare de ordin 2:
u(j−1)k − 2u(j)
k + u(j+1)k
(∆t)2 = v2 u(j)k−1 − 2u(j)
k + u(j)k+1
(∆x)2 , (92)
⇒
u(j+1)k =
(
v∆t∆x
)2(
u(j)k−1 − 2u(j)
k + u(j)k+1
)
+2u(j)k −uj−1
k , k = 2, . . . ,N−
(93)de unde rezulta ca marimea u, într-un anumit punct k si la un
anumit moment de timp j + 1 depinde explicit de valoarea înacel punct si în cele învecinate la momentul de timp anterior j side valoarea în acel punct la momentul j − 1:
Formulare Formule de der. Algoritmi Der.ord.sup Der.part. MDF-1 MDF-2 MDF-3
Exemplul 3 - rezolvarea unei ecuatii cu derivatepartiale de tip hiperbolic
0 0.2 0.4 0.6 0.8 1
−2
−1.5
−1
−0.5
0
0.5
1
1.5
2
tmax = 1.00e+000, not = 1000, idx
t = 400
0 0.2 0.4 0.6 0.8 1
−2
−1.5
−1
−0.5
0
0.5
1
1.5
2
tmax = 1.00e+000, not = 1000, idx
t = 600
Propagarea unui impuls Gauss (stânga) si rezultatul reflexiei dupa atingerea unei frontiere pe care s-a impus
conditie Dirichlet nula (dreapta).
Atunci când o problema necesita atât o discretizare spatiala câtsi una temporala, pentru ca solutia numerica sa fie stabila, estenecesar sa fie îndeplinita conditia lui Courant
|v |∆t ≤ ∆x . (94)
Formulare Formule de der. Algoritmi Der.ord.sup Der.part. MDF-1 MDF-2 MDF-3
Concluzii - MDF
1. În rezolvarea ecuatiilor diferentiale cu metoda diferentelorfinite, trebuie avute în vedere aspectul stabilitatiialgoritmului si aspectul instabilitatii numerice datorateintrarii în zona în care predomina erorile de rotunjire.
2. Stabilitatea se analizeaza diferit daca ecuatia este cuderivate ordinare (ODE) sau cu derivate partiale (PDE).
3. Formulat pe scurt, trebuie ca pasii de discretizare trebuiesa fie suficient de mici pentru ca algoritmul sa fie stabil, darnu prea mici ca sa nu se intre în zona erorilor de rotunjire.
4. Un algoritm eficient ar trebui sa poata sa foloseasca paside discretizare neuniformi, adaptati tipului de variatie asolutiei.