El mètode de les línies per a la resolució numèrica d'equacions en derivades parcials. The method of lines for numerical solutions of partial differential equations
C. Dalfó M. A. Fiol












El mètode de les línies per a
la resolució numèrica d'equacions en derivades parcials

The method of lines for numerical solutions of partial differential equations
Aquesta recerca té el suport del projecte 2017SGR1087 de l'Agència de Gestió d'Ajuts Universitaris i de Recerca (AGAUR) del Govern de Catalunya. 
    C. Dalfó

	Departament  de Matemàtica, Universitat de Lleida, Igualada (Barcelona), Catalunya

	

	<cristina.dalfo@matematica.udl.cat>

		M. A. Fiol

	Deptartament de Matemàtiques, Barcelona Graduate School of Mathematics,

		Universitat Politècnica de Catalunya, Barcelona, Catalunya

		<miguel.angel.fiol@upc.edu>

    
================================================================================================================================================================================================================================================================================================================================



[l]0.3 
    < g r a p h i c s >
    [l][1cm]0.79
		The first author has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 734922.



En aquest treball es descriu un mètode numèric semi-discret per a la resolució d'un tipus d'equacions en derivades parcials. És conegut com el mètode de les línies (MOL, per les sigles en anglès),
i es basa en la discretització de totes les  variables involucrades, excepte una. Ił.lustrem l'aplicació del MOL amb la resolució de l'equació de Laplace en coordenades cartesianes. Els conceptes matemàtics que hi apareixen es comparen amb els corresponents al mètode analític de separació de variables. A més, mostrem que els resultats obtinguts amb el MOL  són aproximacions molt bones de les solucions analítiques.

.5cm

In this paper, we describe a semi-discrete method for a numerical resolution of a type of partial differential equations, called the method of lines (MOL). This method is based on the discretization of all but one of the variables of the problem. We illustrate this method by solving the Laplace equation in Cartesian coordinates. We compare the concepts used by the MOL with respect to the analytical method of variable separation. We show that the results obtained with the MOL are very good approximations of the analytical solutions.


MSC2010: 65M20. 

Paraules clau: Equacions en derivades parcials, discretització d'una variable contínua, anàlisi numèric.

Keywords: Partial differential equations, discretization of a continuous variable, numerical analysis.





§ INTRODUCCIÓ


El mètode de les línies, també anomenat mètode de semi-discretització per Zuazua <cit.> és un mètode per a la resolució d'equacions en derivades parcials (EDP's) amb condicions de contorn i en regions amb una certa simetria. Es basa en discretitzar totes les variables, excepte una (si el problema ho permet, sovint es discretitzen les variables espacials i es manté contínua la variable temporal). Això condueix a un sistema d'equacions diferencials ordinàries (EDO's), les quals es poden resoldre pels mètodes numèrics habituals per a EDO's amb condicions inicials. Aquest mètode va ser introduït a principis dels anys 60 i, des de llavors, la seva exactitud i estabilitat han fet que es generessin una bona quantitat d'articles (vegeu, per exemple, Zafarullah <cit.> i Verwer i Sanz-Serna <cit.>).

Per tant, el MOL és un enfocament diferencial i de diferències finites per resoldre equacions en derivades parcials. Aquest mètode té moltes aplicacions en la resolució de problemes en Física i Enginyeria. Per exemple, aquest mètode s'ha usat àmpliament per experts en tècniques computacionals en problemes d'electromagnetisme, vegeu, per exemple, Berardia i Vurro <cit.>,  Diestel <cit.>, Sadiku i Obiozor <cit.>, o  Shakeri i Dehghan <cit.>.

En aquest treball, il·lustrem el mètode de les línies amb l'equació de Laplace en el pla, per simplicitat. El procés a seguir consta dels punts següents:


  * Discretització de l'equació diferencial en una direcció: Contracció de la regió de solucions a un conjunt de línies, i aproximació de les derivades corresponents.

  * Obtenció d'un sistema d'equacions diferencials ordinàries acoblades.

  * Transformació per obtenir un
sistema d'equacions diferencials ordinàries desacoblades.

  * Resolució de les equacions diferencials.

  * Transformació inversa i introducció de les condicions de contorn.

  * Solució de les equacions.


