Revisión del 11:45 19 abr 2024

1 Introducción

En el trabajo que se presenta a continuación, vamos a estudiar la solución de la ecuación de Laplace, analizando su comportamiento mediante la representación numérica y los errores que conlleva. Además, vamos a estudiar la solución de la ecuación de Poisson mediante la solución fundamental.

Cabe destacar que los códigos y gráficas que se presentan a continuación han sido realizados en Matlab.

2 Conceptos previos

En esta sección, vamos a presentar algunos conceptos esenciales para comprender el trabajo que se presenta.

Funciones armónicas: Sea [math]f: \Omega \rightarrow \mathbb{R} [/math], siendo [math] \Omega[/math] un subconjunto [math]\mathbb{R }^n[/math] se dice armónica en [math] \Omega[/math] si cumple que sobre [math] \Omega[/math] tiene derivadas parciales continuas de primer y segundo orden y satisface la ecuación de Laplace:

Esto es equivalente a cumplir la propiedad del valor medio:

Fórmula de Poisson: La fórmula de Poisson determina la única solución de la ecuación de Laplace con [math] \Omega = B_r (\textbf{x}_0) [/math] y condición frontera [math] u(\textbf{x}) = g(\textbf{x}), \textbf{x} \in \partial B_r (\textbf{x}_0) [/math]. Su expresión es:

Siendo [math] \omega_n [/math] la medida de la frontera de la bola.

Desigualdad de Harnack: Sea [math]u[/math] armónica y [math]u\geq 0 [/math] en [math] \Omega \subset \mathbb{R}^n [/math]. Supongamos [math] \overline{B_R (\textbf{z})} \subset \Omega [/math]. Entonces [math] \forall \textbf{x} \in B_R(\textbf{z}), ~ r=|\textbf{z}-\textbf{x}|[/math] se tiene que

Principio del máximo: Sea [math]\Omega \subset \mathbb{R}^n [/math] y sea [math] u \in C( \Omega ) [/math]. Si [math] u [/math] satisface la propiedad del valor medio 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 [math] u(\textbf{x}) \lt \max \limits_{\partial ~\Omega} u [/math] y [math] u(\textbf{x}) \gt \min \limits_{\partial ~\Omega} u [/math], [math] \forall \textbf{x} \in \Omega [/math].

Además, en este trabajo vamos a emplear series de Fourier.

3 Solución de la ecuación de Laplace

En esta sección, vamos a estudiar la solución de la ecuación de Laplace considerando un problema concreto.

3.1 Solución mediante la fórmula de Poisson

Para obtener la solución de la ecuación de Laplace, vamos a considerar la bola unidad [math]B_1 \subset R^2[/math] centrada en [math] (0, 0) [/math] y el siguiente problema:

Siendo [math] g(\theta)= \max \left\{0, 1-\frac{2}{\pi} |\theta - \pi | \right\} [/math], en coordenadas polares.

Para obtener la solución vamos a emplear la fórmula de Poisson, que en este caso es de la forma:

Nótese que en este caso [math] R=1 [/math] y estamos trabajando en dimensión [math]n=2[/math], luego [math]\omega_n = \omega_2=2 \pi [/math].

Nuestro objetivo es representar esta solución. El problema reside en que necesitamos aplicar integración numérica en Matlab, y esto requiere una traducción de la fórmula a una expresión que podamos introducir en el programa. Por ello, debemos hacer un cambio a coordenadas polares. De esta manera, podremos resolver numéricamente la integral mediante el método del trapecio.

El cambio a polares en la fórmula de Poisson es el siguiente:

Además, observamos que cuando [math] r=R[/math] el denominador del integrando se anula. Esto va a conllevar un error numérico que trataremos más adelante. Para evitar esto a la hora de realizar la representación, como sabemos que la solución en la frontera es [math] g[/math], vamos a considerar [math] g[/math] en la frontera en lugar de los valores obtenidos numéricamente.

Teniendo esto en cuenta, la gráfica de la solución obtenida mediante la fórmula de Poisson es la siguiente:

Gráfica de la solución mediante la fórmula de Poisson en la bola unidad con condición frontera la función [math] g[/math].

Como podemos observar, la solución es continua y su valor en la frontera viene dado por la función [math] g[/math]. Además, se cumple que tanto el máximo como el mínimo, cuyos valores son [math] 1[/math] y [math] 0[/math] respectivamente, se alcanzan en la frontera. Este resultado es lógico de acuerdo al principio del máximo, pues se tiene que la solución no es constante y el conjunto sobre el que trabajamos es acotado.

3.1.1 Programa

Código para representar la solución mediante la fórmula de Poisson.

Nota: En este programa y en los siguientes la solución se obtiene de forma numérica. Discretizaremos el espacio mediante una malla en la cual obtenemos un valor aproximado de la solución. En los siguientes programas los comentarios indicarán los valores que se pueden modificar y qué cambian.

