Ecuación de Laplace y de Poisson. Grupo Eau De Parfum (EDP))

Trabajo realizado por estudiantes
Título	Ecuación de Laplace y de Poisson. Grupo Eau De Parfum (EDP)
Asignatura	EDP
Curso	2023-24
Autores	Lestau Torres, Pablo López Rojo, Celia Muñoz Guijarro, Sofía
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

1 Introducción

En este trabajo, exploraremos las soluciones a las ecuaciones de Laplace y Poisson, aplicando la fórmula de Poisson y técnicas numéricas para abordar problemas con condiciones de frontera de Dirichlet. Investigaremos sus limitaciones, aplicaremos la desigualdad de Harnack y exploraremos el comportamiento de estas ecuaciones en dominios no acotados como [math]\mathbb{R}^2[/math]

2 Contexto histórico

Pierre-Simon Laplace

Los estudios de fenómenos eléctricos y magnéticos tomaron forma a finales del siglo XVIII, estableciendo la electrostática y la magnetostática como disciplinas matemáticas al estilo newtoniano. Dentro de esta rama, Poisson destacó como uno de sus pioneros, ya que entre otras contribuciones amplió y profundizó los estudios sobre el potencial gravitatorio realizados previamente por Euler, Lagrange y Laplace. En 1785, Laplace había propuesto que la variación del potencial en cualquier punto, ya sea dentro o fuera de un cuerpo que ejerce una atracción gravitatoria, se rige por la ecuación que hoy lleva su nombre.

[math] \Delta u={\partial^2 u\over \partial x^2 } + {\partial^2 u\over \partial y^2 } + {\partial^2 u\over \partial z^2 } = 0. [/math]

siendo u una función real. Sin embargo, Poisson identificó en 1812 que esta ecuación no era aplicable a puntos (x,y,z) dentro del cuerpo atrayente. Así, reformuló la ecuación, conocida ahora como la ecuación de Poisson

[math] \Delta u={\partial^2 u\over \partial x^2 } + {\partial^2 u\over \partial y^2 } + {\partial^2 u\over \partial z^2 } = f(x,y,z). [/math]

3 Conocimientos previos

Siméon-Denis Poisson

Antes de comenzar será de gran ayuda introducir una serie de conceptos que utilizaremos a lo largo de este trabajo:

Función armónica:Una función [math]u \in C^{2}(\Omega)[/math]es armónica si [math]\Delta u =0[/math], es decir, si [math]u [/math] es solución de la ecuación de Laplace.

Desigualdad de Harnack:Sea u armónica y [math]u \geq 0[/math] en [math]\Omega \subset \mathbb{R}^n [/math]. Supongamos [math]\overline{B_{R}}(z) \subset \Omega[/math]. Entonces [math]\forall x \in \overline{B_{R}}(z) \text{ con } |z-x|=r \lt R[/math]

[math]{R^{n-2} (R-r)\over (R+r)^{n-1}} u(z)\le u(x)\le {R^{n-2}(R+r)\over (R-r)^{n-1}}u(z).[/math]

Solución fundamental:La solución fundamental de la ecuación de Poisson es una función [math]\Phi[/math] tal que [math]\Delta \Phi =\delta[/math], donde [math]\delta[/math] representa la función de Dirichlet. Esta solución es fundamental en el sentido de que sirve como generadora de soluciones más generales para la ecuación de Poisson bajo diversas condiciones de frontera.