A més de les seves múltiples aplicacions, per exemple, com ja s'ha dit, en electromagnetisme, aquest mètode involucra moltes idees, conceptes i tècniques de la matemàtica discreta, juntament amb els mètodes bàsics de l'anàlisi diferencial. De fet, en el desenvolupament del mètode ens trobem amb els conceptes i els resultats següents:
les equacions diferencials en derivades parcials (eł.líptiques, parabòliques i hiperbòliques), amb les seves aplicacions;
l'estudi dels 
problemes de contorn (amb condicions de Dirichlet o de Neumann);
les tècniques numèriques d'interpolació polinòmica i l'aproximació de derivades;
la teoria de
matrius tridiagonals semidefinides positives (les quals també apareixen en els grafs distància-regulars);
la teoria de
matrius circulants i la seva relació amb els grafs de Cayley sobre grups cíclics;
la diagonalizació de matrius i els algorismes del càlcul dels valors propis i vectors propis;
el teorema dels cercles de Gershgorin (per localitzar els valor propis);
els polinomis ortogonals de variable discreta (també importants en els grafs distància-regulars);
l'estudi de condicions d'ortogonalitat en el cas discret respecte del cas continu;
les transformacions discretes (canvis de base) com, per exemple, la Transformada Discreta de Fourier (DFT);
la resolució de les recurrències lineals i les diferències finites;
la resolució d'equacions diferencials ordinàries (EDO's) lineals;
l'aproximació i càlcul d'errors; i, finalment,
la representació gràfica i la interpretació de les solucions.




§ EL MÈTODE DE LES LÍNIES (MOL): L'EQUACIÓ DE LAPLACE


	La transferència de calor a l'estat estacionari en sòlids segueix l'equació de Laplace i ha estat resolta directament per diversos mètodes numèrics, com el mètode successiu de sobre-relaxació, el mètode implícit de direcció alternativa i el mètode de transitoris falsos.
	Aquests mètodes són alternatius al mètode de les línies (MOL), el qual  



	
	és el que fem servir a continuació. El MOL consisteix en discretitzar totes les variables
	 
	 excepte una. De vegades, el MOL també és anomenat mètode semi-analític (vegeu Subramanian i White <cit.>).



 §.§ L'equació de Laplace en un domini rectangular


Considerem l'equació de Laplace per a la funció potencial V(x,y) en una regió rectangular a× b.

    ∇^2V    =∂^2 V/∂ x^2+∂^2 V/∂ y^2=0,      x∈ (0,a), y∈(0,b),

amb condicions de frontera que, per simplicitat, establim com

    V(0,y)      =V(a,y)=0,          y∈(0,b),
    
    V(x,0)      =0, V(x,b)=f(x),      x∈(0,a).

El primer pas és discretitzar la variable x∈(0,a) en els N punts x_j=jΔ x, per a j=1,2,…,N, on el “pas” h=Δ x=a/N+1 correspon a la distància entre dues de les N línies resultants de la discretització, segons mostrem a la Figura <ref>.



Ara, iterem l'aproximació lineal de la derivada en dos punts successius,

    V'(x_j)≃V(x_j+1)-V(x_j)/h,

i obtenim l'anomenada “aproximació centrada” de la segona derivada

    V”(x_j)≃V'(x_j)-V'(x_j-1)/h≃V(x_j+1)-2V(x_j)
    +V(x_j-1)/h^2,

amb error |ϵ_j|≤ Ch^2, per a una C constant, que correspon a un mètode de segon ordre.
Amb notació simplificada,

    ∂^2 V/∂ x^2≃V_j+1-2V_j+V_j-1/h^2,

on V_j denota la funció V(x_j,y). Si substituïm (<ref>) en (<ref>), obtenim les N equacions diferencials

    ∂^2 V/∂ y^2+1/h^2(V_j+1-2V_j+V_j-1)=0,   ji=1,…,N.

Aleshores, si reemplacem V(x,y) pel vector columna =(V_1,…,V_N)^⊤, podem escriure (<ref>) en forma matricial com

    ∂^2 /∂ y^2-1/h^2=0⃗,

on  és una matriu tridiagonal, és a dir,

    (
         [   ∂^2 V_1/∂ y^2;   ∂^2 V_2/∂ y^2;               ⋮; ∂^2 V_N-1/∂ y^2;   ∂^2 V_N/∂ y^2;                 ])-1/h^2(
         [ p_l  -1   0   ⋯   0   0   0;  -1   2  -1   ⋯   0   0   0;   0  -1   2   ⋯   0   0   0;               ⋱            ;   0   0   0   ⋯  -1   2  -1;   0   0   0   ⋯   0  -1 p_r;     ])
       (
         [   V_1;   V_2;     ⋮; V_N-1;   V_N;       ])=0⃗,

