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

De MateWiki
Saltar a: navegación, buscar
(Introducción)
(Análisis del error de la solución por la fórmula de Poisson)
Línea 158: Línea 158:
 
==== Análisis del error de la solución por la fórmula de Poisson====
 
==== Análisis del error de la solución por la fórmula de Poisson====
  
Para analizar el error, tomaremos el siguiente caso específico:
+
Para analizar el error, tomaremos el siguiente problema específico:
 
+
  
 
<center><math> \begin{cases}
 
<center><math> \begin{cases}
Línea 166: Línea 165:
 
\end{cases}</math></center>
 
\end{cases}</math></center>
  
donde <math> g(x,y)=xy</math>
+
donde <math> g(x,y)=xy</math>.
 +
 
 +
En este caso, una solución es la propia función <math>g</math>, ya que es armónica. Aún así, se usará la solución por la fórmula de Poisson para hallar el error que se produce.
 +
Además, se debe mencionar el error producido por integrar con el método del trapecio. Este viene dado por la fórmula:
  
 
==== Resolución mediante Serie de Fourier====
 
==== Resolución mediante Serie de Fourier====

Revisión del 15:22 17 abr 2024

Trabajo realizado por estudiantes
Título Ecuación de Laplace y ecuación de Poisson.
Asignatura EDP
Curso 2023-24
Autores Alfredo de Lorenzo, Hugo Sanz, Manuel Fdez.
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

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á

A continuación, nos adentraremos en el estudio del potencial Newtoniano en dos dimensiones, específicamente el potencial logarítmico. Exploraremos un ejemplo de su aplicación y analizaremos su comportamiento asintótico, que se caracteriza por una disminución logarítmica a medida que nos alejamos indefinidamente de la fuente de carga. Este comportamiento revela la manera en que las influencias se propagan a larga distancia en sistemas físicos, ofreciendo una perspectiva única sobre la distribución de carga y el flujo de campos en el espacio.

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:

[math] \frac{R^{n-2} (R-r)}{(R+r)^{n-1}} u(x_0) \leq u(x) \leq \frac{R^{n-2} (R+r)}{(R-r)^{n-1}}u(x_0)[/math]

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:

[math] \int_{a}^{b} f(x) \,dx \approx \frac{h}{2} \left[ f(a) + 2f(x_1) + 2f(x_2) + \ldots + 2f(x_{n-1}) + f(b) \right] [/math]

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]:

[math] \begin{cases} \Delta u = 0, & \text{x} \in \Omega \\ u = g, & \text{x} \in \partial \Omega \end{cases}[/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:

[math] U(r, \theta) = \frac{\alpha_{0}}{2} + \sum_{k=1}^{n} \left( \alpha_{k} \left(\frac{r}{R} \right)^k cos(k\theta) + \beta_{k} \left(\frac{r}{R} \right)^k sin(k\theta)\right) [/math]

donde los coeficientes se obtienen de la siguiente forma:

[math]\alpha_{0} =\frac{1}{\pi}\int_{-\pi}^{\pi} \cos\Bigl(k \theta \Bigl) dx [/math]


[math]\alpha_{k} =\frac{1}{\pi}\int_{-\pi}^{\pi} G(\theta) \cos\Bigl(k \theta \Bigl) dx [/math]


[math]\beta_{k} =\frac{1}{\pi}\int_{-\pi}^{\pi} G(\theta) \sin\Bigl(k \theta \Bigl) dx [/math]

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.

[math] \begin{cases} \Delta u = 0, & \text{x} \in B_1 \\ u = g, & \text{x} \in \partial B_1 \end{cases}[/math]

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

4.2.1 Resolución mediante fórmula de Poisson

En este caso, la solución dada por la fórmula 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] y [math]\theta-y=\pi[/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.

Representación gráfica de la solución
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.2 Análisis del error de la solución por la fórmula de Poisson

Para analizar el error, tomaremos el siguiente problema específico:

[math] \begin{cases} \Delta u = 0, & \text{x} \in B_1 \\ u = g, & \text{x} \in \partial B_1 \end{cases}[/math]

donde [math] g(x,y)=xy[/math].

En este caso, una solución es la propia función [math]g[/math], ya que es armónica. Aún así, se usará la solución por la fórmula de Poisson para hallar el error que se produce. Además, se debe mencionar el error producido por integrar con el método del trapecio. Este viene dado por la fórmula:

4.2.3 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

[math] U(r, \theta) = \frac{\alpha_{0}}{2} + \sum_{k=1}^{n} \left( \alpha_{k} \left(\frac{r}{R} \right)^k cos(k\theta) + \beta_{k} \left(\frac{r}{R} \right)^k sin(k\theta)\right) [/math]

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

[math] g(x,y)= xy \longrightarrow g(r, \theta) = r^2 \cos(\theta) \sin(\theta) \hspace{0.1cm} \text{ cuando} \hspace{0.15cm}\begin{cases} x= rcos(\theta) & \text{r} \in [0,1] \\ y= rsin(\theta)& \theta \in [0,2\pi] \end{cases} [/math]

ya que se trabajará en la bola unidad. Entonces, en la frontera, donde [math]r=1[/math], la condición frontera será:

[math] G(\theta)=g(1, \theta) = \cos(\theta) \sin(\theta) \hspace{0.5cm} \theta \in [0,2\pi] [/math]

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

[math]\alpha_{0} =\frac{1}{\pi}\int_{-\pi}^{\pi} \cos\Bigl(k \theta \Bigl) dx [/math]


[math]\alpha_{k} =\frac{1}{\pi}\int_{-\pi}^{\pi} cos(\theta) \sin(\theta) \cos\Bigl(k \theta \Bigl) dx [/math]


[math]\beta_{k} =\frac{1}{\pi}\int_{-\pi}^{\pi} cos(\theta) \sin(\theta) \sin\Bigl(k \theta \Bigl) dx [/math]

La resolución de este problema se ha realizado por métodos numéricos, a continuación se mostrará el código empleado y los resultadoss obtenidos.

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:

[math] \phi(x) = \begin{cases} -\frac{1}{2\pi} log(|x|) & \text{si n=2} \\ \frac{1}{4\pi |x| } & \text{si n=3} \end{cases}[/math]

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:

[math] u(\mathbf{x}) = \int_{\mathbb{R}^n} f(\mathbf{y}) \phi(\mathbf{x} - \mathbf{y}) d\mathbf{y} [/math]


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:

[math] \begin{cases} \Delta u = -f \\ u(x) \rightarrow 0 & \text{si } |x| \rightarrow \infty \end{cases}[/math]

su solución viene dada por el potencial newtoniano:

[math] u(\mathbf{x}) = \int_{\mathbb{R}^3} f(\mathbf{y}) \phi(\mathbf{x} - \mathbf{y}) d\mathbf{y} [/math]


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:

[math] \begin{cases} \Delta u = -f \\ u(x)=-\frac{M}{2\pi} log |x| + O\left(\frac{1}{|x|} \right) & \text{si } |x| \rightarrow \infty \end{cases}[/math]
donde
[math] M= \int _{\mathbb{R}^2} f(\mathbf y) d\mathbf y[/math]

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

[math] u(\mathbf{x}) = \int_{\mathbb{R}^3} f(\mathbf{y}) \phi(\mathbf{x} - \mathbf{y}) d\mathbf{y} [/math]

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 :

[math] u(\mathbf{x}) = -\frac{1}{2\pi}\int_{\mathbb{R}^2} f(\mathbf{y}) log|\mathbf{x} - \mathbf{y}| d\mathbf{y} [/math]

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:

[math] \Delta u = f ;x \in {\mathbb{R}^2} [/math]

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:

[math] u(x)=-\frac{M}{2\pi} log |x| + O(\frac{1}{|x|})[/math]

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

[math] M= \int _{\mathbb{R}^2} f(\mathbf y) d\mathbf y[/math]

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

[math] u(x)=-\frac{1}{2} log |x|[/math]
, para ver si se cumple, se coge un radio de la simetría radial, y se ven sus valores a medida que el radio crece, y se comparan sus valores en una gráfica.

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

De hecho, el potencial logarítmico es la única solución de
[math] \Delta u = -f ;x \in {\mathbb{R}^2} [/math]
que satisface dicho límite asintótico. Para hacer esta sección se ha utilizado el libro "Partial Differentialm Equations in Action From Modelling to Theory" de Sandro Salsa

6 Referencias