Ecuatii diferentiale

2.49.9 Rezolvarea numerica a EDP, pg. 551

Deşi există un mare număr de probleme în care apar EDP pentru care poate fi găsită o soluţie exactă (de obicei sub forma unei serii infinite), există totuşi şi probleme cu EDP pentru care nu poate fi găsită o soluţie exactă. Pentru aceasta se folosesc metode numerice care găsesc o soluţie aproximativă. Metoda de rezolvare depinde de tipul EDP: eliptică, parabolică sau hiperbolică.

Ca şi în cazul ecuaţiilor diferenţiale ordinare (EDO) care trebuie rezolvate prin metode numerice, scopul este de a obţine valori aproximative ale soluţiei intr-un set convenabil de puncte. De obicei aceasta înseamnă aflarea valorilor soluţiei în punctele unui caroiaj rectangular care acceptă o porţiune a domeniului soluţiei. În figură sunt prezentate câteva posibilităţi:

1. 
   caroiaj rectangular ce acoperă o regiune rectangulară în planul xy pentru care putem încerca să rezolvăm ecuaţia Laplace.
   
2. 
   caroiaj rectangular suprapus peste o regiune oarecare
   
3. 
   caroiaj rectangular suprapus peste o regiune rectangulară infinită a planului xt reprezentând domeniul pe care încercam să rezolvam de pildă ecuaţia unidimensională a căldurii
   
4. 
   caroiaj rectangular ce acoperă o regiune rectangulară în spaţiu pentru care putem încerca să rezolvăm ecuaţia Laplace tridimensională.
   

În figurile a, c şi d punctele extreme ale caroiajului se afla pe marginea regiunii, acolo unde valorile soluţiei sunt date ca mărimi de intrare ale problemei.

Pentru regiunea rectangulară din figura b lucrurile nu stau la fel. Formulele generale pe care le vom obţine trebuie modificate pentru acele puncte ale caroiajului care sunt vecine marginii dar nu se află chiar pe aceasta.

Ideea fundamentală pe care se bazează rezolvarea numerică a EDP este următoarea: fiecare dintre derivatele parţiale ce apar în ecuaţie este înlocuită printr-o aproximare cu diferenţe finite. Când aceste diferenţe sunt evaluate în fiecare din punctele caroiajului se obţine un sistem de ecuaţii care poate fi rezolvat direct sau prin metode iterative.

Mai exact, într-un caroiaj plan coordonatele punctelor sunt p_i = p_o + ih şi q_j = q₀ + jk, atunci din definiţia derivatei avem:

(2.90a)

(2.90b)

mai mult dacă se derivează de două ori formula de interpolare a lui Stirling (EC. 18,Sec 5.2))

se pune r>1 şi se neglijează toate diferenţele cu gradul >2, pentru p_i = x rezulta:

(2.91a)

(2.91b)

Aceste observaţii vor fi folosite în continuare pentru rezolvarea ecuaţiilor eliptice, parabolice şi hiperbolice.

2.4.1Cazul 1. Ecuaţii eliptice. Ecuaţii Laplace în coordonate bidimensionale

Dacă se folosesc ecuaţiile 2 pentru a aproxima fiecare dintre derivatele parţiale din ecuaţiile Laplace în coordonate bidimensionale adică se obţine ca o ecuaţie cu diferenţe finite aproximând ecuaţia iniţiala relaţia;

Dacă se presupune că h = k se scoate f_i,j sub forma:

(3) (2.92)

