Revisión actual del 22:19 19 abr 2024

1 Introducción

En el siguiente documento estudiaremos la fórmula de Poisson a partir de varios ejemplos de problemas de Laplace para condiciones frontera de tipo Dirichlet. Analizaremos cuáles son sus limitaciones y los problemas que surgen al intentar solventarlos mediante métodos numéricos. Además, trataremos de entender en profundidad el significado de la desigualdad de Harnack tanto analíticamente como gráficamente. Como consecuencia demostraremos el teorema de Liouville. Por último, quitaremos las condiciones frontera y estudiaremos el carácter de las soluciones de la ecuación de Laplace en un dominio no acotado, en este caso [math]\mathbb{R}^{2}[/math].

2 Preliminares

Antes de comenzar, es necesario presentar la notación que emplearemos y sobre todo, el problema que estudiaremos en el documento. Además, enunciaremos algunos de los resultados más importantes relacionados con dicho problema. Este es de la forma: [math] \left\{ \begin{aligned} \Delta u &= f, & x &\in \Omega \\ u &= g, & x &\in \partial \Omega \end{aligned} \right. [/math] donde [math]\Omega \subset \mathbb{R}^n[/math] es un abierto. Nótese que en caso de que [math]f=0[/math] la ecuación es homogénea y nos referiremos a ella como ecuación de Laplace. Si es no homogénea, se la denomina ecuación de Poisson. Otras variaciones posibles del problema consisten en modificar las condiciones frontera. Aquellas que indican el valor de la solución que buscamos en la frontera son conocidas como condiciones de tipo Dirichlet. Otros tipos de condiciones son:

• Condiciones Neumann: [math]\nabla u \cdot \vec{n}= \partial_{n}u=g \in \partial\Omega [/math]

• Condiciones Robin: [math] u +k \partial_{n}u=g \in \partial\Omega[/math]

• Condiciones Mixtas: [math] \left\{ \begin{aligned} u &= g, & x &\in \Gamma \\ \partial_n u &= h, & x &\in \partial\Omega\setminus\Gamma \end{aligned} \right. [/math]

De ahora en adelante nos referiremos al par ecuación-condición frontera como problema de Laplace o Poisson según corresponda. A continuación enunciaremos dos resultados fundamentales en el estudio de este tipo de problemas.

[math]\textbf{Teorema: (Existencia de soluciones del problema de Laplace en la bola)}.[/math] Sea el problema [math] \left\{ \begin{aligned} \Delta u &= 0, & x &\in B_R(0) \\ u &= g, & x &\in \partial B_R(0) \end{aligned} \right. [/math] donde [math]B_R(x_0) \subset \mathbb{R}^n[/math] es la bola de radio [math]R[/math] centrada en [math]x_0 \in \mathbb{R}^n[/math], [math]g[/math] una función continua y [math]n\geq2[/math]. Entonces la solución [math]u \in \mathcal{C}^2(B_R(0))\cap \mathcal{C}(\overline{B_R(0)})[/math] viene dada por la fórmula de Poisson:

donde [math]\omega_n= 2\frac{\pi^{\frac{n}{2}}}{\Gamma(\frac{n}{2})}[/math]. También se puede expresar en desarrollo en serie de Fourier:

[math]\textbf{Teorema: (Unicidad de soluciones del problema de Laplace)}.[/math] Sea el problema de Laplace o Poisson con cualquier tipo de condición frontera de las mencionadas anteriormente excepto las de Neumann. Entonces existe una única solución [math]u \in C^2(B_R(0))\cap C(\overline{B_R(0)}) [/math] del problema. En el caso de condiciones de tipo Neumann la solución es única salvo constante.

[math]\textbf{Definición: (Función Armónica)}.[/math] Diremos que una función [math]u \in C^2(\Omega)[/math] es armónica si [math]\Delta u =0[/math].

[math]\textbf{Teorema: (Desigualdad de Harnack)}.[/math] Sea u una función armónica y [math]u\geq 0[/math] en [math]\Omega \subset \mathbb{R}^n[/math]. Sea [math]\overline{B_R(z)} \subset \Omega[/math]. Entonces se verifica la desigualdad de Harnack:

Como veremos en interpretación desigualdad de Harnack, este resultado se puede interpretar como que todas las funciones armónicas que pasan por un punto están "encerradas" en una región muy concreta del espacio determinada por dos funciones armónicas "límite".

[math]\textbf{Teorema: (Principio del máximo para funciones armónicas)}[/math] Sean [math]\Omega \subset \mathbb{R}^n[/math] un abierto y [math]u \in \mathcal{C}(\Omega)[/math]. Si [math]u[/math] es armónica y alcanza un máximo o un mínimo en [math]\Omega[/math] entonces [math]u[/math] es constante. Como consecuencia si [math]\Omega[/math] es acotado y [math]u[/math] no es constante entonces.

Un último resultado importante es el cálculo del error para el método del trapecio de integración numérica. Se puede demostrar que el error en este método para la integral [math]\int_a^bf(x)dx[/math] es:

donde [math]\xi\in (a,b)[/math]. Se puede ver una demostración de esto en Error regla del trapecio compuesta

3 Limitaciones de la fórmula de Poisson

Lo primero que hay que estudiar una vez llegamos a un resultado es su ámbito de aplicación. Por ejemplo, debemos cuestionarnos si las hipótesis planteadas son más fuertes de las necesarias. En caso afirmativo implicaría que la conclusión es en realidad más general y la podemos aplicar en más casos. Otra vía de estudio es, dicho de manera informal, cómo de buena es nuestra solución. Un ejemplo sencillo de esto podría ser el teorema de la convergencia de las series de Fourier. Este nos dice que la aproximación en serie de Fourier converge puntualmente allí donde la función a aproximar es continua y al punto medio en las discontinuidades ^[1]. Además, en caso de ser periódica y continua la serie converge uniformemente. Por tanto, nos clasifica cómo de buena es la aproximación, es decir nuestra solución, en diferentes situaciones en las que las hipótesis previas se verifican.

Veamos entonces cómo de buena es la fórmula de Poisson que presentamos en preliminares. A la vista de sus hipótesis y conclusiones, se podría pensar que al haber obtenido una fórmula explicita que podemos aplicar en cualquier caso, el problema en la bola quedaría totalmente resuelto. Sin embargo, esta presenta algunos inconvenientes que se analizarán en las secciones posteriores. Como bien dice el teorema, la solución calculada por esta vía es [math]\mathcal{C}^2(B_R(0))\cap \mathcal{C}(\overline{B_R(0)})[/math]. Sin embargo, salta a la vista que aunque la fórmula efectivamente nos da una solución con la regularidad que buscamos, en muchos casos es muy difícil resolver la integral. Por ello, en la práctica recurrimos a aproximaciones mediante métodos numéricos. Es aquí donde aparece una limitación de la fórmula de Poisson y que estudiaremos en esta sección.
Planteemos el siguiente problema de Laplace en [math]\mathbb{R}^2[/math]:

Como ya adelantamos podríamos encontrarnos integrales que no fueran fácilmente resolubles analíticamente por lo que recurrimos a métodos numéricos. Aproximando la integral por el método del trapecio obtenemos la gráfica de la solución aproximada. Se muestra además el código con el que se ha calculado y dibujado la solución:

Solución del Problema

clear all
close all
Npuntos=50;                            %Número de puntos para la aproximación del  
Nr=30; Nt=60;                           %Número de puntos para el graficado
tint=linspace(0,2*pi,Npuntos);          %Vector de puntos para evaluar en la regla del trapecio
g=@(t)(max(0,1-2/pi*abs(t-pi)));        %Condición frontera
u=@(r,t)((1-r^2)/(2*pi)*trapz(tint,g(tint)./(r^2+1-2*r*cos(tint-t))));
                                        %Expresión de la solución por la fórmula de Poisson
rr=linspace(0,0.95,Nr);                 %Vector de puntos para evaluar en coordenada radial
tt=linspace(0,2*pi,Nt);                 %Vector de puntos para evaluar en coordenada angular
U=zeros(length(rr),length(tt));         %Inicialización de la matriz de evaluación
for i=1:length(rr)
    for j=1:length(tt)
        U(i,j)=u(rr(i),tt(j));          %Evaluación de la solución por fórmula de Poisson para cada punto
    end
end
[RR,TT]=meshgrid(rr,tt);                %Formación de la malla de puntos
X=RR.*cos(TT);                          %Cambio a cartesianas
Y=RR.*sin(TT);
mesh(X,Y,U')
colormap turbo                          %Graficado de la función
xlabel('$x_1$','Interpreter','latex')
ylabel('$x_2$','Interpreter','latex')

En la gráfica se aprecia cómo la aproximación es cada vez peor según nos acercamos a la frontera. La razón de esto se encuentra en que el integrando tiene una singularidad en la frontera de la bola unidad. Como se puede ver, esto provoca que la solución obtenida numéricamente diverja según se acerca a la frontera. Por lo que los problemas que aparecen en la integral exacta únicamente en la frontera, al hacer una aproximación numérica los encontramos también en el interior. Una forma de mitigar este problema es haciendo la malla más fina, sin embargo esto no evita que eventualmente la solución termine divergiendo, aspecto que se estudiará a continuación. En conclusión, en muchos casos no tenemos más remedio que usar estas aproximaciones y debemos imponer directamente la condición frontera en la solución.

Veamos con ayuda de un problema del que conocemos la solución explícita, como haciendo la malla más fina conseguimos reducir los problemas generados al aproximar numéricamente.

Es fácil ver que la ecuación de Laplace implica que la solución debe ser armónica. Por otro lado, podemos comprobar que la función [math]u(x,y)=xy[/math] es armónica al ser un monomio multivariante de grado [math]1[/math]. Es decir, esta es solución del problema al ser igual a la expresión de la condición frontera. Como hemos visto en preliminares esta solución es única, por lo que podemos hallar una aproximación numéricamente y al conocer la solución exacta, conocer el error que estamos cometiendo.

Para estudiar el error, fijemos el punto [math](r,\theta)=(0.9,\frac{\pi}{4})[/math]. Al ser un punto que está lejos de la frontera podemos intuir que no tendrá demasiados problemas en la convergencia de la solución en su entorno; la malla empleada no deberá ser excesivamente fina para evitar dicho problema. Veamos este aspecto graficando el error en dicho punto en escala logarítmica en función del número de puntos empleados en la fórmula del trapecio.

Logaritmo del error en un punto fijo con [math]10^n[/math] puntos

clear all
close all
listaNpuntos=10.^[1:8];                     %Número de puntos de la discretización para aproximar la integral por la regla del trapecio
err=zeros(1,length(listaNpuntos));          %Inicialización del vector error (con logaritmo)
for i=1:length(listaNpuntos)
    Npuntos=listaNpuntos(i);                %Fijamos el número de puntos de la discretización a la componente i-ésima del vector listaNpuntos
    tint=linspace(0,2*pi,Npuntos);          %Partición del intervalo [0,2pi] con Npuntos puntos para la aproximación de la integral
    uexacta=@(r,t)(sin(t).*cos(t).*r.^2);   %Solución exacta del problema (en polares)
    g=@(t)(cos(t).*sin(t));                 %Condición frontera del problema (en polares)
    u=@(r,t)((1-r^2)/(2*pi)*trapz(tint,g(tint)./(r^2+1-2*r*cos(tint-t)))); 
                                            %Función u solución del problema calculada con la fórmula de Poisson
    err(i)=log10(abs(uexacta(0.9,pi/4)-u(0.9,pi/4)));
                                            %Cálculo del logaritmo del error puntual en valor absoluto 
end
plot(log10(listaNpuntos),err,'b')           %Gráfico del error en función del número de puntos
xlabel('$n$','Interpreter','latex')
ylabel('$\log_{10}(error(10^n))$','Interpreter','latex')

Podemos ver cómo el error se comporta como es de esperar. Para un punto lejano a la frontera basta aproximar la integral con [math]10^{-3}[/math] puntos para que la aproximación sea prácticamente exacta. Por lo que para puntos no muy cercanos a la frontera, como hemos mencionado, basta aumentar el número de puntos involucrados en la fórmula del trapecio.

Veamos ahora el problema del aumento error de la aproximación según nos acercamos a la frontera de la bola. Para ello fijamos el número de puntos empleados en la fórmula del trapecio. Para los puntos fijemos el ángulo a [math]\frac{\pi}{4}[/math] y calculemos el error variando [math]r[/math] para cada punto de la forma [math]\lbrace(1-10^{-n},\frac{\pi}{4})\rbrace_{n\in\mathbb{N}\setminus \lbrace 1 \rbrace}[/math].

Logaritmo del error según nos acercamos a la frontera

clear all
close all
Npuntos=1000;                               %Número de puntos de la discretización para aproximar la integral por la regla del trapecio
tint=linspace(0,2*pi,Npuntos);              %Partición del intervalo [0,2pi] con Npuntos puntos para la aproximación de la integral
uexacta=@(r,t)(sin(t).*cos(t).*r.^2);       %Solución exacta del problema (en polares)
g=@(t)(cos(t).*sin(t));                     %Condición frontera del problema (en polares)
u=@(r,t)((1-r^2)/(2*pi)*trapz(tint,g(tint)./(r^2+1-2*r*cos(tint-t))));
                                            %Función u solución del problema calculada con la fórmula de Poisson
nn=10.^-[2:8];                              %Vector de lo que le vamos a restar a 1 para ir acercandonos a la frontera
err=zeros(1,length(nn));                    %Inicialización del vector de error
for i=1:length(nn)
    err(i)=log10(abs(uexacta(1-nn(i),pi/4)-u(1-nn(i),pi/4))); 
                                            %Cálculo del logaritmo del error puntual en valor absoluto
end
plot(-[log10(nn)],err,'b')                  %Gráfico del error en función de n con x=1-nn
xlabel('$n$','Interpreter','latex')
ylabel('$\log_{10}(errorpunt(1-10^{-n}))$','Interpreter','latex')

Se puede observar también que según nos vamos acercando a la frontera los errores son cada vez mayores, como ya adelantamos.

Por otro lado, como hemos visto en preliminares, podemos estimar el error cometido pues conocemos la fórmula del error en el método del trapecio. Calculemos dicha cota.

Particularizando la fórmula de Poisson para este problema tenemos que :

Es fácil ver que expresando la integral en coordenadas polares obtenemos:

Por tanto, basta con estimar el error en esta última integral siguiendo la fórmula del error del método del trapecio Error regla del trapecio compuesta Obteniendo

Para algún [math]\xi \in (0,2\pi)[/math]. Ahora, derivando tenemos que

Finalmente, podemos acotar el error total multiplicando por la constante que falta:

Podemos emplear esta cota en los errores calculados anteriormente

Comparación error real con la cota

clear all
close all
listaNpuntos=10.^[1:8];                     %Número de puntos de la discretización para aproximar la integral por la regla del trapecio
err=zeros(1,length(listaNpuntos));          %Inicialización del vector error (con logaritmo)
for i=1:length(listaNpuntos)
    Npuntos=listaNpuntos(i);                %Fijamos el número de puntos de la discretización a la componente i-ésima del vector listaNpuntos
    tint=linspace(0,2*pi,Npuntos);          %Partición del intervalo [0,2pi] con Npuntos puntos para la aproximación de la integral
    uexacta=@(r,t)(sin(t).*cos(t).*r.^2);   %Solución exacta del problema (en polares)
    g=@(t)(cos(t).*sin(t));                 %Condición frontera del problema (en polares)
    u=@(r,t)((1-r^2)/(2*pi)*trapz(tint,g(tint)./(r^2+1-2*r*cos(tint-t)))); 
                                            %Función u solución del problema calculada con la fórmula de Poisson
    err(i)=log10(abs(uexacta(0.9,pi/4)-u(0.9,pi/4)));
                                            %Cálculo del logaritmo del error puntual en valor absoluto 
end
errcota=@(r,n)((20*pi^3.*(1-r.^2))./(2*pi*3.*n.^2.*(r-1).^6));
                                            %Cota de error de la fórmula del trapecio
plot(log10(listaNpuntos),err,'b')           %Gráfico del error en función del número de puntos

hold on
plot(log10(listaNpuntos),log10(errcota(0.9,listaNpuntos)),'r') 
hold off
xlabel('$n$','Interpreter','latex')
ylabel('$\log_{10}(error(10^n))$','Interpreter','latex')
legend('Error calculado','Cota del error','Interpreter','latex')

Comparación error real con la cota

clear all
close all
Npuntos=1000;                               %Número de puntos de la discretización para aproximar la integral por la regla del trapecio
tint=linspace(0,2*pi,Npuntos);              %Partición del intervalo [0,2pi] con Npuntos puntos para la aproximación de la integral
uexacta=@(r,t)(sin(t).*cos(t).*r.^2);       %Solución exacta del problema (en polares)
g=@(t)(cos(t).*sin(t));                     %Condición frontera del problema (en polares)
u=@(r,t)((1-r^2)/(2*pi)*trapz(tint,g(tint)./(r^2+1-2*r*cos(tint-t))));
                                            %Función u solución del problema calculada con la fórmula de Poisson
nn=10.^-[2:8];                              %Vector de lo que le vamos a restar a 1 para ir acercandonos a la frontera
err=zeros(1,length(nn));                    %Inicialización del vector de error
for i=1:length(nn)
    err(i)=log10(abs(uexacta(1-nn(i),pi/4)-u(1-nn(i),pi/4))); 
                                            %Cálculo del logaritmo del error puntual en valor absoluto
end
errcota=@(r,n)((20*pi^3.*(1-r.^2))./(2*pi*3.*n.^2.*(r-1).^6));
                                            %Cota de error de la fórmula del trapecio
plot(-[log10(nn)],err,'b')                  %Gráfico del error en función de n con x=1-nn
hold on
plot(-[log10(nn)],log10(errcota(1-nn,Npuntos)),'r') 
hold off
xlabel('$n$','Interpreter','latex')
ylabel('$\log_{10}(errorpunt(1-10^{-n}))$','Interpreter','latex')
legend('Error calculado','Cota del error','Interpreter','latex','Location','northwest')

3.1 Solución por serie de Fourier

Una vez mostrados los errores que conlleva el uso de la fórmula de Poisson sobre todo en la frontera del dominio, calcularemos la solución del problema [math] \left\{ \begin{aligned} \Delta u &= 0, & x &\in B_1(0) \\ u &= xy , & x &\in \partial B_1(0) \end{aligned} \right. [/math] esta vez sin usar la fórmula. La manera de calcularla ahora será mediante la serie de Fourier. Como se ve en preliminares, la solución de este problema en polares viene dado por la expresión

donde [math]a_{0},\lbrace a_{k}\rbrace_{k=1}^{\infty},\lbrace b_{k}\rbrace_{k=1}^{\infty}[/math] son los coeficientes de Fourier de la función [math]\mathcal{G}(\theta)=\frac{1}{2}\sin(2\theta)[/math], que es la condición frontera en coordenadas polares. En este caso se puede ver fácilmente que si la expresión en serie de Fourier de [math]\mathcal{G}[/math] es

por unicidad de los coeficientes, se tendría que [math]a_{0}=0,\lbrace a_{k}\rbrace_{k=1}^{\infty}=0,\lbrace b_{k}\rbrace_{k=1,k\neq 2}^{\infty}=0, a_{2}=\frac{1}{2}[/math]. Por tanto, se tiene de manera exacta que la solución del problema en polares es

Nótese cómo la solución ahora no tiene ningún problema en la frontera del dominio, en contrapartida con la dada por la fórmula de Poisson.

3.2 Interpretación desigualdad de Harnack

Como se adeltantó en introducción uno de los objetivos del trabajo era darle una interpretación a la desigualdad de Harnack. Además, obtendremos como cosecuencia inmediata el teorema de Liouville. Para hacerlo, estudiaremos esta desigualdad en el ejemplo anterior.

En primer lugar, hay que recordar que Harnack es cierta si [math]u\geq 0[/math]. Es fácil comprobar que la solución que hemos obtenido no verifica esto, [math]\mathcal{U}(1,\frac{3\pi}{4})=-\frac{1}{2}\lt0[/math]. Sin embargo, lo que podemos hacer es hallar el mínimo de la función y trasladarla de forma que la función obtenida es positiva. Nótese que al estar trabajando en un compacto y con una función continua siempre va a existir dicho mínimo. Vamos a hallarlo:

Recordemos que en el ejemplo anterior, habíamos visto que [math]u(x,y)=xy[/math] era la solución al problema, por lo que es una función armónica. Por el principio del máximo en funciones armónicas preliminares, el máximo en [math]B_1(0)[/math]debe estar en su frontera. Por tanto basta escribir la función en coordenadas polares [math]\mathcal{U}(r,\theta)=\frac{1}{2}r^{2}\sin(2\theta)[/math] y buscar los extremos en [math]\partial B_1(0)[/math]. Derivando:

Tomando [math]\theta\in [0,2\pi)[/math] vemos que hay un mínimo en [math]\theta =\frac{3\pi}{4}[/math] donde [math]\mathcal{U}(1,\frac{3\pi}{4})=-\frac{1}{2}[/math]. Así, podemos definir la siguiente función

Como mencionamos, esta función es la trasladada de forma que [math]\mathcal{V}_{1}[/math] es positiva, por lo que sí podemos aplicar Harnack en el punto [math](0)[/math]:

A continuación se muestra la gráfica en la que se pueden ver las funciones que acotan [math]\mathcal{V}_{1}[/math] en escala logarítmica.

Funciones límite desigualdad de Harnack en [math]B_1(0)[/math]

clear all
close all
M=-1/2; N=100;                          %Definición del valor máximo y del número de puntos 
u=@(r,t)(r.^2)*(cos(t).*sin(t));        %Definición de la solución
v=@(r,t)(r.^2)*(cos(t).*sin(t))-M;      %Definición de v:=u-M
vmax=@(r) v(0,0)*(1+r)./(1-r);          %Función cota superior de la desigualdad de Harnack
vmin=@(r) v(0,0)*(1-r)./(1+r);          %Función cota inferior de la desigualdad de Harnack
                                        %Graficado de las funciones
plot(linspace(0,0.99,N),log(vmax(linspace(0,0.99,N))),'b')
hold on
plot(linspace(0,0.99,N),log(vmin(linspace(0,0.99,N))),'b')
hold off
xlabel('$r$','Interpreter','latex')
legend('Cota con $R=1$','Location','Northwest','Interpreter','latex')

Ahora es fácil darle una interpretación a la desigualdad. Esta dice que toda función armónica que pase por, en este caso, [math]0[/math] está acotada superior e inferiormente por dichas funciones en [math]\overline{B_r(0)} \subset \Omega[/math]. Estas son también funciones armónicas, por tanto, nos indica que todas las funciones armónicas que pasan por dicho punto están encerradas en el área delimitada entre estas dos funciones armónicas "límite" que se ve en la gráfica anterior.

Una vez entendido esto, podemos dar un paso más. Veamos qué pasa si hacemos el radio de la bola compacta más grande. Tomemos por ejemplo [math]B_{2}(0)[/math] y [math]B_{10}(0)[/math]. Al igual que antes debemos trasladar la función para que sea positiva; hallando su mínimo que estará en las fronteras de cada caso. Como la función es la misma que antes el cálculo del mínimo es igual, se encuentra en [math]\frac{3\pi}{4}[/math]. Sin embargo, recordemos que [math]\mathcal{U}(r,\theta)=\frac{1}{2}r^{2}\sin(2\theta)[/math], por lo que dicho valor disminuirá según el radio. Concretamente, [math]\mathcal{U}(r,\theta)= -2[/math] y [math]\mathcal{U}(r,\theta)= -50[/math] respectivamente. Ahora, podemos definir las funciones a las que aplicaremos Harnack:

Las cotas en estos dos casos quedarán:

donde la primera corresponde a [math]B_2(0)[/math] y la segunda a [math]B_{10}(0)[/math]. A continuación se muestran las funciones "límite" en los tres casos:

Funciones límite desigualdad de Harnack en [math]B_1(0),B_5(0)[/math] y [math]B_{10}(0)[/math]

clear all
close all
n=2;                                                            %Dimensión del problema
Rs=[1,5,10];                                                    %Vector de radios
colores=['r','b','g'];                                          %Vector de colores (para graficar)
for i=1:length(Rs)
    R=Rs(i);                                                    %Fijamos el radio a la componente i de Rs
    u=@(r,t)(r.^2)*(cos(t).*sin(t));
    M=@(R)-R.^2/2; N=100;                                       %Calculamos el mínimo de la función u
    v=@(r,t)(r.^2)*(cos(t).*sin(t))-M(R);                       %Definimos v:=u-M
    harnmax=@(R,n,r) v(0,0)*R.^(n-2).*(R-r)./((R+r).^(n-1));    %Función cota superior de la desigualdad de Harnack
    harnmin=@(R,n,r) v(0,0)*R.^(n-2).*(R+r)./((R-r).^(n-1));    %Función cota inferior de la desigualdad de Harnack
                                                                %Graficado de las cotas 
    plot(linspace(0,R-0.01,N),log(harnmax(R,n,linspace(0,R-0.01,N))),colores(i))
    hold on
    plot(linspace(0,R-0.01,N),log(harnmin(R,n,linspace(0,R-0.01,N))),colores(i))
end
hold off
xlabel('$r$','Interpreter','latex')
legend({'Cota con $R=1$','','Cota con $R=5$','','Cota con $R=10$',''},'Location','northwest','Interpreter','latex')

Para poder ver mejor su comportamiento vamos a desplazar las gráficas de tal manera que quedan de la siguiente forma:

Funciones límite desigualdad de Harnack trasladadas en [math]B_1(0),B_5(0)[/math] y [math]B_{10}(0)[/math]

clear all
close all
n=2;                                                            %Dimensión del problema
Rs=[1,5,10];                                                    %Vector de radios
colores=['r','b','g'];                                          %Vector de colores (para graficar)
for i=1:length(Rs)
    R=Rs(i);                                                    %Fijamos el radio a la componente i de Rs
    u=@(r,t)(r.^2)*(cos(t).*sin(t));
    M=@(R)-R.^2/2; N=100;                                       %Calculamos el mínimo de la función u
    v=@(r,t)(r.^2)*(cos(t).*sin(t))-M(R);                       %Definimos v:=u-M
    harnmax=@(R,n,r) v(0,0)*R.^(n-2).*(R-r)./((R+r).^(n-1));    %Función cota superior de la desigualdad de Harnack
    harnmin=@(R,n,r) v(0,0)*R.^(n-2).*(R+r)./((R-r).^(n-1));    %Función cota inferior de la desigualdad de Harnack
                                                                %Graficado de las cotas 
    plot(linspace(0,R-0.01,N),-log(v(0,0))*ones(1,length(linspace(0,R-0.01,N)))+log(harnmax(R,n,linspace(0,R-0.01,N))),colores(i))
    hold on
    plot(linspace(0,R-0.01,N),-log(v(0,0))*ones(1,length(linspace(0,R-0.01,N)))+log(harnmin(R,n,linspace(0,R-0.01,N))),colores(i))
end
hold off
xlabel('$r$','Interpreter','latex')
legend({'Cota con $R=1$','','Cota con $R=5$','','Cota con $R=10$',''},'Interpreter','latex')

Es fácil apreciar que según tomamos bolas de radio más grande el área que encierran las funciones límite se va haciendo más pequeña. Esto nos puede dar la idea de que si tomamos el límite [math]R \to \infty[/math] el área se hará tan pequeña que tendremos la función completamente determinada. Es más, viendo el comportamiento de la gráfica es razonable suponer que dicha función sería una constante. Para ver esto podemos tomar el límite [math]R \to \infty[/math] a ambos lados de la expresión. Sabemos que

Tomando límites $R \to \infty$ a los dos lados :

Como intuíamos llegamos a que la función [math]\mathcal{V}(x)[/math] es constante y por tanto [math]\mathcal{U}(x)[/math] también. Este es un resultado famoso que se conoce como teorema de Liouville.

[math]\textbf{Teorema de Liouville}:[/math] Sea [math]u[/math] una función acotada y armónica en [math]\mathbb{R}^n[/math] entonces [math]u[/math] es constante.

Funciones límite desigualdad de Harnack trasladadas en [math]B_1(0),B_5(0)[/math] y [math]B_{10}(0)[/math] en [math]\mathbb{R}^3[/math]

clear all
close all
n=3;                                                            %Dimensión del problema
Rs=[1,5,10];                                                    %Vector de radios
colores=['r','b','g'];                                          %Vector de colores (para graficar)
for i=1:length(Rs)
    R=Rs(i);                                                    %Fijamos el radio a la componente i de Rs
    u=@(r,t)(r.^2)*(cos(t).*sin(t));
    M=@(R)-R.^2/2; N=100;                                       %Calculamos el mínimo de la función u
    v=@(r,t)(r.^2)*(cos(t).*sin(t))-M(R);                       %Definimos v:=u-M
    harnmax=@(R,n,r) R.^(n-2).*(R-r)./((R+r).^(n-1));           %Función cota superior de la desigualdad de Harnack
    harnmin=@(R,n,r) R.^(n-2).*(R+r)./((R-r).^(n-1));           %Función cota inferior de la desigualdad de Harnack
                                                                %Graficado de las cotas 
    plot(linspace(0,R-0.01,N),log(harnmax(R,n,linspace(0,R-0.01,N))),colores(i))
    hold on
    plot(linspace(0,R-0.01,N),log(harnmin(R,n,linspace(0,R-0.01,N))),colores(i))
end
hold off
xlabel('$r$','Interpreter','latex')
legend({'Cota con $R=1$','','Cota con $R=5$','','Cota con $R=10$',''},'Interpreter','latex')

Lo más llamativo de estas gráficas es como a pesar de que el comportamiento es parecido al de dimensión [math]2[/math], el eje [math]y[/math] se ve modificado. Es decir que las gráficas cuanto mayor es la dimensión más se abren, peores parecen las cotas dadas por la desigualdad de Harnack. A pesar de esto, el Teorema de Liouville se sigue cumpliendo en dimensiones superiores.

4 Ecuación de Poisson en [math]\mathbb{R}^2[/math]

Hasta ahora, siempre se ha estudiado el problema de Laplace en regiones acotadas. A diferencia de las secciones anteriores, ahora se estudiará el caso en un dominio no acotado, concretamente [math]\mathbb{R}^{2}[/math]. Planteemos por tanto el siguiente problema

Cuya solución viene dada por la convolución con la solución fundamental. En dimensión [math]2[/math] la solución fundamental tiene la expresión [math]\Phi(x)=-\frac{1}{2\pi}\log(\left| x \right|)[/math], por lo que concluimos que la solución del problema viene dado por el llamado potencial logarítmico:

Pongamos ahora un ejemplo para ver la apariencia de las soluciones en [math]\mathbb{R}^2[/math], supongamos el siguiente problema

donde [math]1_{B_{1}}[/math] es la función característica en [math]B_{1}[/math], es decir, la función que vale [math]1[/math] en [math]x\in \mathbb{B}_{1}[/math] y [math]0[/math] en [math]x\notin B_{1}[/math]. De manera que la solución vendrá dada por

Haciendo un cambio a coordenadas polares resultaría en la integral

Integrándola numéricamente con Matlab podemos obtener su gráfica.

Solución del problema en [math]\mathbb{R}^2[/math] dada por el potencial logarítmico

clear all
close all
Nr=200; Nt=200; Nr0=100; Nt0=100;                   %Número de puntos de las discretizaciones
F=@(r0,t0,rr,tt)(-1/(4*pi))*log(r0^2+rr.^2-2*r0*rr.*cos(tt-t0)).*rr;
                                                    %Fórmula del potencial logaritmico en polares
rr=linspace(0,0.99999,Nr); tt=linspace(0,2*pi,Nt);  %Vectores discretización radio y angulo para la integral
[RR,TT]=meshgrid(rr,tt);                            %Construcción de la malla para la integral
rr0=linspace(0,10,Nr0); tt0=linspace(0,2*pi,Nt0);   %Vectores discretización radio y angulo para la solución
[RR0,TT0]=meshgrid(rr0,tt0);                        %Construcción de la malla para la solución
U=zeros(Nr0,Nt0);                                   %Inicialización de la matriz de la función solución discretizada
for i=1:Nr0
    for j=1:Nt0
         U(i,j)=trapz(tt,trapz(rr,F(RR0(i,j),TT0(i,j),RR,TT),2));
                                                    %Evaluación de la matriz de la función solución discretizada
    end
end
X0=RR0.*cos(TT0); Y0=RR0.*sin(TT0);                 %Cambio a coordenadas cartesianas
mesh(X0,Y0,U)                                       %Graficado de la función
colormap turbo
xlabel('$x_1$','Interpreter','latex')
ylabel('$x_2$','Interpreter','latex')

Observamos en esta dos zonas diferenciadas. En primer lugar la zona central, donde nos encontramos un máximo. El comportamiento de la solución en esta región viene determinado por la función [math]f[/math], que dicta que en [math]B_{1}(0)[/math], [math]\Delta u=1[/math]. Geométricamente esto dice que para todo punto de esa región, la función vale más en el punto que en su entorno en promedio. Es decir, la función es superarmónica. Este comportamiento se ve más claramente en el punto máximo. La otra zona que destaca es el exterior de [math]B_{1}(0)[/math]. Aquí el comportamiento por su parte lo dicta que [math]\Delta u=0[/math], es decir, que la función es armónica. Geométricamente esto se puede ver en que para todo punto de la región, la función toma el valor promedio de los puntos de su entorno.

Veamos por otro lado el comportamiento de la solución en un punto lejano al origen, cuando [math]\lvert x\rvert\to\infty[/math]. Por teoría se tiene ^[2] que cuando [math]\lvert x\rvert\to\infty[/math], la solución del problema se comporta de manera que

donde [math]M=\int_{\mathbb{R}^2}f(x)\thinspace dx=\lvert B_{1}(0)\rvert=\pi[/math]. Graficando tanto la expresión anterior como la obtenida con el potencial logarítmico en función de [math]r[/math] y para un [math]\theta=0[/math] fijo, se tiene que las expresiones cada vez se parecen más, como es de esperar. Observamos también cómo el valor de la función cuanto mayor es el radio cada vez es menor.

Solución del problema en función de [math]r[/math] para un [math]\theta[/math] fijo y comparación con la estimación asintótica

clear all
close all
Nr=200; Nt=200; Nr0=200;                            %Número de puntos de las discretizaciones
F=@(r0,t0,rr,tt)(-1/(4*pi))*log(r0^2+rr.^2-2*r0*rr.*cos(tt-t0)).*rr;
                                                    %Fórmula del potencial logaritmico en polares
rr=linspace(0,0.99999,Nr); tt=linspace(0,2*pi,Nt);  %Vectores discretización radio y angulo para la integral
[RR,TT]=meshgrid(rr,tt);                            %Construcción de la malla para la integral
rr0=linspace(0,Nr0,Nr0);                            %Vector discretización radio
U=zeros(1,Nr0);                                     %Inicialización del vector de la función solución discretizada
for i=1:Nr0
    U(i)=trapz(tt,trapz(rr,F(rr0(i),0,RR,TT),2));   %Evaluación del vector de la función solución discretizada
end
uinfty=@(r)(-1/2)*log(abs(r))+1./r;                 %Definición de la solución asintótica calculada
plot(rr0,U,'b')                                     %Graficado de las funciones
hold on
plot(rr0,uinfty(rr0),'r')
hold off
xlabel('$r$','Interpreter','latex')
legend('$\mathcal{U}(r,0)$','$u_{\infty}$','Cota del error','Interpreter','latex')

4.1 Referencias

1. ↑ Salsa, S., Verzini, G. (2008). Partial Differential Equations in Action: From Modelling to Theory. Pg 535. Alemania: Springer International Publishing.
2. ↑ Salsa, S., Verzini, G. (2022). Partial Differential Equations in Action: From Modelling to Theory. Pg 150. Alemania: Springer International Publishing.