Clase 15: Ecuaciones diferenciales parciales

Bibliografía de la Clase:

  • Introduction to Numerical Analysis - J.Stoer, R.Bulirsch. Chapter 7: Ordinary Differential Equations.
  • Análisis Numérico - L. Burden, J Faires. Capítulo 12: Soluciones numéricas para las ecuaciones diferenciales parciales
  • Guía profesor Luis Salinas 396

Introducción a EDP

Una ecuación diferencial parcial (EDP) es una ecuación que contiene una función multivariable desconocida (ej. u(x,y)) y sus derivadas parciales. Este tipo de ecuaciones describen varios fenómenos físicos como el calor, el sonido, dinámica de fluidos, etc.

El orden de la EDP está determinado por el grado de la derivada parcial más alto.

Una EDP de segundo orden genérica tiene la forma:

A(x,y)\frac{\partial^2 u}{\partial x^2}(x,y) +
2B(x,y)\frac{\partial^2 u}{\partial x \partial y}(x,y)+
C(x,y)\frac{\partial^2 u}{\partial y^2}(x,y) + \\
a(x,y)\frac{\partial u}{\partial x}(x,y) +
b(x,y)\frac{\partial u}{\partial y}(x,y) +
c(x,y)u(x,y)
= f(x,y)

en una cierta región \Omega \subset \mathbb{R}^2. Se dice que la EDP es de tipo:

  • Hiperbólica: si \Delta=B^2 -AC > 0 en \Omega.
  • Parabólica: si \Delta=B^2-AC=0 en \Omega.
  • Elíptica: si \Delta=B^2-AC < 0 en \Omega.

Ecuaciones diferenciales Elípticas

Ejemplos de EDP elípticas son la ecuación de Poisson y Laplace:

Ecuación de Poisson

Usada en el estudio de problemas físicos como distribución de calor, energía potencial sobre un punto y dinámica de fluidos.

(1)\frac{\partial^2 u}{\partial x^2}(x,y) + \frac{\partial^2
     u}{\partial y^2}(x,y) = f(x,y)

donde f(x,y) es una función definida sobre una región R con borde S.

La ecuación de Poison es del tipo elíptica ya que A=1,B=0 y C=1, por lo tanto \Delta=0-1\times1=-1<0.

Esta ecuación es equivalente a:

\Delta u(x,y) = f(x,y)

donde \Delta es el operador de Laplace o Laplaciano:

\Delta = \frac{\partial^2}{\partial x^2} +
\frac{\partial^2}{\partial y^2}

Ecuación de Laplace

es un caso particular de la ecuación de Poisson donde f(x,y)=0 que significa que la región R es plana.

\Delta u(x,y) = 0

Si el valor de u(x,y) en la región R es determinada por el valor de la función en el borde S se les denomina condiciones de borde de Dirichlet que viene dado por:

u(x,y) = g(x,y), \qquad \forall(x,y) \in S

caption

En la figura, la región R está definida como:

R = \{(x,y) \in \mathbb{R}^2: a \leq x \leq b, c \leq y \leq
d\}
\subset \mathbb{R}^2, \qquad a > 0

Métodos numéricos para EDP

Para resolver numéricamente este tipo de ecuaciones se utiliza una adaptación del método de diferencias finitas. El método consiste en generar una grilla equiespaciada en la región R para obtener aproximaciones de la función en los puntos de la grilla u(x_i,y_i) iterativamente.

El número de puntos de la grilla para x e y son n y m respectivamente, por lo tanto el espaciado entre los puntos es:

h = \frac{b-a}{n}  \qquad \text{equiespaciado para x} \\
k = \frac{d-c}{m}  \qquad \text{equiespaciado para y}

caption

Para expresar numéricamente y'(x_i) e y''(x_i) se utiliza la serie de Taylor alrededor de x_i. De manera aproximada se tiene:

y(x_{i+1}) = y(x_i+h) \approx
y(x_i)+hy'(x_i)+\frac{h^2}{2}y''(x_i)+
\frac{h^3}{6}y'''(x_i)
\\
y(x_{i-1}) = y(x_i-h) \approx
y(x_i)-hy'(x_i)+\frac{h^2}{2}y''(x_i)-
\frac{h^3}{6}y'''(x_i)