on els valors p_l i p_r depenen de les condicions de contorn a x=0 i x=a. Si la condició de contorn és de Dirichlet per l'esquerra [respectivament, per la dreta], aleshores p_l=2 [respectivament, p_r=2]; si la condició de contorn és de Neumann per l'esquerra [respectivament, per la dreta], aleshores p_l=1 [respectivament, p_r=1].

Per resoldre el sistema (<ref>) d'equacions diferencials ordinàries (EDOs) lineals de segon ordre acoblades, necessitem transformar-lo en un sistema d'EDOs desacoblades.
Aleshores, la solució és diagonalitzar la matriu .
En aquest sentit, notem que  és simètrica i, per tant, els seus valors propis λ_1,…,λ_N són reals. A més, el teorema de Gershgorin[Recordem que el Teorema de Gershgorin afirma que, donada una matriu =(a_jk)∈ M_n(ℂ), si es defineixen els cercles D_1,…,D_n amb centre a_jj i radi r_j=∑ _k≠ j|a_jk|, els valors propis de la matriu A es troben en la unió dels D_j.] ens assegura la no negativitat, ja que han de complir |λ_k-2|≤ 2, de manera que utilitzarem la notació λ_k=ω_k^2 per a k=1,…,N. De fet, com s'explica a la subsecció <ref>, la matriu tridiagonal  té tots els valors propis diferents, i dona lloc a una successió de Sturm de polinomis ortogonals, els quals estan relacionats estretament amb els corresponents vectors propis (vegeu, per exemple, Godsil <cit.>).

Per tant, tenim que
^⊤=, on  és la matriu ortogonal (amb columnes ortonormals corresponents als vectors propis de ), és a dir,  ^-1=^⊤, i  és la matriu diagonal dels valors propis de . Aleshores, =^⊤ que, si substituïm a (<ref>), ens dona

    ∂^2 /∂ y^2-1/h^2^⊤=0⃗,

i, si multipliquem per ^⊤(=^-1),

    ∂^2 ^⊤/∂ y^2-1/h^2^⊤=0⃗,

i, si utilitzem el potencial transformat :=^⊤, obtenim

    (∂^2/∂ y^2-1/h^2)=0⃗.

Per tant, obtenim N equacions diferencials ordinàries (en la variable y), del tipus

    V_j”=1/h^2ω_j^2,     j=1,…,N,

amb  solucions (el polinomi característic és
s^2=ω_j^2/h^2,
amb arrels s_1,2=±ω_j/h)

    V_j=α_j e^ω_j/hy+β_j e^-ω_j/hy=A_j coshω_j y/h+B_j sinhω_j y/h,   j=1,…,N.

Per trobar el vector de potencials per a cada j=1,…,N, només cal realitzar la transformada inversa de , o sigui,

    =.

Finalment, imposem les condicions de contorn i resolem les equacions resultants per determinar les constants α_j,β_j, o A_j,B_j.



§ DIAGONALITZACIÓ DE LA MATRIU  I ORTOGONALITAT




 §.§ Valors i vectors propis


Abans de resoldre un exemple concret, en aquest apartat discutim la diagonalització de la matriu  i les propietats d'ortogonalitat involucrades.
Els resultats obtinguts, els valors propis λ_k=ω_k^2 i les matrius de vectors propis =(T_jk), es mostren a la Taula <ref>, segons les condicions de contorn a x=0 i x=a.



Per al càlcul procedim de la manera següent. L'equació a resoldre és

    (
         [ p_l  -1   0   ⋯   0   0   0;  -1   2  -1   ⋯   0   0   0;   0  -1   2   ⋯   0   0   0;               ⋱            ;   0   0   0   ⋯  -1   2  -1;   0   0   0   ⋯   0  -1 p_r;     ])
       (
         [   v_1;   v_2;     ⋮; v_N-1;   v_N ])= λ(
         [   v_1;   v_2;     ⋮; v_N-1;   v_N ]),

d'on

    (p_l-λ)v_1    -v_2 =0, 
    
    (2-λ)v_j    -v_j-1 =v_j+1,   j=2,…,N-1 
    
    (p_r-λ)v_N    -v_N-1 =0,

on, per als valors usats de p_l i p_r, el teorema de Gershgorin i el fet que la matriu és simètrica impliquen que els valors propis estan a l'interval [0,4].

Per resoldre l'equació en diferències (<ref>), podríem usar la solució v_j=s^j, però s'obté un resultat més compacte amb v_j=e^ijα (és a dir, agafant s=e^iα=cosα+isinα). Si substituïm aquest valor a (<ref>), i dividim per e^i(j-1)α s'obté

    e^iα=(2-λ)-e^-iα  ⇒   2-λ=2cosα.