Fórmula de Poisson :La única solución [math] u \in C^2(B_R(x_0))\cap C(\overline {B_R})[/math] del problema de Laplace [math] \left \{ \begin{array}{ll} \Delta u=0, x\in B_R\\ u=g, x\in \partial B_R \\ \end{array} \right. [/math], se conoce como fórmula de Poisson, siendo su expresión :

[math] u(x) = \frac{R^2 - |x - x_0|^2}{\omega_n R} \int_{\partial B_R(x_0)} \frac{g(\sigma)}{| (x - x_0)-\sigma|^n} d\sigma [/math]

donde [math] \omega_n [/math] representa la medida de la superficie de la esfera.

Cabe destacar que vamos a utilizar Series de Furier

4 Ecuación de Laplace

Para iniciar el análisis de la ecuación de Laplace, primero es necesario plantear el problema que se desea resolver.

Consideremos [math]B_1 \subset \mathbb{R}^2[/math] como la bola unidad centrada en el origen [math](0, 0)[/math]. Analizaremos el siguiente problema:

[math]\left \{ \begin{array}{ll} \Delta u = 0, & x \in B_1 \\ u = g, & x \in \partial B_1 \end{array} \right. [/math]

Nuestro objetivo es visualizar la solución de la ecuación de Laplace descrita arriba, empleando la fórmula de Poisson y las series de Fourier. Además, experimentaremos con cambios de la función frontera [math]g[/math] para observar cómo afectan la solución.

4.1 Solución de la ecuación de Laplace usando la fórmula de Poisson

Tomaremos nuestro sistema en coordenadas polares [math](r,\theta)[/math], siendo la función frontera [math] g(\theta) = \max\left\{0, 1 - \frac{2}{\pi} |\theta - \pi|\right\}[/math].

Nos centraremos en resolver el problema descrito anteriormente utilizando la fórmula de Poisson en la bola de radio 1 centrada en el origen, que en coordenadas polares tiene la siguiente expresión:

[math] u(r,\theta)=\frac{1-r^2}{2\pi}\int_{\partial B_1} \frac{g(\phi)}{1+r^2-2rcos(\theta-\phi)} d\phi [/math]

Para aproximar la solución, emplearemos la fórmula del trapecio en la integración numérica, debido a su eficacia en la aproximación de integrales definidas.

import numpy as np
import matplotlib.pyplot as plt
 
# Definición de la función g(theta)
def g(theta):
    return np.maximum(0, 1 - 2/np.pi * np.abs(theta - np.pi))
 
# Función para calcular la integral de Laplace-Poisson
def laplace_poisson(r, theta):
    # Valores de phi para integración
    phi_values = np.linspace(0, 2*np.pi, 1000)
    # Definición del integrando
    integrand = g(phi_values) * (1 - r**2) / (1 + r**2 - 2*r*np.cos(theta - phi_values))
    # Calcular la integral usando el método trapezoidal
    return np.trapz(integrand, phi_values) / (2 * np.pi)
 
# Función para calcular la solución de Laplace
def laplace_solution(R, Theta):
    # Crear un array vacío para almacenar los valores de la solución de Laplace
    laplace_values = np.zeros_like(R)
   
    # Iterar sobre cada punto en la malla
    for i in range(len(R)):
        for j in range(len(R[0])):
            r = R[i, j]
            theta = Theta[i, j]
            laplace_values[i, j] = laplace_poisson(r, theta)
    return laplace_values
 
# Generar malla de valores de r y theta
num_points = 100
r_values = np.linspace(0, 1, num_points)
theta_values = np.linspace(0, 2*np.pi, num_points)
R, Theta = np.meshgrid(r_values, theta_values)
 
# Calcular la solución de Laplace para cada punto en la malla
laplace_values = laplace_solution(R, Theta)
 
# Graficar
plt.figure(figsize=(8, 6))
plt.contourf(R * np.cos(Theta), R * np.sin(Theta), laplace_values, levels=50, cmap='viridis')
plt.colorbar(label='Solución de Laplace')
plt.title('Solución de Laplace en el disco unidad sin imponer g en la frontera')
plt.xlabel('x')
plt.ylabel('y')
plt.axis('equal')
plt.grid(True)
plt.show()

Observamos que en la frontera la solución vale 0, cuando por continuidad esto no debería ser así. La razón principal por la que en la frontera no se puede utilizar directamente la fórmula de Poisson es porque la fórmula se deriva asumiendo que el punto en el que estamos calculando la solución está en el interior del dominio. Cuando estamos en la frontera, esta suposición ya no es válida, y necesitamos utilizar directamente la condición de contorno.

Para abordar adecuadamente este problema, es imprescindible aplicar directamente la condición de contorno en la frontera. Esto implica que, en lugar de utilizar la fórmula de Poisson para esos puntos, debemos asignar los valores de la función de contorno [math]g[/math] especificados para la frontera.

import numpy as np
import matplotlib.pyplot as plt

# Definición de la función g(theta)
def g(theta):
    return np.maximum(0, 1 - 2/np.pi * np.abs(theta - np.pi))

# Función para calcular la integral de Laplace-Poisson
def laplace_poisson(r, theta):
    # Valores de phi para integración
    phi_values = np.linspace(0, 2*np.pi, 1000)
    # Definición del integrando
    integrand = g(phi_values) * (1 - r**2) / (1 + r**2 - 2*r*np.cos(theta - phi_values))
    # Calcular la integral usando el método trapezoidal
    return np.trapz(integrand, phi_values) / (2 * np.pi)

# Función para calcular la solución de Laplace
def laplace_solution(R, Theta):
    # Crear un array vacío para almacenar los valores de la solución de Laplace
    laplace_values = np.zeros_like(R)
    
    # Iterar sobre cada punto en la malla
    for i in range(len(R)):
        for j in range(len(R[0])):
            r = R[i, j]
            theta = Theta[i, j]
            
            if r < 1:  # Dentro de la bola unitaria la calculada
                laplace_values[i, j] = laplace_poisson(r, theta)
            else:  # En el límite exigimos la funcion g
                laplace_values[i, j] = g(theta)
    
    return laplace_values

# Generar malla de valores de r y theta
num_points = 100
r_values = np.linspace(0, 1, num_points)
theta_values = np.linspace(0, 2*np.pi, num_points)
R, Theta = np.meshgrid(r_values, theta_values)

# Calcular la solución de Laplace para cada punto en la malla
laplace_values = laplace_solution(R, Theta)

# Graficar
plt.figure(figsize=(8, 6))
plt.contourf(R * np.cos(Theta), R * np.sin(Theta), laplace_values, levels=50, cmap='viridis')
plt.colorbar(label='Solución de Laplace')
plt.title('Solución de Laplace en el disco unidad')
plt.xlabel('x')
plt.ylabel('y')
plt.axis('equal')
plt.grid(True)
plt.show()

4.2 Limitaciones de la fórmula de Poisson

La fórmula de Poisson enfrenta limitaciones cerca de las fronteras del dominio debido a la singularidad de la integral. Esto se manifiesta en errores de aproximación, especialmente cuando se utiliza la fórmula del trapecio para la integración. Para investigar estos errores, analizaremos la solución exacta [math]u(x,y)=xy[/math], que es armónica y mantiene su valor en la frontera [math]g(x,y)=xy[/math]. Mediante diferentes discretizaciones de la fórmula del trapecio, evaluaremos la precisión de la fórmula de Poisson y el impacto de la discretización en la aproximación de la solución.

Por comodidad y por lo ya realizado anteriormente, llevaremos a cabo un cambio a coordenadas polares.Siendo [math]u(r,\theta)=r^2 cos(\theta)\sin(\theta)[/math] y dado que tomaremos la bola unidad entonces la condición frontera se expresa como [math]g(\theta)= cos(\theta)\sin(\theta)[/math].

Para evaluar la precisión de estas aproximaciones, analizaremos el error en un punto específico, alejado de la frontera. El punto elegido para este análisis es [math] (r, \theta) = (0.9, \frac{\pi}{4})[/math].

El error será cuantificado en una escala logarítmica definida como [math] f(n) = \log_{10}(\text{error}(10^n))[/math], donde [math]\text{error}(10^n)[/math] representa el error computado al utilizar [math] 10^n [/math] puntos para la evaluación integral. Posteriormente, estos errores serán representados gráficamente para proporcionar una visión clara del comportamiento del error conforme aumenta el número de puntos de la discretización.

#Definimos Solución Exacta
def exact_solution(r, theta):
    return r**2 * np.sin(theta) * np.cos(theta)
#Definimos Condición Frontera
def g(theta):
    return np.sin(theta) * np.cos(theta)
#Solución de la ecuación laplace
def laplace_poisson(r, theta, num_points_trapz):
    phi_values = np.linspace(0, 2*np.pi, num_points_trapz)
    integrand = g(phi_values) * (1 - r**2) / (1 + r**2 - 2*r*np.cos(theta - phi_values))
    return np.trapz(integrand, phi_values) / (2 * np.pi)
#Calculamos la solución en todos los puntos
def laplace_solution(R, Theta, num_points_trapz):
    laplace_values = np.zeros_like(R)
    for i in range(len(R)):
        for j in range(len(R[0])):
            r = R[i, j]
            theta = Theta[i, j]
            laplace_values[i, j] = laplace_poisson(r, theta, num_points_trapz)
    return laplace_values
import numpy as np
import matplotlib.pyplot as plt

# Función para calcular el error de la solución de Poisson
def calculate_error(num_points_trapz):
    r = 0.9
    theta = np.pi / 4
    
    # Calcula el valor exacto de la solución en el punto dado
    exact_value = exact_solution(r, theta)
    
    # Aproxima el valor de la solución utilizando el método de la integral trapezoidal
    approx_value = laplace_poisson(r, theta, num_points_trapz)
    
    # Calcula el error absoluto entre la solución exacta y la aproximada
    error = np.abs(exact_value - approx_value)
    
    # Si el error es muy pequeño, asigna un valor pequeño para evitar errores en el logaritmo
    if error == 0:
        error = 1e-16 
    
    # Retorna el logaritmo en base 10 del error multiplicado por 10 elevado a la potencia de num_points_trapz
    return np.log10(error)

# Lista de valores para el número de puntos en la integral trapezoidal
num_points_trapz_values = [10**n for n in range(1, 8)]

# Calcula el error para cada valor de número de puntos
errors = [calculate_error(num_points) for num_points in num_points_trapz_values]

# Grafica el error en función del logaritmo del número de puntos
plt.figure(figsize=(8, 6))
plt.plot(np.log10(num_points_trapz_values), errors, color='blue',  marker='o', linestyle='-')
plt.title('Error vs. Numero de puntos para el trapecio')
plt.xlabel('$\log_{10}$(Numero de puntos para el trapecio)')
plt.ylabel('$\log_{10}$(error $(10^n))$')
plt.ylim(-16,5)
plt.grid(True)
plt.show()

En este punto queremos observar la relación que hay entre el error teórico de la integral al resolverla numéricamente por la regla del trapecio con el error obtenido, el fin de este estudio es ver cuan buena es la solución que reporta la fórmula de Poisson. Considerando la siguiente integral

[math] I = \int_a^bf(x)dx [/math]

Si la resolvemos mediante el método del trapecio, presentará el siguiente error

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

Donde n es el número de puntos de la discretización. En nuestro caso la función f sería:

[math] \psi(r,\theta, \phi)=\frac{1-r^2}{2\pi} \frac{ cos(\phi)\sin(\phi)}{1+r^2-2rcos(\theta-\phi)} = \frac{1-r^2}{4\pi} \frac{ \sin(2\phi)}{1+r^2-2rcos(\theta-\phi)} [/math]

Donde, por ejemplo, nos interesa calcular el error de la fórmula del trapecio en el punto [math] (r, \theta) = (0.9, \frac{\pi}{4})[/math]. Para ello primero sustituimos sendas variables en la función [math]\psi[/math]. Además nuestro caso se integra en el intervalo [math][0,2\pi][/math] por tanto [math]a=0,b=2\pi[/math].

Derivamos la expresión [math] \psi(0.9, \frac{\pi}{4}, \phi) [/math] dos veces respecto de [math] \phi [/math] y calculamos el máximo para de esta froma acotar el error como sigue,

[math] |E_t| = |-\frac{1}{12n^2}\psi''(0.9, \frac{\pi}{4}, \phi)(2\pi)^3| \leq \max_{\phi \in [0, 2\pi]} \left| -\frac{1}{12n^2} \psi''(0.9, \frac{\pi}{4}, \phi) (2\pi)^3 \right| [/math]

El máximo del valor absoluto de la segunda derivada lo calculamos de forma numérica pues, debido a la expresión, era difícil de calcularlo de manera simbólica; además se obtiene una cota más precisa del error que la que se obtendría acotando cada término. Dicho valor se alcanza aproximadamente en [math]\phi = \frac{\pi}{4}[/math], lo cual además tiene sentido por la simetría en el denominador de [math] \psi(r,\theta, \phi) [/math] que alcanza el mínimo cuando [math] \theta = \phi [/math]. Mostramos un código a continuación de como se calcula de manera simbólica la segunda derivada así como una representación del valor absoluto de la misma.

# Define la expresión
expr := ((1 - r^2)/(4*Pi))*sin(2*φ) / (1 + r^2 - 2*r*cos(θ - φ));

# Segunda derivada respecto a φ
expr_phi_2 := diff(expr_phi_1, φ);

# Simplificar la segunda derivada respecto a φ
a := simplify(expr_phi_2);

# Evaluar la segunda derivada simplificada para r = 0.9 y θ = 0.785398
r := 0.9;
θ := 0.785398;
b := abs(evalf(a));

# Graficar la segunda derivada evaluada para φ en el rango de [0, 2π]
plot([b], φ = 0 .. 2*Pi, legend = "Valor absoluto de la segunda derivada respecto a φ");

No obstante esto a priori no sirve para hallar el valor buscado, pues no sabemos si tomando una subdivisión mayor del intervalo la función presenta anomalías cerca de [math] \frac{\pi}{4}[/math] resultando así una valor distinto, no obstante es improbable pues [math]\psi(0.9, \frac{\pi}{4}, \phi)[/math] es diferenciable en [math][0,2\pi][/math]. Para hallar el máximo de esta función, implementaremos el siguiente código que divide una malla en distintos números de subintervalos y calcula los valores de la segunda derivada de una función con respecto a [math]\phi[/math] en cada subdivisión. Luego, registra el máximo en valor absoluto de estos valores y lo representa gráficamente en función del número de subintervalos en la malla.Cuando observamos que la curva que representa el máximo de la segunda derivada parece estancarse, sugiere que la función no presenta anomalías y por tanto el valor límite obtenido está cerca del valor real.

import numpy as np
import matplotlib.pyplot as plt

# Definimos la segunda derivada para acotar el error del método del trapecio
def double_prime(phi, r, theta):
    return np.abs(((2*(r**2 - 1)*(-r*(4*r*np.sin(phi - theta)**2/(r**2 - 2*r*np.cos(phi - theta) + 1) - np.cos(phi - theta))*np.sin(phi)*np.cos(phi)/(r**2 - 2*r*np.cos(phi - theta) + 1) -
    2*r*np.sin(phi)**2*np.sin(phi - theta)/(r**2 - 2*r*np.cos(phi - theta) + 1) + 2*r*np.sin(phi - theta)*np.cos(phi)**2/(r**2 - 2*r*np.cos(phi - theta) + 1) + 2*np.sin(phi)*np.cos(phi))/(r**2 - 2*r*np.cos(phi - theta) + 1)))/ (2 * np.pi))

# Función para encontrar el máximo de la segunda derivada en función del número de divisiones de la malla
def find_maximum(r, theta, num_divisions):
    phi_values = np.linspace(0, 2*np.pi, num_divisions)
    double_prime_values = double_prime(phi_values, r, theta)
    max_value = np.max(double_prime_values)
    return max_value

# Parámetros
r = 0.9
theta = np.pi/4
max_values = []
num_divisions_values = range(10**3, 10**5, 100)

# Calcular el máximo para cada número de divisiones
for num_divisions in num_divisions_values:
    max_value = find_maximum(r, theta, num_divisions)
    max_values.append(max_value)

# Graficar
plt.figure(figsize=(8, 6))
plt.xscale('log')
plt.plot(num_divisions_values, max_values, color='blue')
plt.xlabel('Número de divisiones de la malla')
plt.ylabel('Máximo de la segunda derivada')
plt.title('Máximo de la segunda derivada en función del número de divisiones de la malla')
plt.grid(True)

# Encontrar el valor que maximiza la segunda derivada
max_value_index = np.argmax(max_values)
optimal_num_divisions = num_divisions_values[max_value_index]
optimal_max_value = max_values[max_value_index]
optimal_phi = np.linspace(0, 2*np.pi, optimal_num_divisions)[np.argmax(double_prime(np.linspace(0, 2*np.pi, optimal_num_divisions), r, theta))]

print("El valor de phi que maximiza la segunda derivada es aproximadamente {:.6f}, con un valor máximo de {:.6f}".format(optimal_phi, optimal_max_value))
plt.show()

Este código reporta que el valor de [math]\phi[/math] que maximiza el valor absoluto de la segunda derivada es aproximadamente [math]0.785406\approx \frac{\pi}{4}[/math], con un valor máximo de 278.202831. Mostramos el error teórico dado por la expresión.

[math] g(n) = log\max_{\phi \in [0, 2\pi]}|E_t(n,\phi)|[/math]

import numpy as np
import matplotlib.pyplot as plt

# Definimos la segunda derivada para acotar el error del método del trapecio
def double_prime(phi, r, theta):
    return ((2*(r**2 - 1)*(-r*(4*r*np.sin(phi - theta)**2/(r**2 - 2*r*np.cos(phi - theta) + 1) - np.cos(phi - theta))*np.sin(phi)*np.cos(phi)/(r**2 - 2*r*np.cos(phi - theta) + 1) -
    2*r*np.sin(phi)**2*np.sin(phi - theta)/(r**2 - 2*r*np.cos(phi - theta) + 1) + 2*r*np.sin(phi - theta)*np.cos(phi)**2/(r**2 - 2*r*np.cos(phi - theta) + 1) + 2*np.sin(phi)*np.cos(phi))/(r**2 - 2*r*np.cos(phi - theta) + 1)))/ (2 * np.pi)

# Calculamos el valor de double_prime(phi, 1, pi/4) para el valor de phi dado
phi = 0.785406 # Calculado en el ejercicio anterior
r = 0.9
theta = np.pi/4
valor_double_prime = double_prime(phi, r, theta)

# Calculamos el error teórico de la regla del trapecio para n puntos
def theoretical_error(n):
    a = 0
    b = 2*np.pi
    return np.log(abs(-(b - a)**3 / (12 * n**2) * valor_double_prime ))

# Generamos distintos valores de n
num_puntos_trapz_valores = [10**n for n in range(1, 8)]
errores_teoricos = [theoretical_error(num_puntos) for num_puntos in num_puntos_trapz_valores]

# Imprimimos el error teórico
plt.figure(figsize=(8, 6))
plt.plot(np.log10(num_puntos_trapz_valores), errores_teoricos, color='red', marker='s', linestyle='--')
plt.title('Error Teórico vs. Número de Puntos en la Regla del Trapecio')
plt.xlabel('Número de Puntos en la Regla del Trapecio')
plt.ylabel('Error Teórico')
plt.ylim(-16,5)
plt.grid(True)
plt.show()

Comparamos el error teórico con el error obtenido entre la solución exacta y la solución calculada.

# Generamos valores para diferentes números de puntos en la regla del trapecio
num_points_trapz_values = [10**n for n in range(1, 8)]

# Calculamos errores para cada número de puntos utilizando la función calculate_error
errors = [calculate_error(num_points) for num_points in num_points_trapz_values]

# Calculamos errores teóricos para cada número de puntos utilizando la función theoretical_error
theoretical_errors = [theoretical_error(num_points) for num_points in num_points_trapz_values]

plt.figure(figsize=(8, 6))

# Graficamos los errores calculados en una escala logarítmica
plt.plot(np.log10(num_points_trapz_values), errors, color='blue', marker='o', linestyle='-', label='Error Práctico')

# Graficamos los errores teóricos en una escala logarítmica
plt.plot(np.log10(num_points_trapz_values), theoretical_errors,  color='red',marker='s', linestyle='--', label='Error Teórico')
plt.title('Error Práctico vs. Teórico para Diferentes N')
plt.xlabel('$\log_{10}$(Numero de puntos para el trapecio)')
plt.ylabel('$\log_{10}$(error $(10^n))$')
plt.ylim(-16,5)
plt.legend()
plt.grid(True)
plt.show()

Observamos que la solución obtenida por la fórmula de Poisson mediante la fórmula del trapecio es buena en un punto en el interior del dominio. Sin embargo, el error en el método del trapecio cerca de la frontera de la región de integración conforme nos acercamos a [math]r=1[/math] crece. Consideramos fijar el número de puntos en la fórmula del trapecio en un valor determinado 100. Luego, trazamos la gráfica del error para puntos de la forma [math](r, \theta) = (1 - 10^{-n}, \pi/4)[/math], donde [math]n[/math] representa la precisión de la aproximación. Se espera un aumento del error.

import numpy as np
import matplotlib.pyplot as plt

num_points_trapz = 100

# Función para calcular el error de la solución de Poisson
def calculate_error(num_points_trapz):
    r = 0.9
    theta = np.pi / 4
    
    # Calcula el valor exacto de la solución en el punto dado
    exact_value = exact_solution(r, theta)
    
    # Aproxima el valor de la solución utilizando el método de la integral trapezoidal
    approx_value = laplace_poisson(r, theta, num_points_trapz)
    
    # Calcula el error absoluto entre la solución exacta y la aproximada
    error = np.abs(exact_value - approx_value)
    
    # Si el error es muy pequeño, asigna un valor pequeño para evitar errores en el logaritmo
    if error == 0:
        error = 1e-16 
    
    # Retorna el logaritmo en base 10 del error multiplicado por 10 elevado a la potencia de num_points_trapz
    return np.log10(error)

# Lista de valores para el número de puntos en la integral trapezoidal
num_points_trapz_values = [10**n for n in range(1, 8)]

# Calcula el error para cada valor de número de puntos
errors = [calculate_error(num_points) for num_points in num_points_trapz_values]

# Función para calcular la diferencia entre la solución exacta y la aproximada mediante el método de Laplace-Poisson
def calculate_error(r, theta, num_points_trapz):
    exact_value = exact_solution(r, theta)  # Calcula la solución exacta
    approx_value = laplace_poisson(r, theta, num_points_trapz)  # Calcula la solución aproximada
    error = np.abs(exact_value - approx_value)  # Calcula el error absoluto entre las dos soluciones
    if error == 0:
        error = 1e-16  # Evita errores en el cálculo del logaritmo si el error es cero
    return np.log10(error)  # Retorna el logaritmo en base 10 del error absoluto

n_values = np.arange(1, 8)  # Genera valores de n desde 1 hasta 10

# Genera valores de r y calcula los errores para diferentes valores de n
num_points_trapz = 100  # Número de puntos en la regla del trapecio
errors = []
for n in n_values:
    r = 1 - 1 / 10**n  # Calcula el valor de r
    theta = np.pi / 4  # Ángulo theta fijo
    errors.append(calculate_error(r, theta, num_points_trapz))  # Calcula el error para el valor de r actual

# Grafica los errores en función de n
plt.figure(figsize=(8, 6))
plt.plot(n_values, errors, color = 'blue', marker='o', linestyle='-')
plt.title('Error conforme nos acercamos a la frontera')
plt.xlabel('$n$ que sitasifacen la ecuación $r=1 - 1/10^n$')  # Etiqueta del eje x
plt.ylabel('$\log_{10}$(error $(10^n))$')  # Etiqueta del eje y
plt.grid(True)
plt.show()

Efectivamente el error aumenta conforme aumenta [math]n[/math], como era de esperar debido al carácter singular de la integral.

4.3 Solución por Series de Fourier

Procedemos ahora a resolver este problema mediante el uso de las series de Fourier,siendo de nuevo la solución [math] u(x,y) =xy [/math], siendo su valor en la frontera [math] g(x,y) =xy [/math]. Tomando nuestro problema en coordenadas polares obtenemos

[math] \begin{cases} \Delta u=u_{rr}+\frac{1}{r}u_r+\frac{1}{r^2}u_{\theta\theta}=0, r\in(0,R),\theta\in[0,2\pi]\\ u(R,\theta)=G(\theta)=R^2 cos(\theta)sin(\theta), \theta\in[0,2\pi]. \end{cases} [/math]

Donde su solución 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]

Para que se cumpla la condición frontera se debe verificar [math] u(R,\theta)=G(\theta) [/math], donde R=1,al estar en la bola unidad. Siendo el desarrollo en series de Fourier de [math]G(\theta) [/math] considerando la base trigonométrica usual ,[math] L^2([-\pi,\pi]) [/math]:

[math] G(\theta) = \frac{a_0}{2} + \sum_{k=1}^\infty (a_k \cos(k\theta) + b_k \sin(k\theta)) [/math]

Donde los coeficientes de Fourier son

[math] a_0=\frac{1}{\pi} \int_{-\pi}^{\pi} G(\theta) \, d\theta [/math]

[math] a_k = \frac{1}{\pi} \int_{-\pi}^{\pi} G(\theta) \cos(k\theta) \, d\theta [/math]

[math] b_k = \frac{1}{\pi} \int_{-\pi}^{\pi} G(\theta)\sin(k\theta) \, d\theta [/math]

Mientras que los coeficientes de Fourier para [math] G(\theta)=\cos(\theta)\sin(\theta)=\frac{1}{2}\sin(2\theta)[/math] son [math] a_0 = 0, a_k = 0 \forall k , b_k = 0 \forall k \neq 2 [/math] y [math] b_2 = \frac{1}{2} [/math].

Esto significa que la serie de Fourier de [math] G(\theta)[/math] se reduce simplemente a [math]u(r,\theta) = \frac{1}{2} r^2\sin(2\theta)[/math]con [math]r\in [0,1][/math] y [math]\theta \in [0,2\pi][/math]

Debido a que la solución por serie de Fourier tiene solo un número finito de términos, no tiene mucho sentido dibujar la gráfica con el error en función del número de términos, ya que el error desaparece completamente al alcanzar el segundo término, y no hay variación posterior del error al aumentar el número de términos.

4.4 Desigualdad de Harnack

En esta sección, examinaremos la desigualdad de Harnack. Para ello, buscaremos la región la que deben estar todas las soluciones armónicas en la bola \( B_R \), que coinciden con \( u \), en \( (0,0) \)

La desigualdad de Harnack se aplica a una función \( f \) no negativa definida en una bola cerrada en \( \mathbb{R}^n \) con radio \( R \) y centro \( x_0 \). Luego, no podemos aplicarla directamente a la función \( u(x,y) = xy \), ya que toma valores negativos. Por lo tanto, generaremos una nueva función, traslación de la anterior, \( v = u - M \), donde \( M \) es el mínimo de la función \( u(x,y) \). Esta función es tanto positiva como armónica, ya que \( \Delta u = \Delta v = 0 \), lo que implica que cumple las hipótesis de Harnack y nos permite aplicar la desigualdad.

Comenzamos hallando el mínimo de la función \( g(x,y) = xy \) definida en \( \partial B_R(0) \). Pues vemos que la función \( u(x,y) = xy \), que expresada en polares responde a\( u(r,\theta) = r^2\sin(\theta)\cos(\theta) = \frac{r^2}{2}\sin(2\theta) \), alcanza su valor más bajo en la frontera cuando \( r \) tiende a \( R \) y el seno de \( \theta \) alcanza su mínimo, que es \( -1 \). De modo que

\(M=\min_{\partial B_R(0)}g(r,\theta) =\min_{\partial B_R(0)} R^2\sin(\theta)\cos(\theta)= \min_{\partial B_R(0)}\frac{R^2}{2}\sin(2\theta) =-\frac{R^2}{2}\)

Luego, la función a la que se le va a aplicar Harnack será \( v(r,\theta) = u(r,\theta) - M=\frac{r^2}{2}\sin(2\theta)+\frac{R^2}{2} \geq 0\)

A continuación, se mostrará la expresión para dimensión 2 en la bola de radio \(R\) centrada en el origen. Teniendo en cuenta que \( v(0,\theta) = u(0,\theta)-M=-M \)

\( \frac{(R-r)}{(R+r)} \frac{R^2}{2} \leq v(r,\theta) \leq \frac{(R+r)}{(R-r)}\frac{R^2}{2}, \forall (r,\theta) \in B_{R} \)

Dado que buscamos las cotas para \(u \), deshacemos el cambio de variable \(u(r,\theta)=v(r,\theta)+M \) y obtenemos, sumando a ambos lados \(-\frac{R^2}{2} \), la siguiente desigualdad

\( -\frac{(R^2r)}{(R+r)} \leq u(r,\theta) \leq \frac{(R^2r)}{(R-r)}, \forall (r,\theta) \in B_{R} \)

Finalmente, buscamos graficar estas dos cotas para representar la región en que deben estar todas las soluciones. Se harán tres gráficas para apreciar las diferencias a mediada que el radio aumenta

import numpy as np
import matplotlib.pyplot as plt

# Radios de las bolas y colores asociados
Radios = [1, 2, 10]
Colores = ['blue', 'red', 'purple']

# Establecer el número de filas y columnas para los subgráficos
num_filas = 1
num_columnas = len(Radios)

# Tamaño de las figuras
figsize = (15, 5)  # Ancho x Alto

# Crear una figura que contenga todas las subgráficas
plt.figure(figsize=figsize)

# Iterar sobre cada radio y color
for indice, (R, color) in enumerate(zip(Radios, Colores)):
    # Vector de la discretización para los radios
    radios_discretizacion = np.linspace(0, R - 0.01, 100)
    
    # Calcular las cotas inferiores y superiores utilizando comprensiones de listas
    cota_inf = [-(R**2 * r) / (R + r) for r in radios_discretizacion]  # Cota inferior
    cota_sup = [(R**2 * r) / (R - r) for r in radios_discretizacion]  # Cota superior
    
    # Desplazar los valores obtenidos para aplicar la escala logarítmica
    m = min(min(cota_inf), min(cota_sup)) - 1  # Encontrar el valor mínimo entre ambas cotas y restamos uno pues log(0) no está definido
    cota_inf = np.array(cota_inf) - m  # Ajustar la cota inferior
    cota_sup = np.array(cota_sup) - m  # Ajustar la cota superior
    
    # Crear subgráfico en la posición adecuada
    plt.subplot(num_filas, num_columnas, indice + 1)
    
    # Trazar las cotas inferior y superior
    plt.plot(radios_discretizacion, cota_inf, '--', color=color, label='Cota inferior (R={})'.format(R))
    plt.plot(radios_discretizacion, cota_sup, color=color, label='Cota superior (R={})'.format(R))
    
    # Rellenar el área entre las cotas
    plt.fill_between(radios_discretizacion, cota_inf, cota_sup, color=color, alpha=0.1)
    
    # Etiquetas y título
    plt.xlabel('Radios(r)')
    plt.ylabel('Cota')
    plt.legend(loc='upper left')
    plt.title("Bola de Radio {}".format(R))
    plt.grid(True)

# Ajustar el espaciado entre subgráficas para evitar superposiciones
plt.tight_layout()

# Mostrar todas las figuras
plt.show()

Observamos como el comportamiento a medida que aumenta el radio es similar. Para poder extraer otras conclusiones más notables tomaremos el eje y en escala logarítmica. De modo que no podremos graficar la cota dada por la función \(u\) pues el logaritmo de un valor negativo no está definido. Tomaremos la función \(v\) que ha sido construida como una función positiva. Mostraremos una gráfica dependiente del radio y en escala logarítmica para diferentes radios de la bola: \(R={1,2,10} \)

import numpy as np
import matplotlib.pyplot as plt

# Radios de las bolas y colores asociados
Radios = [1, 2, 10]
Colores = ['blue', 'red', 'purple']

# Establecer el número de filas y columnas para los subgráficos
num_filas = 1
num_columnas = len(Radios)

# Tamaño de las figuras
figsize = (15, 5)  # Ancho x Alto

# Crear una figura que contenga todas las subgráficas
plt.figure(figsize=figsize)

# Iterar sobre cada radio y color
for indice, (R, color) in enumerate(zip(Radios, Colores)):
    # Vector de la discretización para los radios
    radios_discretizacion = np.linspace(0, R - 0.01, 100)
    
    # Calcular las cotas inferiores y superiores utilizando comprensiones de listas
    cota_inf = [R**2*(R - r)/(2*(R + r)) for r in radios_discretizacion]  # Cota inferior
    cota_sup = [R**2*(R + r)/(2*(R - r)) for r in radios_discretizacion]  # Cota superior
    
    # Crear subgráfico en la posición adecuada
    plt.subplot(num_filas, num_columnas, indice + 1)
    
    # Trazar las cotas inferior y superior
    plt.plot(radios_discretizacion, cota_inf, '--', color=color, label='Cota inferior (R={})'.format(R))
    plt.plot(radios_discretizacion, cota_sup, color=color, label='Cota superior (R={})'.format(R))
    
    # Rellenar el área entre las cotas
    plt.fill_between(radios_discretizacion, cota_inf, cota_sup, color=color, alpha=0.1)
    
    # Etiquetas y título
    plt.xlabel('Radios(r)')
    plt.ylabel('Cota')
    plt.yscale('log')
    plt.legend(loc='upper left')
    plt.title("Bola de Radio {}".format(R))
    plt.grid(True)

# Ajustar el espaciado entre subgráficas para evitar superposiciones
plt.tight_layout()

# Mostrar todas las figuras
plt.show()

Dado que las escalas son diferentes, mostraremos una foto solapada con todas ellas para que sea más visible este resultado.

Ahora podemos apreciar con mayor exactitud que aquellas funciones armónicas que pasen por la intersección entre la cota superior e inferior para un radio preciso, se van a quedar encerradas entre sendas funciones.

Estudiemos ahora qué pasaría si cambiamos la dimensión del espacio, en vez de en una circunferencia, en una esfera. Empleamos de nuevo la función \(v(r,\theta, \phi) \) que posee la no negatividad para aplicar Harnack donde el mínimo responde a la misma forma, \(\frac{R^2}{2} \), pues su función no ha cambiado. Análogamente, obtenemos

\( R\frac{(R-r)}{(R+r)^2} \frac{R^2}{2} \leq v(r,\theta, \phi) \leq R\frac{(R+r)}{(R-r)^2}\frac{R^2}{2}, \forall (r,\theta, \phi) \in B_{R} \)

Y dado que buscamos las cotas de \(u(r,\theta, \phi)=v(r,\theta, \phi)-M \) sumando a ambos lados de las desigualdades M llegamos a la siguiente acotación:

\( \frac{(R-r)}{(R+r)^2}\frac{R^3}{2}- \frac{R^2}{2} \leq v(r,\theta, \phi) \leq R\frac{(R+r)}{(R-r)^2}\frac{R^3}{2}- \frac{R^2}{2}, \forall (r,\theta, \phi) \in B_{R} \)

Representado de nuevo las funciones que definen la cota superior e inferior en función del radio para valores \(R={1,2,10} \), obtenemos las siguientes imágenes:

import numpy as np
import matplotlib.pyplot as plt

# Radios de las bolas y colores asociados
Radios = [1, 2, 10]
Colores = ['blue', 'red', 'purple']

# Establecer el número de filas y columnas para los subgráficos
num_filas = 1
num_columnas = len(Radios)

# Tamaño de las figuras
figsize = (15, 5)  # Ancho x Alto

# Crear una figura que contenga todas las subgráficas
plt.figure(figsize=figsize)

# Iterar sobre cada radio y color
for indice, (R, color) in enumerate(zip(Radios, Colores)):
    # Vector de la discretización para los radios
    radios_discretizacion = np.linspace(0, R - 0.01, 100)

    # Calcular las cotas inferiores y superiores utilizando comprensiones de listas
    cota_inf = [R**2*(R*(R - r) - (R + r)**2)/(2*(R + r)**2) for r in radios_discretizacion]  # Cota inferior
    cota_sup = [R**2*(R*(R + r) - (R - r)**2)/(2*(R - r)**2)  for r in radios_discretizacion]  # Cota superior
    
    # Crear subgráfico en la posición adecuada
    plt.subplot(num_filas, num_columnas, indice + 1)
    
    # Trazar las cotas inferior y superior
    plt.plot(radios_discretizacion, cota_inf, '--', color=color, label='Cota inferior (R={})'.format(R))
    plt.plot(radios_discretizacion, cota_sup, color=color, label='Cota superior (R={})'.format(R))
    
    # Rellenar el área entre las cotas
    plt.fill_between(radios_discretizacion, cota_inf, cota_sup, color=color, alpha=0.1)
    
    # Etiquetas y título
    plt.xlabel('Radios(r)')
    plt.ylabel('Cota')
    plt.legend(loc='upper left')
    plt.title("Bola de Radio {}".format(R))
    plt.grid(True)

# Ajustar el espaciado entre subgráficas para evitar superposiciones
plt.tight_layout()

# Mostrar todas las figuras
plt.show()

De nuevo, realizaremos el cambio a escala logarítmica del eje y para apreciar mejor los resultados, de modo que la gráfica corresponderá a la función positiva \(v\).

import numpy as np
import matplotlib.pyplot as plt

# Radios de las bolas y colores asociados
Radios = [1, 2, 10]
Colores = ['blue', 'red', 'purple']

# Establecer el número de filas y columnas para los subgráficos
num_filas = 1
num_columnas = len(Radios)

# Tamaño de las figuras
figsize = (15, 5)  # Ancho x Alto

# Crear una figura que contenga todas las subgráficas
plt.figure(figsize=figsize)

# Iterar sobre cada radio y color
for indice, (R, color) in enumerate(zip(Radios, Colores)):
    # Vector de la discretización para los radios
    radios_discretizacion = np.linspace(0, R - 0.01, 100)

    # Calcular las cotas inferiores y superiores utilizando comprensiones de listas
    cota_inf = [R**3*(R - r)/(2*(R + r)**2) for r in radios_discretizacion]  # Cota inferior
    cota_sup = [R**3*(R + r)/(2*(R - r)**2)  for r in radios_discretizacion]  # Cota superior
    
    # Crear subgráfico en la posición adecuada
    plt.subplot(num_filas, num_columnas, indice + 1)
    
    # Trazar las cotas inferior y superior
    plt.plot(radios_discretizacion, cota_inf, '--', color=color, label='Cota inferior (R={})'.format(R))
    plt.plot(radios_discretizacion, cota_sup, color=color, label='Cota superior (R={})'.format(R))
    
    # Rellenar el área entre las cotas
    plt.fill_between(radios_discretizacion, cota_inf, cota_sup, color=color, alpha=0.1)
    
    # Etiquetas y título
    plt.xlabel('Radios(r)')
    plt.ylabel('Cota')
    plt.yscale('log')  # Escala logarítmica en el eje y
    plt.legend(loc='upper left')
    plt.title("Bola de Radio {}".format(R))
    plt.grid(True)

# Ajustar el espaciado entre subgráficas para evitar superposiciones
plt.tight_layout()

# Mostrar todas las figuras
plt.show()

Siendo la gráfica conjunta la siguiente

Donde los resultados que apreciamos no difieren en gran medida a los anteriores.

Finalmente, resulta significativo examinar las diferencias que surgen al cambiar la dimensión de análisis. Con el propósito de facilitar una comparación, hemos integrado las dos representaciones gráficas previamente presentadas, en las cuales se empleaba una escala logarítmica, en una sola imagen. Al emplear la escala logarítmica, se resaltan las diferencias de manera más equilibrada.

Podemos apreciar como las cotas van aumentando a medida que la dimensión crece y viceversa. Esto implica que la región de existencia de soluciones armónicas en una bola aumenta a medida que lo hace la dimensión.

5 Ecuación de Poisson bidimensional

Para calcular la solución de la ecuación de Poisson [math] \Delta u = f [/math], donde [math]\textbf{x} \in \mathbb{R}^2[/math] y [math]f[/math] es la función característica de la bola de radio 1, podemos utilizar el potencial logarítmico. A partir del potencial logarítmico podemos utilizar la solución fundamental para obtener la solución de la ecuación.

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

Por comodidad y la simetría de la función característica utilizaremos coordenadas polares. Por tanto las expresiones equivalentes son

Potencial logarítmico en polares [math] (x_1,x_2)=(r\cos(\theta), r\sin(\theta))[/math].

[math] \Phi(\textbf{x}) = \Phi(r,\theta)= -\frac{1}{2\pi} \log |(r\cos(\theta), r\sin(\theta))| = -\frac{1}{2\pi} \log r [/math]

Si además tomamos otro valor [math] \textbf{y} = (s\cos(\phi), s\sin(\phi)) [/math]. Entonces podemos expresar el potencial logarítmico en función de tres variables

[math] \Phi(\textbf{x}-\textbf{y}) = \Phi(r,\theta, s, \phi)=-\frac{1}{2\pi} \log |(r\cos(\theta), r\sin(\theta))-(s\cos(\phi), s\sin(\phi))| = -\frac{1}{2\pi} \log \left( \sqrt{r^2 - 2rs \cos(\theta - \phi) + s^2} \right) [/math]

Por otro lado función [math] f [/math] la escribimos como

[math] f(\textbf{y})=f(s)= \begin{cases} 1 & \text{si } s \leq 1 \\ 0 & \text{si } s \gt 1 \end{cases} [/math]

Por tanto la integral en coordenadas polares tiene la siguiente expresión:

[math] u(\textbf{x}) = u(r,\theta)= -\frac{1}{2\pi} \int_{0}^{2\pi} \int_{0}^{1} \log \left( \sqrt{r^2 - 2rs \cos(\theta - \phi) + s^2} \right) s \, ds \, d\phi [/math]

Resolvemos la integral de manera numérica mediante la integración por el método del trapecio dos veces (para la varible [math] s [/math] y [math] \phi [/math]), mostramos además la solución de forma tridimensional [math] (r\cos(\theta),r\sin(\theta), u(r,\theta)) [/math] para así ver la propiedad asintótica de la misma.

import numpy as np
import matplotlib.pyplot as plt
 
# Definimos la función característica
def characteristic_function(s, phi):
    # Devuelve 1 si s <= 1, 0 en caso contrario
    if s <= 1:
        return 1
    else:
        return 0
 
# Definimos el potencial logarítmico
def log_potential(r, theta, s, phi):
    # Calcula el potencial logarítmico para una posición (r, theta) y un punto (s, phi)
    return -1/(2*np.pi) * np.log(np.sqrt(r**2 - 2*r*s*np.cos(theta - phi) + s**2))*s
 
# Definimos la función para aproximar la integral doble
def double_integral_approximation(r, theta, num_segments_s, num_segments_phi):
    # Generamos los valores de s y phi
    s_values = np.linspace(0, 1, num_segments_s)
    phi_values = np.linspace(0, 2*np.pi, num_segments_phi)
   
    # Inicializamos una matriz para almacenar los valores de la función integrando
    integrand_values = np.zeros((num_segments_s, num_segments_phi))
   
    # Calculamos los valores de la función integrando para cada punto (s, phi)
    for i, s in enumerate(s_values):
        for j, phi in enumerate(phi_values):
            integrand_values[i, j] =  log_potential(r, theta, s, phi)*characteristic_function(s, phi)
   
    # Aproximamos la integral doble utilizando la regla del trapecio
    approximation = np.trapz(np.trapz(integrand_values, axis=1, dx=(2*np.pi)/num_segments_phi), axis=0, dx=1/num_segments_s)
    return approximation
 
# Generar una malla de puntos (r, theta)
num_points_r = 50
num_points_theta = 50
r_values = np.linspace(0.01, 10**3, num_points_r)
theta_values = np.linspace(0, 2*np.pi, num_points_theta)
R, Theta = np.meshgrid(r_values, theta_values)
 
# Calcular la aproximación de la integral doble en cada punto de la malla
num_segments_s = 50
num_segments_phi = 50
approximations = np.zeros_like(R)
for i in range(num_points_r):
    for j in range(num_points_theta):
        approximations[i, j] = double_integral_approximation(R[i, j], Theta[i, j], num_segments_s, num_segments_phi)
 
# Graficar
fig = plt.figure(figsize=(10, 6))
ax = fig.add_subplot(111, projection='3d')
surf = ax.plot_surface(R * np.cos(Theta), R * np.sin(Theta), approximations, cmap='viridis', edgecolor='none')
ax.set_title('Solución ecuación de Poisson')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('Aproximación')
plt.colorbar(surf, ax=ax, label='Aproximación')
plt.show()

Para verificar que la solución tiene el comportamiento esperado en el infinito, podemos estudiar cómo se comporta [math]u(\mathbf{x})[/math] cuando [math]|\mathbf{x}| \to \infty[/math]. O dicho de otra forma estudiar [math]\lim_{{r \to \infty}} u(r, \theta)[/math]. Asintóticamente la solución debería comportarse como

[math] u(\mathbf{x}) = -\frac{M}{2\pi} \log |\mathbf{x}| + O\left(\frac{1}{|\mathbf{x}|}\right) \quad |\mathbf{x}| \to +\infty [/math]

[math] \text{Donde} \quad M = \int_{\mathbb{R}^2} f(\mathbf{y}) \, d\mathbf{y} . [/math]

Que escrito en coordenadas polares y con nuestros datos,

[math] u(r,\theta) = -\frac{M}{2\pi} \log r + O\left(\frac{1}{r}\right) \quad r \to +\infty [/math]

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

Para realizar esta verificación, evaluaremos la solución en valores de [math]r[/math] cada vez más grandes dado un [math]\theta[/math] fijo, y compararemos estos valores con la función [math]u_{\infty}(r,\theta)=-\frac{1}{2} \log(r)[/math]. Si la solución se comporta correctamente, esperaríamos que los valores calculados se aproximen a los valores de la función [math]u_{\infty}(r,\theta)[/math] a medida que [math]r[/math] aumenta. Para ello calcularemos el error entre la solución obtenida y la función.

def expected_behavior(r):
    return -1/2 * np.log(r)-10/r
 
# Genera un rango de valores para r
r_values = np.linspace(0.1, 1000, 1000)  # Valores de r desde 0.1 hasta 10^3
 
# Calcula los valores de la solución aproximada y el comportamiento esperado en el infinito
approx_values = []  # Aquí guardarás los valores aproximados de la solución
expected_values = []  # Aquí guardarás los valores del comportamiento esperado en el infinito
for r in r_values:
    # Calcula la aproximación de la solución en el punto r y theta arbitrario
    approx_value = double_integral_approximation(r, np.pi/4, 100, 100)
    approx_values.append(approx_value)
   
    # Calcula el comportamiento esperado en el infinito
    expected_value = expected_behavior(r)
    expected_values.append(expected_value)
 
# Calcula la diferencia entre la aproximación y el comportamiento esperado en el infinito
differences = np.log10(np.abs(np.array(approx_values) - np.array(expected_values)))
 
# Grafica la diferencia en función de r
plt.plot(r_values, differences, color='blue')
plt.xlabel('r')
plt.ylabel('|Aprox - Comport. en el infinito|')
plt.title('Diferencia entre la aproximación y el comportamiento esperado en el infinito')
plt.grid(True)
plt.show()