y luego se suman para obtener una expresión para y''(x_i):

y(x_{i+1})+y(x_{i-1}) & \approx 2y(x_i)+h^2 y''(x_i) \\
y''(x_i) & \approx  \frac{y(x_{i+1})-2y(x_i)+y(x_{i-1})}{h^2}

Aproximaciones usando series de Taylor

Podemos usar las expansiones de Taylor mostradas anteriormente para resolver la ecuación de Poisson expresando:

\frac{\partial^2 u}{\partial x^2}(x_i,y_j) \approx
\frac{u(x_{i+1},y_j)-2u(x_i,y_j)+u(x_{i-1},y_j)}{h^2}
\\
\frac{\partial^2 u}{\partial y^2}(x_i,y_j) \approx
\frac{u(x_i,y_{j+1})-2u(x_i,y_j)+u(x_i,y_{j-1})}{k^2}

ahora podemos expresar la ecuación (1) en el punto (x_i,y_j) aproximadamente como:

\frac{u(x_{i+1},y_j)-2u(x_i,y_j)+u(x_{i-1},y_j)}{h^2} &+ \\
\frac{u(x_i,y_{j+1})-2u(x_i,y_j)+u(x_i,y_{j-1})}{k^2} & \approx
f(x_i,y_j)

Agrupando términos, y expresando u_{i,j} \approx u(x_i,y_j) se tiene:

(2)\boxed{
     2\Big[ \Big( \frac{h}{k}\Big)^2+1\Big] u_{i,j} -
     u_{i+1,j}-u_{i-1,j} -
     \Big( \frac{h}{k}\Big)^2 \big( u_{i,j+1} + u_{i,j-1} \big)
     = -h^2 f(x_i,y_j)}

Para definir las condiciones de borde, se debe especificar que el valor de u(x,y) en la frontera está definido por la función g(x,y).

(3)u_{0,j} = g(x_0,y_j) \quad \text{and} \quad u_{n,j} =
     g(x_n,y_j) \qquad \forall j=0,1,\dots,m \\
     u_{i,0} = g(x_i,y_0) \quad \text{and} \quad u_{i,m} =
     g(x_i,y_m) \qquad \forall i=1,2,\dots,n-1 \\

Las ecuaciones (2) y (3) definen un sistema de ecuaciones que puede ser representado matricialmente.

Ejercicio en clases

Determine la distribución del calor en estado estable en una superficie metálica de dimensiones 0.5m \times 0.5m utilizando una grilla de tamaño n=m=4. Las condiciones de borde en los ejes x e y son 0 y el calor en los otros dos bordes crece linealmente desde 0 a 100 grados.

Solución

Dado que n=m=4 se tiene que:

h=k=\frac{0.5}{4}=\frac{1}{8}

y que f(x,y)=0 ya que es una placa 2D:

(4)4 u_{i,j}-u_{i+1,j}-u_{i-1,j} -
     u_{i,j+1}-u_{i,j-1}
     = 0 \qquad \forall i=1,2,3 \quad , j=1,2,3

Sus condiciones de frontera son:

u_{0,j} = 0 \qquad \forall j=0,1,2,3 \\
u_{i,0} = 0 \qquad \forall i=0,1,2,3 \\
u_{i,4} = 200 x_i \qquad \forall i=1,2,3,4 \\
u_{4,j} = 200 y_j \qquad \forall j=1,2,3,4

Por lo tanto se tiene que

u_{0,1}=u_{0,2}=u_{0,3}=u_{0,4} = 0 \\
u_{0,0}=u_{1,0}=u_{2,0}=u_{3,0}=u_{4,0}=0 \\
u_{1,4}=u_{4,1}=25 \\
u_{2,4}=u_{4,2}=50\\
u_{3,4}=u_{4,3}=75\\
u_{4,4}=100

Para el resto de los puntos se cumple el sistema de ecuaciones (4).

Para i=1,j=1

4u_{1,1}-u_{2,1}-u_{0,1}-u_{1,2}-u_{1,0} = 0 \\
4u_{1,1}-u_{2,1}-u_{1,2}=u_{0,1}-u_{1,0} \\
4u_{1,1}-u_{2,1}-u_{1,2}=0 \\