%%% EJERCICIO 1 %%%
 
clear all
close all
 
% Definimos los intervalos de integración y sus discretizaciones
 
th0=0; thf=2*pi; %Valor inicial y final de theta
r0=0; rf=1;      %Valor inicial y final del radio
 
divth=10^-2; % División del intervalo de theta
divr=10^-2;  % División del intervalo del radio
 
Th=th0:divth:thf; %Discretización del intervalo de theta
Rr=r0:divr:rf;    %Discretización del intervalo del radio
 
[TH,RR]=meshgrid(Th,Rr);
 
% Definimos la condición frontera
 
g=@(s)(1-2*abs(s-pi)/pi).*(1-2*abs(s-pi)/pi>0);  % Función en la frontera
G=g(Th);  % Valor de la función en la frontera (r=rf)
 
% Definimos la solución con la fórmula de Poisson
 
IntP=ones(size(TH)); % Iniciamos los datos que guardarán los valores de la integral 
for i=1:length(Th)
    for j=1:length(Rr)
        Aux=G./(rf^2+Rr(j)^2-2*rf*Rr(j)*cos(Th(i)-Th));
        IntP(j,i)=trapz(Th,Aux); % Evaluamos la integral por cada ángulo distinto
    end
end
 
U=(rf^2-RR.^2)/(2*pi).*IntP; % Obtenemos los valores de la función.
U(end,:)=G;                  % Valores en la frontera
 
% Cambio de variable
 
XX=RR.*cos(TH);
YY=RR.*sin(TH);
YY(:,end)=zeros(length(Rr),1); % sin(2pi) no es 0 por error numérico de matlab
 
% Dibujamos
 
surf(XX,YY,U)
title("Solución de la Ecuación de Laplace")
shading interp

3.2 Error numérico

Como ya se ha mencionado, hay un error de aproximación en la frontera de la bola debido al carácter singular del integrando de la fórmula de Poisson por la naturaleza numérica del método del trapecio. Debido a esto, tendremos que en la frontera el valor de la solución obtenida numéricamente es infinito.

Para estudiar el error en detalle, vamos a elegir la función armónica [math]u(x, y) = xy [/math], cuyo valor en la frontera también [math] g(x,y)=xy [/math].

Comenzamos calculando su solución mediante la fórmula de Poisson. Al igual que antes, empleando un cambio a polares se tiene que la solución es

con [math]g(\theta)= \frac{1}{2} sen(2 \theta)[/math], que es [math]g(x,y)[/math] en polares y en la frontera.

Ahora obtenemos el error numérico.

Los dos principales factores que van a determinar el error mediante el método del trapecio el tamaño de la discretización y la cercanía a la singularidad, que como ya hemos visto es en la frontera.

Para visualizar esto vamos a representar el error de la función tomando diferentes discretizaciones, considerando una sucesión de puntos que tiende a la frontera.

Para ello, vamos a obtener el error en escala logarítmica [math] f(n)=log_{10} error(10^n) [/math], siendo [math] error(10^n) [/math] el error para una discretización de tamaño [math] 10^n [/math]. Esto lo haremos para cada uno de los puntos de la sucesión.

Concretamente vamos a considerar la sucesión de puntos en coordenadas polares [math] \left\{ \left( 1-10^k, \frac{\pi}{4} \right) \right\}_{k=1}^6[/math].

La gráfica obtenida considerando lo anterior es:

Gráfica del error en escala logarítmica [math] f(n)=log_{10} error(10^n) [/math], siendo [math] error(10^n) [/math] el error en el punto [math] \left\{ \left( 1-10^k, \frac{\pi}{4} \right) \right\}_{k=1}^6 [/math] y considerando [math] 10^n [/math] puntos en la fórmula del trapecio.

Como vemos, en todos los casos el error se estabiliza en torno a un valor, que es mayor a medida que nos acercamos a la frontera. A su vez, el tiempo de convergencia a dicho valor también aumenta con el radio. Esto coincide con el comportamiento que esperábamos, ya que habíamos comentado que la singularidad del integrando en la frontera sería el principal motivo de error numérico.

Además, podemos apreciar que el error numérico disminuye a medida que el valor de [math]n[/math] crece, pues con ello aumenta el tamaño de la discretización.

Nota: En la gráfica se puede observar como en los últimos instantes de la [math]n[/math] el error numérico crece ligeramente, esto se debe a un error de cálculo de Matlab al tratar con valores muy pequeños.

Además, podemos estimar el error máximo del método del trapecio mediante la siguiente fórmula:

siendo [math] x_i [/math] un valor perteneciente al intervalo [math] [b-a] [/math].

En nuestro caso, el intervalo es [math][0, 2\pi][/math] y [math] f(\theta)=\frac{ sen(2 \theta)}{1+r^2-2rcos(\alpha - \theta)}[/math].

Calculando la segunda derivada de [math] f[/math] obtenemos:

A continuación, la acotamos maximizando el numerador y minimizando el denominador, de modo que la expresión resultante es:

Teniendo en cuenta que fuera de la integral tenemos [math] \frac{(1-r^2)}{4 \pi}[/math] y la fórmula del error máximo mencionada anteriormente, la cota del error es:

A continuación, comparamos los errores que hemos obtenido con la cota anterior. De modo que la gráfica es:

Gráfica de comparación de errores obtenidos con la cota de la fórmula del error. Se considera la [math] n [/math] de la fórmula como [math] 10^n [/math] puntos.

Se observa como los errores obtenidos son inferiores a la cota. Sin embargo, no parece ser la más óptima pues la diferencia entre los valores es considerablemente grande.

Por ello, vamos a maximizar la función del error en [math] \theta [/math] para cada [math]r[/math].

Observando el integrando en la fórmula de Poisson, se puede deducir que la única singularidad se encuentra en la frontera, concretamente cuando [math] \textbf{x}= \textbf{y}[/math]. Esta singularidad se mantiene en la segunda derivada respecto de [math] \theta [/math] del integrando.

Considerando el cambio a polares, [math] \textbf{x}=(r , \alpha)[/math] y [math] \textbf{y}= (1, \theta)[/math], se tiene que dicha singularidad se alcanza cuando [math] r=1[/math] y [math] \theta= \alpha[/math]. Por ello, vamos a considerar [math] 1-r \sim 0 [/math] y [math] \alpha - \theta \sim 0 [/math] y trabajaremos con infinitésimos.

Ahora vamos a fijar el valor de [math] \alpha [/math] y vamos a estudiar cómo se comporta [math] f(\theta)[/math] considerando estos infinitésimos. Sea [math]h(\theta) = 1+r^2-2rcos(\alpha - \theta) = (1-r)^2 + 2r(1-cos(\alpha - \theta)) \sim (1-r)^2 + r (\alpha - \theta)^2= o((1-r)^2) + o(r (\alpha - \theta)^2) [/math].

Obtenemos también la expresión en función de los infinitésimos de los siguientes términos de la fórmula de la derivada:

Si denotamos [math]num_i[/math], con [math]i=1,...,4[/math] a los numeradores del a segunda derivada de [math]f(\theta)[/math] y además consideramos denominador común, tenemos que

donde se tiene que:

[math] 1) ~~~num_1 \cdot h \sim \left \{ \begin{array}{ll} o((\alpha - \theta)^4) + o ( (1-r)^2 (\alpha - \theta)) & \quad \alpha \neq \left\{ \frac{\pi}{4} + k \frac{\pi}{2} \right\}_{k \in \mathbb{Z}} \\ o((\alpha - \theta)^5) + o ((1-r)^2 (\alpha - \theta)^2) & \quad \alpha = \left\{ \frac{\pi}{4} + k \frac{\pi}{2} \right\}_{k \in \mathbb{Z}} \end{array} \right. [/math]

[math] 2) ~~~(num_2+ num_3) \cdot h^2 \sim \left \{ \begin{array}{ll} o((\alpha - \theta)^4) + o ( (1-r)^4) + o ((\alpha - \theta)^2 (1-r)^2) & \quad \alpha \neq \left\{k \frac{\pi}{2} \right\}_{k \in \mathbb{Z}} \\ o((\alpha - \theta)^5) + o ((\alpha - \theta) (1-r)^4) + o ((\alpha - \theta)^3 (1-r)^2) & \quad \alpha \neq \left\{ k \frac{\pi}{2} \right\}_{k \in \mathbb{Z}} \end{array} \right. [/math]

[math] 3) ~~~num_4 \sim \left \{ \begin{array}{ll} o(\alpha - \theta) \end{array} \right. [/math]

Teniendo en cuenta que la suma de infinitésimos es el máximo de esos infinitésimos, tenemos que:

Finalmente, si fijamos [math]r[/math], el máximo se obtiene cuando [math] \theta= \alpha[/math], es decir, [math]máx |f''|= |f''(\alpha)|[/math]. Como en la sucesión de puntos hacia la frontera que hemos considerado el ángulo no varía, para cada [math]r[/math] la cota es [math]f(\alpha, r)[/math] (nótese que la función [math]f[/math] también depende de [math]r[/math]). Entonces la cota del error es:

Para visualizar esto, vamos a representar el error junto a las nuevas cotas:

Gráfica de comparación de errores obtenidos con la segunda cota del error.

Como podemos observar, el orden de la cota de error es considerablemente inferior a la que habíamos tomado antes y el error del método del trapecio sigue siendo menor

3.2.1 Programa

Código para representar el error del método del trapecio sin considerar la cota:

%%% EJERCICIO 2 SIN COTA %%%
 
clear all
close all
format long
 
% Definimos los intervalos de integración y sus discretizaciones
 
th0=0; thf=2*pi; %Valor inicial y final de theta
r0=0; rf=1;      %Valor inicial y final del radio
 
% Definimos la condición frontera
 
g=@(s,r)sin(2*s)/2;  % Función en la frontera
 
% Definimos la solución con la fórmula de Poisson
 
N=8; % N máx
K=6; % Máx dif
 
dif=10.^-(1:K);
Rs=1-dif; Ths=pi/4*ones(size(Rs)); % Puntos donde comprobaremos las diferencias
 
error=zeros(length(Rs),N);
 
for n=1:min(N,8) % Con 9 puede dar problemas!!
    Dth=linspace(th0,thf,10^n);
    G=g(Dth);  % Discretización para la integral
    for i=1:length(Rs)
        Aux=G./(rf^2+Rs(i)^2-2*rf*Rs(i)*cos(Ths(i)-Dth));
        U=(rf^2-Rs(i)^2)/(2*pi).*trapz(Dth,Aux); % Evaluamos la integral en el punto concreto
        error(i,n)=abs(U-Rs(i)^2*sin(2*Ths(i))/2);
    end
end
 
colors=["r","b","g","m","k","c"]; % Cambiar si se ponen más puntos
 
for i=1:K
    semilogy(1:N,error(i,:),colors(i),"LineWidth",2)
    hold on
    leyenda(i)="(1-10^{-"+num2str(i)+"}, \pi/4)" ;
end
title("Error del Método del Trapecio")
legend(leyenda)
xlabel("Discretización de tamaño 10^{n}")
ylabel("Error")

Código para comparar el error del método del trapecio con la primera cota:

%%% EJERCICIO 2 y 1ª Cota%%%
 
clear all
close all
format long
 
% Definimos los intervalos de integración y sus discretizaciones
 
th0=0; thf=2*pi; %Valor inicial y final de theta
r0=0; rf=1;      %Valor inicial y final del radio
 
% Definimos la condición frontera
 
g=@(s,r)sin(2*s)/2;  % Función en la frontera
 
% Definimos la solución con la fórmula de Poisson
 
N=8; % N máx
K=6;  % Máx dif
 
dif=10.^-(1:K);
Rs=1-dif; Ths=pi/4*ones(size(Rs)); %Puntos donde comprobaremos las diferencias
 
error=zeros(length(Rs),N);
 
for n=1:min(N,8)
    Dth=linspace(th0,thf,10^n);
    G=g(Dth);  % Discretización para la integral
    for i=1:length(Rs)
        Aux=G./(rf^2+Rs(i)^2-2*rf*Rs(i)*cos(Ths(i)-Dth));
        U=(rf^2-Rs(i)^2)/(2*pi).*trapz(Dth,Aux); % Evaluamos la integral en el punto concreto
        error(i,n)=abs(U-Rs(i)^2*sin(2*Ths(i))/2);
    end
end
 
k=1:K;
 
% Establecemos la cota
 
C=(2*pi)^3/12;
 
n=1:N;
 
[NN,RR]=meshgrid(n,Rs);
Acot=@(r)(r.^6-r.^5+6*r.^3+r.^2-r+2)./(2*pi*(1-r).^5);
f=@(r,n)(C*Acot(r).*(10.^(-2*n)));
 
Cota=f(RR,NN);
colors=["r","b","g","m","k","c"]; % Cambiar si se ponen más puntos
for i=k
    semilogy(n,error(i,:),colors(i),"LineWidth",2) 
    hold on
    leyenda(i)="(1-10^{-"+num2str(i)+"}, \pi/4) y su Cota" ;
end
for i=k
    semilogy(n,Cota(i,:),"--"+colors(i),"LineWidth",2)
    hold on
end
title("Error del Método del Trapecio y 1ª Cota")
legend(leyenda)
xlabel("Discretización de tamaño 10^{n}")
ylabel("Error")

Código para comparar el error del método del trapecio con la segunda cota:

%%% EJERCICIO 2 y Cota del Máximo %%%
 
clear all
close all
format long
 
% Definimos los intervalos de integración y sus discretizaciones
 
th0=0; thf=2*pi; %Valor inicial y final de theta
r0=0; rf=1;      %Valor inicial y final del radio
 
% Definimos la solución con la fórmula de Poisson
 
N=8; % N máx
K=6;
dif=10.^-(1:K);
Rs=1-dif; alpha=pi/4; %Puntos donde comprobaremos las diferencias
 
U=zeros(length(dif),N);
error=zeros(length(dif),N);
for n=1:N
    Th=linspace(0,2*pi,10^n);
    for i=1:K
        Int=sin(2*Th)./(1+Rs(i)^2-2*Rs(i)*cos(Th-alpha));
        Aux=trapz(Th,Int);
        U=(1-Rs(i)^2).*Aux/(4*pi);
        error(i,n)=abs(U-Rs(i)^2*sin(2*alpha)/2);
    end
end
 
syms r t
 
f=sin(2*t)./(1+r.^2-2*r.*cos(t-alpha)); % Integrando simbólico
 
dft=diff(f,t); % 1ª derivada
 
dft2=diff(dft,t); % 2ª derivada
 
Dft2=matlabFunction(dft2); % La pasamos a función handle para evaluar más fácil
 
A=Dft2(Rs,alpha); % La cota de f'' para cada rádio
 
n=1:8;
[NN,AA]=meshgrid(n,A);
 
RR=(ones(N,1)*Rs)';
 
Cota=abs(2*pi^2/12*(1-RR.^2).*AA.*10.^(-2*NN)); % Le añadimos los factores restantes y obtenemos la cota
 
colors=["r","b","g","m","k","c"]; % Cambiar si se ponen más puntos
for i=1:K
    semilogy(n,error(i,:),colors(i),"LineWidth",2) 
    hold on
    leyenda(i)="(1-10^{-"+num2str(i)+"}, \pi/4) y su Cota" ;
end
for i=1:K
    semilogy(n,Cota(i,:),"-."+colors(i),"LineWidth",2)
    hold on
end
title("Error del Método del Trapecio y Cota del Máximo")
legend(leyenda)
xlabel("Discretización de tamaño 10^{n}")
ylabel("Error")

3.3 Error numérico de la solución por serie de Fourier

Ya hemos visto que al obtener la solución mediante la fórmula de Poisson, obtenemos un error debido a la naturaleza numérica del método del trapecio. Ahora vamos a comprobar que, considerando de nuevo la solución [math] u(x,y)=xy [/math] con valor en la frontera [math] g(x,y)=xy [/math], si empleamos series de Fourier, no vamos a obtener error.

Con las cuentas pertinentes se obtiene que la familia de soluciones es [math]\left \{ 1/2, r^k cos(k \theta), r^k sin(k \theta) \right\}_{k \in \mathbb{N}} [/math]. Por el principio de superposición, al ser una base de [math]L^2[/math], sabemos [math] u(\theta,r)= \frac{\alpha_0 }{2} + \sum_{k=1}^{\infty} (\alpha_{k} r^k cos(k \theta) + \beta_k r^k sin(k \theta)) [/math] también es solución en el interior de la bola.

Sin embargo, para que sea solución en la frontera, se debe cumplir que [math] u(\theta,r)=G(\theta)[/math], siendo [math] G(\theta)=g(x,y)[/math] con [math] (x,y) \in \partial B_1(0) [/math]. Si consideramos la base trigonométrica usual, se tiene que [math] G(\theta)= \frac{a_0 }{2} + \sum_{k=1}^{\infty} (a_{k} cos(k \theta) + b_k sin(k \theta)) [/math]

De esta manera se tiene que los coeficientes de la serie de Fourier son:

En este caso [math] G(\theta)=g(x,y)=g(1,\theta) = cos(\theta) sin(\theta)= \frac{ sin(2 \theta)}{2} [/math]. Por tanto, se tiene que [math]a_0=0 [/math], [math]a_k=0 ~ \forall k ~ \in \mathbb{N} [/math] y [math]b_k=0 ~ \forall k ~\in \mathbb{N} [/math] excepto [math] b_2 = \frac{1}{2}[/math]. Luego, teniendo en cuenta que [math] R=1 [/math], la solución es [math]u(r,\theta)= \frac{1}{2} r^2 sin(2\theta)[/math].

Como podemos observar, la solución por series de Fourier no va a dar errores numéricos puesto que la solución es un múltiplo del segundo término de la serie. Por ello, en este caso, es más conveniente trabajar con este método.

3.4 Aplicación de la desigualdad de Harnack

En esta sección, vamos a aplicar la desigualdad de Harnack ^[1] sobre el problema sobre el que estamos trabajando. Al igual que en secciones anteriores vamos a considerar [math] g(x,y)=xy [/math].

Sea [math] M [/math] el mínimo de [math] g [/math] en la frontera [math] \partial B_1(0) [/math] y [math] v=u-M \geq 0[/math] armónica. Nuestro objetivo es encontrar la región en la que deben estar todas las soluciones armónicas en la [math]B_1[/math]. En primer lugar, vamos a resolverlo para el caso general, es decir, para [math]B_R[/math].

Para calcular [math] M [/math] partimos de [math] g(r, \theta )= \frac{r^2}{2} \sin{(2\theta)} [/math]. Dado que [math]r[/math] es siempre positivo y el seno puede tomar valores negativos, vemos que el mínimo de esta función se alcanza cuando [math]r[/math] alcanza su máximo (que en este caso es [math] R [/math]) y [math] \sin{(2\theta)} [/math] alcanza su mínimo, (que en este caso es [math] -1 [/math] cuando [math]\theta=\frac{3\pi}{4}, \frac{7\pi}{4} )[/math]. Por tanto, llegamos a que [math] M=-\frac{R^2}{2} [/math].

Aplicamos ahora la desigualdad de Harnack tomando [math]n=2[/math], [math]R=1[/math] y [math] r=|x-z|=|x| [/math] (ya que tomamos [math] z=(0,0) [/math]) para [math]v[/math]:

Dado que [math] v(r,\theta)=u(r,\theta)+\frac{R^2}{2}[/math], tenemos que la región en la que deben estar todas las soluciones armónicas en la [math]B_R[/math], que en [math](0,0)[/math] valen [math]0[/math], es la siguiente :

Ahora, si particularizamos para [math]B_1[/math], el resultado es el siguiente:

A continuación, vamos a visualizar cómo varía la región donde se encuentran las funciones armónicas a medida que aumenta el radio ([math] R [/math]). Para ello, vamos a representar las cotas que determinan esa región en el siguiente vídeo:

Representación de las cotas superior e inferior para dimensión 2 que determinan la región donde se encuentran las funciones armónicas a medida que aumenta el radio [math] R [/math] entre [math] 1 [/math] y [math] 20 [/math], considerando [math] r \in [0,1) [/math].

Como podemos observar, el comportamiento a medida que varía el radio es similiar, únicamente aumenta su tamaño. Lo más detacable es que la cota superior en la frontera tiende a infinito. Además, la cota inferior alcanza su valor mínimo también en la frontera.

Para poder apreciar mejor esto, vamos a aplicar el comando "semilogy". Sin embargo, como [math]u [/math] toma valores negativos, debemos considerar la función [math]v[/math], que por definición es siempre positiva.

En el siguiente vídeo se ilustra lo anterior:

Representación de las cotas superior e inferior para dimensión 2 considerando la función [math] v [/math] y aplicando el comando "semilogy".

En este vídeo se puede apreciar mejor cómo la cota inferior alcanza el mínimo en [math] R [/math]. Este resultado es justo lo que esperamos viendo la expresión analítica de la cota inferior, se observa que el numerador de la expresión se anula cuando [math]r=R[/math]. También se puede observar que la cota superior tiende a infinto en [math]R[/math].

Veamos ahora el resultado que se obtiene para dimensión [math]3[/math]. Dado que [math]g(x,y,z)=xy[/math] no depende de [math]z[/math], el mínimo sigue siendo [math] M=-\frac{R^2}{2} [/math]. Aplicando la desigualdad de Harnack tenemos que:

Por tanto, la región en la que deben estar las soluciones armónicas en la [math]B_R[/math] en dimensión [math]3[/math] es:

Al igual que para dimensión [math]2[/math], vamos a representar las cotas que determinan la región donde se encuentran las funciones armónicas a medida que aumenta el radio [math] R [/math] entre [math] 1 [/math] y [math] 20 [/math]:

Representación de las cotas superior e inferior para dimensión 3 que determinan la región donde se encuentran las funciones armónicas a medida que aumenta el radio [math] R [/math] entre [math] 1 [/math] y [math] 20 [/math], considerando [math] r \in [0,1) [/math].

Como podemos ver, el comportamiento de las cotas superior e inferior es muy similar al obtenido para dimensión 2. La principal diferencia es que, en este caso, la pendiente es más pronunciada. Al igual que antes, vamos a aplicar el comando "semilogy" para poder apreciar mejor el comportamiento de ambas cotas:

Representación de las cotas superior e inferior para dimensión 3 considerando la función [math] v [/math] y aplicando el comando "semilogy".

En esta gráfica se puede apreciar que la pendiente de la cota inferior es menos pronunciada que en el caso de dimensión 2.

3.4.1 Programa

Códigos para los vídeos en 2 dimensiones:

clear all
close all
 
R0=1; Rf=20;  % Definimos el intervalo de valores de R
divR=10^-1;    % División de la discretización
 
RR=R0:divR:Rf; % Discretización del intervalo de R
 
Hmin=@(r,R)R.^2.*(R-r)./(R+r)/2 - R.^2/2;  % Cota inferior
Hmax=@(r,R)R.^2.*(R+r)./(R-r)/2 - R.^2/2;  % Cota superior
 
pelicula=VideoWriter('Video_Ej4_1_P3'); % Creamos el vídeo
pelicula.FrameRate=10;
open(pelicula); % Lo abrimos
for j=1:length(RR)
    figura = figure(1);
    rr=0:divR:RR(j);
    HMIN=Hmin(rr,RR(j));
    HMAX=Hmax(rr,RR(j));
    plot(rr,HMIN,"b","LineWidth",1.5)
    hold on
    plot(rr,HMAX,"r","LineWidth",1.5)
    hold off
    xlim([0 RR(j)])
    ylim([min(HMIN) max(HMAX)])
    xlabel("radio")
    legend("cota inferior","cota superior")
    title("Desigualdad de Harnack")
    subtitle("R= "+num2str(RR(j)))
    imagen=getframe(figura);
    writeVideo(pelicula,imagen); % Inserto la imagen
 
end
close(pelicula);

Aplicando el comando "semilogy" a la función [math] v [/math]:

clear all
close all
 
R0=1; Rf=20;  % Definimos el intervalo de valores de R
divR=10^-1;    % División de la discretización
 
RR=R0:divR:Rf; % Discretización del intervalo de R
 
Hmin=@(r,R)R.^2.*(R-r)./(R+r)/2; %  Cota inferior
Hmax=@(r,R)R.^2.*(R+r)./(R-r)/2 ;%  Cota superior
 
pelicula=VideoWriter('Video_Ej4_log_P3'); % Creamos el vídeo
pelicula.FrameRate=10;
open(pelicula); % Lo abrimos
for j=1:length(RR)
    figura = figure(1);
    rr=0:divR:RR(j);
    HMIN=Hmin(rr,RR(j));
    HMAX=Hmax(rr,RR(j));
    semilogy(rr,HMIN,"b","LineWidth",1.5)
    hold on
    semilogy(rr,HMAX,"r","LineWidth",1.5)
    hold off
    xlim([0 RR(j)])
    ylim([min(HMIN) max(HMAX)])
    xlabel("radio")
    legend("cota inferior","cota superior")
    title("Desigualdad de Harnack")
    subtitle("R= "+num2str(RR(j)))
    imagen=getframe(figura);
    writeVideo(pelicula,imagen); % Inserto la imagen
 
end
close(pelicula);

Códigos para los vídeos en 3 dimensiones:

R0=1; Rf=20;  % Definimos el intervalo de valores de R
divR=10^-1;    % División de la discretización
 
RR=R0:divR:Rf; % Discretización del intervalo de R
 
Hmin=@(r,R)R.^3.*(R-r)./((R+r).^2)/2 - R.^2/2; 
Hmax=@(r,R)R.^3.*(R+r)./((R-r).^2)/2 - R.^2/2;
 
pelicula=VideoWriter('Video_Ej4_2_P3'); % Creamos el vídeo
pelicula.FrameRate=10;
open(pelicula); % Lo abrimos
for j=1:length(RR)
    figura = figure(1);
    rr=0:divR:RR(j);
    HMIN=Hmin(rr,RR(j));
    HMAX=Hmax(rr,RR(j));
    plot(rr,HMIN,"b","LineWidth",1.5)
    hold on
    plot(rr,HMAX,"r","LineWidth",1.5)
    hold off
    xlim([0 RR(j)])
    ylim([min(HMIN) max(HMAX)])
    xlabel("radio")
    legend("cota inferior","cota superior")
    title("Desigualdad de Harnack")
    subtitle("R= "+num2str(RR(j)))
    imagen=getframe(figura);
    writeVideo(pelicula,imagen); % Inserto la imagen
 
end
close(pelicula);

Aplicando el comando "semilogy" a la función [math] v [/math]:

R0=1; Rf=20;  % Definimos el intervalo de valores de R
divR=10^-1;    % División de la discretización
 
RR=R0:divR:Rf; % Discretización del intervalo de R
 
Hmin=@(r,R)R.^3.*(R-r)./((R+r).^2)/2; 
Hmax=@(r,R)R.^3.*(R+r)./((R-r).^2)/2;
 
pelicula=VideoWriter('Video_Ej4_2_log_P3'); % Creamos el vídeo
pelicula.FrameRate=10;
open(pelicula); % Lo abrimos
for j=1:length(RR)
    figura = figure(1);
    rr=0:divR:RR(j);
    HMIN=Hmin(rr,RR(j));
    HMAX=Hmax(rr,RR(j));
    semilogy(rr,HMIN,"b","LineWidth",1.5)
    hold on
    semilogy(rr,HMAX,"r","LineWidth",1.5)
    hold off
    xlim([0 RR(j)])
    ylim([min(HMIN) max(HMAX)])
    xlabel("radio")
    legend("cota inferior","cota superior")
    title("Desigualdad de Harnack")
    subtitle("R= "+num2str(RR(j)))
    imagen=getframe(figura);
    writeVideo(pelicula,imagen); % Inserto la imagen
 
end
close(pelicula);

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

Ahora vamos a estudiar la ecuación de Poisson en [math] \mathbb{R}^2 [/math]. Concretamente vamos a usar el potencial logarítmico para, a partir de la solución fundamental, aproximar la ecuación de Poisson.

4.1 Representación de una solución particular

En esta sección, vamos a aproximar la ecuación de Poisson:

Cuando [math] f [/math] es la función característica de la bola de radio [math] 1[/math].

Para ellos vamos a emplear la solución fundamental cuya expresión en [math] \mathbb{R}^2 [/math] es la siguiente:

Sabemos que la solución a la ecuación de Poisson viene dada por la convolución

Por tanto, la solución en este caso, vendrá dada por la siguiente expresión:

Haciendo el cambio de variable [math] \textbf{y}= ( r cos (\theta), r sin( \theta)) [/math] y [math] \textbf{x} = (l cos(\alpha), l sin (\alpha))[/math] se tiene que

Por tanto la solución es:

Sin embargo, esta solución no depende de [math] \alpha [/math], esta propiedad nos permite reducir la complejidad del programa con el cual representaremos la solución. Para verlo, vamos a derivar la solución anterior con respecto de [math] \alpha [/math] y observaremos que es cero.

Nótese que la última igualdad se tiene porque [math] cos(2 \pi- \alpha) = cos(- \alpha)[/math].

Luego, la solución no depende de [math] \alpha [/math], siendo entonces la solución:

Representamos esta solución:

Gráfica de la solución de la ecuación de Poisson con [math] f [/math] la función característica de la bola de radio [math] 1[/math] y [math] l \in [0,10] [/math] y [math] \alpha \in [0,2 \pi] [/math].

Como podemos observar, la solución es positiva únicamente en la bola unidad y alcanza el valor máximo en el punto [math] (0,0) [/math], ya que es decreciente en [math] l[/math].

4.1.1 Programa

Código para representar la gráfica de la solución:

clear all
close all
 
% Definimos los intervalos de integración y sus discretizaciones
 
th0=0; thf=2*pi; %Valor inicial y final de theta
r0=0; rf=1;      %Valor inicial y final del radio
 
divth=10000; % División del intervalo de theta
divr=10000;  % División del intervalo del radio
 
Th=linspace(th0,thf); %Discretización del intervalo de theta
Rr=linspace(r0,rf);   %Discretización del intervalo del radio
 
[TH,RR]=meshgrid(Th,Rr);
 
U=zeros(size(TH)); % Inicializamos los valores de u
 
f=@(a,l,r)l.*log(abs(l.^2+r.^2-2*r.*l.*cos(a))); % Función integrando
 
for i=1:length(Rr)
    U(i,:)=integral2(@(a,l)f(a,l,Rr(i)),0,2*pi,0,1)*ones(1,length(Th));  % Evaluamos la integral 
end
 
% Cambio de variable
 
XX=RR.*cos(TH);
YY=RR.*sin(TH);
 
% Dibujamos
 
surf(XX,YY,-U/(4*pi))
shading interp
title("Solución de la ecuación de Poisson para r \in [0,"+num2str(rf)+"]")

4.2 Comportamiento de la solución en el infinito

A continuación, vamos a estudiar su comportamiento en el infinito.

Comenzamos acotando el integrando de la solución inferior y superiormente. Conociendo que [math] r \in [0,1] [/math] y [math] \theta \in [0,2 \pi] [/math] se tiene que:

Observamos que ambas acotaciones, cuando [math] l \rightarrow \infty [/math], son del orden de [math] 2 r log(l)[/math]. Por tanto, para valores grandes de [math] l [/math], se tiene que la solución es

Por tanto, la solución se comporta como el logaritmo del radio en el infinito. Para visualizar esto , como la solución no depende del ángulo [math] \alpha [/math], la hemos representado como gráfica en función de [math]l[/math] y hemos aplicado el comando "semilogx" para eliminar la dependencia logarítmica en [math]l[/math] . Teniendo esto en cuenta, se debería visualizar una recta con pendiente [math]-\frac{1}{2} [/math] cuando consideramos un valor elevado de [math]l[/math].

Gráfica de comparación entre el comportamiento en el infinito de la solución y el logaritmo, representada en base logarítmica.

Como podemos observar, efectivamente, las gráficas de la solución y el logaritmo coinciden, siendo la recta con pendiente [math]-\frac{1}{2} [/math]. Por tanto, hemos obtenido los resultados esperados según el estudio analítico previo.

4.2.1 Programa

Código para representar la gráfica de comparación entre la solución y el logaritmo en base logarítmica:

clear all
close all
 
% Definimos los intervalos de integración y sus discretizaciones
 
r0=10^10; rf=10^15;      %Valor inicial y final del radio
 
divr=10^4;  % División del intervalo del radio
 
Rr=linspace(r0,rf,divr);   %Discretización del intervalo del radio
 
U=zeros(size(Rr)); % Inicializamos los valores de u
 
f=@(a,l,r)l.*log(l.^2+r.^2-2*r.*l.*cos(a));
 
for i=1:length(Rr)
    U(i)=integral2(@(a,l)f(a,l,Rr(i)),0,2*pi,0,1);  % Evaluamos la integral 
end
 
% Dibujamos
 
semilogx(Rr,-U/(4*pi),"r","LineWidth",3)
hold on 
semilogx(Rr,-log(Rr)/2,"b","LineWidth",1)
title("Comportamiento de la función en el infinito")
xlabel("l")
legend("u(l)","-log(l)/2")

5 Referencias

• Error de la regla del trapecio
• Series de Fourier