4.5. Práctica Grupal: Problemas de Contorno#

Un problema de contorno es una ecuación diferencial donde se imponen unos valores que debe verificar la función solución en puntos situados en la frontera (contorno) del intervalo de definición de la solución.

Por ejemplo: Determinar que función de una variable \(U(x)\) verifica que su derivada segunda más la propia función vale cero, dentro del intervalo \([0, \frac{3\pi}{2}]\).

Es decir, estamos buscando \(U(x)\) tal que \(U''(x) + U(x) = 0\) o lo que es lo mismo \(U''(x) = -U(x)\) en \([0, \frac{3\pi}{2}]\).

¿Se os ocurre alguna función que cumpla lo anterior?

Por ejemplo \(U(x) = cos(x)\) ó \(U(x) = sin(x)\), \(U(x) = cos(x) + sin(x) \)… hay infinitas soluciones. Si añadimos las condiciones de contorno \(U(0) = 2\); \(U(\frac{3\pi}{2})= 2\), tenemos el siguiente problema de contorno:

Hallar \(U(x) \in \mathcal{C^2[0, \frac{3\pi}{2}]}\) tal que:

(4.25)#\[\begin{align} &U''(x) + U(x) = 0 \text{ con }\\ &U(0) = 2 \\ &U(\dfrac{3\pi}{2})= 2 \end{align}\]

Este problema tiene solución única que es: \(U(x) = 2cos(x)-2sin(x)\)

4.5.1. Resolución Numérica. Operador diferencial#

Hemos visto que, a partir de los valores de una función en unos puntos \(\left[ a=x_0, x_1, x_2,...,x_N = b \right]\), podemos calcular el valor de la derivada de dicha función en esos puntos multiplicando por una matriz de derivación. En forma vectorial lo hemos escrito como: \(U' = DU\) siendo \(U\) un vector que contiene el valor de \(U(x)\) en el conjunto de puntos equiespaciados del intervalo \(\left[ a,b \right]\) y \(U'\) el vector de las derivadas en esos mismos puntos.

Si quisieramos hallar el valor de \(U''(x) + U(x)\) en los puntos de discretización podriamos construir el operador \(M = DD + I\), con \(I\) la matriz identidad, de manera que al calcular \(MU\) tendríamos:

(4.26)#\[\begin{equation} MU = (DD+I)U = DDU + IU = U'' + U \end{equation}\]

Pero este no es el caso que queremos resolver. En nuestro ejemplo no conocemos el valor de \(U(x)\), de hecho es precisamente eso lo que buscamos, una función que cumpla \(U''(x) + U(x)= 0\).

Vamos a plantear el problema de manera que las incógnitas sean los valores de la función en \(\left[ x_0, x_1, x_2,...,x_N \right]\), es decir el vector de incognitas será:

(4.27)#\[\begin{equation} U = \begin{pmatrix} U(x_0)\\ U(x_1)\\ U(x_2)\\ \vdots \\ U(x_N) \end{pmatrix} \end{equation}\]

y la información que tenemos es que \(U''(x) + U(x)= 0\) que en forma matricial podemos escribir como:

(4.28)#\[\begin{align} \begin{pmatrix} & & \\ & D & \\ & & \end{pmatrix} \begin{pmatrix} & & \\ & D & \\ & & \end{pmatrix} \begin{pmatrix} \\ U \\ \quad \end{pmatrix} + \begin{pmatrix} & & \\ & I & \\ & & \end{pmatrix} \begin{pmatrix} \\ U \\ \quad \end{pmatrix}= \begin{pmatrix} \\ 0 \\ \quad \end{pmatrix} \end{align}\]

O

(4.29)#\[\begin{align} \begin{pmatrix} & & \\ & M & \\ & & \end{pmatrix} \begin{pmatrix} \\ U \\ \quad \end{pmatrix} = \begin{pmatrix} \\ 0 \\ \quad \end{pmatrix} \end{align}\]

siendo:

(4.30)#\[\begin{align} \begin{pmatrix} & & \\ & M & \\ & & \end{pmatrix} = \begin{pmatrix} & & \\ & D & \\ & & \end{pmatrix} \begin{pmatrix} & & \\ & D & \\ & & \end{pmatrix} + \begin{pmatrix} & & \\ & I & \\ & & \end{pmatrix} \end{align}\]

