Ecuación de Laplace y ecuación de Poisson(Grupo 2 1/2)

Trabajo realizado por estudiantes
Título	Ecuación de Laplace y ecuación de Poisson.(Grupo 2 1/2)
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 teóricos.

En primer lugar, se resolverá un problema específico de Dirichlet en dos dimensiones por distintos métodos y comparándolos a través del error que se produce en cada uno de ellos. Además se investigará en el motivo de este error.

A continuación, se profundizará en el estudio del potencial Newtoniano en dos dimensiones, específicamente el potencial logarítmico. Explorando un ejemplo de su aplicación y analizando su comportamiento asintótico, que se caracteriza por una disminución logarítmica a medida que la distancia aumenta 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 bidimensional, es decir, [math]n=2[/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}[/math]

Coordenadas polares. [math]\Delta f = \frac{\partial^2f}{\partial r^2}+ \frac{1}{r} \frac{\partial f}{\partial r} + \frac{1}{r^2} \frac{\partial^2 f}{\partial \theta^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] \overline{B_R (x_0)} \subset \mathbb {R^{n}}[/math], entonces para cualquier [math]x \in \overline{B_R(x_0)} [/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 [math]h[/math] 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]).

[math] E = -\frac{1}{12 n^{2}}f''(\xi)(b-a)^3 [/math]

Donde:

• [math] E [/math] es el error de aproximación,
• [math] a[/math] y [math] b [/math] son los límites de integración,
• [math] n [/math] es el número de subintervalos utilizados en la aproximación,
• [math]f''(\xi) [/math] es la segunda derivada de la función integrada evaluada en algún punto [math] \xi [/math] dentro del intervalo [math] [a, b][/math].

Esta fórmula explica que el error de aproximación es proporcional al cuadrado del ancho del intervalo de integración [math] (b-a)^2 [/math] y disminuye proporcionalmente al cuadrado del número de subintervalos utilizados [math] n^2 [/math]. Además, está multiplicado por el valor máximo de la segunda derivada de la función integrada en el intervalo [math] [a, b] [/math], lo que significa que el error aumentará si la función es más "curva" en ese intervalo.

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. Laplace 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ó e 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 y fórmula de Poisson a su solución. 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, profesional y académica.

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}^2 \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 :

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

con [math]g(\theta) [/math] 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 las distintas formas de hallar la solución y por el método del trapecio. También se estudiará el región 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{R-r^2}{2\pi}\int_{\partial B_1} \frac{g(y)}{R+r^2-2Rrcos(\theta-y)} dy [/math]