Per tant, l'equació característica és

    λ=2(1-cosα)=4sin^2(α/2),

cosa que significa que, per a cada valor propi λ, tenim dos possibles valors de α:

    α_1,2=± 2arcsin(√(λ)/2).

Així, amb α=α_1, la solució general de la recurrència és

    v_j=Ae^ijα+Be^-ijα.

Per trobar els possibles valors de les constants A, B, i α, imposem les “condicions inicials” (<ref>) i (<ref>).
Suposem que tenim condicions de Dirichlet per la dreta i l'esquerra, és a dir, p_l=p_r=2. (Els altres casos es resolen de manera similar).
De la primera equació, s'obté que

    (2-λ)v_1-v_2    = 2cosα(Ae^iα+Be^-iα)
    -Ae^i2α-Be^-i2α=0
       ⇒   A+B=0,

on hem aplicat cosα=e^iα+e^-iα/2 i sinα=e^iα-e^-iα/2i.
Per tant, B=-A i, llevat d'una constant multiplicativa, podem suposar que

    v_j=sin(jα),     j=1,…,N.

D'altra banda, de (<ref>) s'obté

    (2-λ)v_1-v_2     =2cosα(Ae^iNα+Be^-iNα)
    -Ae^i(N-1)α-Be^-i(N-1)α)=0
       ⇒  sin((N+1)α)=0.

Per tant, els possibles valors de α són α_k=kπ/N+1, per a k=1,…, N, amb el que obtenim els valors propis

    λ_k=4sin^2(kπ/2(N+1)),   k=1,…,N,

els quals compleixen que λ_1<λ_2<⋯<λ_N,
amb els respectius vectors propis

    _̌k=(v_jk)=(sin(kπ/N+1), sin(2kπ/N+1),…,sin(Nkπ/N+1))^⊤.

Finalment, com que les columnes de  han de ser ortonormals, tenim

    _̌k^2=∑_j=1^N sin^2(jkπ/N+1)=1/2∑_j=1^N (1-cos(2jkπ/N+1))=⋯ = N+1/2,

on hem usat que sin^2 x=1/2(1-cos(2x)) i la suma d'una sèrie geomètrica  ∑_j=1^N r^j= r(1-r^N)/1-r.
Concloem que els elements de la matriu  són

    T_jk=(v_jk)=√(2/N+1)sin(jkπ/N+1),    j,k=1,…,N,

com s'indica a la Taula <ref>.



 §.§ Ortogonalitat



Considerem ara les propietats d'ortogonalitat que apareixen en el MOL, i les comparem amb les corresponents del mètode analític.
En primer lloc, notem que, per ser una matriu simètrica, els elements de les subdiagonals de la matriu =(p_jk) satisfan p_j,j+1p_j+1,j>0 per a j=1,…,N i, per tant,  donen lloc a una família de polinomis ortogonals p_0,p_1,…, p_N-1 de variable discreta[És a dir, polinomis ortogonals respecte d'un producte escalar discret en els punts d'una malla x_1,…,x_N∈: ⟨ f, g⟩=∑_j=1^N w_jf(x_j)g(x_j) per a uns certs “pesos” w_j>0.]. Aquest és també el cas, per exemple, de la matriu quocient d'un graf distància-regular, els polinomis ortogonals del qual corresponen als anomenats “polinomis distància”. Els detalls d'aquesta aplicació, i altres de tipus combinatori, poden trobar-se a l'article de Cámara, Fàbrega, Fiol, i Garriga <cit.>. Aquests polinomis constitueixen una successió de Sturm (vegeu Godsil <cit.> or Chihara <cit.>) i satisfan una relació de recurrència de tres termes (com a l'expressió (<ref>)) 

    p_k+1(x)=(2-x)p_k(x)-p_k-1(x),   k=1,…,N-2,

inicialitzada amb p_0(x)=1 i p_1(x)=2-x. A més, p_N(x)=(2-x)p_N-1(x)-p_N-2(x)
correspon al polinomi característic de la matriu  i, per tant, les seves arrels són els valors propis λ_k, amb vectors propis

    _̌k=(p_0(λ_k),p_1(λ_k),…,p_N-1(λ_k))^⊤,     k=1,…, N,