Es decir, resolver la ecuación diferencial se convierte en resolver un sistema lineal donde la matriz del sistema es el operador diferencial asociado a la ecuación. Ahora solo falta por ver como imponer las condiciones de contorno.

4.5.2. Resolución Numérica. Condiciones de Contorno.#

En el ejemplo inicial hemos visto que existen infinitas funciones que cumplen \(U''(x) = -U(x)\) en \([0, \frac{3\pi}{2}]\). En terminos del problema matricial \(Ax = b\) que hemos planteado para resolverlo diríamos que el problema no tiene solución única (Teorema de Rouche-Frobenius). Necesitamos incluir las condiciones de contorno para que el sistema lineal sea compatible determinado y tenga, por tanto, solución única.

De momento el sistema se puede escribir como:

(4.31)#\[\begin{align} \begin{pmatrix} M_{00} & M_{01} & \dots & M_{0N} \\ M_{10} & M_{11} & \dots & M_{1N} \\ \vdots & & \ddots & \vdots \\ M_{N0} & M_{N1} & \dots & M_{NN} \end{pmatrix} \begin{pmatrix} U(x_0)\\ U(x_1)\\ \vdots \\ U(x_N) \end{pmatrix} = \begin{pmatrix} \\ 0 \\ \vdots \\ 0 \end{pmatrix} \end{align}\]

Las condiciones de contorno dicen que \(U(0) = U(x_0) = 2\) y que \(U(\frac{3\pi}{2})= U(x_N) = 2\). Vamos a forzar al sistema linbeal \(Ax = b\) para que verifique esas ecuaciones. El modo de hacerlo va a ser sustituir la primera y la última ecuación, -la primera y la última fila de la matriz \(M\) y del vector \(b\)-, para imponer las condiciones de contorno:

(4.32)#\[\begin{align} \begin{pmatrix} {\color{red} 1 }& {\color{red} 0 } & \color{red} {\dots} & {\color{red} 0 } & {\color{red} 0} \\ M_{10} & M_{11} & \dots & M_{1N{\tiny -1}} & M_{1N} \\ \vdots & & \ddots & & \vdots \\ M_{N{\tiny -1} 0} & M_{N{\tiny -1}1} & \dots & M_{N{\tiny -1}N{\tiny -1}} & M_{N{\tiny -1}N} \\ \color{green}{ 0} & \color{green}{0 }& \color{green}{\dots} & \color{green}{0} & \color{green}{1 } \end{pmatrix} \begin{pmatrix} \color{red} {U(x_0) }\\ U(x_1)\\ \vdots \\ U(x_{N-1}) \\ \color{green}{U(x_N)} \end{pmatrix} = \begin{pmatrix} {\color{red} 2 }\\ 0 \\ \vdots \\ 0 \\ \color{green}{2} \end{pmatrix} \end{align}\]

4.5.2.1. Practica 1.#

Resolver el problema de contorno:

(4.33)#\[\begin{align} &U''(x) + U(x) = 0 \text{ con }\\ &U(0) = 2 \\ &U(\frac{3\pi}{2})= 2 \end{align}\]

Y comparar graficamente la solución numérica obtenida con la verdadera \(U(x) = 2cos(x)-2sin(x)\).

4.5.2.2. Práctica 2.#

Resolver el problema de contorno:

(4.34)#\[\begin{align} &U''(x) + 9U(x) = cos(x) \text{ con }\\ &U'(0) = 5 \\ &U(\pi)= \frac{-9}{8} \end{align}\]

Fijaos que una de las condiciones de contorno está puesta en la derivada de la función. ¿Cómo cambia eso la manera de implementarla?

# Práctica 1 

n = 81
a = 0.
b = 3.*np.pi/2.
x_i = np.linspace(a, b, n)
dx = (b-a)/(n-1)

# Defino la matriz de derivación D





# Condiciones de contorno




### Representación gráfica
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[1], line 5
      3 n = 81
      4 a = 0.
----> 5 b = 3.*np.pi/2.
      6 x_i = np.linspace(a, b, n)
      7 dx = (b-a)/(n-1)

NameError: name 'np' is not defined
# Práctica 1 

n = 81
a = 0.
b = np.pi
x_i = np.linspace(a, b, n)
dx = (b-a)/(n-1)

# Defino la matriz de derivación D





# Condiciones de contorno




### Representación gráfica