donde [math]R=1[/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. Además esto se ve reflejado en la teoría dado que la fórmula de Poisson asegura estar definida únicamente en el interior del dominio.

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');

Se puede apreciar que la función cumple la proposición enunciada anteriormente. Ya que el máximo y el mínimo de la función se ubica en la frontera, y esta es estable, es decir, la superficie es suave.

4.2.2 Análisis del error de la solución por la fórmula de Poisson

Para analizar el error, se tomará el siguiente problema específico:

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

donde [math] g_{2}(x,y)=xy[/math].

En este caso una solución es la propia función [math]g_{2}[/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.

En primer lugar, se mide el error producido por la fórmula de Poisson en el punto [math](r,\theta)=(0.9,\frac{\pi}{4})[/math] variando la cantidad de puntos [math]10^n[/math] del mallado usado para evaluar la integral mediante el método del trapecio.

Se representará la función [math]Error(n) = log_{10}(|g_{2}(r,\theta)-u(r,\theta)|)[/math], donde [math] u [/math] es la solución por la fórmula de Poisson.

Representación gráfica del error

% Definimos el dominio de la superficie
R = 1; % Radio de la bola,
nmax = 7; % Número maximo de mallado, 10^namx
r = 0.9;
error = zeros(length(nmax));
for i = 1:nmax
    Xomega = linspace(0, 2*pi, 10^i); % Valores de omega
    % Definimos la solución dada y la solución por la fómula de Poisson
    g = @(r,omega) r.^2.*cos(omega).*sin(omega);
    u = @(r,omega) (R^2 - r^2) / (2 * pi) * trapz(Xomega, g(r,Xomega) ./ (R^2 + r^2 - 2.*R .* r .* cos(omega - Xomega)));
    error(i) = log10(g(r,pi/4)-u(r,pi/4));
end
% Graficamos el resultado
plot(linspace(0,nmax,length(error)),error)
xlabel('n')
ylabel('Error')

En la imagen se aprecia como el error decae bruscamente para los primeros valores de n y para los siguientes valores de n, decrece de manera más lenta. Esto se debe a que el error que se produce es debido a la forma de integrar con la regla del trapecio. Para mostrarlo, se estimará el máximo del error producido por la regla del trapecio para cada valor de n, en la malla de puntos en la que se integra en la fórmula de Poisson. De esta forma se obtendrá una cota superior.

Se debe mencionar el error producido por integrar con el método del trapecio. Este viene dado por la fórmula:

[math] E = -\frac{1}{12(10^n)^2} f''(y)(2\pi)^3 [/math]

donde [math]f(y) = \frac{g_{2}(r,y)}{1+r^2-2rcos(\theta-y)}[/math] y [math]g_{2}[/math] ha sido parametrizada en coordenadas polares, resultando [math]g_{2}(r,\theta)=r^2cos(\theta)sin(\theta)[/math]

Representación gráfica de la comparación de errores

% Definimos el dominio de la superficie
R = 1; % Radio de la bola
nmax = 5; % Número máximo de n, 10^nmax
G = @(omega) R.^2.*cos(omega).*sin(omega); % Definimos la función en la frontera
% Definimos la función a integrar para hallar el error máximo que se
% produce en el método del trapecio
syms y r theta 
f = r^2*cos(y)*sin(y)/(1+r^2-2*r*cos(theta-y)); % Parte de la integral de la fórmula de Poisson
d = diff(f,y); % Derivamos una vez
df = diff(d,y); % Derivamos dos veces
r_value = 0.9; % Radio del punto que evaluamos
theta_value = pi/4; % Ángulo del punto que evaluamos

Error_trapecio_total = zeros(length(nmax));
Error_funcion = zeros(length(nmax));
for i = 1:nmax
    Error_trapecio = zeros(10^i,1); 
    Y = linspace(0,2*pi,10^i); % Espacio de valores donde se hace el máximo del error del trapecio
    Xomega = linspace(0, 2*pi, 10^i); % Valores de omega
    % Definimos la solución dada y la solución por la fórmula de Poisson
    g = @(r,omega) r.^2.*cos(omega).*sin(omega);
    u = @(r,omega) (R^2 - r^2) / (2 * pi) * trapz(Xomega, G(Xomega) ./ (R^2 + r^2 - 2.*R .* r .* cos(omega - Xomega)));
    Error_funcion(i) = log10(abs(g(r_value,pi/4)-u(r_value,theta_value)));
    for j=1:10^i
        y_value= Y(j);
        Error_trapecio(j)= log10(abs(-1/(12*(10^i)^2) * (2*pi)^3*subs(df, [y,r,theta], [y_value,r_value,theta_value])));
    end
    Error_trapecio_total(i) = max(Error_trapecio); %Seleccionamos el máximo del error del trapecio en la malla
end

% Graficamos el resultado
plot(linspace(1,nmax,nmax),Error_funcion)
hold on
plot(linspace(1,nmax,nmax),Error_trapecio_total)
xlabel('n')
ylabel('Errores') 
legend(["Error(n)","Error trapecio"])

Se puede comprobar que, en efecto, el error de la regla del trapecio proporciona una cota que va disminuyendo a la vez que el error producido por la fórmula de Poisson. Por lo que queda reflejado la evidente relación entre ambos.

Ahora se fijan el número de puntos en los que se divide la malla, y variará el punto en el que se evalúa [math](r,\theta)=(1-10^{-n},\frac{\pi}{4})[/math] acercándose a la frontera de la bola.

Representación gráfica de la comparación de errores

% Definimos el dominio de la superficie
R = 1; % Radio de la bola
nmax = 10; % Número máximo de n para el radio
n = 3; % Número de n puntos
G = @(omega) R.^2.*cos(omega).*sin(omega); % Definimos la función en la frontera
% Definimos la función a integrar para hallar el error máximo que se
% produce en el método del trapecio
syms y r theta 
f = r^2*cos(y)*sin(y)/(1+r^2-2*r*cos(theta-y));
df = diff(f,y,2); % Derivamos dos veces
theta_value = pi/4; % Ángulo del punto que evaluamos

Error_trapecio_total = zeros(length(nmax));
Error_funcion = zeros(length(n));
for i = 1:nmax
    Error_trapecio = zeros(10^n,1); 
    Y = linspace(0,2*pi,10^n); % Espacio de valores donde se hace el máximo del error del trapecio
    Xomega = linspace(0, 2*pi, 10^n); % Valores de omega
    % Definimos la solución dada y la solución por la fómula de Poisson
    g = @(r,omega) r.^2.*cos(omega).*sin(omega);
    u = @(r,omega) (R^2 - r^2) / (2 * pi) * trapz(Xomega, G(Xomega) ./ (R^2 + r^2 - 2.*R .* r .* cos(omega - Xomega)));
    r_value = 1-10^(-i);
    Error_funcion(i) = log10(abs(g(r_value,pi/4)-u(r_value,theta_value)));
    for j=1:10^n
        y_value= Y(j);
        Error_trapecio(j)= log10(abs(-1/12 * (2*pi)^3*subs(df, [y,r,theta], [y_value,r_value,theta_value])));
    end
    Error_trapecio_total(i) = max(Error_trapecio);
end

% Graficamos el resultado
plot(linspace(1,nmax,nmax),Error_funcion)
hold on
plot(linspace(1,nmax,nmax),Error_trapecio_total)
xlabel('n')
ylabel('Errores') 
legend(["Error(n)","Error trapecio"])

Como se puede apreciar, a medida que el punto se acerca a la frontera, el error crece. Esto se debe a que, como se ha visto anteriormente, la fórmula de Poisson en la frontera no está determinada. Respecto al error producido por la regla del trapecio, sucede lo mismo.

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 resultados obtenidos.

Representación de las solución en función de los términos

% Definimos el radio de la bola unidad. 
R=1;
% Definimos la función frontera en coordenadas polares.
g = @(y) R.^2.*cos(y).*sin(y);
% Discretizar el intervalo de integración.
y = linspace(-pi, pi, 1000); 
% Discretizar el dominio de la solución en polares.
X=linspace(0, 2*pi, 100);
Radio=linspace(0,1, 50);
% Definimos el número de coeficientes.
terminos=[5,10,100];


%Creamos un bucle para que recorra la lista "terminos", y nos resuelva
%el problema en función del número de los términos de la serie. 
for j=1:length(terminos)
    for i=terminos(j)
        % Calculamos el coeficiente a0 de la serie de fourier.
        a_0= (1/pi)*trapz(y, g(y));
        
        % Definimos el vector para meter los coeficientes de Fourier.
        a_k = zeros(i, 1);
        b_k = zeros(i, 1);
        % Calcular las integrales el método del trapecio 
        for n = 1:i
            h1= @(y) g(y).*sin(n.*y);
            h2= @(y) g(y).*cos(n.*y);
            resultado_seno = (1/pi)*trapz(y,g(y).*sin(n*y));
            resultado_cos = (1/pi)*trapz(y,g(y).*cos(n*y));
        
            % Almacenar el resultado
            a_k(n) = resultado_cos; 
            b_k(n) = resultado_seno;
        end
        a_k
        b_k
        % Obtenemos la solución final para el problema por serie de Fourier.
        fn=@(r,x) a_0/2;
        for n=1:i
            fn = @(r,x) a_k(n) .* (r/R).^n.* cos(n.*x) + b_k(n) .* (r/R).^n.* sin(n.*x)+  fn(r,x);
        end
        
        % Definimos la solución exacta.
        h= @(r,x) r.^2.*cos(x).*sin(x);

        %Mallado para la representación
        [r, x] = meshgrid(Radio, X);
        z=fn(r,x); % evaluamos la sol por serie de Fourier en la malla de puntos.
        Z= h(r,x); % evaluamos la sol exacta en la malla de puntos.

        % Representamos ambas en una misma gráfica. 
            subplot(2,2,j+1)
            surf(r .* cos(x), r .* sin(x), z)
            grid on; % Activar la cuadrícula
            title(['Gráfico de la solución con ',num2str(terminos(j))], ' términos');
            subplot(2,2,1)
            surf(r .* cos(x), r .* sin(x), Z)
            grid on; % Activar la cuadrícula
            title(['Gráfico de la solución exacta']);
    end
end

Se puede observar en estas gráficas como el número de términos de la serie de Fourier no afecta a la representación de la solución. Esto se debe a que los coeficientes de la serie de Fourier son todos nulos excepto [math]\beta_2[/math] ya que [math]\beta_2=\frac{1}{2}[/math], entonces el número de términos no afecta al desarrollo de la solución por Serie de Fourier.

En la siguiente sección se verá esto de forma más clara.

4.2.4 Análisis del error de la solución por Serie de Fourier

A continuación, se va a estimar el error que se obtiene al aproximar la serie de Fourier pintando el máximo del error frente al número de términos en escala logarítmica. Se va a estudiar error como el supremo del error entre la solución obtenida por serie de Fourier y la solución exacta.

Una vez calculado los errores de forma numérica, obtenemos la siguiente gráfica en función del número de términos.

Representación del error con serie de Fourier

Se ha realizado en función del número de términos. A pesar de que se ha visto en la sección anterior que el único término de la serie de Fourier que afectaba era [math] \beta_{2}[/math], en la gráfica se ve como según se aumentan el número de términos el [math] \log_{10}(error(n)) [/math] aumenta. Esto se puede deber a que, como se han empleado métodos numéricos, los coeficientes no son exactamente 0 cuando n suficientemente grande y estamos trabajando en escala logarítmica, por lo que el cambio en el orden del error se aprecia mas significativo.

A pesar de esto, se está hablando de errores de [math] 10^{-15} - 10^{-16} [/math], por ello, se está hablando de errores suficientemente bajos como para considerarlos significativos. Entonces se ve que la solución es suficientemente próxima a la exacta.

Por otro lado, destacar que para [math] n=1 [/math], la solución presenta un error considerable, aproximadamente [math] \frac{1}{2} [/math], ya que el único término de la serie de Fourier que se considerará para que la solución coincida con la exacta, no se ha calculado para esa solución.

Finalmente destacar, que en detrimento de la solución obtenida mediante la fórmula de Poisson el error es muy pequeño en todo punto de la superficie; por lo tanto, se "arregla" el problema que se tenía para la solución anterior en los puntos próximos a la frontera.

Se muestra el código empleado para el cálculo de errores y la representación de este.

% Definimos el radio de la bola unidad. 
R=1;
% Define la solución exacta.
h= @(r,x) r.^2.*cos(x).*sin(x);

% Definimos la función frontera en coordenadas polares.
g = @(y) R.^2.*cos(y).*sin(y);

% Discretizar el intervalo de integración
y = linspace(-pi, pi, 1000); 

% Discretizar el dominio de la solución en polares.
THETA=linspace(0, 2*pi, 100);
Radio=linspace(0,1, 100);

% Definimos el número de terminos de la serie.
terminos = 100;

%Calculamos el coeficiente a0 de la serie de fourier.
a_0= (1/pi)*trapz(y, g(y));

%definimos la matriz donde se almacenarán los errores.
error = zeros(length(terminos));
for i = 1:terminos
    % Definimos el vector para almacenar los coeficientes de Fourier.
    a_k = zeros(i, 1);
    b_k = zeros(i, 1);
    % Calcular las integrales usando el método del trapecio 
    for n = 1:i
        resultado_seno = (1/pi)*trapz(y,g(y).*sin(n*y));
        resultado_cos = (1/pi)*trapz(y,g(y).*cos(n*y)); 
    
        % Almacenar el resultado
        a_k(n) = resultado_cos; 
        b_k(n) = resultado_seno;
    end
    % Obtenemos la solución final para el problema por serie de Fourier.
    fn=@(r,x) a_0/2;
    for n=1:i
        fn = @(r,x) a_k(n) .* (r/R).^n.* cos(n.*x) + b_k(n) .* (r/R).^n.* sin(n.*x)+  fn(r,x);
    end

    % Genera una malla de puntos en el dominio de interés
    [x, r] = meshgrid(Radio, THETA);
    
    % Evalúa ambas funciones en cada punto de la malla
    Z_v = h(x, r);
    Z_e = fn(x, r);
    
    % Encuentra el máximo valor entre las dos funciones en cada punto
    Z_sup = max(Z_v- Z_e);
    
    % Encuentra el supremo de los valores máximos
    supremo = max(Z_sup, [], 'all');
    
    % Muestra el resultado
    disp(['El supremo entre las dos funciones es: ', num2str(supremo)]);

    % Encuentra el supremo de los valores máximos
    error(i) = log10(max(Z_sup, [], 'all'));
end
%Representamos la cota superior del error para cada n.
plot(linspace(0,terminos,length(error)),error)

4.2.5 Análisis de la Desigualdad de Harnack

A continuación, mediante la desigualdad de Harnack se tratará de encontrar la región en la que deben estar todas las soluciones armónicas en la bola unidad que valgan lo mismo que nuestra solución del problema en el punto [math] (0,0)[/math], donde está centrada la bola.

En la sección 2.2 se ha definido la desigualdad, por tanto a continuación, se aplicará para [math] n=2 [/math] ya que [math] B_R=B_R (0) \subset \mathbb {R^{2}}[/math].

Para una fácil comprensión de esta desigualdad se ha realizado un código donde se representan ambos extremos de la desigualdad:

[math] \frac{1-r}{1+r} v(0,0) \leq v(x) \leq \frac{1+r}{1-r}v(0,0)[/math]

donde va a depender del radio, con [math] r \in [0,1) [/math] y [math] v= u - M[/math] donde [math] M= min(g) \text{ en } B_1[/math].

Se calcula [math] M[/math] en función del radio de la bola del dominio definido, tal que:

[math]M := \min_{(x, y) \in B_R} \{g(x,y)\} = \min_{\theta \in [0,2\pi)} \{g(\theta)\} = \min_{\theta \in [0,2\pi)} \{R^2 cos(\theta)sen(\theta)\}= \min_{\theta \in [0,2\pi)} \{R^2 \frac{\sin(2\theta)}{2}\} =\frac{R^2}{2} \min_{\theta \in [0,2\pi)} \{\sin(2\theta)\}= - \frac{R^2}{2} .[/math]

Para [math] B_1[/math], [math] M=\frac{1}{2}[/math], se puede observar la región en la que deben estar todas las soluciones armónicas en la [math] B_1[/math], que valen lo mismo que v en el punto [math] (0,0)[/math].

Representación de la desigualdad de Harnack para [math] B_1[/math]

Para una mejor apreciación se utilizará la escala logarítmica, de esta manera obtenemos la siguiente región:

Representación de la desigualdad de Harnack para [math] B_1[/math] en escala logarítmica

Por otro lado, se va a realizar este mismo proceso, en escala logarítmica, para los dominios [math] B_2 \text{ y } B_{10} [/math], entonces:

[math] \frac{R-r}{R+r} v(0,0) \leq v(x) \leq \frac{R+r}{R-r}v(0,0)[/math]

y se compararán con la anterior.

Representación de la desigualdad de Harnack para los diferentes radios en escala logarítmica

Se puede observar como las cotas presentan la misma estructuras independientes del radio de la bola, sin embargo, varía el valor de [math] v(0,0)[/math], es decir, se establece una traslación en la región donde están todas las soluciones armónicas en la respectiva bola y sean iguales a [math] v[/math] en [math] (0,0)[/math].

Se debe destacar que en estas gráficas se observa como la región de estas gráficas disminuye considerablemente según se aumenta el radio, ya que si se comparan las gráficas para [math] R=1[/math] y [math] R=2[/math] y se restringe el intervalo observado al [math] [0,1] [/math], la región de las soluciones se ve reducida considerablemente. Esto ocurre de forma más significativa para el caso de [math] R=10[/math].

Finalmente, vamos a comparar estas cotas con las que se obtendrían para [math] n=3 [/math], donde la cota va a estar definida como:

[math] \frac{R(R-r)}{(R+r)^{2}} v(0,0,0) \leq v(x) \leq \frac{R(R+r)}{(R-r)^{2}}v(0,0,0) [/math]

Representación de la desigualdad de Harnack para los diferentes radios en escala logarítmica

En estas gráficas, se puede observar como la región para [math] n=3 [/math] aumenta considerablemente, sobre todo en la cota superior.

Por tanto, aumenta el dominio total donde se definen las soluciones armónicas que satisfagan la condición de [math] v(0,0,0) [/math] cumple lo mencionado para las bolas en [math] n=2 [/math].

A continuación, se define el código empleado a lo largo de esta sección.

% Definimos los diferentes radios con los que vamos a hacer el estudio.
Radios=[1,2,10];

% Creamos un bucle que varía en función de la dimensión (está definido para
% los casos n=2 y n=3). 
for gr=2:3
        % Una vez definido el grado, vamos a aplicar la desigualdad de 
        % Harnock para todos los radios.
    for j=1:length(Radios) 
        for R=Radios(j)

            % Definimos la función frontera en coordenadas polares.
            g = @(y) R.^2.*cos(y).*sin(y);
            % Discretizar el intervalo de integración
            y = linspace(-pi, pi, 1000); 

            %Calculamos el coeficiente a0 de la serie de fourier.
            a_0= (1/pi)*trapz(y, g(y));
            % Definimos el número de coeficientes.
            i=100;
            % Definimos el vector para meter los coeficientes de Fourier.
            a_k = zeros(i, 1);
            b_k = zeros(i, 1);
            
            % Calcular las integrales usando el método del trapecio 
            for n = 1:i
                resultado_seno = (1/pi)*trapz(y,g(y).*sin(n*y));
                resultado_cos = (1/pi)*trapz(y,g(y).*cos(n*y)); 
            
                % Almacenar el resultado
                a_k(n) = resultado_cos; 
                b_k(n) = resultado_seno;
            end
            % Obtenemos la solución final para el problema por serie de Fourier.
            fn=@(r,x) a_0/2;
            for n=1:i
                fn = @(r,x) a_k(n) .* (r/R).^n.* cos(n.*x) + b_k(n) .* (r/R).^n.* sin(n.*x)+  fn(r,x);
            end
            
            %Calculamos el mínimo de la función- condición frontera.
            x=fminbnd(g,0,2*pi);
            M=g(x) %mínimo evaluado en la función.
            
            % Se define v=u-M tq v>=0.
            v= @(r,x) fn(r,x) - M;
            %Se evalúa esa función en x_0 = 0.
            p=v(0,0)
            
            % Se define el extremo inferior de la desigualdad
            f1=@(r) ((R).^(gr-2)*(R-r)./((R+r).^(gr-1)))*p;
            % Se define el extremo superior de la desigualdad
            f2=@(r) ((R).^(gr-2)*(R+r)./((R-r).^(gr-1)))*p;

            % Defininimos el dominio del radio
            r = linspace(0, R, 100);
            
            %Represenación de ambas cotas en escala logarítmica.
            figure(j)
            hold on
            plot(r, log10(f1(r)));
            plot(r, log10(f2(r)));
            xlabel('r');
            ylabel('f(r)');
            title('Gráfico de f(r)');
            legend('Desigualdad inferior', 'Desigualdad superior')
        end
    end
end

Se ha realizado el estudio sobre [math] v(x)[/math], sin embargo, se busca estudiar el comportamiento relacionado con las soluciones de [math] u(x)[/math]. Por tanto, se debe deshacer el cambio de tal forma que:

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

Finalmente, como se conoce el valor de [math] M [/math], se llega a:

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

generalizando a cualquier grado y radio.

Representando entonces las regiones donde están definidas las soluciones armónicas para [math] n=2 [/math] y [math] n=3 [/math] en función de los radios se obtiene:

Representación de la desigualdad de Harnack para los diferentes radios en escala logarítmica

Y en estas gráficas se pueden observar las mismas condiciones que se han definido para el caso de [math] v(x)[/math].

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]. Se va a trabajar con un ejemplo práctico, donde se verá el comportamiento del potencial logarítmico, y se verá cómo se vomporta cuando la norma del punto tiende a infinito, viendo la diferencia con dicha asíntota. En primer lugar se debe definir la solución fundamental del laplaciano.

5.1 Solución fundamental del Laplaciano.

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 [math] u(\mathbf{x}) [/math] , que pertenece a [math]C^2(\mathbb{R}^2)[/math], viene dada por el potencial newtoniano.

5.3 Potencial logarítmico

El teorema visto en el caso anterior, se trata de un teorema general, en dimensión [math]n=2 [/math] 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 \hspace{0.4cm} 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. Otra forma, debido a que la integral no se puede resolver de forma analítica y que se hace de forma numérica, se puede resolver sin realizar dicho cambio. Una de las opciones puede ser definiendo la función característica de la bola, e integrando sobre un cuadrado, pero esto se puede omitir teniendo en cuenta que por la simetría de la función, la integral es proporcional al área en el que se integra, por lo que por facilidad a la hora de resolver el problema se queda:

[math] -\frac{1}{2\pi}\int_{\partial B_1 0} log|\mathbf{x} - \mathbf{y}| d\mathbf{y} = -\frac{\pi}{4}\frac{1}{2\pi}\int_{-1} ^1 \int_{-1} ^1 log|\mathbf{x} - \mathbf{y}| dx dy [/math]

Ya que [math] \frac{\pi}{4} [/math] es la razón entre el área de la bola unidad y el cuadrado de lado 1. Finalmente con ello se obtiene 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). 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.