trobats anteriorment a (<ref>) i (<ref>) (aquests últims llevat d'una constant multiplicativa).
L'ortogonalitat esmentada s'obté respecte d'un producte escalar discret en els punts de la malla
λ_1,λ_2,…, λ_N, amb λ_1<λ_2<⋯ <λ_N. 
Per tant, un cop que normalitzem
els vectors _̌k, és a dir _k=_̌k/_̌k, obtenim de nou els elements de la matriu  donats a (<ref>):

    T_j,k=(_k)_j,     j,k=1,…,N,

és a dir,

    T_j,k=√(2/N+1)sin(jkπ/N+1)
    =p_j-1(λ_k)/(∑_j=0^N-1p_j^2(λ_k))^1/2,
         j,k=1,…,N.

(Compareu amb les funcions pròpies que apareixen en el cas analític, vegeu la secció <ref>).
En aquest context,
observem que el potencial transformat =^⊤ no és més que el vector de coeficients de Fourier obtingut a l'expressar V en termes d'uns polinomis ortogonals (que, en els punts de la malla, tenen els valors donats per cada fila de ).

Per exemple, en el cas de la matriu  a (<ref>)
amb N=5 i p_l=p_r=2, s'obtenen els polinomis

    p_0(x)=    1,
    
    p_1(x)=    2-x,
    
    p_2(x)=    x^2-4x+3,
    
    p_3(x)=    -x^3+6x^2-10x+4, 
    
    p_4(x)=    x^4-8x^3+21x^2-20x+5, 
    
    p_5(x)=    -x^5+10x^4-36x^3+56x^2-35x+6.

Els valors propis de  són (arrels de p_5(x)):

    λ_1=2-√(3), λ_2=1, λ_3=2, λ_4=3, λ_5=2+√(3).

A la Figura <ref> es mostren les gràfiques de p_j(x), per a j=0,…,4. Noteu les paritats (parells i senars) d'aquests polinomis respecte de x=λ_3=2, induïdes per les corresponents simetries de la matriu




    =(
    [  √(3)/6     1/2  √(3)/3     1/2  √(3)/6;     1/2     1/2       0     1/2     1/2;  √(3)/3       0 -√(3)/3       0  √(3)/3;     1/2    -1/2       0     1/2    -1/2;  √(3)/6    -1/2  √(3)/3    -1/2  √(3)/6;         ]).




§ UN EXEMPLE




Ara anem a veure com apliquem el MOL a la resolució d'un problema concret.

Volem resoldre l'equació de Laplace (potencial)

    ∂^2 V/∂ x^2+∂^2 V/∂ y^2=0,

amb condicions de contorn de tipus Dirichlet V(0,y)=V(a,y)=V(x,0)=0 i V(x,b)=100, amb a=b=1 a la regió representada a la Figura <ref>. També volem calcular el valor del potencial en els punts indicats en aquesta figura.



La matriu  té valors propis (VAPS) ω_k^2=4sin^2kπ/2N+2, per a k=1,2,…,N, i la matriu de vectors propis (VEPS)  té components T_jk=√(2/N+1)sinjkπ/N+1 (vegeu la Taula <ref>).

Si escollim N=15, resulta que h= x=1/16, i les abscisses dels punts que volem calcular són x_4=1/4, x_8=1/2 i x_12=3/4.

De V_j=∑_k=1^N T_jkV_k, per a i=1,2,…,N, obtenim la solució

    (
         [ V_1; V_2;   ⋮; V_N;     ])=
       (
         [ T_11 T_12    ⋯ T_1N; T_21 T_22    ⋯ T_2N;              ⋱     ; T_N1 T_N2    ⋯ T_NN;      ])
       (
         [ A_1 coshω_1 y/h+B_1 sinhω_1 y/h; A_2 coshω_2 y/h+B_2 sinhω_2 y/h;                               ⋮; A_N coshω_N y/h+B_N sinhω_N y/h;                                 ]).

Ara, per determinar les constants A_j i B_j, imposem les condicions de contorn:


  * V(x,0)=0 ⇒


    (
         [ V_1; V_2;   ⋮; V_N;     ])=
       (
         [ T_11 T_12    ⋯ T_1N; T_21 T_22    ⋯ T_2N;              ⋱     ; T_N1 T_N2    ⋯ T_NN;      ])
       (
         [ A_1; A_2;   ⋮; A_N;     ])=0

d'on
A_j=0, per a j=1,2,…,N.


  * V(x,1)=100 ⇒





    (
         [ V_1; V_2;   ⋮; V_N;     ])=
       (
         [ T_11 T_12    ⋯ T_1N; T_21 T_22    ⋯ T_2N;              ⋱     ; T_N1 T_N2    ⋯ T_NN;      ])
       (
         [ B_1 sinhω_1/h; B_2 sinhω_2/h;             ⋮; B_N sinhω_N/h;               ])=
       (
         [ 100; 100;   ⋮; 100;     ]),

d'on
  
    (
         [ B_1 sinhω_1/h; B_2 sinhω_2/h;             ⋮; B_N sinhω_N/h;               ])=
       (
         [ T_11 T_12    ⋯ T_1N; T_21 T_22    ⋯ T_2N;              ⋱     ; T_N1 T_N2    ⋯ T_NN;      ])^⊤(
         [ 100; 100;   ⋮; 100;     ]):=
       (
         [ c_1; c_2;   ⋮; c_N;     ]).

Per tant,
  
    B_j=c_j/sinhω_j y/h,     j=1,2,…,N.

  
Així, arribem a la solució

    V_j(y)     =     V(x_j,y)=∑_k=1^N T_jkB_ksinh(ω_j y/h)
        =    ∑_k=1^N √(2/N+1)sin(jkπ/N+1)
           c_k/sinh(ω_k/h)sinh(ω_k y/h),

la qual està representada a la Figura <ref>.

En els punts demanats els valors que s'obtenen (calculats amb el programa de la Secció <ref>) són

    V(1/4,3/4)=43.1008,   V(1/2,1/2)=24.9644,   V(3/4,1/4)=6.7984

que, comparats amb els valors calculats aplicant els resultats de la Secció <ref> també amb el programa de la Secció <ref>,

    V(1/4,3/4)=43.2028,    V(1/2,1/2)=24.9999,   V(3/4,1/4)=6.7971,

demostren ser molt bones aproximacions.



§ EL MÈTODE ANALÍTIC



Ara volem resoldre el problema anterior amb el mètode analític, per poder-lo comparar amb el resultat obtingut amb el MOL. 
És a dir, ja que podem conèixer la solució analítica, emprem l'equació de Laplace com a test pel mètode de les línies.  
	
Considerem doncs de nou l'equació

    ∇^2V=∂^2 V/∂ x^2+∂^2 V/∂ y^2=0,

per a x∈(0,a) i y∈(0,b), amb condicions de contorn

    V(0,y)   = V(a,y)=V(x,0)=0,
    
    V(x,b)   =f(x),

com es veu a la Figura <ref>.



La solució analítica s'obté per separació de variables, és a dir, suposant que el potencial és de la forma
V(x,y)=X(x)Y(y), d'on trobem que X i Y satisfan les equacions

    X”(x)+ω^2 X(x)     =0,
    
    Y”(y)-ω^2 Y(y)    =0.

Imposant les condicions de contorn homogènies, resulta que ω=ω_k=kπ/a per a k=1,2,… i obtenim  les funcions

    X(x)     =Bsin(kπ/ax),
    
    Y(y)     =Csinh(kπ/ay),

amb B i C constants.
Per tant, les solucions del problema homogeni són:

    V_k(x,y)=A_ksinh(kπ/ay)sin(kπ/ax).

Considerem ara una solució en forma de sèrie

    V(x,y)=∑_k=1^∞A_ksinh(kπ/ay)
    sin(kπ/ax).




Si imposem la condició de contorn no homogènia

    V(x,b)=f(x)=∑_k=1^∞A_ksinh(kπ/ab)
    sin(kπ/ax),

i considerem que {sin(kπ/ax):k=1,2,…} és una base ortogonal a (0,a), obtenim el corresponent desenvolupament en sèrie de Fourier de f(x), amb coeficients

    α_k=A_ksinh(kπ/ab)=2/a∫_0^a f(x)sin(kπ/ax)  dx.

Per tant, la solució analítica és

    V(x,y)=∑_k=1^∞α_ksinh(kπ/ay)/sinh(kπ/ab)sin(kπ/ax).

Observem que les condicions de contorn iguals a zero a x=0 i x=a donen solucions periòdiques en sinus, mentre que les condicions de contorn diferents a y=0 i y=b donen solucions no periòdiques en sinus hiperbòlics.

En l'exemple anterior utilitzant el MOL, consideràvem el cas particular a=b=1 i f(x)=100. En aquest cas, ara obtenim que α_n=400/nπ per a n senar i  α_n=0 per a n parell. Aleshores, la solució analítica és

    V(x,y)=400/π∑_k=1^∞1/2k-1sinh((2k-1)π y)/sinh((2k-1)π)sin((2k-1)π x),

la representació gràfica de la qual es mostra a la Figura <ref>.



§ SOBRE L'ERROR



Podem calcular l'error comès quan utilitzem el mètode de les línies, simplement amb la diferència entre els resultats obtinguts amb el MOL i el mètode analític. 
De fet, se sap que, sota condicions generals de regularitat de les solucions, l'error del mètode de les línies és de l'ordre de O(h^2) (vegeu, per exemple, Schieser <cit.>, o Zuazua <cit.>).

Nosaltres estudiem primer el comportament de l'error
 en el cas de l'exemple anterior (Secció <ref>). Vegeu la Figura <ref>. Observem que tenim 8 funcions d'error (per a x_j=j/16, j=1,…,8) malgrat que N=15 perquè, per simetria en les condicions de contorn, hi ha 7 funcions d'error que es superposen a d'altres 7 (i en queda una sense superposar). Com era d'esperar a causa de la discontinuïtat en les condicions de contorn, els errors són més petits quan som més a prop del centre de la gràfica (és a dir, per a x_8=1/2).







Per obtenir resultats més precisos, 
considerem ara un altre exemple més senzill, el qual ens permet donar la solució numèrica amb una fórmula tancada. 
 Així, suposem ara que les condicions de contorn són
V(0,y)=V(1,y)=V(x,0)=0, com abans, i que V(x,1)=sinh(π)sin(π x). Aleshores, per trobar la solució analítica no cal fer l'últim pas del desenvolupament en sèrie de Fourier, i s'obté 

    V(x,y)=sinh(π y) sin (π x),

	la gràfica de la qual es mostra a la Figura <ref>.
Pel que fa al mètode de línies, comentem primer que aquest cas va ser resolt computacionalment per Subramanian i White <cit.> utilitzant matrius exponencials, encara que sense arribar a una expressió tancada.
En canvi, les tècniques descrites aquí ens permeten obtenir la següent solució compacta

    V(x_j,y)=sinh[π(h) y]sin(π x_j),

on x_j=jh, j=1,…,N, i (x)=sin(π x)/π x. Vegeu la Figura <ref>.

Notem que, com lim_h→ 0(h)=1,  quan h→ 0, j→∞, i jh→ x, 
	i aleshores 
V(x_j,y) → V(x,y), com era d'esperar. De fet, en aquest cas és fàcil comprovar que, per a cada x_j=x i y,  l'error ϵ=V(x,y)-V(x_j,y) és de l'ordre de O(h^2), com ja s'ha comentat.















§ PROGRAMA EN MATLAB



A continuació, escrivim en el llenguatge del programa MATLAB les instruccions per calcular el valor del potencial en els punts que es demanen a l'exemple anterior.





>> AA=1;

>> BB=1;

>> N=15;

>> %Determinem el vector alpha

>> H=AA/(N+1);

>> LAMBDA=2*sin((1:N)*pi*0.5/(N+1));

>> ALPHA=LAMBDA/H;

>> %Calculem la matriu de transformació

>> S=sqrt(2/(N+1));

>> T=zeros(N,N);

>> for J=(1:N)

 for K=(1:N)

        T(J,K)=S*sin(J*K*pi/(N+1));

      end

  end

>> V=100*ones(N,1);

>> C=inv(T)*V;

>> A=ALPHA';

>> B=C./sinh(BB*A);

>> %Calculem V en els punts donats

>> V1=0; V2=0; V3=0;

>> for K=1:N

 V1=V1+T(4,K)*B(K)*sinh(ALPHA(K)*0.75);

 V2=V2+T(8,K)*B(K)*sinh(ALPHA(K)*0.5);

 V3=V3+T(12,K)*B(K)*sinh(ALPHA(K)*0.25);

 end

>> diary

>> V1, V2, V3

 V1 = 43.1008

 V2 = 24.9644

 V3 = 6.7984

>> diary off

.5cm







































§ ALTRES EXEMPLES RESOLUBLES AMB EL MOL


Finalment, detallem alguns altres exemples que es poden resoldre amb el mètode de les línies.



  1. 
L'equació de Laplace en la regió quadrada 1×1 amb condicions de contorn:
  
    
  (a) Dirichlet-Neumann.
    
  (b) Neumann-Dirichlet.
    
  (c) Neumann-Neumann.
  
  (Els valors a utilitzar de la matriu  i els valors propis ω_k^2 estan a la Taula <ref>).
  
  2.
  
  
L'equació de Laplace en coordenades cilíndriques
  ρ^2∂^2 V/∂ρ^2+ρ∂ V/∂ρ+∂^2 ϕ/∂ϕ^2=0, amb les  condicions de contorn V(ρ,ϕ)=0 per a 0≤ρ≤ a i
  V(ρ,ϕ)=V_0 per a a<ρ≤ b, amb la discretització de la coordenada angular ϕ.
    

La solució analítica és V(ρ)=V_0ln(ρ/a)/ln(b/a).

  

  3.
  
    
L'equació d'ona ∂^2 v/∂ t^2-∂^2 v/∂ x^2=q(x,t), amb 0≤ x≤ L i 0≤ t≤ T, les condicions inicials v(x,0)=f_1(x) i ∂ v/∂ t(x,0)=f_2(x) per a 0≤ x≤ L, la condició de contorn de Dirichlet v(0,t)=g_1(x) i la condició no local ∫_0^L v(x,t)dx=g_2(x) per a 0≤ t≤ T.
    
A la taula següent hi ha alguns casos particulars amb la corresponent solució analítica.
    







§ CONCLUSIONS

	El mètode de les línies ha demostrat ser una bona eina per resoldre alguns tipus d'equacions en derivades parcials amb condicions de contorn. Com que es basa en discretitzar totes les variables excepte una, es pot utilitzar directament en l'equació de Laplace, la qual hem fet servir per il·lustrar aquest mètode. Una alternativa és aplicar el mètode de les falses transicions, afegint una derivada temporal i discretitzant les variables espacials.
	Si estem al pla xy, l'error d'aquest mètode alternatiu és de l'ordre de Δ_x^2+Δ_y^2. Per tant, és
	
	 millor aplicar el mètode semi-discret (o semi-analític), que és el que hem aplicat, i en el qual hem vist que l'error és de l'ordre de Δ_x^2.
	
	Per una altra banda, 
	es pot 
	estudiar com canvia l'error si, en comptes de fer línies equiespaiades, fem línies situades als zeros dels polinomis de 
	Txebixov 
	
	(els quals tenen un millor comportament respecte de l'error). En aquest context, es pot consultar el treball de Youssef i Shukur <cit.>. Finalment, proposem el següent problema obert: Estudiar  les semblances i les diferències entre el MOL i la integració al llarg de les rectes característiques de l'equació en derivades parcials de primer ordre ∂ V/∂ t+c∂ V/∂ x=f(x,t), amb c>0.
	


§ AGRAÏMENTS


Els autors volen agrair els comentaris i suggeriments del revisor anònim, els quals han contribuït considerablement a millorar aquest treball.






9
BeVu16
M. Berardia, M. Vurro,
The numerical solution of Richards' equation by means of method of lines and ensemble Kalman filter,
Math. Comput. Simulation 125 (2016) 38–47.

cffg09
M. Cámara, J. Fàbrega, M. A. Fiol, E. Garriga,
Some families of orthogonal polynomials of a discrete variable and their applications to graphs and codes,
Electron. J. Combin. 16 (2009) #R83.

c78
T. S. Chihara,
An Introduction to Orthogonal Polynomials,

Gordon and Breach, New York, 1978.

d84
H. Diestel, Analysis of planar multiconductor transmission-line systems with the method of lines,
Arch. Elektron. Uebertragungstech 41 (1987) 169–175.

g93
C. D. Godsil,
Algebraic Combinatorics,
Chapman and Hall, New York, 1993.






SaOb
M. N. O. Sadiku, C. N. Obiozor,
A simple introduction to the method of lines,
Int. J. Elec. Eng. Educ. 37/3 (2000) 282–296.

s91
W. E. Schiesser, 
	The Numerical Method of Lines, 
	Academic Press, New York, 1991.

ShDe08
F. Shakeri, M. Dehghan,
The method of lines for solution of the one-dimensional wave equation subject to an integral conservation condition,
Comput. Math. Appl. 56 (2008) 2175–2188.

sw04
V. R. Subramanian, R. E. White,
Semianalytical method of lines for solving elliptic partial differential equations, Chem. Eng. Sci. 59 (2004) 781–788.

vs84
J. G. Verwer, J. M. Sanz-Serna, 
Convergence of method of lines approximations to partial differential equations,
Computing 33 (1984), no. 3-4, 297–313.


yss14
I. K. Youssef, A. M. Shukur, The line method combined with spectral Chebyshev for space-time fractional diffusion equation, Appl. Comput. Math. 3 (2014), no. 6,  330–336.


Za70
	A. Zafarullah,
Applications of the method of lines to parabolic differential equations with error estimates,
J. Assoc. Comput. Mach. 17 (1970), no. 2, 294–302.

Zu09a
E. Zuazua,
Introducción al análisis numérico de ecuaciones en
	derivadas parciales de evolución,
Universidad Autónoma de Madrid, 2003.

Zu09b
E. Zuazua,
Métodos numéricos de resolución de ecuaciones en derivadas parciales,
Basque Center for Applied Mathematics (BCAM), 2009.











