Diferencia entre revisiones de «Ecuación de Laplace y ecuación de Poisson»

Revisión del 14:46 17 abr 2024

1 Introducción

En este artículo se trabajará en la ecuación de Laplace y la ecuación de Poisson, calculando distintas soluciones en distintos escenarios, graficándolas y modificando parámetros con el objetivo de poder alcanzar una mayor comprensión sobre la teoría de estas ecuaciones en derivadas parciales y relacionarla directamente con los resultados.

En primer lugar se tratará...

2 Preliminares

2.1 Laplaciano de una función

Sea una función [math]f:\mathbb{R}^n \rightarrow \mathbb{R}[/math] (se tomará en todo el artículo el espacio tridimensional, es decir, [math]n=3[/math]), se define el operador diferencial de segundo orden [math]\Delta f[/math], denominado laplaciano de [math]f[/math] como: Coordenadas cartesianas. [math]\Delta f = \frac{\partial^2 f}{\partial x^2} + \frac{\partial^2 f}{\partial y^2} + \frac{\partial^2 f}{\partial z^2}[/math]

Coordenadas cilíndricas. [math]\Delta f = \frac{1}{r} \frac{\partial}{\partial r} \left( r \frac{\partial f}{\partial r} \right) + \frac{1}{r^2} \frac{\partial^2 f}{\partial \theta^2} + \frac{\partial^2 f}{\partial z^2} [/math]

2.2 Desigualdad de Harnack

Sea u una función armónica (es decir, que [math] \Delta u = 0[/math]) y no negativa en la bola [math] B_R=B_R (0) \subset \mathbb {R^{n}}[/math], entonces para cualquier [math]x \in B_R[/math] la desigualdad de Harnack establece que:

2.3 método del trapecio

El método del trapecio es una técnica de integración numérica utilizada para aproximar el valor de una integral definida. Consiste en dividir el área bajo una curva en múltiples trapecios y sumar las áreas de estos trapecios para obtener una aproximación de la integral.

Para una función [math] f(x)[/math] en el intervalo [math][a, b][/math], el método del trapecio aproxima la integral definida como:

donde $h$ es el ancho de cada subintervalo [math](b - a) / n[/math] y [math]x_i = a + i \cdot h[/math] son los puntos donde se evalúa la función.

Este método es una aproximación simple pero efectiva de la integral y su precisión aumenta con el número de subintervalos utilizados ([math]n[/math]).

3 Contexto histórico

Para intentar entender las ecuaciones de la mejor forma posible, se debe conocer las raíces de donde surgieron. Por ello, se explicará brevemente.

En primer lugar, Pierre-Simon Laplace fue un astrónomo, físico y matemático francés. Este trabajó en la teoría del calor y en el estudio de los campos gravitatorios. Fue al estudiar el potencial gravitacional y el potencial eléctrico, con la hipótesis de que estos no tuvieran fuentes ni sumideros, cuando ideó y investigó las propiedades de lo que se conoce hoy como ecuación de Laplace. Además posteriormente se utilizó en más campos como la hidrodinámica y otros aspectos de la física.

Posteriormente, el matemático y físico francés Siméon Denis Poisson, realizó una generalización significativa de la ecuación de Laplace y su correspondiente fórmula para resolverla. De esta forma, se le dio el nombre de ecuación de Poisson. Que ambos investigaran temas semejantes no es coincidencia, ya que Laplace fue profesor de Poisson en una escuela de París y, gracias a ello, desarrollaron una gran relación de amistad.

ref: http://automata.cps.unizar.es/Biografias/Laplace.htm , https://en.wikipedia.org/wiki/Poisson%27s_equation , https://en.wikipedia.org/wiki/Sim%C3%A9on_Denis_Poisson , https://en.wikipedia.org/wiki/Laplace%27s_equation , https://en.wikipedia.org/wiki/Pierre-Simon_Laplace

4 Ecuación de Laplace

4.1 Planteamiento

Una vez definido el laplaciano de una función, la ecuación de Laplace es:

[math] \Delta u = 0 [/math]

donde [math]u:\mathbb{R}^3 \rightarrow \mathbb{R}[/math]

El problema de Dirichlet para la ecuación de Laplace, surge de manera inmediata al planteamiento de esta misma. De esta forma, se trata de encontrar una función que cumpla la ecuación de Laplace en un abierto [math]\Omega \subset \mathbb{R}^3[/math] y además fijar lo que vale la solución en la frontera del abierto [math]\partial \Omega[/math]:

De este problema se conoce el siguiente resultado:

Si [math] \Omega[/math] es un dominio acotado y g es una función continua en su frontera, es decir [math]g \in C(\partial \Omega)[/math]. Entonces, el problema tiene una única solución [math] u \in C^2(\Omega) \cap C(\overline{\Omega)} [/math]. Además, sean [math] u_{g_1} \text{ y } u_{g_2} [/math] las soluciones de los problemas correspondientes con [math] g_1, g_2 \in C(\partial\Omega) [/math]. Entonces se tiene:

[math] \text{(a) (Comparación). Si } g_1 \geq g_2 \text{ en } \partial\Omega \text{ y } g_1 \neq g_2 \text{ , entonces } u_{g_1} \gt u_{g_2} \text{ en } \Omega.[/math]

[math] \text{(b) (Estabilidad). } |u_{g_1}(x) - u_{g_2}(x)| \leq \max_{\partial\Omega} |g_1 - g_2| \text{ para todo } x \in \Omega. [/math]

Por otro lado, si [math] \Omega [/math] es la bola de centro [math] 0 [/math] y radio [math] R [/math], [math] B_R(0) [/math]. Entonces la solución [math] u \in C^2(B_R(0))[/math] viene dada por la fórmula de Poisson:

En coordenadas cartesianas.

[math] u(x)=\frac{R^2-|x|^2}{2\pi}\int_{\partial B_R} \frac{g(y)}{|x-y|^2} dy [/math]

En coordenadas polares.

[math] u(r,\theta)=\frac{R^2-r^2}{2\pi}\int_{\partial B_R} \frac{g(y)}{R^2+r^2-2Rrcos(\theta-y)} dy [/math]

No obstante, otra forma de obtener una solución para este problema es mediante series de Fourier. Tras realizar separación de variables sobre el problema modificado a coordenadas polares y aplicar el principio de superposición, se llega a la siguiente expresión de la solución:

donde los coeficientes se obtienen de la siguiente forma:

donde [math]G(\theta) [/math] es la condición frontera en coordenadas polares.

4.2 Ejemplo

A continuación, se planteará un problema de Laplace y se estudiará la solución por los métodos de la fórmula de Poisson y desarrollo por series de Fourier. Además se analizará el error que se produce por los distintos métodos y el cálculo numérico. También se estudiará el espacio en el que están las soluciones, utilizando la desigualdad de Harnack.

Sea [math] B_1 ⊂ R^2 [/math] la bola unidad centrada en [math](0, 0)[/math]. Se considera el siguiente problema.

donde la función [math]g[/math] viene descrita en coordenadas polares y se define como [math] g(\theta)=max\{0, 1-\frac{2}{\pi} |\theta - \pi|\} [/math].

En este caso, la solución dada por la ecuación de Poisson es:

[math] u(r,\theta)=\frac{1-r^2}{2\pi}\int_{\partial B_1} \frac{g(y)}{1+r^2-2rcos(\theta-y)} dy [/math]

Se debe apreciar que en la frontera de la bola [math]\partial B_1[/math] el denominador de la integral se hace cero cuando el coseno se anula. Es decir, cuando [math]\theta=y[/math]. Es por este motivo, que se obtienen los siguientes gráficos de la solución si se introduce la función.

Representación gráfica de la solución

Se puede observar como en la frontera, la superficie tiene singularidades por el motivo mencionado. Definiendo la función en la frontera, directamente como [math]g[/math] y representando la solución en el disco unidad, mediante el código adjunto a continuación, se obtienen las siguientes imágenes.

g = @(y) max(0, 1 - (2/pi) * abs(y - pi)); % Definimos la función g

% Definimos el dominio de la superficie
R=1; % Radio de la bola
n = 2; % Número de puntos, 10^n
Xr = linspace(0, R, 10^n); % Valores de r
Xomega = linspace(0, 2*pi, 10^n); % Valores de omega
[r, omega] = meshgrid(Xr, Xomega);

% Calculamos los valores de u(r, omega) para cada par (r, omega)
u = zeros(size(r));
for i = 1:numel(Xr)
    for j = 1:numel(Xomega)
        if r(i,j)==R     % Definimos la función en la frontera
            u(i,j) = g(omega(i,j));
        else     % Definimos la función en el resto del dominio
            y = linspace(0, 2*pi, 1000); % Discretizar el intervalo de integración
            integral_func = g(y) ./ (R^2 + r(i, j)^2 - 2*R * r(i, j) * cos(omega(i, j) - y));
            u(i, j) = (R^2 - r(i, j)^2) / (2 * pi) * trapz(y, integral_func);
        end
    end
end

% Representamos la superficie
surf(r .* cos(omega), r .* sin(omega), u);
xlabel('x');
ylabel('y');
zlabel('z');
title('Superficie de u');