Para este procedimiento se ha utilizado el siguiente código:

% Definir la función f que es el potencial logarítmico, y la vamos a
% multiplicar por -pi/4, ya que es la razón entre el área del cuadrado de
% lado 2 y la bola de radio 1, y se necesita para corregir el dominio de integración 
f = @(x1, y1, x2, y2) -pi/4*(1/(2*pi))*log(sqrt((x1-x2).^2 + (y1-y2).^2));

% Definir los límites del intervalo
%Definimos a y b los cuales están definidos entre -1 y 1 para optimizar, ya
%que cuando estos valores estén fuera de la bola unidad, la función será 0,
%por lo que podemos limitarnos a este caso.
a = linspace(-1,1,100);
b = linspace(-1,1,100);
%La función es en R^2, por lo que vamos a coger un dominio lo
%suficientemente grande como para ver como se comporta la función cuando
% la norma de x tiende a infinito. Con tener vectores de 200 puntos nos
% sirve para ver cómo es el comportamiento asintótico.
c = linspace(-10^3,10^3,100);
d = linspace(-10^3,10^3,100);

% Generar las mallas de puntos en las que vamos a representar la función.
% La primera malla se define para establecer el dominio en el que se va a
% integrar, en este caso un cuadrado 
%lado 2 centrado en el (0,0)
[X, Y] = meshgrid(a, b);
%La segunda malla de puntos es en la que se representará el dominio de la
%función
[X2, Y2] = meshgrid(c,d);
% Inicializar la matriz de valores en la que se guardarán los valores de
% la función
values = zeros(size(X)); 
%Crear un bucle para evaluar la función de forma numérica. Como la función
%tiene 4 variables y se integra en 2, el bucle se emplea para que dos de
%las variables tengan un valor numérico, y los vamos cambiando con las
%iteraciones.
for i = 1:length(a)
    for j = 1:length(b)
        % Definir fnum para evaluar f en c(i) y d(j) las mallas de puntos 
        % definidas, de tal forma que se quede una función de 2 variables
        % la cual ya se podrá integrar.
        fnum = @(x, y) f(x, y, c(i), d(j)); 
        %Una vez se han sustituido los valores, se integra de forma
        %numérica en dos dimensiones en la malla de puntos [X,Y]
        values(i, j) = trapz(b, trapz(a, fnum(X, Y), 2)); 
    end
