Ecuación de Laplace. Otelo, Yan y Mika

1 Introducción

En este articulo estudiaremos la ecuación de Laplace en la bola unidad de dimensión dos con diferentes condiciones frontera. En concreto, analizaremos la fórmula de Poisson y los errores cometidos al aproximar la integral de esta mediante la regla del trapecio con varias discretizaciones del dominio. Además, estudiaremos el error cometido al fijar una discretización e ir tomando puntos que se acercan a la frontera. También estudiaremos la solución por serie de Fourier y la desigualdad de Harnack. Por último, aproximaremos la solución de la Ecuación de Poisson con el potencial logarítmico de dos dimensiones.

2 Solución de la Ecuación de Laplace con la fórmula de Poisson y la regla del trapecio

Consideramos el problema

[math] \left\{\Delta u =0, \hspace{5mm} x \in B_1\atop \hspace{5mm} u = g \hspace{5mm} x\in\partial B_1 \right. [/math]

donde g es la función de una variable [math] g(\theta) = \max\{0, 1-\frac{2}{\pi}|\theta - \pi|\} [/math]

Como la función está en polares, para hallar la solución del problema vamos a usar la fórmula de Poisson en polares:

[math] u(r, \theta) = \frac{R^2 - r^2}{2\pi} \int_{0}^{2\pi}\frac{G(\varphi)}{R^2 + r^2 - 2Rr\cos(\theta - \varphi)}d\varphi[/math]

Para resolver numéricamente la integral de la fórmula en MatLab usaremos la regla del trapecio. Esta regla consiste en dividir [math] [0,2\pi] [/math] en N subintervalos para después realizar la media de los valores del integrando en los dos extremos de cada subintervalo y multiplicarla por la longtitud de estos.

Si aplicamos la fórmula sobre la adherencia de la bola obtenemos discontinuidad en la frontera. Esto se muestra en la gráfica a continuación:

Solución discontinua

Para solucionar esto, imponemos directamente la condición frontera y aplicamos la fórmula en el interior de la bola. De esta manera, eliminamos la discontinuidad y conseguimos la siguiente gráfica:

Solución tras imponer la condición frontera

2.1 Código

%Resumen de código: 
%Crearemos un malla de puntos en coordenadas polares y en cada punto
%evaluaremos usando la fórmula de Poisson, realizando la integral con la
%regla del trapecio. Sin embargo, en la frontera (r=1) evaluaremos
%directamente la función g definida a continuación

%condicion frontera
g =@(x) max(0,1-(2/pi).*abs(x-pi));

%dimensión de la malla
Q = 100;

%coordenadas polares para construir la malla
R = 1;
D = linspace(0,R,Q);%radios
Z = linspace(0, 2*pi,Q);%ángulos
evaluaciones = zeros(Q,Q);%matriz malla

%para la regla del trapecio
N=1000;                         %número de puntos
a=0; b=2*pi;                    %extremos del intervalo
h=(b-a)/N;                      %incremento
u=a:h:b; 
w=ones(N+1,1);                  %vector de pesos
w(1)=1/2; w(N+1)=1/2;

%imponemos la condición frontera (r = 1) evaluando en g
for j = 1:Q
    evaluaciones(Q,j) = g(Z(j));
end

%fórmula de Poisson en el interior
for i =1:Q-1 %hasta Q-1 para no aplicarla en la frontera
    r = D(i);
    for j = 1:Q
        z = Z(j).*ones(1,N+1);  %tenemos que convertirlo en vector para poder operar con u
        f = (g(u)./(R^2 + r^2 - (2*R*r).*cos(z - u)))';  %funcion que integramos
        int = h*w'*f;      %resultado de la integral                                  
        evaluaciones(i,j) = ((R^2 - r^2)/(2*pi))*int;   %formula de poisson en coordenadas polares
    end
end

%pasamos a cartesianas para hacer la gráfica
X = D'*cos(Z);
Y = D'*sin(Z);

%representación 3D
surf(X,Y,evaluaciones)
title('Imponiendo la condición frontera')
xlabel('X')
ylabel('Y')
zlabel('u')

3 Limitaciones de la fórmula de Poisson relacionadas con la regla del trapecio

Como hemos comentado en la introducción, las limitaciones de la fórmula de Poisson provienen principalmente de la aproximación de la integral. En concreto, la fórmula de trapecio que utilizamos incorpora un error en la aproximación y en este apartado vamos a analizar cómo varía en función de la discretización, es decir, en función del número de subintervalos que tomemos para dividir [math] [0,2\pi][/math]. Para ello, vamos a considerar el mismo problema pero esta vez con la condición frontera

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

Además, vamos a considerar como solución exacta la misma función [math] u(x,y) = xy [/math]

Solución exacta del problema

Con ella calcularemos el error que obtenemos al usar la fórmula de trapecio en un punto lejos de la frontera, en concreto, en [math] (r, \theta) = (0.9, \pi/4) [/math]. Para poder apreciarlo mejor, lo calculamos en escala logarítmica, es decir, si el número de subintervalos es [math] N = 10^n [/math] y denotamos [math] u_N [/math] la solución obtenida para esta discretización, el error será

[math] f(n) = \log_{10} error(10^n) = \log_{10} |u - u_N|[/math]

Para asegurarnos de que los errores cometidos en la aproximación no son muy grandes y todo va como se espera, vamos a usar una estimación del error de la fórmula del trapecio de Wikipedia. Esta nos dice que, si utilizamos la regla del trapecio para integrar una función [math] f[/math] en el intervalo [math] [a,b][/math] con una división en [math] N[/math] subintervalos, el error estimado para un valor [math] c \in [a,b][/math] es

[math] error(c) = \left|-\frac{(b-a)^3}{12N^2}f''(c)\right|[/math]

Por tanto, para nuestro caso tomaremos el valor absoluto y tendremos

[math] error(c) = \frac{(2\pi)^3}{12N^2}\left|\frac{\partial^2}{\partial \varphi^2}\left(\frac{G(\varphi)}{R^2 + r^2 - 2Rr\cos(\theta^* - \varphi)}\right)_{\varphi = c}\right|[/math]

donde [math] \theta^*[/math] es una constante.

Hemos hecho una acotación del valor absoluto de la doble derivada

[math] \left|\frac{\partial^2}{\partial \varphi^2}\left(\frac{G(\varphi)}{R^2 + r^2 - 2Rr\cos(\theta^* - \varphi)}\right)_{\varphi = c}\right| \leq \frac{116}{(1-r)^6} [/math]

Obteniendo así

[math] error(c) \leq \frac{116(2\pi)^3}{12N^2(1-r)^6} = \frac{232\pi^310^6}{3N^2}[/math]

A continuación se muestra la representación gráfica del error calculado y el error estimado en el punto [math] (r, \theta) = (0.9, \pi/4) [/math].

Error calculado y cota del error estimado seún la discretizaión de la fórmula del trapecio

3.1 Código

%Resumen de código: 
%Para este apartado evaluamos soloamente en el punto (r, theta) = (0.9, pi/4)
%de nuevo con la fórmula de Poisson y esta vz comparamos con la solución
%exacta xy también evaluda en el punto

%condición frontera en polares
g =@(x) sin(x).*cos(x); 
R = 1;
imagen_u = 0.9^2*sin(pi/4)*cos(pi/4);
max_exp =8;
errores = zeros(1,max_exp); cota_error_teorico = zeros(1,max_exp);
for n=1:max_exp
    %para la regla del trapecio
    N=10^n;                         %número de puntos
    a=0; b=2*pi;                    %extremos del intervalo
    h=(b-a)/N;                      %incremento
    u=a:h:b; 
    w=ones(N+1,1);                  %vector de pesos
    w(1)=1/2; w(N+1)=1/2;
   
    %fórmula de Poisson en el punto
    r = 0.9;
    z = (pi/4).*ones(1,N+1);
    f = (g(u)./(R^2 + r^2 - (2*R*r).*cos(z - u)))';%funcion que integramos 
    evaluacion = ((R^2 - r^2)/(2*pi))*h*w'*f ;%formula de poisson para coordenadas polares
    errores(1,n) = log10(abs(evaluacion - imagen_u));

    %acotación de la estimación teórica del error
    cota_error_teorico(1,n) = log10(232*(pi)^3*10^6/(3*N^2));
end

figure(1)
hold on
plot(1:max_exp, errores, 'b')
plot(1:max_exp, cota_error_teorico, 'r')
hold off
title('Errores según el exponente del número de puntos de la discretización')
legend("error calculado", "cota error teórico")
xlim([1,max_exp])
ylim([-20,10])
xlabel('exponente del número de puntos de la discretización (10^n)')
ylabel('errores')

Ahora, vamos a calcular el error en distintos puntos. Para ello, fijamos [math] N =100 [/math] para la fórmula del trapecio y hacemos que los puntos en los que aplicamos la fórmula de Poisson se vayan acercando a la frontera definiéndolos en coordenadas polares como [math] (r, \theta) = (1 - 10^{-n}, \pi/4) [/math]. En consecuancia, la cota estimada del error pasa tener esta forma:

[math] error(c) \leq \frac{116(2\pi)^3}{12N^2(1-r)^6} = \frac{232\pi^310^{6n - 4}}{3}[/math]

Por tanto, a medida que [math] r [/math] se acerca a 1, la cota de la estimación teórica tiende a infinito y la gráfica anteriormente realizada queda así:

Error calculado y cota del error estimado cuando r tiende a 1

3.2 Código

%Resumen de código
% Fijamos ahora el número de puntos en la fórmula del trapecio en 100 y
% dibujamos la gráfica del error cuando nos acercamos a la frontera, es decir, para puntos
% de la forma (r, θ) = (1 − 10^−n, π/4). 

%condicion frontera (R =1) en polares
g =@(x) sin(x).*cos(x);
R = 1;

%para la regla del trapecio
N=100;                         %número de puntos
a=0; b=2*pi;                    %extremos del intervalo
h=(b-a)/N;                      %incremento
u=a:h:b; 
w=ones(N+1,1);                  %vector de pesos
w(1)=1/2; w(N+1)=1/2;

%comparación de las evaluaciones
max_exp =10;
errores = zeros(1,max_exp);cota_error_teorico2 = zeros(1,max_exp);
for n=1:max_exp
    %fórmula de Poisson para cada r en el punto
    r = 1 - 10^(-n);
    z = (pi/4).*ones(1,N+1);  %de nuevo tien que ser vector
    imagen_u = r^2*g(pi/4);
    f = (g(u)./(R^2 + r^2 - (2*R*r).*cos(z - u)))';%funcion que integramos
    evaluacion = ((R^2 - r^2)/(2*pi))*h*w'*f  ;%formula de poisson para coordenadas polares
    errores(1,n) = log(abs(evaluacion - imagen_u));
    cota_error_teorico2(1,n) = log10(232*(pi)^3*10^(6*n-4)/3);
end
%la cota del error teorico no nos sirve cuando r se aproxima a
%1 poeque el denominador de la cota f'' tiende a 0
figure(1)
hold on
plot(1:max_exp, errores, 'b')
plot(1:max_exp, cota_error_teorico2, 'r')
hold off
title('Error y cota estimada del error a medida que r se acerca 1')
legend("error calculado", "cota error teórico")
xlim([1,max_exp])
xlabel('n que hace que r tienda a 1')
ylabel('errores')

5 Desigualdad de Harnack

En esta sección vamos a analizar lo que nos dice las desigualdades de Harnack. Esta es la siguiente: Sea [math]u(x)[/math] una función armónica y no negativa, [math]u(x)\geq 0[/math] en un dominio [math]\Omega \in \mathbb{R}^n[/math] y [math]\mathbb{B}_R(z) \in \Omega[/math]. Entonces, para todo [math]x \in \overline{\mathbb{B}_R(z)}[/math],

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

donde [math] r = |x-z| [/math].

Como la función [math]u(x,y)=xy[/math] es claramente armónica falta comprobar que es no negativa. Para asegurar que es no negativa, vamos a considerar la función armónica y no negativa [math]v(x,y) = u(x,y) - M [/math], donde M es el mínimo de [math]g(x,y) = xy[/math] en la frontera [math]\partial \mathbb{B}_1(0)[/math]. Se calcula el mínimo M de la función [math]g[/math] para cualquier radio R y se obtiene lo siguiente:

[math]M = -\frac{R^2}{2} [/math]

A continuación, aplicamos las desigualdades de Harnack a la función [math]v[/math] en la bola [math]\mathbb{b}_R(0)[/math] para acotarla de la siguiente manera:

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

Deshaciendo el cambio para encontrar las cotas de las soluciones armónicas y desarrollando los cálculos, se obtiene:

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

Ahora, sustituyendo el valor de la M, la región donde deben estar las soluciones armónicas en la bola [math]\mathbb{B}_R(0)[/math] es:

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

Podemos observar que las cotas solo depende del radio [math]r \in [0,R)[/math].Además, se tiene que la cota inferior dada por la desigualdad de Harnack es siempre negativa.

Vamos a dibujar estas regiones con el siguiente código de matlab para las bolas [math] \mathbb{B}_1(0) [/math],[math] \mathbb{B}_2(0)[/math] y [math] \mathbb{B}_{10}(0) [/math]:

%En este código vamos a representar la región donde deben estar las soluciones a partir de las desigualdades de Harnack para las bolas de radio 1, 2 y 10. 

g =@(x,y) y^2.*sin(x).*cos(x);   % Condición frontera 
R = [1,2,10];    % Radios de las bolas

for k=1:3
    M = -R(k)^2/2;  % Mínimo de la función de la condición frontera
    rmax = R(k)-0.01; num = 100; % rmax es el valor máximo del radio cercano a la frontera y num el número de puntos de la discretización
    D = linspace(0, rmax, num); % Vector de la discretización
    valores_inferiores = zeros(1,num); valores_superiores = zeros(1,num);   % Vector donde se almacena los valores de las cotas inferiores y superiores
    for i=1:num
        valores_inferiores(i) = -(R(k)^2*D(i))/(R(k)+D(i));  % Cota inferior dada por la desigualdad de Harnack
        valores_superiores(i) = (R(k)^2*D(i))/(R(k)-D(i));  % Cota superior dada por la desigualdad de Harnack
    end

    % Desplazamos los valores de obtenidos para poder aplicar la escala logarítmica
    m = min(valores_inferiores);
    valores_inferiores = log10(valores_inferiores - (m-1).*ones(1,num));
    valores_superiores = log10(valores_superiores - (m-1).*ones(1,num));
    
    % Representación gráfica
    figure(k)
    hold on
    plot(D,valores_inferiores,'r')
    plot(D, valores_superiores, 'g')
    xlabel('Radios(r)')
    ylabel('Cota')
    legend('Cota inferior','Cota superior','Location','northwest')
    title("Región para la bola de Radio " + num2str(R(k)))
    hold off
end

Región para la bola de radio 1

Región para la bola de radio 2

Región para la bola de radio 10

Nótese que como nuestro objetivo es representarlo en escala logarítmica, ya que se aprecia mejor, y la cota inferior siempre toma valores negativo; desplazamos los valores para hacerlos todos positivos y poder reescalar de forma logarítmica. A los valores calculados por la desigualdad de Harnack les restamos el mínimo de ellos, de esta forma son todos positivos.

La región que queda entre las cotas superiores e inferiores proporcionadas por las desigualdades de Harnack es donde deben estar las soluciones armónicas. Se observa que a medida que aumenta el radio de la bola, la cota inferior se va ampliando.

A continuación, vamos a repetir lo anterior pero para dimensión 3. Tomando las mismas funciones [math]u(x,y,z) = xy[/math] y [math]g(x,y,z)=xy[/math] anteriores pero en dimensión 3, se obtiene el mismo mínimo anterior:

[math] M = -\frac{R^2}{2} [/math]

Con este mínimo y en dimensión 3, las desigualdades de Harnack son:

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

Sustituyendo el valor del mínimo M, nos queda lo siguiente:

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

Podemos ver que las cotas solo dependen del radio [math] r \in [0,R)[/math] y que la cota inferior toma valores negativos.

Modificamos el código anterior para representar la región dada por las desigualdades de Harnack para las bolas de radio 1, 2 y 10 pero esta vez en dimensión 3.

%En este código vamos a representar la región donde deben estar las soluciones a partir de las desigualdades de Harnack para las bolas de radio 1, 2 y 10. 

g =@(x,y) y^2.*sin(x).*cos(x);   % Condición frontera 
R = [1,2,10];    % Radios de las bolas

for k=1:3
    M = -R(k)^2/2;  % Mínimo de la función de la condición frontera
    rmax = R(k)-0.01; num = 100; % rmax es el valor máximo del radio cercano a la frontera y num el número de puntos de la discretización
    D = linspace(0, rmax, num); % Vector de la discretización
    matriz_inferiores = zeros(1,num); matriz_superiores = zeros(1,num);   % Matrices donde se almacena los valores de las cotas inferiores y seperiores
    for i=1:num
        matriz_inferiores(i) = (-(R(k)^2*D(i))*(D(i)+3*R(k)))/(2*(R(k)+D(i))^2);  % Cota inferior dada por la desigualdad de Harnack
        matriz_superiores(i) = ((R(k)^2*D(i))*(3*R(k)-D(i)))/(2*(R(k)-D(i))^2);  % Cota superior dada por la desigualdad de Harnack
    end

    % Desplazamos los valores de obtenidos para poder aplicar la escala logarítmica
    m = min(matriz_inferiores);
    matriz_inferiores = log10(matriz_inferiores - (m-1).*ones(1,num));
    matriz_superiores = log10(matriz_superiores - (m-1).*ones(1,num));
    
    % Representación gráfica
    figure(k)
    hold on
    plot(D,matriz_inferiores,'r')
    plot(D, matriz_superiores, 'g')
    xlabel('Radios(r)')
    ylabel('Cota')
    legend('Cota inferior','Cota superior','Location','northwest')
    title("Región para la bola de Radio " + num2str(R(k)))
    hold off
end

Región para la bola de radio 1 en dimensión 3

Región para la bola de radio 2 en dimensión 3

Región para la bola de radio 10 en dimensión 3

Volvemos a desplazar los valores para poder graficar en escala logarítmica ya que los valores de la cota inferior son negativos. En las gráficas podemos ver que tienen la misma forma que en dos dimensiones. Sin embargo, esto es engañoso ya que la solución depende de tres variables, el radio, el ángulo azimutal y el ángulo longitudinal. Entonces en dos dimensiones, la región sería rotar las curvas de las cotas alrededor del eje z; pero en 3 dimensiones, habría que rotarla también alrededor del eje x hasta 180º.

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

En este apartado daremos solución a un caso concreto de la Ecuación de Poisson. Previamente, mostramos algunas nociones teóricas que serán utilizadas en la resolución del problema.

LLamamos solución fundamental del Laplaciano en [math] \mathbb{R}^2 [/math] a la función que viene dada por la expresión [math] \phi (x)=\frac{-1}{2\pi}log(|x|)[/math], donde [math] log(|x|) [/math] recibe el nombre de potencial logarítmico. En clase hemos visto que, utilizando esta expresión, podemos encontrar una única solución para el problema

[math] \left\{ -\Delta u=f \hspace{4mm} x\in\mathbb{R}^2 \atop \lim_{|x|\to\infty} u(x)=\frac{-M}{2\pi}log(|x|) + o (\frac{1}{(|x|)}), \hspace{6mm} \text{ donde} \hspace{3mm} M= \int_{\mathbb{R}^2} f(y) dy \right. [/math]

que será de la forma

[math] u(x)=\int_{\mathbb{R}^2} \phi(x-y)f(y) dy. [/math]

6.1 El problema

En este apartado, intentaremos encontrar una solución para la siguiente ecuación:

[math] \{\hspace{2mm} \Delta u=\mathbb{1}_{\mathbb{B}(0,1)}, \hspace{4mm} x\in\mathbb{R}^2. [/math]

Esta ecuación es prácticamente idéntica a la descrita arriba, salvo por un signo negativo. Pero, al ser una ecuación lineal, si cambiamos el signo de la solución habremos obtenido la solución con f cambiada de signo, y puesto que nuestra función f es la característica queda igual. La función solución resultante es:

[math] u(x)=\int_{\mathbb{R}^2} \phi(x-y)f(y) dy = \int_{\mathbb{R}^2} \frac{-1}{2\pi}log(|x-y|) \cdot \mathbb{1}_{\mathbb{B}(0,1)} (y) dy = \int_{\mathbb{B}(0,1)} \frac{-1}{2\pi}log(|x-y|)dy. [/math]

Puesto que esta es una integral irresoluble analíticamente, hallaremos una solución numérica. Para ello, primero haremos un cambio de variable para [math] x [/math] y pasaremos a coordenadas polares la integral, que después de simplificar resulta:

[math] \left\{ x=(r_xcos(\theta_x),r_xsen(\theta_x)) \atop y=(r_ycos(\theta_y),r_ysen(\theta_y)) \right. \longrightarrow \int_{\mathbb{B}(0,1)} \frac{-1}{2\pi}log(|x-y|)dy= \frac{-1}{4\pi}\int_{0}^{2\pi} \int_{0}^1 r_ylog(r_x^2+r_y^2+cos(\theta_y-\theta_x) dr_y d\theta_y. [/math]

Obteniendo así la función solución en polares

[math] U(r_x,\theta_x)= \frac{-1}{4\pi}\int_{0}^{2\pi} \int_{0}^1 r_ylog(r_x^2+r_y^2+cos(\theta_y-\theta_x) dr_y d\theta_y, [/math]

cuya representación gráfica hecha por Matlab es

Representación gráfica de la solución.

%Creamos los vectores donde evaluaremos, tanto para x como para y.

Nx=100; Ny=100;
rxi=0; rxf=10; Oxi=0; Oxf=2*pi; ryi=0; ryf=0.9999; Oyi=0; Oyf=2*pi;
rx=linspace(rxi,rxf,Nx); Ox=linspace(Oxi,Oxf,Nx); ry=linspace(ryi,ryf,Ny); Oy=linspace(Oyi,Oyf,Ny);
[RRy,OOy]=meshgrid(ry,Oy);

soluciones=zeros(length(rx),length(Oy));

%fijámos las x, y para cada valor resolvemos la integral bidimensional
%resultante
for i=1:length(rx)
    for j=1:length(Ox)
        soluciones(i,j)=(-1/(4*pi))*integral2( @(rry,OOy) rry.*log(rx(i).^2+rry.^2-2.*rx(i).*rry.*cos(OOy-Ox(j))) ,ryi,ryf,Oyi,Oyf); 
    end
end
%Pasamos a cartesianas y ploteamos.
x1=rx'*cos(Ox);
x2=ry'*sin(Ox);
surf(x1,x2,soluciones)

En esta imagen, podemos ver como se cumple el principio del máximo. Si tomamos bolas que no intersequen a [math] B_1(0) [/math] donde el laplaciano de nuestra función (por la función característica) es cero, la función solución será armónica. En esta imagen se puede apreciar como los máximos y mínimos se toman siempre en las fronteras de dichas bolas, puesto que la función no es constante. Sin embargo, si cogemos cualquier bola que contenga a [math] B_1(0) [/math], el laplaciano no será nulo, y de hecho se puede ver como el máximo se tomar en el interior de dichas bolas, siendo esto posible por que en esos casos no se cumple el principio del máximo al no ser la solución armónica.

Por último, podemos comprobar que la solución obtenida tiene el comportamiento esperado en el infinito. Teniendo en cuenta que

[math] M= \int_{\mathbb{R}^2} \mathbb{1}_{\mathbb{B}(0,1)} (y) dy = \pi [/math]

se tiene que cuando [math] |x|\to\infty [/math] nuestra función debería tener valores prácticamente iguales a [math] \frac{1}{2}log(|x|) [/math]. Tomando valores de [math] |x| [/math] muy altos y utilizando Matlab obtenemos los siguientes resultados

Comparación de la solución con el comportamiento esperado.

clear all

%Igual que el programa anterior.
%Variando rxf podemos ver la solución en valores mayores de x
Nx=100; Ny=100;
rxi=0; rxf=10^20; Oxi=0; Oxf=2*pi; ryi=0; ryf=0.9999; Oyi=0; Oyf=2*pi;
rx=linspace(rxi,rxf,Nx); Ox=linspace(Oxi,Oxf,Nx); ry=linspace(ryi,ryf,Ny); Oy=linspace(Oyi,Oyf,Ny);
[RRy,OOy]=meshgrid(ry,Oy);

soluciones=zeros(length(rx),length(Oy));


for i=1:length(rx)
    for j=1:length(Ox)
        soluciones(i,j)=(-1/(4*pi))*integral2( @(rry,OOy) rry.*log(rx(i).^2+rry.^2-2.*rx(i).*rry.*cos(OOy-Ox(j))) ,ryi,ryf,Oyi,Oyf); 
    end
end

%Definimos un vector con las imagenes de la función esperada en el infinito
soluciones1=zeros(1,length(rx));
for i=1:length(rx)
    soluciones1(i)=(-1/2)*log(rx(i));
end

%mostramos el último valor.
soluciones(100,100)'
soluciones1(100)

verificando así que el comportamiento es el esperado.