Valoarea lui f în orice punct din caroiaj este egală cu media aritmetică a valorilor lui f din cele 4 puncte alăturate. Configuraţia din figură se numeşte "stea". Dacă relaţia 3 se calculează pentru fiecare punct al caroiajului care nu este pe contur pentru care valorile iniţiale ale lui f sunt cunoscute, se obţine un sistem simultan de ecuaţii liniare cu necunoscutele f_ij. Numărul de ecuaţii este egal cu numărul de puncte din reţea pentru care trebuie calculată valoarea lui f; cel puţin pentru regiunile rectangulare se poate arăta că acest sistem de ecuaţii are întotdeauna o soluţie unică nebanală. De obicei numărul punctelor şi deci şi cel al ecuaţiilor este cuprins intre câteva sute şi câteva mii; cu toate acestea fiecare ecuaţie este relativ simplă deoarece nu conţine mai mult de 5 necunoscute. Pentru a ilustra formularea şi rezolvarea unui astfel de sistem să determinăm distribuţia staţionară a temperaturii într-o regiune pătrată ca cea din figură folosind un caroiaj obţinut prin împărţirea fiecărei laturi în 4 parţi egale. Necunoscutele problemei sunt temperaturile în cele 9 puncte ale caroiajului care nu sunt pe contur şi în care temperatura nu este indicată prin condiţii la limită. De la început se observă că datorită simetriei problemei f₁₁ = f₃₁ , f₁₂ = f₃₂ şi f₁₃ = f ₃₃ astfel încât problema presupune de fapt rezolvarea a numai şase ecuaţii cu şase necunoscute: f_11, f_12, f_13, f_21, f_22, f_23. dacă se aplică relaţia 3 în fiecare din cele 6 puncte ale reţelei P_11, P_12, P_13, P_21, P_22, P₂₃ şi ţinând cont de simetrie şi de valorile pe contur se obţine pentru P₁₁:

4f₁₁- f₀₁- f₁₀- f₂₁- f₁₂ = 0 dar f₀₁ = f₁₀ = 0 rezultă

4f₁₁- f₂₁- f₁₂ = 0 (4) (2.93)

Similar pentru punctele P₁₂, P₁₃, P₂₁, P₂₂, P₂₃ se obţine respectiv:

4f₁₂- f₁₁- f₂₂ – f₁₃ = 0 (5) (2.94)

4f₁₃- f₁₂- f₂₃ = 3/16 (6) (2.95)

4f₂₁- 2f₁₁- f₂₂ = 0 (7) (2.96)

4f₂₂- f₂₁- 2f₁₂ – f₂₃ = 0 (8) (2.97)

4f₂₃- f₂₂- 2f₁₃ = 1/4 (9) (2.98)

Dacă se folosesc ecuaţiile 4, 5, 6, pentru a elimina f₂₁, f₂₂, f₂₃ din ecuaţiile 7, 8, 9 se obţine sistemul :

15f₁₁- 8f₁₂+ f₁₃ = 0

-8f₁₁+16f₁₂-8f₁₃ = -3/16

f₁₁- 8f₁₂+15 f₁₃ = 1

Rezultă apoi că

f₁₁ = 0.0151 = f₃₁ ; f₂₁ = 0.0212

f₁₂ = 0.0391 = f₃₂ ; f₂₂ = 0.0547

f₁₃ = 0.0865 = f₃₃ ; f₂₃ = 0.1194

Valorile exacte determinate prin soluţia serie obţinuta prin metoda separării variabilelor sunt:

f₁₁ = f₃₁ = 0.0137; f₁₂ = f₃₂ = 0.0364; f₁₃ = f₃₃ = 0.0833

f₂₁ = 0.0194; f₂₂ = 0.0153; f₂₃ = 0.1159

Prin această metodă pot fi rezolvate şi probleme în care de-a lungul unei porţiuni sau a întregului contur este folosită derivata în locul funcţiei însăşi.

Să presupunem că de-a lungul conturului AB al regiunii din figură, valoarea derivatei este de  ori valoarea funcţiei aproximând derivată printr-o diferenţă, condiţia la limită = f devine = f

Astfel vom avea:sau. Dacă se aplică acum ecuaţia 3 de exemplu în P₃₃ rezultă:

4f₃₃- f₄₃- f₃₄ – f₂₃ – f₃₂ = 0 Û 4f₃₃- f₃₃/(1- lh)– 3/16 – f₂₃ –f₃₂ = 0 Þ

[(3-4lh)/( 1-lh)]^. f₃₃- f₂₃- f₃₂ = 3/16

De-a lungul unei margini/laturi izolate , l = 0 şi în particular ultima relaţie devine:

3f₃₃- f₂₃- f₃₂ = 3/16

