El m\`etode de les l\'inies per a la resoluci\'o num\`erica d'equacions en derivades parcials. The method of lines for numerical solutions of partial differential equations

TL;DR

This paper introduces the method of lines (MOL), a semi-discrete numerical technique for solving partial differential equations, demonstrated through the Laplace equation, showing high accuracy compared to analytical solutions.

Contribution

It presents a detailed description of the MOL for PDEs and compares its effectiveness with the classical variable separation method.

Findings

MOL provides accurate approximations of analytical solutions.

The method discretizes all but one variable, simplifying PDE solving.

Results show MOL's effectiveness for Laplace equation in Cartesian coordinates.

Abstract

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.

Tables

Table 1. Taula 1: Components de la matriu T 𝑇 T i valors propis λ k = ω k 2 subscript 𝜆 𝑘 superscript subscript 𝜔 𝑘 2 \lambda_{k}=\omega_{k}^{2} segons les condicions de contorn a x = 0 𝑥 0 x=0 i x = a 𝑥 𝑎 x=a .

Esquerra	Dreta	$T_{jk}$	${\lambda }_k={\omega }_k^2$
Dirichlet	Dirichlet	$\sqrt{\frac{2}{N+1}}\sin \left(\frac{jk\pi }{N+1}\right)$	$4{\sin }^2\left(\frac{k\pi }{2N+2}\right)$
Dirichlet	Neumann	$\sqrt{\frac{2}{N+0.5}}\sin \left(\frac{j\left(k-0.5\right)\pi }{N+0.5}\right)$	$4{\sin }^2\left(\frac{\left(k-0.5\right)\pi }{2N+1}\right)$
Neumann	Dirichlet	$\sqrt{\frac{2}{N+0.5}}\cos \left(\frac{\left(j-0.5\right)\left(k-0.5\right)\pi }{N+0.5}\right)$	$4{\sin }^2\left(\frac{\left(k-0.5\right)\pi }{2N+1}\right)$
Neumann	Neumann	$\sqrt{\frac{2}{N}}\cos \left(\frac{\left(j-0.5\right)\left(k-1\right)\pi }{N}\right)$ si $k>1$	$4{\sin }^2\left(\frac{\left(k-1\right)\pi }{2N}\right)$
		$\frac{1}{\sqrt{N}}$ si $k=1$

L	T	$q(x,t)$	$f_1\left(x\right)$	$f_2\left(x\right)$	$g_1\left(t\right)$	$g_2\left(t\right)$	$v(x,t)$ analítica
1	4	$0$	$0$	$\pi \cos \left(\pi x\right)$	$\sin \left(\pi t\right)$	$0$	$\cos \left(\pi x\right)\sin \left(\pi t\right)$
1	5	$-2\left(x-t\right)e^{-x-t}$	$0$	$x{e}^{-x}$	$0$	$-2t{e}^{-t-1}$	$xt{e}^{-x-t}$
			$0$		$0$	$+t{e}^{-t}$	$xt{e}^{-x-t}$
1	4	$0$	$\cos \left(\pi x\right)$	$0$	$\cos \left(\pi t\right)$	$0$	$\frac{1}{2}\cos \left(\pi \left(x+t\right)\right)$
							$+\frac{1}{2}\cos \left(\pi \left(x-t\right)\right)$
1	5	$2x^5+2x^3-2x^2$
		$-(20x^3+6x-2)\cdot$	$0$	$-x^5-x^3+x^2$	$0$	$\frac{t\left(t-1\right)}{12}$	$(x^5+x^3-x^2)\cdot$
		$\left(t^2-t\right)$					$\left(t^2-t\right)$

Equations

$${\nabla }^{2}V$$

$$V(0,y)$$

$$V(x,0)$$

$$V^{\prime }(x_j)\simeq \frac{V(x_{j+1})-V(x_j)}{h},$$

$$V^{\prime \prime }(x_j)\simeq \frac{V^{\prime }(x_j)-V^{\prime }(x_{j-1})}{h}\simeq \frac{V(x_{j+1})-2V(x_j)+V(x_{j-1})}{h^2},$$

$$\frac{{\partial }^{2}V}{\partial x^2}\simeq \frac{V_{j+1}-2V_j+V_{j-1}}{h^2},$$

$$\frac{{\partial }^{2}V}{\partial y^2}+\frac{1}{h^2}(V_{j+1}-2V_j+V_{j-1})=0,ji=1,\dots ,N.$$

$$\frac{\partial^{2}\mbox{\boldmath$V$}}{\partial y^{2}}-\frac{1}{h^{2}}\mbox{\boldmath$P$}\mbox{\boldmath$V$}=\mbox{\boldmath$0$},$$

$${\normalsize{\left(\begin{array}[]{c}\frac{\partial^{2}V_{1}}{\partial y^{2}}\\ \frac{\partial^{2}V_{2}}{\partial y^{2}}\\ \vdots\\ \frac{\partial^{2}V_{N-1}}{\partial y^{2}}\\ \frac{\partial^{2}V_{N}}{\partial y^{2}}\\ \end{array}\right)-\frac{1}{h^{2}}\left(\begin{array}[]{ccccccc}p_{l}&-1&0&\cdots&0&0&0\\ -1&2&-1&\cdots&0&0&0\\ 0&-1&2&\cdots&0&0&0\\ &&&\ddots&&&\\ 0&0&0&\cdots&-1&2&-1\\ 0&0&0&\cdots&0&-1&p_{r}\\ \end{array}\right)\left(\begin{array}[]{c}V_{1}\\ V_{2}\\ \vdots\\ V_{N-1}\\ V_{N}\\ \end{array}\right)=\mbox{\boldmath$0$},}}$$

$$\frac{\partial^{2}\mbox{\boldmath$V$}}{\partial y^{2}}-\frac{1}{h^{2}}\mbox{\boldmath$T$}\mbox{\boldmath$D$}\mbox{\boldmath$T$}^{\top}\mbox{\boldmath$V$}=\mbox{\boldmath$0$},$$

$$\frac{\partial^{2}\mbox{\boldmath$T$}^{\top}\mbox{\boldmath$V$}}{\partial y^{2}}-\frac{1}{h^{2}}\mbox{\boldmath$D$}\mbox{\boldmath$T$}^{\top}\mbox{\boldmath$V$}=\mbox{\boldmath$0$},$$

$$\left(\frac{\partial^{2}}{\partial y^{2}}-\frac{1}{h^{2}}\mbox{\boldmath$D$}\right)\overline{\mbox{\boldmath$V$}}=\mbox{\boldmath$0$}.$$

$${\overline{V}}_j^{\prime \prime }=\frac{1}{h^2}{\omega }_j^2,j=1,\dots ,N,$$

$${\overline{V}}_j={\alpha }_j{e}^{\frac{{\omega }_j}{h}y}+{\beta }_j{e}^{-\frac{{\omega }_j}{h}y}=A_j\cosh \frac{{\omega }_{j}y}{h}+B_j\sinh \frac{{\omega }_{j}y}{h},j=1,\dots ,N.$$

$$\mbox{\boldmath$V$}=\mbox{\boldmath$T$}\overline{\mbox{\boldmath$V$}}.$$

$$\left(\begin{array}[]{ccccccc}p_{l}&-1&0&\cdots&0&0&0\\ -1&2&-1&\cdots&0&0&0\\ 0&-1&2&\cdots&0&0&0\\ &&&\ddots&&&\\ 0&0&0&\cdots&-1&2&-1\\ 0&0&0&\cdots&0&-1&p_{r}\\ \end{array}\right)\left(\begin{array}[]{c}v_{1}\\ v_{2}\\ \vdots\\ v_{N-1}\\ v_{N}\end{array}\right)=\lambda\left(\begin{array}[]{c}v_{1}\\ v_{2}\\ \vdots\\ v_{N-1}\\ v_{N}\end{array}\right),$$

$$(p_l-\lambda )v_1$$

$$(2-\lambda )v_j$$

$$(p_r-\lambda )v_N$$

$$e^{i\alpha }=(2-\lambda )-e^{-i\alpha }\Rightarrow 2-\lambda =2\cos \alpha .$$

$$\lambda =2(1-\cos \alpha )=4{\sin }^2\left(\frac{\alpha }{2}\right),$$

$${\alpha }_{1,2}=\pm 2\arcsin \left(\frac{\sqrt{\lambda }}{2}\right).$$

$$v_j=A{e}^{ij\alpha }+B{e}^{-ij\alpha }.$$

$$(2-\lambda )v_1-v_2$$

$$\Rightarrow A+B=0,$$

$$v_j=\sin (j\alpha ),j=1,\dots ,N.$$

$$(2-\lambda )v_1-v_2$$

$$\Rightarrow \sin ((N+1)\alpha )=0.$$

$${\lambda }_k=4{\sin }^2\left(\frac{k\pi }{2(N+1)}\right),k=1,\dots ,N,$$

$$\mbox{\boldmath$v$}_{k}=(v_{jk})=\left(\sin\left(\frac{k\pi}{N+1}\right),\sin\left(\frac{2k\pi}{N+1}\right),\ldots,\sin\left(\frac{Nk\pi}{N+1}\right)\right)^{\top}.$$

$$\|\mbox{\boldmath$v$}_{k}\|^{2}=\sum_{j=1}^{N}\sin^{2}\left(\frac{jk\pi}{N+1}\right)=\frac{1}{2}\sum_{j=1}^{N}\left(1-\cos\left(\frac{2jk\pi}{N+1}\right)\right)=\cdots=\frac{N+1}{2},$$

$$T_{jk}=(v_{jk})=\sqrt{\frac{2}{N+1}}\sin \left(\frac{jk\pi }{N+1}\right),j,k=1,\dots ,N,$$

$$p_{k+1}(x)=(2-x)p_k(x)-p_{k-1}(x),k=1,\dots ,N-2,$$

Full text

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.

††

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.

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.

1 Introducció

El mètode de les línies, també anomenat mètode de semi-discretització per Zuazua [13] é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 [12] i Verwer i Sanz-Serna [10]).

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 [1], Diestel [4], Sadiku i Obiozor [6], o Shakeri i Dehghan [8].

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. 2. 2.

Obtenció d’un sistema d’equacions diferencials ordinàries acoblades. 3. 3.

Transformació per obtenir un sistema d’equacions diferencials ordinàries desacoblades. 4. 4.

Resolució de les equacions diferencials. 5. 5.

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

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 (el.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.

2 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 [9]).

2.1 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\times b$.

[TABLE]

amb condicions de frontera que, per simplicitat, establim com

[TABLE]

El primer pas és discretitzar la variable $x\in(0,a)$ en els $N$ punts $x_{j}=j\Delta x$, per a $j=1,2,\ldots,N$, on el “pas” $h=\Delta x=\frac{a}{N+1}$ correspon a la distància entre dues de les $N$ línies resultants de la discretització, segons mostrem a la Figura 1.

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

[TABLE]

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

[TABLE]

amb error $|\epsilon_{j}|\leq Ch^{2}$, per a una $C$ constant, que correspon a un mètode de segon ordre. Amb notació simplificada,

[TABLE]

on $V_{j}$ denota la funció $V(x_{j},y)$. Si substituïm (4) en (1), obtenim les $N$ equacions diferencials

[TABLE]

Aleshores, si reemplacem $V(x,y)$ pel vector columna \mbox{\boldmathV}=(V_{1},\ldots,V_{N})^{\top}, podem escriure (5) en forma matricial com

[TABLE]

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

[TABLE]

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 (6) 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 $P$. En aquest sentit, notem que $P$ és simètrica i, per tant, els seus valors propis $\lambda_{1},\ldots,\lambda_{N}$ són reals. A més, el teorema de Gershgorin111Recordem que el Teorema de Gershgorin afirma que, donada una matriu \mbox{\boldmathA}=(a_{jk})\in M_{n}(\mathbb{C}), si es defineixen els cercles $D_{1},\ldots,D_{n}$ amb centre $a_{jj}$ i radi $r_{j}=\sum_{k\neq 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 $|\lambda_{k}-2|\leq 2$, de manera que utilitzarem la notació $\lambda_{k}=\omega_{k}^{2}$ per a $k=1,\ldots,N$. De fet, com s’explica a la subsecció 3.2, la matriu tridiagonal $P$ 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 [5, §8,5]).

Per tant, tenim que \mbox{\boldmathT}^{\top}\mbox{\boldmathP}\mbox{\boldmathT}=\mbox{\boldmathD}, on $T$ és la matriu ortogonal (amb columnes ortonormals corresponents als vectors propis de $P$), és a dir, \mbox{\boldmathT}^{-1}=\mbox{\boldmathT}^{\top}, i $D$ és la matriu diagonal dels valors propis de $P$. Aleshores, \mbox{\boldmathP}=\mbox{\boldmathT}\mbox{\boldmathD}\mbox{\boldmathT}^{\top} que, si substituïm a (6), ens dona

[TABLE]

i, si multipliquem per \mbox{\boldmathT}^{\top}(=\mbox{\boldmathT}^{-1}),

[TABLE]

i, si utilitzem el potencial transformat \overline{\mbox{\boldmathV}}:=\mbox{\boldmathT}^{\top}\mbox{\boldmathV}, obtenim

[TABLE]

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

[TABLE]

amb solucions (el polinomi característic és $s^{2}=\frac{\omega_{j}^{2}}{h^{2}}$, amb arrels $s_{1,2}=\pm\frac{\omega_{j}}{h}$)

[TABLE]

Per trobar el vector de potencials per a cada $j=1,\ldots,N$, només cal realitzar la transformada inversa de \overline{\mbox{\boldmathV}}, o sigui,

[TABLE]

Finalment, imposem les condicions de contorn i resolem les equacions resultants per determinar les constants $\alpha_{j},\beta_{j}$, o $A_{j},B_{j}$.

3 Diagonalització de la matriu $P$ i ortogonalitat

3.1 Valors i vectors propis

Abans de resoldre un exemple concret, en aquest apartat discutim la diagonalització de la matriu $P$ i les propietats d’ortogonalitat involucrades. Els resultats obtinguts, els valors propis $\lambda_{k}=\omega_{k}^{2}$ i les matrius de vectors propis \mbox{\boldmathT}=(T_{jk}), es mostren a la Taula 1, 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

[TABLE]

d’on

[TABLE]

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 (12), podríem usar la solució $v_{j}=s^{j}$, però s’obté un resultat més compacte amb $v_{j}=e^{ij\alpha}$ (és a dir, agafant $s=e^{i\alpha}=\cos\alpha+i\sin\alpha$). Si substituïm aquest valor a (12), i dividim per $e^{i(j-1)\alpha}$ s’obté

[TABLE]

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

[TABLE]

cosa que significa que, per a cada valor propi $\lambda$, tenim dos possibles valors de $\alpha$:

[TABLE]

Així, amb $\alpha=\alpha_{1}$, la solució general de la recurrència és

[TABLE]

Per trobar els possibles valors de les constants $A$, $B$, i $\alpha$, imposem les “condicions inicials” (11) i (13). 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

[TABLE]

on hem aplicat $\cos\alpha=\frac{e^{i\alpha}+e^{-i\alpha}}{2}$ i $\sin\alpha=\frac{e^{i\alpha}-e^{-i\alpha}}{2i}$. Per tant, $B=-A$ i, llevat d’una constant multiplicativa, podem suposar que

[TABLE]

D’altra banda, de (13) s’obté

[TABLE]

Per tant, els possibles valors de $\alpha$ són $\alpha_{k}=\frac{k\pi}{N+1}$, per a $k=1,\ldots,N$, amb el que obtenim els valors propis

[TABLE]

els quals compleixen que $\lambda_{1}<\lambda_{2}<\cdots<\lambda_{N}$, amb els respectius vectors propis

[TABLE]

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

[TABLE]

on hem usat que $\sin^{2}x=\frac{1}{2}(1-\cos(2x))$ i la suma d’una sèrie geomètrica $\sum_{j=1}^{N}r^{j}=\frac{r(1-r^{N})}{1-r}$. Concloem que els elements de la matriu $T$ són

[TABLE]

com s’indica a la Taula 1.

3.2 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 \mbox{\boldmathP}=(p_{jk}) satisfan $p_{j,j+1}p_{j+1,j}>0$ per a $j=1,\ldots,N$ i, per tant, donen lloc a una família de polinomis ortogonals $p_{0},p_{1},\ldots,p_{N-1}$ de variable discreta222És a dir, polinomis ortogonals respecte d’un producte escalar discret en els punts d’una malla $x_{1},\ldots,x_{N}\in\mathbb{R}$: $\langle f,g\rangle=\sum_{j=1}^{N}w_{j}f(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 [2]. Aquests polinomis constitueixen una successió de Sturm (vegeu Godsil [5, §8.5] or Chihara [3]) i satisfan una relació de recurrència de tres termes (com a l’expressió (12))

[TABLE]

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 $P$ i, per tant, les seves arrels són els valors propis $\lambda_{k}$, amb vectors propis

[TABLE]

trobats anteriorment a (17) i (18) (aquests últims llevat d’una constant multiplicativa). L’ortogonalitat esmentada s’obté respecte d’un producte escalar discret en els punts de la malla $\lambda_{1},\lambda_{2},\ldots,\lambda_{N}$, amb $\lambda_{1}<\lambda_{2}<\cdots<\lambda_{N}$. Per tant, un cop que normalitzem els vectors \mbox{\boldmathv}_{k}, és a dir \overline{\mbox{\boldmathv}}_{k}=\mbox{\boldmathv}_{k}/\|\mbox{\boldmathv}_{k}\|, obtenim de nou els elements de la matriu $T$ donats a (19):

[TABLE]

és a dir,

[TABLE]

(Compareu amb les funcions pròpies que apareixen en el cas analític, vegeu la secció 5). En aquest context, observem que el potencial transformat \overline{\mbox{\boldmathV}}=\mbox{\boldmathT}^{\top}\mbox{\boldmathV} 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 $T$).

Per exemple, en el cas de la matriu $P$ a (10) amb $N=5$ i $p_{l}=p_{r}=2$, s’obtenen els polinomis

[TABLE]

Els valors propis de $P$ són (arrels de $p_{5}(x)$):

[TABLE]

A la Figura 2 es mostren les gràfiques de $p_{j}(x)$, per a $j=0,\ldots,4$. Noteu les paritats (parells i senars) d’aquests polinomis respecte de $x=\lambda_{3}=2$, induïdes per les corresponents simetries de la matriu

[TABLE]

4 Un exemple

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

Volem resoldre l’equació de Laplace (potencial)

[TABLE]

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 3. També volem calcular el valor del potencial en els punts indicats en aquesta figura.

La matriu $P$ té valors propis (VAPS) $\omega_{k}^{2}=4\sin^{2}\frac{k\pi}{2N+2}$, per a $k=1,2,\ldots,N$, i la matriu de vectors propis (VEPS) $T$ té components $T_{jk}=\sqrt{\frac{2}{N+1}}\sin\frac{jk\pi}{N+1}$ (vegeu la Taula 1).

Si escollim $N=15$, resulta que $h=\triangle x=\frac{1}{16}$, i les abscisses dels punts que volem calcular són $x_{4}=\frac{1}{4}$, $x_{8}=\frac{1}{2}$ i $x_{12}=\frac{3}{4}$.

De $V_{j}=\sum_{k=1}^{N}T_{jk}\overline{V}_{k}$, per a $i=1,2,\ldots,N$, obtenim la solució

[TABLE]

Ara, per determinar les constants $A_{j}$ i $B_{j}$, imposem les condicions de contorn:

• •

$V(x,0)=0\Rightarrow$

[TABLE]

d’on $A_{j}=0$, per a $j=1,2,\ldots,N$.

• •

$V(x,1)=100\Rightarrow$

[TABLE]

d’on

[TABLE]

Per tant,

[TABLE]

Així, arribem a la solució

[TABLE]

la qual està representada a la Figura 4.

En els punts demanats els valors que s’obtenen (calculats amb el programa de la Secció 7) són

[TABLE]

que, comparats amb els valors calculats aplicant els resultats de la Secció 5 també amb el programa de la Secció 7,

[TABLE]

demostren ser molt bones aproximacions.

5 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ó

[TABLE]

per a $x\in(0,a)$ i $y\in(0,b)$, amb condicions de contorn

[TABLE]

com es veu a la Figura 5.

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

[TABLE]

Imposant les condicions de contorn homogènies, resulta que $\omega=\omega_{k}=\frac{k\pi}{a}$ per a $k=1,2,\ldots$ i obtenim les funcions

[TABLE]

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

[TABLE]

Considerem ara una solució en forma de sèrie

[TABLE]

Si imposem la condició de contorn no homogènia

[TABLE]

i considerem que $\{\sin\left(\frac{k\pi}{a}x\right):k=1,2,\ldots\}$ és una base ortogonal a $(0,a)$, obtenim el corresponent desenvolupament en sèrie de Fourier de $f(x)$, amb coeficients

[TABLE]

Per tant, la solució analítica és

[TABLE]

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 $\alpha_{n}=\frac{400}{n\pi}$ per a $n$ senar i $\alpha_{n}=0$ per a $n$ parell. Aleshores, la solució analítica és

[TABLE]

la representació gràfica de la qual es mostra a la Figura 6.

6 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 [7], o Zuazua [13, 14]).

Nosaltres estudiem primer el comportament de l’error en el cas de l’exemple anterior (Secció 4). Vegeu la Figura 7. Observem que tenim 8 funcions d’error (per a $x_{j}=j/16$, $j=1,\ldots,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(\pi)\sin(\pi x)$. Aleshores, per trobar la solució analítica no cal fer l’últim pas del desenvolupament en sèrie de Fourier, i s’obté

[TABLE]

la gràfica de la qual es mostra a la Figura 8. Pel que fa al mètode de línies, comentem primer que aquest cas va ser resolt computacionalment per Subramanian i White [9] 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

[TABLE]

on $x_{j}=jh$, $j=1,\ldots,N$, i $\operatorname{sinc}(x)=\frac{\sin(\pi x)}{\pi x}$. Vegeu la Figura 9. Notem que, com $\lim_{h\rightarrow 0}\operatorname{sinc}(h)=1$, quan $h\rightarrow 0$, $j\rightarrow\infty$, i $jh\rightarrow x$, i aleshores $V(x_{j},y)\rightarrow 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 $\epsilon=V(x,y)-V(x_{j},y)$ és de l’ordre de $O(h^{2})$, com ja s’ha comentat.

7 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)pi0.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)=Ssin(JK*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

8 Altres exemples resolubles amb el MOL

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

L’equació de Laplace en la regió quadrada $1\times 1$ amb condicions de contorn:

• $(a)$

Dirichlet-Neumann.

• $(b)$

Neumann-Dirichlet.

• $(c)$

Neumann-Neumann.

(Els valors a utilitzar de la matriu $T$ i els valors propis $\omega_{k}^{2}$ estan a la Taula 1).

• 2.

L’equació de Laplace en coordenades cilíndriques $\rho^{2}\frac{\partial^{2}V}{\partial\rho^{2}}+\rho\frac{\partial V}{\partial\rho}+\frac{\partial^{2}\phi}{\partial\phi^{2}}=0$, amb les condicions de contorn $V(\rho,\phi)=0$ per a $0\leq\rho\leq a$ i $V(\rho,\phi)=V_{0}$ per a $a<\rho\leq b$, amb la discretització de la coordenada angular $\phi$. La solució analítica és $V(\rho)=V_{0}\frac{\ln\left(\frac{\rho}{a}\right)}{\ln\left(\frac{b}{a}\right)}$.

• 3.

L’equació d’ona $\frac{\partial^{2}v}{\partial t^{2}}-\frac{\partial^{2}v}{\partial x^{2}}=q(x,t)$, amb $0\leq x\leq L$ i $0\leq t\leq T$, les condicions inicials $v(x,0)=f_{1}(x)$ i $\frac{\partial v}{\partial t}(x,0)=f_{2}(x)$ per a $0\leq x\leq L$, la condició de contorn de Dirichlet $v(0,t)=g_{1}(x)$ i la condició no local $\int_{0}^{L}v(x,t)dx=g_{2}(x)$ per a $0\leq t\leq T$. A la taula següent hi ha alguns casos particulars amb la corresponent solució analítica.

9 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 $\Delta_{x}^{2}+\Delta_{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 $\Delta_{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 [11]. 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 $\frac{\partial V}{\partial t}+c\frac{\partial V}{\partial 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.

Referències

• [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.

• [2]

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$.

• [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ón al análisis numérico de ecuaciones en derivadas parciales de evolución, Universidad Autónoma de Madrid, 2003.

• [14]

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