4.2.1 Resolución mediante Serie de Fourier

A continuación, se va a resolver este problema mediante Serie de Fourier.

Como se ha mencionado anteriormente, la solución para los problemas con la ecuación de Laplace viene definida por

entonces se deben pasar los datos del problema a coordenadas polares, es decir,

ya que se trabajará en la bola unidad.

Por tanto, para este problema concreto se tendrá que los coeficientes se obtienen de la siguiente forma:

5 Ecuación de Poisson

A continuación, se definirá la ecuación de Poisson: [math] \Delta u = f [/math]

donde [math]u:\mathbb{R}^n \rightarrow \mathbb{R}[/math] y [math] f \in C^2(\mathbb{R})[/math].

En primer lugar se debe definir la solución fundamental del laplaciano.

5.1 Solución fundamental del Laplaciano para dimensiones pequeñas.

Se define como la solución fundamental a [math] \phi(x): \mathbb{R}^n \rightarrow \mathbb{R}[/math] a:

5.2 Potencial newtoniano.

El potencial newtoniano, también conocido como potencial de Newton, es un operador en cálculo vectorial que sirve como el inverso del Laplaciano negativo. Se aplica a funciones suaves que disminuyen lo suficientemente rápido hacia cero en el infinito. El potencial newtoniano es un operador integral singular que se define mediante la convolución con una función que presenta una singularidad en el origen. Esta función, es la solución fundamental del laplaciano, [math] \phi(x)[/math].

Se define el potencial newtoniano de la siguiente manera:

Los resultados conocidos para la resolución de problemas con la ecuación de Poisson en función de la dimensión son los siguientes:

Teorema. Sea [math]f \in C^2(\mathbb{R}^3)[/math] con soporte compacto. Sea [math]u[/math] el potencial newtoniano de [math]f[/math], definido por el potencial Newtoniano. Entonces, [math]u[/math] es la única solución en [math]\mathbb{R}^3[/math] de [math]\Delta u = -f[/math] que pertenece a [math]C^2(\mathbb{R}^3)[/math] y se anula en el infinito. Es decir, dado:

su solución viene dada por el potencial newtoniano:

Teorema. Sea [math]f \in C^2(\mathbb{R}^2)[/math] con soporte compacto. Sea [math]u[/math] el potencial newtoniano de [math]f[/math], definido por el potencial Newtoniano. Entonces, [math]u[/math] es la única solución en [math]\mathbb{R}^2[/math] de:

entonces, su solución viene dada por el potencial newtoniano:

y que pertenece a [math]C^2(\mathbb{R}^2)[/math].

5.3 Potencial logarítmico

El teorema visto en el caso anterior, se trata de un teorema general, en dimensión n=2 este potencial es el logarítmico, por lo que la solución en este caso se obtendría con :

Donde x e y son puntos de [math]{R}^{2}[/math]. Para ver como se emplea el potencial Newtoniano en este caso, se ilustra mediante un ejemplo, en el cual se tiene:

Donde f es la función característica de la bola de radio 1. Para resolver este ejemplo hay varias formas, en primer lugar uno al ver que la función f es la característica de la bola, lo puede plantear como una integral sobre dicha bola en vez de sobre [math]{\mathbb{R}^2}[/math]. Para ello, al estar en una simetría radial es interesante realizar un cambio a coordenadas polares , ya que de dicha forma los límite de integración quedan de forma más sencilla. De todas formas, debido a que la integral no se puede resolver de forma analítica, hay que hacerlo de forma numérica con un ordenador, por ello debido a la sencillez de la función característica, no se realizará dicho cambio.

Definimos la función característica, y con ello la solución [math]u(x)[/math], y para resolver el ejercicio se plantean dos mallas de puntos, una que represente [math]{\mathbb{R}^2}[/math], y otra que sea un cuadrado de lado 2 centrado en el (0,0), con el fin de optimizar el coste computacional y tener más valores en la bola de radio 1. Una vez creadas dichas mallas se evalúa la función en los puntos de la primera malla, integrándolos de forma numérica con los puntos de la segunda. El código que se ha empleado para este procedimiento es

Como se puede observar, la función no tiende a 0 a medida que la norma de la x aumenta, esto es porque el comportamiento asintótico del potencial logarítmico es:

Cuando [math]|x| \rightarrow \infty[/math], y podemos ver cómo eso se corresponde con la gráfica, donde

Al estar integrando sobre la bola unidad, M es [math]\pi[/math], por lo que finalmente la función asintóticamente será

De esta forma, se puede comprobar que la función se comporta de forma correcta en el infinito. siendo en el infinito como la función asintótica

6 Referencias