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

El m`etode de les l´ınies per a la resoluci´o num`erica d’equacions en derivades parcials The method of lines for numerical solutions of partial differential equations C. Dalf´o Departament de Matem`atica, Universitat de Lleida, Igualada (Barcelona), Catalunya [email protected] M. A. Fiol Deptartament de Matem`atiques, Barcelona Graduate School of Mathematics, Universitat Polit`ecnica de Catalunya, Barcelona, Catalunya [email protected] Resum En aquest treball es descriu un m`etode num`eric semi-discret per a la resoluci´o d’un tipus d’equacions en derivades parcials. ´ Es conegut com el m`etode de les l´ınies (MOL, per les sigles en angl`es), i es basa en la discretitzaci´o de totes les variables involucrades, excepte una. Il.lustrem l’aplicaci´o del MOL amb la resoluci´o de l’equaci´o de Laplace en coordenades cartesianes. Els conceptes matem`atics que hi apareixen es comparen amb els corresponents al m`etode anal´ıtic de separaci´o de variables. A m´es, mostrem que els resultats obtinguts amb el MOL s´on aproximacions molt bones de les solucions anal´ıtiques. 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´o d’una variable cont´ınua, an`alisi num`eric. Keywords: Partial differential equations, discretization of a continuous variable, numerical analysis. Per apar`eixer al Butllet´ı de la Societat Catalana de Matem`atiques (2019). 1 2 1 Introducci´o El m`etode de les l´ınies, tamb´e anomenat m`etode de semi-discretitzaci´o per Zuazua [13] ´es un m`etode per a la resoluci´o 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´e cont´ınua la variable temporal). Aix`o condueix a un sistema d’equacions diferencials ordin`aries (EDO’s), les quals es poden resoldre pels m`etodes num`erics habituals per a EDO’s amb condicions inicials. Aquest m`etode 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 [12] i Verwer i Sanz-Serna [10]). Per tant, el MOL ´es un enfocament diferencial i de difer`encies finites per resoldre equacions en derivades parcials. Aquest m`etode t´e moltes aplicacions en la resoluci´o de problemes en F´ısica i Enginyeria. Per exemple, aquest m`etode s’ha usat `ampliament per experts en t`ecniques computacionals en problemes d’electromagnetisme, vegeu, per exemple, Berardia i Vurro [1], Diestel [4], Sadiku i Obiozor [6], o Shakeri i Dehghan [8]. En aquest treball, il·lustrem el m`etode de les l´ınies amb l’equaci´o de Laplace en el pla, per simplicitat. El proc´es a seguir consta dels punts seg¨uents: 1. Discretitzaci´o de l’equaci´o diferencial en una direcci´o: Contracci´o de la regi´o de solucions a un conjunt de l´ınies, i aproximaci´o de les derivades corresponents. 2. Obtenci´o d’un sistema d’equacions diferencials ordin`aries acoblades. 3. Transformaci´o per obtenir un sistema d’equacions diferencials ordin`aries desacoblades. 4. Resoluci´o de les equacions diferencials. 5. Transformaci´o inversa i introducci´o de les condicions de contorn. 6. Soluci´o de les equacions. A m´es de les seves m´ultiples aplicacions, per exemple, com ja s’ha dit, en electromagnetisme, aquest m`etode involucra moltes idees, conceptes i t`ecniques de la matem`atica discreta, juntament amb els m`etodes b`asics de l’an`alisi diferencial. De fet, en el desenvolupament del m`etode ens trobem amb els conceptes i els resultats seg¨uents: les equacions diferencials en derivades parcials (el.l´ıptiques, parab`oliques i hiperb`oliques), amb les seves aplicacions; l’estudi dels problemes de contorn (amb condicions de Dirichlet o de Neumann); les t`ecniques num`eriques d’interpolaci´o polin`omica i l’aproximaci´o de derivades; la teoria de matrius tridiagonals semidefinides positives (les quals tamb´e apareixen en els grafs dist`ancia-regulars); la teoria de matrius circulants i la seva relaci´o amb els grafs de Cayley sobre grups c´ıclics; la diagonalizaci´o de matrius i els algorismes del c`alcul dels 3 valors propis i vectors propis; el teorema dels cercles de Gershgorin (per localitzar els valor propis); els polinomis ortogonals de variable discreta (tamb´e importants en els grafs dist`ancia-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´o de les recurr`encies lineals i les difer`encies finites; la resoluci´o d’equacions diferencials ordin`aries (EDO’s) lineals; l’aproximaci´o i c`alcul d’errors; i, finalment, la representaci´o gr`afica i la interpretaci´o de les solucions. 2 El m`etode de les l´ınies (MOL): L’equaci´o de Laplace La transfer`encia de calor a l’estat estacionari en s`olids segueix l’equaci´o de Laplace i ha estat resolta directament per diversos m`etodes num`erics, com el m`etode successiu de sobre-relaxaci´o, el m`etode impl´ıcit de direcci´o alternativa i el m`etode de transitoris falsos. Aquests m`etodes s´on alternatius al m`etode de les l´ınies (MOL), el qual ´es el que fem servir a continuaci´o. El MOL consisteix en discretitzar totes les variables excepte una. De vegades, el MOL tamb´e ´es anomenat m`etode semi-anal´ıtic (vegeu Subramanian i White [9]). 2.1 L’equaci´o de Laplace en un domini rectangular Considerem l’equaci´o de Laplace per a la funci´o potencial V(x, y) en una regi´o rectangular a×b. ∇2V=∂2V ∂x2+∂2V ∂y2= 0, x ∈(0, a), y ∈(0, b),(1) amb condicions de frontera que, per simplicitat, establim com V(0, y) = V(a, y)=0, y ∈(0, b),(2) V(x, 0) = 0, V (x, b) = f(x), x ∈(0, a).(3) El primer pas ´es discretitzar la variable x∈(0, a) en els Npunts xj=j∆x, per a j= 1,2, . . . , N, on el “pas” h= ∆x=a N+1 correspon a la dist`ancia entre dues de les N l´ınies resultants de la discretitzaci´o, segons mostrem a la Figura 1. Ara, iterem l’aproximaci´o lineal de la derivada en dos punts successius, V0(xj)≃V(xj+1)−V(xj) h, i obtenim l’anomenada “aproximaci´o centrada” de la segona derivada V00(xj)≃V0(xj)−V0(xj−1) h≃V(xj+1)−2V(xj) + V(xj−1) h2, 4 0 x1x2a h x y V1VN b xN Figura 1: Discretitzaci´o de la variable x. amb error |j| ≤ Ch2, per a una Cconstant, que correspon a un m`etode de segon ordre. Amb notaci´o simplificada, ∂2V ∂x2≃Vj+1 −2Vj+Vj−1 h2,(4) on Vjdenota la funci´o V(xj, y). Si substitu¨ım (4) en (1), obtenim les Nequacions diferencials ∂2V ∂y2+1 h2(Vj+1 −2Vj+Vj−1)=0, ji = 1, . . . , N. (5) Aleshores, si reemplacem V(x, y) pel vector columna V= (V1, . . . , VN)>, podem escriure (5) en forma matricial com ∂2V ∂y2−1 h2P V =0,(6) on P´es una matriu tridiagonal, ´es a dir,          ∂2V1 ∂y2 ∂2V2 ∂y2 . . . ∂2VN−1 ∂y2 ∂2VN ∂y2          −1 h2          pl−1 0 ··· 000 −1 2 −1··· 000 0−1 2 ··· 000 ... 000··· −1 2 −1 000··· 0−1pr                 V1 V2 . . . VN−1 VN        =0,(7) on els valors pliprdepenen de les condicions de contorn a x= 0 i x=a. Si la condici´o de contorn ´es de Dirichlet per l’esquerra [respectivament, per la dreta], aleshores pl= 2 [respectivament, pr= 2]; si la condici´o de contorn ´es de Neumann per l’esquerra [respectivament, per la dreta], aleshores pl= 1 [respectivament, pr= 1]. Per resoldre el sistema (6) d’equacions diferencials ordin`aries (EDOs) lineals de segon ordre acoblades, necessitem transformar-lo en un sistema d’EDOs desacoblades. Aleshores, la soluci´o ´es diagonalitzar la matriu P. En aquest sentit, notem que P´es sim`etrica i, per 5 tant, els seus valors propis λ1, . . . , λNs´on reals. A m´es, el teorema de Gershgorin1ens assegura la no negativitat, ja que han de complir |λk−2| ≤ 2, de manera que utilitzarem la notaci´o λk=ω2 kper a k= 1, . . . , N. De fet, com s’explica a la subsecci´o 3.2, la matriu tridiagonal Pt´e tots els valors propis diferents, i dona lloc a una successi´o de Sturm de polinomis ortogonals, els quals estan relacionats estretament amb els corresponents vectors propis (vegeu, per exemple, Godsil [5, §8,5]). Per tant, tenim que T>P T =D, on T´es la matriu ortogonal (amb columnes ortonormals corresponents als vectors propis de P), ´es a dir, T−1=T>, i D´es la matriu diagonal dels valors propis de P. Aleshores, P=T DT >que, si substitu¨ım a (6), ens dona ∂2V ∂y2−1 h2T DT >V=0, i, si multipliquem per T>(= T−1), ∂2T>V ∂y2−1 h2DT >V=0, i, si utilitzem el potencial transformat V:= T>V, obtenim ∂2 ∂y2−1 h2DV=0.(8) Per tant, obtenim Nequacions diferencials ordin`aries (en la variable y), del tipus V00 j=1 h2ω2 j, j = 1, . . . , N, amb solucions (el polinomi caracter´ıstic ´es s2=ω2 j h2, amb arrels s1,2=±ωj h) Vj=αje ωj hy+βje−ωj hy=Ajcosh ωjy h+Bjsinh ωjy h, j = 1, . . . , N. (9) Per trobar el vector de potencials per a cada j= 1, . . . , N, nom´es cal realitzar la transformada inversa de V, o sigui, V=T V . Finalment, imposem les condicions de contorn i resolem les equacions resultants per determinar les constants αj, βj, o Aj, Bj. 1Recordem que el Teorema de Gershgorin afirma que, donada una matriu A= (ajk)∈Mn(C), si es defineixen els cercles D1,...,Dnamb centre ajj i radi rj=Pk6=j|ajk|, els valors propis de la matriu Aes troben en la uni´o dels Dj. 6 Esquerra Dreta Tjk λk=ω2 k Dirichlet Dirichlet q2 N+1 sin jkπ N+14 sin2kπ 2N+2 Dirichlet Neumann q2 N+0.5sin j(k−0.5)π N+0.54 sin2(k−0.5)π 2N+1  Neumann Dirichlet q2 N+0.5cos (j−0.5)(k−0.5)π N+0.54 sin2(k−0.5)π 2N+1  Neumann Neumann q2 Ncos (j−0.5)(k−1)π Nsi k > 1 4 sin2(k−1)π 2N 1 √Nsi k= 1 Taula 1: Components de la matriu Ti valors propis λk=ω2 ksegons les condicions de contorn a x= 0 i x=a. 3 Diagonalitzaci´o de la matriu Pi ortogonalitat 3.1 Valors i vectors propis Abans de resoldre un exemple concret, en aquest apartat discutim la diagonalitzaci´o de la matriu Pi les propietats d’ortogonalitat involucrades. Els resultats obtinguts, els valors propis λk=ω2 ki les matrius de vectors propis T= (Tjk), es mostren a la Taula 1, segons les condicions de contorn a x= 0 i x=a. Per al c`alcul procedim de la manera seg¨uent. L’equaci´o a resoldre ´es          pl−1 0 ··· 000 −1 2 −1··· 000 0−1 2 ··· 000 ... 000··· −1 2 −1 000··· 0−1pr                 v1 v2 . . . vN−1 vN        =λ       v1 v2 . . . vN−1 vN        ,(10) d’on (pl−λ)v1−v2= 0,(11) (2 −λ)vj−vj−1=vj+1, j = 2, . . . , N −1 (12) (pr−λ)vN−vN−1= 0,(13) on, per als valors usats de plipr, el teorema de Gershgorin i el fet que la matriu ´es sim`etrica impliquen que els valors propis estan a l’interval [0,4]. Per resoldre l’equaci´o en difer`encies (12), podr´ıem usar la soluci´o vj=sj, per`o s’obt´e un resultat m´es compacte amb vj=eijα (´es a dir, agafant s=eiα = cos α+isin α). Si substitu¨ım aquest valor a (12), i dividim per ei(j−1)αs’obt´e eiα = (2 −λ)−e−iα ⇒2−λ= 2 cos α. (14) 7 Per tant, l’equaci´o caracter´ıstica ´es λ= 2(1 −cos α) = 4 sin2α 2,(15) cosa que significa que, per a cada valor propi λ, tenim dos possibles valors de α: α1,2=±2 arcsin √λ 2!. Aix´ı, amb α=α1, la soluci´o general de la recurr`encia ´es vj=Aeijα +Be−ijα.(16) Per trobar els possibles valors de les constants A,B, i α, imposem les “condicions inicials” (11) i (13). Suposem que tenim condicions de Dirichlet per la dreta i l’esquerra, ´es a dir, pl=pr= 2. (Els altres casos es resolen de manera similar). De la primera equaci´o, s’obt´e que (2 −λ)v1−v2= 2 cos α(Aeiα +Be−iα)−Aei2α−Be−i2α= 0 ⇒A+B= 0, on hem aplicat cos α=eiα+e−iα 2i sin α=eiα−e−iα 2i. Per tant, B=−Ai, llevat d’una constant multiplicativa, podem suposar que vj= sin(jα), j = 1, . . . , N. D’altra banda, de (13) s’obt´e (2 −λ)v1−v2= 2 cos α(AeiNα +Be−iNα)−Aei(N−1)α−Be−i(N−1)α)= 0 ⇒sin((N+ 1)α) = 0. Per tant, els possibles valors de αs´on αk=kπ N+1 , per a k= 1, . . . , N, amb el que obtenim els valors propis λk= 4 sin2kπ 2(N+ 1), k = 1, . . . , N, (17) els quals compleixen que λ1< λ2<··· < λN, amb els respectius vectors propis vk= (vjk) = sin kπ N+ 1,sin 2kπ N+ 1,...,sin Nkπ N+ 1>.(18) Finalment, com que les columnes de Than de ser ortonormals, tenim kvkk2= N X j=1 sin2jkπ N+ 1=1 2 N X j=1 1−cos 2jkπ N+ 1=··· =N+ 1 2, on hem usat que sin2x=1 2(1 −cos(2x)) i la suma d’una s`erie geom`etrica PN j=1 rj= r(1−rN) 1−r. Concloem que els elements de la matriu Ts´on Tjk = (vjk) = r2 N+ 1 sin jkπ N+ 1, j, k = 1, . . . , N, (19) com s’indica a la Taula 1. 8 3.2 Ortogonalitat Considerem ara les propietats d’ortogonalitat que apareixen en el MOL, i les comparem amb les corresponents del m`etode anal´ıtic. En primer lloc, notem que, per ser una matriu sim`etrica, els elements de les subdiagonals de la matriu P= (pjk) satisfan pj,j+1pj+1,j >0 per a j= 1, . . . , N i, per tant, donen lloc a una fam´ılia de polinomis ortogonals p0, p1, . . . , pN−1de variable discreta2. Aquest ´es tamb´e el cas, per exemple, de la matriu quocient d’un graf dist`ancia-regular, els polinomis ortogonals del qual corresponen als anomenats “polinomis dist`ancia”. Els detalls d’aquesta aplicaci´o, i altres de tipus combinatori, poden trobar-se a l’article de C´amara, F`abrega, Fiol, i Garriga [2]. Aquests polinomis constitueixen una successi´o de Sturm (vegeu Godsil [5, §8.5] or Chihara [3]) i satisfan una relaci´o de recurr`encia de tres termes (com a l’expressi´o (12)) pk+1(x) = (2 −x)pk(x)−pk−1(x), k = 1, . . . , N −2,(20) inicialitzada amb p0(x) = 1 i p1(x)=2−x. A m´es, pN(x) = (2 −x)pN−1(x)−pN−2(x) correspon al polinomi caracter´ıstic de la matriu Pi, per tant, les seves arrels s´on els valors propis λk, amb vectors propis vk= (p0(λk), p1(λk), . . . , pN−1(λk))>, k = 1, . . . , N, trobats anteriorment a (17) i (18) (aquests ´ultims llevat d’una constant multiplicativa). L’ortogonalitat esmentada s’obt´e 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 vk, ´es a dir vk=vk/kvkk, obtenim de nou els elements de la matriu Tdonats a (19): Tj,k = (vk)j, j, k = 1, . . . , N, ´es a dir, Tj,k =r2 N+ 1 sin jkπ N+ 1=pj−1(λk) PN−1 j=0 p2 j(λk)1 2 , j, k = 1, . . . , N. (Compareu amb les funcions pr`opies que apareixen en el cas anal´ıtic, vegeu la secci´o 5). En aquest context, observem que el potencial transformat V=T>Vno ´es m´es que el vector de coeficients de Fourier obtingut a l’expressar Ven termes d’uns polinomis ortogonals (que, en els punts de la malla, tenen els valors donats per cada fila de T). Per exemple, en el cas de la matriu Pa (10) amb N= 5 i pl=pr= 2, s’obtenen els 2´ Es a dir, polinomis ortogonals respecte d’un producte escalar discret en els punts d’una malla x1,...,xN∈R:hf, gi=PN j=1 wjf(xj)g(xj) per a uns certs “pesos” wj>0. 9 Figura 2: Els polinomis ortogonals pj(x), per a j= 0,...,4. polinomis p0(x) = 1, p1(x) = 2 −x, p2(x) = x2−4x+ 3, p3(x) = −x3+ 6x2−10x+ 4, p4(x) = x4−8x3+ 21x2−20x+ 5, p5(x) = −x5+ 10x4−36x3+ 56x2−35x+ 6. Els valors propis de Ps´on (arrels de p5(x)): λ1= 2 −√3, λ2= 1, λ3= 2, λ4= 3, λ5= 2 + √3. A la Figura 2 es mostren les gr`afiques de pj(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 T=      √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       . 16 Figura 8: Gr`afica del potencial V(x, y) en el cas a=b= 1, V(0, y) = V(a, y) = V(x, 0) = 0 iV(x, b) = sinh(π) sin(πx). litzant matrius exponencials, encara que sense arribar a una expressi´o tancada. En canvi, les t`ecniques descrites aqu´ı ens permeten obtenir la seg¨uent soluci´o compacta V(xj, y) = sinh[πsinc(h)y] sin(πxj),(23) on xj=jh,j= 1, . . . , N, i sinc(x) = sin(πx) πx . Vegeu la Figura 9. Notem que, com limh→0sinc(h) = 1, quan h→0, j→ ∞, i jh →x, i aleshores V(xj, y)→V(x, y), com era d’esperar. De fet, en aquest cas ´es f`acil comprovar que, per a cada xj=xiy, l’error =V(x, y)−V(xj, y) ´es de l’ordre de O(h2), com ja s’ha comentat. 7 Programa en MATLAB A continuaci´o, 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. 17 Figura 9: Gr`afica del potencial Vj(y) = V(xj, y), xj=jh,j= 1,...,15, per a V(x, b) = sinh(π) sin(πx). >> 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´o >> 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); 18 >> C=inv(T)*V; >> A=ALPHA’; >> B=C./sinh(BB*A); >> %Calculem Ven 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 8 Altres exemples resolubles amb el MOL Finalment, detallem alguns altres exemples que es poden resoldre amb el m`etode de les l´ınies. 1. L’equaci´o de Laplace en la regi´o quadrada 1 ×1 amb condicions de contorn: (a) Dirichlet-Neumann. (b) Neumann-Dirichlet. (c) Neumann-Neumann. (Els valors a utilitzar de la matriu Ti els valors propis ω2 kestan a la Taula 1). 2. L’equaci´o de Laplace en coordenades cil´ındriques ρ2∂2V ∂ρ2+ρ∂V ∂ρ +∂2φ ∂φ2= 0, amb les condicions de contorn V(ρ, φ) = 0 per a 0 ≤ρ≤aiV(ρ, φ) = V0per a a<ρ≤b, amb la discretitzaci´o de la coordenada angular φ. La soluci´o anal´ıtica ´es V(ρ) = V0 ln(ρ a) ln(b a). 3. L’equaci´o d’ona ∂2v ∂t2−∂2v ∂x2=q(x, t), amb 0 ≤x≤Li 0 ≤t≤T, les condicions inicials v(x, 0) = f1(x) i ∂v ∂t (x, 0) = f2(x) per a 0 ≤x≤L, la condici´o de contorn de Dirichlet v(0, t) = g1(x) i la condici´o no local RL 0v(x, t)dx =g2(x) per a 0 ≤t≤T. A la taula seg¨uent hi ha alguns casos particulars amb la corresponent soluci´o anal´ıtica. 19 L T q(x, t)f1(x)f2(x)g1(t)g2(t)v(x, t) anal´ıtica 1 4 0 0 πcos(πx) sin(πt) 0 cos(πx) sin(πt) 1 5 −2(x−t)e−x−t0xe−x0−2te−t−1xte−x−t 0 0 +te−txte−x−t 1 4 0 cos(πx) 0 cos(πt) 0 1 2cos(π(x+t)) +1 2cos(π(x−t)) 1 5 2x5+ 2x3−2x2 −(20x3+ 6x−2)·0−x5−x3+x20t(t−1) 12 (x5+x3−x2)· (t2−t) (t2−t) 9 Conclusions El m`etode 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´o de Laplace, la qual hem fet servir per il·lustrar aquest m`etode. Una alternativa ´es aplicar el m`etode de les falses transicions, afegint una derivada temporal i discretitzant les variables espacials. Si estem al pla xy, l’error d’aquest m`etode alternatiu ´es de l’ordre de ∆2 x+ ∆2 y. Per tant, ´es millor aplicar el m`etode semi-discret (o semi-anal´ıtic), que ´es el que hem aplicat, i en el qual hem vist que l’error ´es de l’ordre de ∆2 x. 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 [11]. Finalment, proposem el seg¨uent problema obert: Estudiar les semblances i les difer`encies entre el MOL i la integraci´o al llarg de les rectes caracter´ıstiques de l’equaci´o 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`onim, els quals han contribu¨ıt considerablement a millorar aquest treball. Refer`encies [1] 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. 20 [2] M. C´amara, J. F`abrega, 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. [3] T. S. Chihara, An Introduction to Orthogonal Polynomials, Gordon and Breach, New York, 1978. [4] H. Diestel, Analysis of planar multiconductor transmission-line systems with the method of lines, Arch. Elektron. Uebertragungstech 41 (1987) 169–175. [5] C. D. Godsil, Algebraic Combinatorics, Chapman and Hall, New York, 1993. [6] 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. [7] W. E. Schiesser, The Numerical Method of Lines, Academic Press, New York, 1991. [8] 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. [9] V. R. Subramanian, R. E. White, Semianalytical method of lines for solving elliptic partial differential equations, Chem. Eng. Sci. 59 (2004) 781–788. [10] 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. [11] 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. [12] 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. [13] E. Zuazua, Introducci´on al an´alisis num´erico de ecuaciones en derivadas parciales de evoluci´on, Universidad Aut´onoma de Madrid, 2003. [14] E. Zuazua, M´etodos num´ericos de resoluci´on de ecuaciones en derivadas parciales, Basque Center for Applied Mathematics (BCAM), 2009.