para i=1,j=2:

4u_{1,2}-u_{2,2}-u_{1,3}-u_{1,1}=u_{0,2} \\
4u_{1,2}-u_{2,2}-u_{1,3}-u_{1,1}=0 \\

para i=1,j=3:

4u_{1,3}-u_{2,3}-u_{1,2} = u_{0,3}+u_{1,4}\\
4u_{1,3}-u_{2,3}-u_{1,2} = 25\\

para i=2,j=1

4u_{2,1}-u_{3,1}-u_{1,1}-u_{2,2}-u_{2,0}=0 \\
4u_{2,1}-u_{3,1}-u_{1,1}-u_{2,2}=0

para i=2,j=2

4u_{2,2}-u_{3,2}-u_{1,2}-u_{2,3}-u_{2,1}=0

para i=2,j=3

4u_{2,3}-u_{3,3}-u_{1,3}-u_{2,4}-u_{2,2}=0 \\
4u_{2,3}-u_{3,3}-u_{1,3}-u_{2,2}=u_{2,4} \\
4u_{2,3}-u_{3,3}-u_{1,3}-u_{2,2}=50 \\

para i=3,j=1

4u_{3,1}-u_{4,1}-u_{2,1}-u_{3,2}-u_{3,0}=0 \\
4u_{3,1}-u_{2,1}-u_{3,2}=u_{4,1}+u_{3,0} \\
4u_{3,1}-u_{2,1}-u_{3,2}=25 \\

para i=3,j=2

4u_{3,2}-u_{4,2}-u_{2,2}-u_{3,3}-u_{3,1}=0 \\
4u_{3,2}-u_{2,2}-u_{3,3}-u_{3,1}=u_{4,2} \\
4u_{3,2}-u_{2,2}-u_{3,3}-u_{3,1}=50

para i=3,j=3

4u_{3,3}-u_{4,3}-u_{2,3}-u_{3,4}-u_{3,2}=0 \\
4u_{3,3}-u_{2,3}-u_{3,2}=u_{4,3} + u_{3,4}\\
4u_{3,3}-u_{2,3}-u_{3,2}=150

Matricialmente se expresa como:

(5)\begin{bmatrix}
     4 & -1 & 0 & -1 & 0 & 0 & 0 & 0 & 0 \\
     -1 & 4 & -1 & 0 & -1 & 0 & 0 & 0 & 0 \\
     0 & -1 & 4 & 0 & 0 & -1 & 0 & 0 & 0 \\
     -1 & 0 & 0 & 4 & -1 & 0 & -1 & 0 & 0 \\
     0 & -1 & 0 & -1 & 4 & -1 & 0 & -1 & 0 \\
     0 & 0 & -1 & 0 & -1 & 4 & 0 & 0 & -1 \\
     0 & 0 & 0 & -1 & 0 & 0 & 4 & -1 & 0 \\
     0 & 0 & 0 & 0 & -1 & 0 & -1 & 4 & -1 \\
     0 & 0 & 0 & 0 & 0 & -1 & 0 & -1 & 4 \\
     \end{bmatrix}
     \begin{bmatrix}
     u_{1,1} \\
     u_{1,2} \\
     u_{1,3} \\
     u_{2,1} \\
     u_{2,2} \\
     u_{2,3} \\
     u_{3,1} \\
     u_{3,2} \\
     u_{3,3}
     \end{bmatrix}
     =
     \begin{bmatrix}
     0\\
     0\\
     25\\
     0\\
     0\\
     50\\
     25\\
     50\\
     150
     \end{bmatrix}

Ejercicio propuesto 1

Estudie cómo se utiliza el método de Gauss Seidel para encontrar los u_(x_i,y_i) y encuentre la solución de la ecuación (5). Encuentre el error asociado al método numérico considerando que la solucíón analítica es u(x,y)=400xy.

Ejercicio propuesto 2

En un recinto cuadrado de lado 1 considere la ecuación de Poisson:

\Delta u = -0.0625

con condición de borde igual a cero en los cuatro costados. Encuentre la solución de este problema usando el método de diferencias finitas, considerando h=k=1/4.