end
%Representamos nuestra función con la malla de puntos de R^2, y los valores
%obtenidos de evaluar de forma numérica la función en dicho punto.

surf(X2,Y2,values)

Representación de la solución

Representación de la solución en un dominio más grande

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 se ve 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.

%Con el mismo planteamiento que antes se elige un radio de la solución para compararlo
repr=zeros(length(a));
ejexrepr=zeros(length(a));
for i = 1:length(a)
        fnum = @(x, y) f(x, y, c(i), d(i));
        
        repr(i) = trapz(b, trapz(a, fnum(X,Y), 2));
        ejexrepr(i)=sqrt(2*(c(i)^2));
end


hold on
%Definimos la asíntota teórica
fasint= @(x) -1/2*log(x);
valores_y = fasint(ejexrepr);
% Crear el gráfico de las líneas y los puntos
hold on
plot(ejexrepr, valores_y, '-','DisplayName','Función asintótica'); % Línea
plot(ejexrepr, repr, '*','DisplayName','Función obtenida'); % Punto
legend('Función asintótica', 'Función obtenida')
xlabel('Eje x')
ylabel('Eje y')
title('Gráfico función asintótica y función obtenida')

Representación de la función teórica y la aproximación frente al radio

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, a continuación se muestra una gráfica del error frente al radio, donde se puede ver mejor cómo se aproxima a la asíntota.

xlabel('radio')
ylabel('logaritmo del error')
title('Gráfico del logaritmo del error de la aproximación')
plot(ejexrepr,log(abs(valores_y-repr)))

Representación del logaritmo del error frente al radio

De hecho, el potencial logarítmico es la única solución de [math] \Delta u = -f \hspace{0.4cm} x \in {\mathbb{R}^2} [/math] que satisface dicho límite asintótico.