În orice situaţie sistemul de ecuaţii obţinut astfel poate fi rezolvat la fel ca în cazul precedent.

Pentru a obţine relaţia corespunzătoare într-un punct al caroiajului apropiat de o margine neregulată, este convenabilă folosirea aproximaţiilor cu diferenţiale finite pentru derivate. Astfel să presupunem că în figură, punctele A şi B de pe contur se găsesc la o distanţă q_Ah şi q_Bh , unde q_B şi q_A sunt fiecare < 1. Folosind diferenţialele divizate de ordinul 2 pe direcţiile x şi y ca aproximări pentru avem conform Ex. 26 sect 5.1 pagina 255, în care s-a demonstrat că:

, ca:

şi

_sau

(10) (2.99)

Aceasta este ecuaţia care va fi folosită în punctele caroiajului ale căror vecini imediaţi cad în afara limitelor regiunii.

Există şi altă metodă pentru determinarea prin aproximaţie cu diferenţe finite a Laplacianului şi poate fi folosita pentru determinarea valorii soluţiei în punctele caroiajului, acest lucru se face printr-o metoda iterativă ce decurge astfel: să ne reamintim că valoarea soluţiei într-un punct este media aritmetică a celor 4 valori învecinate. Astfel, după o estimare iniţială valorile soluţiei pentru fiecare punct al caroiajului, ele pot fi corectate şi îmbunătăţite prin deplasarea sistematică prin caroiaj şi înlocuirea fiecărei valori conform ecuaţiei 3. Când facem aceasta, valoarea corectată trebuie folosită imediat în calculele ulterioare. Desigur că în zonele cu margini neregulate trebuie folosită ecuaţia 10 pentru a corecta valorile soluţiei în acele puncte ale caroiajului ai căror vecini sunt în afara limitei. Ca o ilustrare a acestei metode, să reconsiderăm problema precedentă începând cu estimările prezentate în figura 9.24a, se obţine la o prima rafinare a valorii f₁₃:

= 0.0919

Continuând deplasarea prin caroiaj după traseul indicat sau după oricare alt traseu folosind valorile corectate imediat de ele sunt disponibile (dar fără a mai ţine cont de simetria problemei) se obţin valorile din figura 9.24b.Valorile de deasupra reprezintă rezultatele după prima iteraţie iar cele de dedesubt reprezintă rezultatele după 5 iteraţii.

2.4.2Cazul 2. Ecuaţii parabolice. Transfer de căldură unidimensional

Pentru ecuaţia unidimensională a căldurii , regiunea planului xt pentru care se caută soluţia este întotdeauna infinită datorită creşterii infinite a timpului. Un caroiaj tipic pentru această situaţie a fost prezentat în figura 2.90c. ca o aproximaţie cu diferenţe finite a ecuaţiei căldurii folosim ecuaţia (1b) şi (2a), avem:

în mod clar se observă că ar fi convenabil să alegem h şi k astfel încât m = 1/2.

Valorile lui f pe contur sunt date iniţiale ale problemei; altfel încât condiţia iniţială impusă f(x,0) ne da valorile f₀₀, f₁₀, f₂₀,..

Similar, condiţia finală de forma f(0,t) = g₁(t) , f(l,t) = g₂(t), unde g₁ şi g₂ sunt de obicei (dar nu obligatoriu) constante, furnizează valorile f₀₁, f₀₂, f₀₃, … şi f_l1, f_l2, f_l3,. Capetele izolate pot fi ele abordate aşa cum reiese din discuţia precedenta despre ecuaţiile Laplace.

Odată ce au fost determinate din condiţiile iniţiale valorile lui f în punctele matricei de pe contur aflarea soluţiei în restul matricei decurge într-un mod clar utilizând modelul de extrapolare din ecuaţia 11 care este prezentată în figura 9.25. Prima oară se calculează valorile f₁₁, f₂₁, …,f_l-1,1 folosind valorile cunoscute f₀₀, f₁₀, f₂₀, …f_l0. Apoi folosind aceste valori şi condiţiile la limită f₀₁ şi f_l1, vom merge mai departe prin calcularea valorilor lui f în punctele celui de-al treilea rând şi aşa mai departe.