Ecuación de Laplace y Poisson (equipo LUA)
| Trabajo realizado por estudiantes | |
|---|---|
| Título | Ecuación de Laplace y Poisson. Grupo LUA |
| Asignatura | EDP |
| Curso | 2023-24 |
| Autores | Luis Carreras Hoyos, Lucía Gil Ruiz y Alejandra Hernández Sieber |
| Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura | |
Contenido
- 1 Introducción
- 2 Conceptos previos
- 3 Ecuación de Laplace
- 4 Ecuación de Poisson en [math] \mathbb{R}^2 [/math]
- 5 Aplicaciones
- 6 Referencias
1 Introducción
Las ecuaciones de Poisson y Laplace son dos conceptos fundamentales en el campo de las ecuaciones diferenciales parciales. Estas ecuaciones tienen una amplia gama de aplicaciones en diversas disciplinas, desde la física y la ingeniería hasta las matemáticas aplicadas y la ciencia de datos; y son herramientas poderosas para comprender y analizar una variedad de problemas del mundo real.
Comencemos con la ecuación de Laplace, que aparece con frecuencia en el estudio de fenómenos de estado estacionario. Sus soluciones se denominan armónicas. La ecuación de Laplace matemáticamente se expresa como:
donde [math] u [/math] es una función escalar que representa el potencial y [math] \Delta [/math] es el operador laplaciano, que es la suma de las segundas derivadas parciales de [math] u [/math] con respecto a las coordenadas espaciales. En otras palabras, la ecuación de Laplace establece que el laplaciano del potencial es igual a cero en un dominio dado.
Por otro lado, la ecuación de Poisson es una generalización de la ecuación de Laplace y es muy importante en la teoría de los campos conservativos, como el campo eléctrico, magnético o gravitatorio. Matemáticamente, se expresa como:
En este trabajo, exploraremos en detalle la teoría detrás de la ecuación de Laplace y la ecuación de Poisson, examinando sus propiedades fundamentales, métodos de resolución y aplicaciones.
2 Conceptos previos
- Funciones armónicas: Diremos que una función [math] u(x) \in \mathcal{C}^2 (\Omega) [/math] es armónica si [math] \Delta u(x)=0 [/math]. Las funciones armónicas verifican la propiedad del valor medio (¡y viceversa!), esto es: Sea [math] u [/math] armónica en [math] \Omega \subset \mathbb{R}^n [/math], con [math] u \in \mathcal{C}^2(\Omega) [/math] y [math] u \in \mathcal{C}(\bar{\Omega}) [/math]. Entonces, para cualquier bola de centro [math] x \in \Omega [/math] y radio [math] R [/math] tal que [math] \mathbb{B}_R(x) \in \Omega [/math] se verifica:
- Unicidad del problema de Laplace y de Poisson en un dominio acotado: Se puede demostrar que siendo [math] \Omega \subset \mathbb{R}^n [/math] un dominio suave y acotado, entonces existe a lo sumo una solución [math] u(x) \in \mathcal{C}^2(\Omega) \cap \mathcal{C}^1(\bar{\Omega}) [/math], que satisface en [math] \partial \Omega [/math] una condición del tipo Dirichlet, Robin o mixta. En caso de la condición Neumann, la solución del problema es única salvo constante. Si se desea conocer más al respecto se recomienda consultar la sección [math] 3.2 [/math] de la página [math] 116 [/math] del libro Partial Differential Equations in Action. Springer. Sandro Salsa.
- Desigualdad de Harnack: Sea [math] u(x) [/math] una función armónica y [math]u \geq 0 [/math] en [math] \Omega \subset R^n[/math]. Supongamos que [math]\overline{ \mathbb{B}_{R}}(z)\subset \Omega [/math], entonces se verifica:
- Principio del máximo: Sea [math] \Omega \subset \mathbb{R}^n [/math] y [math] u \in \mathcal{C} (\Omega) [/math]. Si [math] u [/math] satisface la propiedad del valor medio y alcanza un máximo o mínimo en [math] \Omega [/math], entonces [math] u [/math] es constante. Como consecuencia, si [math] \Omega [/math] es acotado y la función [math] u [/math] no es constante, entonces:
3 Ecuación de Laplace
En esta primera sección, vamos a estudiar la ecuación de Laplace. Para comenzar, planteamos el siguiente sistema en dimensión [math] n [/math]:
Con [math] R \in \mathbb{R} [/math], [math] p [/math] un punto de [math] \mathbb{R}^n [/math] y [math] w_n [/math] la medida de la esfera de radio [math] 1 [/math]. La solución [math] u(x) \in \mathcal{C}^2(\mathbb{B}_{R}(p)) \cap \mathcal{C}(\overline{\mathbb{B}_{R}(p)})[/math] del sistema anterior viene dada por la fórmula de Poisson. Esta es:
Cabe destacar que esta fórmula es válida cuando la función [math] g [/math] es continua y para una dimensión [math] n \geq 2 [/math]. La demostración de ello puede encontrarse en la sección [math] 3.3.5 [/math] de la página [math] 127 [/math] del libro Partial Differential Equations in Action. Springer. Sandro Salsa.
3.1 Solución de la ecuación de Laplace por la fórmula de Poisson en un caso particular
Particularizamos el sistema anterior para el caso de la esfera unitaria con [math] R=1 [/math], centrada en el punto [math] p= (0,0) [/math], en dimensión [math] n=2 [/math] y la función [math] g(\theta) =max \left\{0, 1- \frac {2} {\pi} |\theta - \pi|\right\} [/math] expresada en coordenadas polares [math] (r, \theta) [/math]. Para dicho caso, tendríamos el siguiente sistema:
Calculamos ahora la solución de la ecuación de Laplace usando la fórmula de Poisson. Notar que, para dimensión [math] n=2 [/math], la expresión [math] w_2 = 2 \pi [/math]. Por tanto, la fórmula de Poisson en este caso será:
Siguiendo la demostración del Teorema [math] 3.13 [/math] en la página [math] 127 [/math] del libro Partial Differential Equations in Action. Springer. Sandro Salsa, al aplicar coordenadas polares, de manera que [math] x_1= r \cdot cos (\theta) [/math] y [math] x_2= r \cdot sen (\theta) [/math] se llega a que:
Volviendo a las coordenadas cartesianas, deshaciendo el cambio, esto equivale a:
Hay que destacar que, en la frontera no se puede usar la fórmula mencionada previamente debido a la singularidad de la integral. Es por ello por lo que hay que imponer directamente la condición frontera. Representamos ahora estos resultados en Matlab con el siguiente código:
clear all
close all
% Añadimos el radio de la bola
radio_bola=1;
% Añadimos la condición frontera
g=@(phi) max(0,1-2/pi.*abs(phi-pi));
% Definimos el número de puntos para la discretización
num_puntos_theta=600;
num_puntos_r=100;
erres=linspace(0,radio_bola,num_puntos_r);
thetas=linspace(0,2*pi,num_puntos_theta);
eses=thetas;
G=g(eses);
[Thetas,Erres]=meshgrid(thetas,erres);
% Calculamos la solución
sol=zeros(size(Thetas));
integral=zeros(size(Thetas));
for j=1:length(erres)
for i=1:length(thetas)
aux=G./(radio_bola^2+erres(j)^2-2*erres(j)*radio_bola*cos(eses-thetas(i)));
integral(j,i)=trapz(eses,aux);
end
end
coef=1/(2*pi)*(radio_bola^2-Erres.^2);
solucion=coef.*integral;
% Sustituimos los valores para r=1 por la condición frontera
solucion(end,:)=g(thetas);
% Dibujamos el resultado
figure(1)
colormap(jet)
surf(Thetas,Erres,solucion,'EdgeColor','none')
sgtitle('Representación de la solución')
colorbar
colormap('jet')
hold on
xlabel('theta')
ylabel('r')
figure(2)
sgtitle('Representación de la solución')
surf(Erres.*cos(Thetas), Erres.*sin(Thetas), solucion,'EdgeColor','none');
colorbar
colormap('jet')
xlabel('x_1')
ylabel('x_2')
En primer lugar, definimos los datos: el radio de la bola y la función establecida en la condición frontera. A continuación, realizamos la discretización del dominio pasando así de un problema continuo a un problema discreto. Posteriormente definimos y calculamos la solución dada por la fórmula de Poisson habiendo aplicado coordenadas polares, aproximando dichas integrales mediante la fórmula del trapecio. Después, como se ha mencionado previamente, se establece la condición frontera directamente como solución en la frontera. Por último, representamos los resultados, siendo la primera figura la solución en coordenadas polares [math] (r,\theta) [/math] y la segunda figura en coordenadas cartesianas [math] (x_1,x_2) [/math] una vez deshecho el cambio.
En la gráfica representada en coordenadas polares, resulta destacable como la solución obtenida es continua pero no derivable en la frontera de la bola de radio [math] r=1 [/math]. Esto también se puede apreciar con claridad en la gráfica representada en coordenadas cartesianas, en concreto en el punto [math] (x_1,x_2)=(-1,0) [/math]. Este comportamiento se debe directamente a la imposición de la condición frontera como solución en la frontera, pues ésta no es derivable y sus propiedades se traspasan a la solución.
Por último, se puede observar cómo se verifica el principio del máximo en ambas gráficas. En este caso, la solución obtenida no es constante, lo cual indica que el máximo y el mínimo de la misma se encuentran en la frontera de su dominio de definición, pues dicho dominio es acotado y, por ello, se puede aplicar el principio del máximo. Es decir, se cumple que:
Siendo [math] \Omega [/math] el dominio de definición de la solución, en este caso, la esfera unitaria [math] \mathbb{B}_1(0) [/math].
3.2 Cálculo de los errores con la fórmula de Poisson
La fórmula de Poisson no es del todo precisa. Cuenta así con algunas limitaciones, como la ya mencionada previamente acerca del carácter singular de la integral cerca de la frontera del dominio. Recordamos que para solucionar este problema imponíamos directamente la condición frontera en la solución.
Por otro lado, con respecto al cálculo de la solución empleando la fórmula de Poisson, también se producen algunos errores. Estos se deben principalmente a la incorporación del método del trapecio para calcular dichas integrales, ya que al fin y al cabo sigue siendo un método de aproximación numérica. Este inconveniente será el objetivo de esta sección.
3.2.1 Errores cometidos con la fórmula de Poisson variando discretizaciones
Para estudiar exactamente la magnitud de dichos errores, vamos a elegir una solución exacta, [math] u(x,y)=x \cdot y [/math]. Es sencillo comprobar que dicha función es armónica dado que
Además, su valor en la frontera de la bola [math] \mathbb{B}_1(0) [/math] también es [math] g(x,y)=x \cdot y [/math]. Es decir, la función [math] u(x,y) [/math] es la solución exacta del problema:
Usando el código de la sección anterior, podemos visualizar esta solución para una mayor comprensión del problema:
Sin embargo, al igual que ocurría en el problema anterior, sabemos que la solución de este problema viene dada por la fórmula de Poisson.
En primer lugar, vamos a expresar la solución exacta [math] u(x,y) [/math] en coordenadas polares por cuestión de comodidad, quedando así que [math] u(r,\theta)= r^2 \cdot cos(\theta) \cdot sin(\theta) [/math] y la fórmula de Poisson en coordenadas polares es:
Esta integral la vamos a resolver de forma numérica usando la fórmula del trapecio y calcularemos el error tomando diferentes discretizaciones para dicha fórmula, por lo que en primer lugar debemos obtener la solución dada por la fórmula de Poisson tomando todas esas discretizaciones.
Por lo tanto, en este estudio podemos fijar un punto concreto de la bola [math] \mathbb{B}_1 (0) [/math]. En particular, hemos tomado el punto de estudio (en coordenadas polares) [math] (r, \theta) = \left(0.9, \frac{\pi}{4}\right) [/math]. Además, hemos estudiado los errores tomando [math] 10^n [/math] puntos al evaluar el integrando aplicando la fórmula del trapecio para [math] n=1,\dots,8 [/math]. Para calcular los errores correspondientes, hemos elaborado un código en Matlab. Sin embargo, este código incluye estudios que realizaremos posteriormente, por lo que lo explicaremos al final de esta subsección. En concreto, hemos obtenido la siguiente gráfica del error:
En primer lugar, es importante mencionar que el método del trapecio es un método de aproximación numérica, y la precisión de la aproximación aumenta a medida que se aumenta el número de subintervalos. En la imagen se puede observar cómo los errores disminuyen drásticamente cuando tomamos una discretización lo suficientemente precisa, como es [math] 10^ {3} [/math] puntos. Aun así, hay que destacar lo pequeños que son estos errores, incluso el obtenido para una discretización de [math] 10^1 [/math] puntos, cuyo valor es [math] 0.21 [/math]. De hecho, el error incluso llega a alcanzar el valor de [math] 6,66 \cdot 10^{-16} [/math] cuando se toman [math] 10^6 [/math] puntos en la discretización. Por último, es llamativo como a pesar de que aumentamos el número de puntos para la discretización en la fórmula del trapecio, los errores aumentan ligeramente en alguno de los puntos. Esto puede observarse en valores como [math] n=7 [/math] u [math] n=8 [/math]. Realmente, no tenemos una explicación fundamentada para este fenómeno. La hipótesis que barajamos es que se trata de valores tan pequeños que es posible que la versión de Matlab de la que disponemos no sea lo suficientemente precisa para manejar dichos datos, dando lugar a estas posibles contradicciones.
Además, podemos comparar los errores obtenidos con la estimación del error máximo usando la fórmula del error teórica que se obtiene para la fórmula del trapecio, que es la siguiente:
donde en nuestro caso [math] \left[a,b\right]=\left[0,2\pi \right] [/math] es el dominio de integración, [math] n[/math] el número de subintervalos, [math]c[/math] un número cualquiera del intervalo [math] \left[0,2\pi \right] [/math] y [math] f [/math] es el integrando de la integral que resolvemos por la fórmula del trapecio. Es decir,
Para más información con respecto a la fórmula del error del trapecio, consultar la primera referencia citada al final del trabajo.
Teniendo esto en cuenta, lo primero que debemos calcular es [math] f’’(s) [/math], que es:
A continuación, vamos a estimar el error máximo a partir de esta fórmula. Para ello, debemos acotar la misma. Teniendo en cuenta las expresiones del seno y el coseno del ángulo doble, y que [math] sin(\phi) \leq 1 [/math] y [math] cos(\phi) \leq 1 [/math] para cualquier ángulo [math] \phi [/math], podemos acotar el valor absoluto de la segunda derivada por:
Por otro lado, es muy importante mencionar que en la fórmula del trapecio hemos optado por evaluar la segunda derivada en el punto crítico que hace máximo su valor absoluto. Sin embargo, hemos calculado dicho máximo de manera numérica. Esta elección implica correr ciertos riegos, pues el cálculo numérico siempre va acompañado de un cierto error de cálculo. En busca de minimizar este error, hemos optado por hacer una discretización con [math] 10^7 [/math] puntos para calcular dicho máximo.
Una vez hemos estimado el error máximo a partir del error teórico de la fórmula del trapecio, hemos representado todos en una gráfica para lograr comparalos.
En la primera gráfica, se ha representado el error numérico, ya estudiado anteriormente, el error teórico y la cota calculada. En primer lugar, recordamos que la fórmula del trapecio es un método de aproximación numérica y, por lo tanto, la precisión de la aproximación ha de aumentar a medida que se aumenta el número de subintervalos. Ya habíamos comentado que esto sí sucedía con el error numérico, pero gracias a esta gráfica también podemos apreciar que esto también ocurre para el error teórico.
Con respecto a la cota, recordamos que la expresión obtenida depende únicamente de la variable [math] r [/math] y como en este estudio, [math] r [/math] toma un valor fijo, es lógico que dicha cota sea constante con respecto al número de puntos de las discretizaciones. Sin embargo, al dividir entre [math] \frac{(2 \pi)^3}{12 \cdot n^2} [/math] presente en la fórmula del error teórica para la fórmula del trapecio se obtiene la gráfica vista previamente.
Otro punto que debemos mencionar es que el error teórico siempre es mayor que el error numérico, pues hemos tomado el máximo de [math] f’’(s) [/math].
Con respecto a la última gráfica, esta representa la diferencia entre ambos errores. En ella, se puede apreciar claramente como dicha diferencia se acerca a cero a medida que aumentamos el número de puntos en la discretización. De hecho, el error numérico apenas disminuye a partir de los [math] 10^3 [/math] puntos en la discretización, mientras que el error teórico no deja de disminuir.
Para terminar esta subsección, vamos a explicar el código elaborado para obtener todos los resultados presentados. Dicho código es el siguiente:
clear all
close all
format long
% Definimos la solución exacta en coordenadas polares
u=@(r,theta) r.^2.*cos(theta).*sin(theta);
% Añadimos el radio de la bola
radio_bola=1;
% Añadimos la condición frontera
g=@(theta) u(radio_bola,theta);
% Definimos el número de puntos para la discretización
num_puntos_theta=600;
num_puntos_r=100;
% Definimos los puntos a estudiar
erres=0.9;
thetas=pi/4;
discretizacion_maximo=linspace(0,2*pi,10^7);
% Definimos las distintas discretizaciones para la fórmula del trapecio
num_puntos_trapz=[1,2,3,4,5,6,7,8];
% Definimos el coeficiente
coef=1/(2*pi)*(radio_bola^2-erres.^2);
% Definimos la segunda derivada para la fórmula del error teórico
f2=@(erre,theta,ese) ((erre.^2-1).*(4.*erre.^2.*cos(ese).*sin(ese).*sin(ese-theta).^2+((4.*erre.^2.*cos(ese).^2-4.*erre.^2.*sin(ese).^2).*cos(ese-theta)+(2.*erre.^3+2.*erre).*sin(ese).^2+(-2.*erre.^3-2.*erre).*cos(ese).^2).*sin(ese-theta)-6.*erre.^2.*cos(ese).*sin(ese).*cos(ese-theta).^2+(7.*erre.^3+7.*erre).*cos(ese).*sin(ese).*cos(ese-theta)+(-2.*erre.^4-4.*erre.^2-2).*cos(ese).*sin(ese)))./(pi.*(2.*erre.*cos(ese-theta)-erre.^2-1).^3);
cotaerror=@(erre) (1-erre.^2).*(erre.^4+11/2.*erre.^3+7.*erre.^2+11/2.*erre+1)/(pi.*(1-erre).^6);
errores_numericos=zeros(1,length(num_puntos_trapz));
errores_teoricos=zeros(1,length(num_puntos_trapz));
coeficientes=[];
% Calculamos la solución
for discr=1:length(num_puntos_trapz)
eses=linspace(0,2*pi,10^num_puntos_trapz(discr));
G=g(eses);
sol=zeros(size(thetas));
integral=zeros(size(thetas));
for j=1:length(erres)
for i=1:length(thetas)
aux=G./(radio_bola^2+erres(j)^2-2*erres(j)*radio_bola*cos(eses-thetas(i)));
integral(j,i)=trapz(eses,aux);
end
end
solucion=coef.*integral;
% Calculamos los errores
errores_numericos(discr)=abs(solucion-u(erres,thetas));
errores_teoricos(discr)=(((2*pi)^3/(12*(length(eses)-1)^2)))*max(abs(f2(erres,thetas,discretizacion_maximo)));
coeficientes=[coeficientes,(((2*pi)^3/(12*(length(eses)-1)^2)))];
end
diferencia_errores=abs(errores_numericos-errores_teoricos);
% Representamos gráficamente los errores numéricos
figure(1)
semilogy(num_puntos_trapz,errores_numericos,'o-','LineWidth',1.5);
title('Error para el punto con r=0.9 y theta=pi/4 para distintas discretizaciones')
xlabel('Valores de n para las 10^n discretizaciones')
ylabel('Error en escala logarítmica')
% Representamos gráficamente los errores numéricos y teóricos
figure(2)
semilogy(num_puntos_trapz,errores_numericos,'o-', 'LineWidth',1.5)
hold on
semilogy(num_puntos_trapz,errores_teoricos,'o-','LineWidth',1.5)
semilogy(num_puntos_trapz,coeficientes.*cotaerror(0.9),'g','LineWidth',1.5)
title('Errores para distintas discretizaciones en escala logarítmica')
xlabel('Valores de n para las 10^n discretizaciones')
ylabel('Errores en escala logarítmica')
legend('Error númerico','Error teórico','Cota calculada')
figure(3)
semilogy(num_puntos_trapz,diferencia_errores,'o-','LineWidth',1.5);
title('Diferencia del error con r=0.9 y theta=pi/4 para distintas discretizaciones')
xlabel('Valores de n para las 10^n discretizaciones')
ylabel('Diferencia de los errores en escala logarítmica')
En primer lugar, definimos los datos necesarios: la solución exacta (expresada en coordenadas polares), el radio de la bola y la condición frontera. A continuación, definimos el número de puntos para la discretización, los puntos a estudiar y un vector para las distintas discretizaciones que usaremos para el método del trapecio. Seguidamente, definimos el coeficiente que acompaña a la integral dado que no depende del parámetro [math] s [/math] y, por tanto, no influirá a la hora de derivar la expresión. Por último, definimos la segunda derivada para la fórmula del error teórico.
Cabe destacar que puesto que la variable [math] \theta [/math] y la variable [math] s [/math] pertenecen al mismo dominio [math] [0, 2\pi] [/math], podríamos haber definido ambas como la misma variable en el código. Sin embargo, el objetivo de este código es estudiar lo que ocurre cuando aumentamos el número de puntos en la discretización de la variable sobre la que integramos, [math] s [/math], llegando a tomar valores muy altos de ese valor. Por ello, tomar [math] s [/math] y [math] \theta [/math] como la misma variable, aumenta mucho la complejidad y el tiempo de ejecución, es por eso por lo que hemos optado por definirlas en variables distintas, con distinto número de puntos para sus respectivas discretizaciones.
Posteriormente, calculamos la solución de manera análoga a la anterior sección, tomando en este caso como solución la integral obtenida multiplicada por el coeficiente definido anteriormente. Además, vamos guardando registro de los errores cometidos (numérica y teóricamente). Por último, representamos los errores numéricos gráficamente en escala logarítmica y los comparamos gráficamente con los errores teóricos.
3.2.2 Errores cometidos con la fórmula de Poisson cerca de la frontera
En último lugar nos interesa estudiar los errores cometidos en los puntos cercanos a la frontera de la bola unidad [math] \mathbb{B}_1(0) [/math]. Para ello, fijaremos un número concreto de puntos en la fórmula del trapecio. Por ejemplo, para [math] n=2 [/math]. Es decir, [math] 10^2=100 [/math] puntos.
Dado que la integral de la fórmula de Poisson cuenta con una singularidad en su frontera, estudiaremos dicho error en puntos lo suficientemente cercanos a esta. Tomamos entonces puntos de la forma [math] (r, \theta) = (1-10^ {-n}, \pi / 4) [/math].
Para ello hemos elaborado el siguiente código en Matlab:
clear all
close all
format long
% Definimos la solución exacta en coordenadas polares
u=@(r,theta) r.^2.*cos(theta).*sin(theta);
% Añadimos el radio de la bola
radio_bola=1;
% Añadimos la condición frontera
g=@(theta) u(radio_bola,theta);
% Definimos el número de puntos para la discretización
num_puntos_theta=600;
num_puntos_r=100;
% Definimos los puntos a estudiar y la matriz de mallado
erres=[1-10^-1,1-10^-2,1-10^-3,1-10^-4,1-10^-5,1-10^-6,1-10^-7];
thetas=pi/4;
[Thetas,Erres]=meshgrid(thetas,erres);
discretizacion_maximo=linspace(0,2*pi,10^7);
% Fijamos el valor de la discretización para la fórmula del trapecio
num_puntos_trapz=10^2;
eses=linspace(0,2*pi,num_puntos_trapz);
% Definimos la segunda derivada para la fórmula del error teórico
f2=@(erre,theta,ese) ((erre.^2-1).*(4.*erre.^2.*cos(ese).*sin(ese).*sin(ese-theta).^2+((4.*erre.^2.*cos(ese).^2-4.*erre.^2.*sin(ese).^2).*cos(ese-theta)+(2.*erre.^3+2.*erre).*sin(ese).^2+(-2.*erre.^3-2.*erre).*cos(ese).^2).*sin(ese-theta)-6.*erre.^2.*cos(ese).*sin(ese).*cos(ese-theta).^2+(7.*erre.^3+7.*erre).*cos(ese).*sin(ese).*cos(ese-theta)+(-2.*erre.^4-4.*erre.^2-2).*cos(ese).*sin(ese)))./(pi.*(2.*erre.*cos(ese-theta)-erre.^2-1).^3);
cotaerror=@(erre) (1-erre.^2).*(erre.^4+11/2.*erre.^3+7.*erre.^2+11/2.*erre+1)./(pi.*(1-erre).^6);
% Calculamos la solución
G=g(eses);
sol=zeros(size(Thetas));
integral=zeros(size(Thetas));
errores_teoricos=zeros(size(Thetas));
for j=1:length(erres)
for i=1:length(thetas)
aux=G./(radio_bola^2+erres(j)^2-2*erres(j)*radio_bola*cos(eses-thetas(i)));
integral(j,i)=trapz(eses,aux);
end
errores_teoricos(j)=(((2*pi)^3/(12*(length(eses)-1)^2)))*max(abs(f2(erres(j),thetas,discretizacion_maximo)));
end
coef=1/(2*pi)*(radio_bola^2-Erres.^2);
solucion=coef.*integral;
% Calculamos los errores
errores_numericos=abs(solucion-u(Erres,Thetas));
diferencia_errores=abs(errores_numericos-errores_teoricos)
% Dibujamos los errores en función de los valores de r
figure(1)
semilogy(erres,errores_numericos,'o-','LineWidth',1.5);
title('Error para distintos valores de r')
xlabel('Valores de r')
ylabel('Error en escala logarítmica')
figure(2)
semilogy(erres,errores_numericos, 'LineWidth',1.5)
hold on
semilogy(erres,errores_teoricos,'LineWidth',1.5)
semilogy(erres,(((2*pi)^3/(12*(length(eses)-1)^2))).*cotaerror(erres),'g--','LineWidth',1.5)
title('Errores para distintos valores de r')
xlabel('Valores de r')
ylabel('Errores en escala logarítmica')
legend('Error númerico','Error teórico','Cota del error')
figure(3)
semilogy(erres,diferencia_errores,'o-','LineWidth',1.5);
title('Diferencia de los errores con theta=pi/4 para distintos valores de r')
xlabel('Valores de r')
ylabel('Diferencia de los errores error en escala logarítmica')
Primero de todo, establecemos los datos necesarios: la solución exacta en coordenadas polares, el radio de la bola y la condición frontera. Seguidamente definimos el número de puntos para la discretización y definimos los puntos cercanos a la frontera de la forma [math] (r, \theta) = (1-10^ {-n}, \pi / 4) [/math] con [math] n=1,2,3,4,5,6,7 [/math]. Posteriormente definimos la discretización elegida para el método del trapecio y calculamos la solución de manera análoga a la realizada en la sección [math] 3.1 [/math]. Restando esta solución a la exacta obtenemos los errores numéricos.
Con respecto al error teórico seguimos el mismo proceso realizado anteriormente. Definimos la segunda derivada y la vamos calculando en cada una de nuestras discretizaciones. Tomamos el valor máximo, y almacenamos estos resultados en una lista. Además, tomamos la misma cota de error teórica del apartado anterior:
Por último, se representa gráficamente el error teórico, el error numérico, la diferencia de ambas y la cota teórica.
En las gráficas de los errores, estos van aumentando a medida que nos acercamos a la frontera. Cuando la función es suave y continua, los métodos de integración numérica son bastante precisos. Sin embargo, cerca de la frontera, la singularidad de la función dificulta la aproximación precisa de la integral.
Además, a pesar de que ambos errores aumentan, el error teórico alcanza valores más significativos. Tanto que hemos optado por representar el error numérico en una sola gráfica para apreciarla mejor. Los valores tan altos del error teórico se deben a que estamos tomando el máximo de la segunda derivada del integrando. Estos valores tienden a infinito cada vez que nos vamos acercando más a la frontera. Asimismo, sucede con la cota del error; al establecer una cota superior para nuestro error teórico, esta diverge a una velocidad mayor.
Con respecto a la gráfica de la diferencia, esta va aumentando a medida que nos acercamos a la frontera. Como se ha explicado, el error teórico toma valores cada vez más cercanos a infinito, demostrando así la dificultad que presenta la fórmula de Poisson cerca de la frontera, el principal inconveniente que describimos al principio de esta sección.
3.3 Solución de la ecuación de Laplace usando series de Fourier
En esta sección, vamos a calcular la solución de la ecuación de Laplace por series de Fourier. Consideramos nuevamente el problema
La resolución de este apartado está basada en la prueba vista en clase para demostrar que la solución de este problema viene dada por la fórmula de Poisson. Nuevamente, esta puede encontrarse como el Teorema [math] 3.13 [/math] de la página [math] 127 [/math] del libro Partial Differential Equations in Action. Springer. Sandro Salsa.
En primer lugar, es importante destacar que vamos a trabajar en coordenadas polares, por lo que reescribimos el problema usando que [math] u(x,y)=U(r,\theta) [/math] y [math] g(x,y)=G(\theta) [/math]:
La primera expresión equivale al laplaciano expresado en coordenadas polares. Esta expresión y otras fórmulas así expresadas pueden consultarse en el Apéndice B.1 en la página [math] 665 [/math] del libro Partial Differential Equations in Action. Springer. Sandro Salsa.
Regresando al problema, puesto que [math] u [/math] debe de ser continua en [math] \overline{\mathbb{B}_1(0)} [/math], entonces [math] U [/math] y [math] G [/math] tienen que ser continuas en [math] [0,1] \times [0, 2 \pi] [/math] y [math] [0,2 \pi] [/math] respectivamente. Por lo tanto, se concluye que [math] G(\theta) [/math] y [math] U(r,\theta) [/math] tienen que ser [math] 2\pi [/math]-periódicas y [math] U(0, \theta) [/math] no puede depender del parámetro [math] \theta [/math].
A continuación, aplicamos separación de variables de manera que [math] U(r,\theta)=R(r) \cdot T(\theta) [/math] y obtenemos un sistema para cada una de las funciones [math] T(\theta) [/math] y [math] R(r) [/math]. Por un lado,
Estudiando los casos en función del signo de [math] \lambda [/math], obtenemos como solución para este sistema la familia de funciones [math] \left\{sin(k \theta),cos(k\theta) \right\}_{k \in \mathbb{N}} [/math], habiendo tomado [math] \sqrt \lambda=k \in \mathbb{N}[/math].
Con respecto al caso de [math] R(r) [/math], tenemos la ecuación
Teniendo en cuenta que en el origen nuestra función solución tiene que ser continua, obtenemos el conjunto de funciones [math] \left\{ \frac{1}{2}, r^{k}cos(k \theta),r^{k}sin(k \theta) \right\}_{k \in \mathbb{N}} [/math]. Notar que, dicha familia forma una base de las funciones [math] L^2([-\pi,\pi]) [/math]. Seguidamente, aplicamos el principio de superposición obteniendo la siguiente expresión:
Finalmente, ajustamos la condición frontera, que recordamos que era [math]U(1,\theta)=\sin(\theta)\cos(\theta)[/math]. A continuación, teniendo en cuenta que [math] U(1,\theta)=sin(\theta)cos(\theta)=\frac{sin(2\theta)}{2}[/math], los coeficientes de Fourier que debemos imponer para que se verifique la condición frontera son [math] a_0=0, ~~ a_k=0 ~~ \forall k \in \mathbb{N}[/math] y [math] b_k=0 ~~ \forall k \neq 2 ~~\mathrm{con}~~ k \in \mathbb{N} [/math] y [math] b_2 = \frac{1}{2} [/math]. Por lo tanto, la solución al sistema original es:
Por último, hay que destacar que la serie de Fourier tiene solo un número finito de términos, y por tanto no tiene mucho sentido comparar, en este caso, los errores cometidos término a término con la solución exacta.
3.4 Desigualdad de Harnack
En esta última subsección vamos a analizar la desigualdad de Harnack para nuestro problema. Como hemos mencionado en los conceptos previos, para que se verifique dicha desigualdad se tiene que cumplir que [math] u [/math] sea una función armónica con [math]u \geq 0 [/math] en [math] \Omega \subset R^n[/math]. La desigualdad se verifica en puntos dentro de una bola [math]\overline{\mathbb{B}_{r}(z)}\subset \Omega [/math].
En nuestro caso, la dimensión es [math]n=2[/math]. Consideramos ahora la bola [math] \mathbb{B}_{1}(0)[/math]. En ella se verifica que la función armónica [math]v :=u-M \geq 0 [/math], siendo [math]u[/math] la solución a nuestro problema y [math]M[/math] el mínimo de [math]u(x,y)=x \cdot y [/math]. Sin embargo, como [math]u[/math] es armónica y por tanto, verifica el Principio del máximo, para buscar el mínimo de [math]u(x,y) [/math] en [math] \mathbb{B}_{1}(0)[/math], basta limitarnos a buscar el mínimo en su frontera, es decir, el mínimo de [math] g(x,y) [/math] en [math] \partial \mathbb{B}_{1}(0)[/math].
Para ello, expresaremos la condición frontera [math] g(x,y) = x \cdot y [/math] en polares resultando así que [math]g(\theta)=cos(\theta)sin(\theta)=\frac{sin(2\theta)}{2}[/math].
Calculamos ahora los puntos críticos de dicha función y evaluamos en su segunda derivada para comprobar cuáles de dichos puntos críticos son mínimos, obteniendo finalmente que el mínimo se encuentra en [math]g\left(\frac{3\pi}{4}\right)=g\left(\frac{7\pi}{4}\right)=-\frac{1}{2} [/math].
Por lo tanto, podemos definir nuestra función [math] v [/math] como [math] v:=u+\frac{1}{2}\geq 0[/math]. Aplicando ahora la desigualdad de Harnack a [math] v [/math] en la bola [math] B_{1}(0)[/math] obtenemos:
Sustituyendo [math]v(0)=u(0)+\frac{1}{2}[/math] y generalizando la desigualdad para cualquier función armónica [math]w(x) [/math] con [math]w(0)=u(0)[/math] en [math] \mathbb{B}_{1}(0)[/math] logramos la desigualdad:
Una vez obtenidas estas desigualdades, podemos representar en Matlab la región en las que están todas las soluciones armónicas [math]w(x)[/math] en [math] \mathbb{B}_{1}(0)[/math] que valen lo mismo que [math] u [/math] en el punto [math] (0,0) [/math]. Cabe destacar que durante toda esta sección vamos a emplear continuamente dos códigos distintos y ambos se explican brevemente al final de esta.
Tomando las dos funciones de los extremos de la desigualdad obtenemos:
En la primera gráfica, la región no se puede apreciar bien pues cuando r es próximo a 1, la función que constituye la frontera superior de la región tiende a infinito. Por ese mismo motivo hemos representado una segunda gráfica tomando valores de [math]r \in [0,0.9][/math]. Como se puede apreciar en esta segunda gráfica, la región sigue el mismo comportamiento al analizado en clase.
Se puede apreciar como, tomando [math]r=0[/math], ambas fronteras coinciden en [math]u(0)=0[/math]. Esto tiene sentido lógico con la teoría, pues estamos estudiando la región en la que se encuentran todas las funciones armónicas que pasan por [math]u(0)=0[/math]. Sin embargo, en la segunda gráfica se puede apreciar como la función que ejerce de frontera inferior de la región toma valores negativos. Esto se debe a que si recordamos la expresión de dicha función:
Esta comienza valiendo [math] u(0)=0 [/math] cuando [math] |x|=0 [/math] y a medida que el módulo de [math] |x| [/math] aumenta, la función disminuye, es de hecho estrictamente decreciente, tomando así siempre valores negativos. Finalmente, cuando [math] |x| =1[/math], alcanza el mínimo [math] M=-\frac{1}{2} [/math].
Para apreciar la región mejor empleamos escala logarítmica, pero para ello, tendremos que desplazar las funciones frontera para que tomen valores positivos. Tal y como hemos explicado en el párrafo anterior, el mínimo de la región se alcanza en [math]M=-\frac{1}{2}[/math], por lo que bastará con desplazar la región dicho valor para así no obtener valores negativos y que el logaritmo quede bien definido. Obtenemos la siguiente gráfica:
Ahora, vamos a comparar las gráficas para distintos dominios: [math] \mathbb{B}_{2}(0)[/math], [math] \mathbb{B}_{10}(0)[/math]. En este caso, nuestra solución queda como [math]u(r,\theta)=r^{2}cos(\theta)sin(\theta)[/math].
Para [math] \mathbb{B}_{2}(0) [/math] con [math] r=2[/math], el mínimo de la función es [math]u\left(\frac{3\pi}{4}\right)=u\left(\frac{7\pi}{4}\right)=-2 [/math]. De manera análoga al caso anterior, obtenemos la siguiente desigualdad:
Obtenemos con ello las siguientes representaciones gráficas:
Para la bola [math] \mathbb{B}_{10}(0) [/math] con [math] r=10[/math], el mínimo de la función es [math]u\left(\frac{3\pi}{4}\right)=u\left(\frac{7\pi}{4}\right)=-50 [/math]. Por lo tanto, en este caso obtenemos la siguiente desigualdad:
Obtenemos con ello las siguientes gráficas como resultado:
Prestando atención a las soluciones obtenidas se puede observar como la curva descrita por la frontera límite superior es más pronunciada a medida que aumentamos los radios. Además, hay que destacar que a simple vista puede parecer que restringiéndonos al dominio [math] x \in [0,1] [/math] las regiones definidas para las bolas de radio mayor estén incluidas en las de radio menor. Sin embargo, esto no es así. Para ello, hemos representado la comparación de dichas regiones en una misma gráfica restringiéndonos al dominio [math] x \in [0,1] [/math]. En ella se observan perfectamente estas diferencias
Por último, vamos a comparar con las cotas que se obtendrían en dimensión [math] n=3 [/math]. La región de la solución [math] w(x) [/math] en [math] \mathbb{B}_{R}(0)[/math] es la siguiente:
En este caso [math] M [/math] es el mínimo de nuestra función frontera [math]g(x,y)=x \cdot y[/math] evaluada sobre los puntos [math] (x_{0},y_{0},z_{0}) [/math] de la esfera [math] S_{1}(0) [/math].
Como la función no depende de la variable [math] z [/math], calcular el mínimo es equivalente a calcularlo sobre el plano [math] (x_{0}, y_{0},0) [/math]. Este plano representa la proyección de la esfera en el espacio tridimensional con respecto al plano [math] z=0 [/math], el cual equivale al disco unidad de dimensión [math] n=2 [/math]. Por consiguiente, el mínimo en la esfera es idéntico al calculado en el disco previamente. Por tanto, tenemos:
Para el caso [math] R=1 [/math]
Tomando [math] R=2 [/math] obtenemos la siguiente desigualdad
Para el caso [math] R=10 [/math] obtenemos la siguiente desigualdad:
Comparando los resultados obtenidos para dimensión [math] n=3 [/math] obtenemos las mismas conclusiones que para [math] n=2 [/math] dado que el comportamiento es análogo. La principal diferencia con respecto a la dimensión anterior se produce en los valores que se alcanzan en el eje [math] y [/math] de definición. De nuevo, si nos restringimos al dominio [math] x \in [0,1] [/math] podemos observar como las regiones definidas para las bolas de radio mayor no están incluidas en las de radio menor. La comparación de dichas regiones es la siguiente:
Analizando ahora las regiones obtenidas para distintas dimensiones, pero mismo radio, en este caso, si se puede observar como las regiones de definición de las funciones armónicas en dimensión [math] n=2 [/math] están incluidas en las de dimensión [math] n=3 [/math]. Realicemos una comparativa de dichas regiones en escala logarítmica para los tres radios considerados:
Hay que destacar que estas comparaciones pueden realizarse en escala logarítmica dado que se alcanzan los mismos mínimos en las regiones, y por tanto son desplazadas por la misma constante [math] M [/math] para conseguir así una correcta definición a la hora de aplicar logaritmos.
A continuación, se incluyen los códigos empleados durante esta sección. Por una parte, el primer código es el empleado para representar las regiones con los distintos radios y dimensiones y, el segundo código para las comparaciones entre estas. Estos son los siguientes:
clear all
close all
% Definimos la función u(x,y), solución exacta
u=@(x,y) x*y;
u_0=u(0,0);
% Definimos los radios de las bolas en las que vamos a analizar la desigualdad de Harnack
Rs=[1,2,10];
% Definimos los valores de n para analizar la desigualdad para distintas dimensiones
enes=[2,3];
for j=1:length(enes)
n=enes(j);
for i=1:length(Rs)
R=Rs(i);
% Definimos la función g(x,y)=x*y en coordenadas polares, G(theta)
G=@(theta) R^2*cos(theta)*sin(theta);
% Calculamos el mínimo de dicha función en la bola de radio R
minimo=fminbnd(G,0,2*pi);
M=G(minimo);
% Definimos la discretización de r=[0,R)
erres=linspace(0,R,1000);
erres=erres(1:end-1);
% Definimos las dos funciones cotas
w1=@(r) (u_0-M)*R^(n-2)*(R-r)./((R+r).^(n-1))+M;
w2=@(r) (u_0-M)*R^(n-2)*(R+r)./((R-r).^(n-1))+M;
% Representamos gráficamente la región sin escala logarítmica
figure((n-2)*2*length(Rs)+(2*i-1))
sgtitle('Representación de la región')
subplot(1,2,1)
w1_r=w1(erres); % Evaluamos las funciones w1 y w2 en los valores de r obtenidos a partir de la discretización
w2_r=w2(erres);
plot(erres,w1_r,'b','linewidth',2)
hold on
plot(erres,w2_r,'r','LineWidth',2)
xlabel('r')
fill([erres,erres(length(erres):-1:1)],[w1_r,w2_r(length(erres):-1:1)],[0.8 0.7 0.8]) % Rellenamos el área entre las dos funciones
title("r=[0,"+num2str(erres(end))+"]")
subplot(1,2,2)
plot(erres,w1_r,'b','linewidth',2)
xlabel('r')
hold on
plot(erres,w2_r,'r','LineWidth',2)
fill([erres,erres(length(erres):-1:1)],[w1_r,w2_r(length(erres):-1:1)],[0.8 0.7 0.8]) % Rellenamos el área entre las dos funciones
title("r=[0,"+num2str(R-0.1)+"]")
xlim([0,R-0.1])
% Representamos la región desplazada con escala logarítmica
figure((n-2)*2*length(Rs)+2*i)
w1_rM=w1_r-M; % Desplazamos las funciones w1 y w2 para que ambas sean positivo y aplicamos escala logarítmica
w2_rM=w2_r-M;
semilogy(erres,w1_rM,'r','linewidth',1.5)
hold on
semilogy(erres,w2_rM,'b','LineWidth',1.5)
fill([erres,erres(length(erres):-1:1)],[w1_rM,w2_rM(length(erres):-1:1)],[0.8 0.7 0.8])
sgtitle('Representación de la región desplazada en escala logarítmica')
end
end
Con respecto a este primer código empleado, lo primero que hacemos es definir los datos necesarios: la solución en la frontera de la bola ([math] g(x,y) [/math]) un vector de radios ([math] R=1,2,10 [/math] en nuestro caso) y otro de dimensiones ([math] n=2,3 [/math]). Además, denotamos como [math]u_0[/math] al punto [math]u(0,0)[/math].
Posteriormente empleamos coordenadas polares en la solución en la frontera y calculamos los valores de [math] M [/math] para cada uno de los radios. A continuación, definimos la desigualdad de Harnack en función de los radios y dimensiones y representamos las regiones obtenidas, incluyendo también para cada una de ellas la región en escala logarítmica.
clear all
close all
% Definimos la función u(x,y), solución exacta
u=@(x,y) x*y;
u_0=u(0,0);
% Definimos los radios de las bolas en las que vamos a analizar la desigualdad de Harnack
Rs=[1,2,10];
% Definimos los valores de n para analizar la desigualdad para distintas dimensiones
enes=[2,3];
% Definimos los colores
colores=['b','r','g'];
for j=1:length(enes)
n=enes(j);
for i=1:length(Rs)
R=Rs(i);
% Definimos la función g(x,y)=x*y en coordenadas polares, G(theta)
G=@(theta) R^2*cos(theta)*sin(theta);
% Calculamos el mínimo de dicha función en la bola de radio R
minimo=fminbnd(G,0,2*pi);
M=G(minimo);
% Definimos la discretización de r=[0,R)
erres=linspace(0,R,1000);
erres=erres(1:end-1);
% Definimos las dos funciones cotas
w1=@(r) (u_0-M)*R^(n-2)*(R-r)./((R+r).^(n-1))+M;
w2=@(r) (u_0-M)*R^(n-2)*(R+r)./((R-r).^(n-1))+M;
w1_r=w1(erres); % Evaluamos las funciones w1 y w2 en los valores de r obtenidos a partir de la discretización
w2_r=w2(erres);
w1_rM=w1_r-M; % Desplazamos las funciones w1 y w2 para que ambas sean positivo y aplicamos escala logarítmica
w2_rM=w2_r-M;
% Representamos la región desplazada con escala logarítmica
figure(i)
semilogy(erres,w1_rM,colores(j),'linewidth',1.5)
hold on
semilogy(erres,w2_rM,colores(j),'LineWidth',1.5)
fill([erres,erres(length(erres):-1:1)],[w1_rM,w2_rM(length(erres):-1:1)],colores(j))
xlabel('r')
sgtitle("Representación de la región para R= "+num2str(R)+" en dimensiones n=2 y n=3")
legend('','','n=2','','','n=3')
alpha(0.5)
% Representamos la región
figure(n+2)
plot(erres,w1_r,colores(i),'linewidth',1.5)
hold on
plot(erres,w2_r,colores(i),'LineWidth',1.5)
fill([erres,erres(length(erres):-1:1)],[w1_r,w2_r(length(erres):-1:1)],colores(i))
xlabel('r')
sgtitle("Representación de la región para n="+num2str(n)+" para los radios R=1,2,10")
alpha(0.5)
xlim([0,0.9])
legend('','','R=1','','','R=2','','','R=10')
end
end
Este último código tiene un funcionamiento similar al citado anteriormente. En primer lugar, definimos los datos necesarios para el problema: la solución exacta, los radios de las bolas y las dimensiones deseadas. Posteriormente, aplicamos coordenadas polares a la solución exacta, calculamos el mínimo de dicha función (definido previamente como [math] M [/math]) y definimos los parámetros necesarios para la discretización requerida. A continuación, definimos las funciones cota de nuestra desigualdad de Harnack, las evaluamos en nuestra discretización y las trasladamos dicho coeficiente [math] M [/math] calculado para así poder realizar posteriormente una escala logarítmica. Por último, representamos las regiones comparativas.
4 Ecuación de Poisson en [math] \mathbb{R}^2 [/math]
En esta sección nos centraremos en el estudio de la solución de la ecuación de Poisson. Para ello, planteamos el siguiente problema:
Donde [math] f [/math] es la función característica de la bola de radio [math] 1 [/math]. Observamos como la función [math] f \in \mathcal{C}^2(\mathbb{R}^3) [/math] que hemos definido tiene soporte compacto.
Conocemos que si [math] f \in \mathcal{C}^2(\mathbb{R}^3) [/math] es una función con soporte compacto, entonces la única solución del problema:
Viene dada por el potencial Newtoniano:
La demostración de esta prueba puede encontrarse como el Teorema [math] 3.33 [/math] en la sección [math] 3.6.2 [/math] de la página [math] 149 [/math] del libro Partial Differential Equations in Action. Springer. Sandro Salsa.
Sin embargo, en dimensión [math] n=2 [/math] existe una versión de dicho teorema sustituyendo el potencial newtoniano por el potencial logarítmico. De esta manera:
Donde en ambos casos, la aplicación [math] \phi(x) [/math] se define como la solución fundamental del laplaciano en dimensiones [math] 3 [/math] y [math] 2 [/math] respectivamente. En dimensión [math] n=3 [/math]: [math] \phi(x)=\frac{1}{4 \pi |x|} [/math], y en dimensión [math] n=2 [/math] : [math] \phi(x)=\frac{-1}{2 \pi} \cdot log |x|[/math].
Estas soluciones verifican que son singulares en [math] x=0 [/math] y [math] \Delta \phi = 0 [/math] fuera de [math] x=0 [/math].
Como la función f es la función característica de la bola de radio [math] 1 [/math] podemos escribir [math] u(x)= \int_{\mathbb{R}^2} \frac{-1}{2 \pi} \cdot log|x-y| \cdot f(y) \cdot dy =\int_{\mathbb{B}_{1}(0)} \frac{-1}{2 \pi} \cdot log|x-y| \cdot dy [/math].
Ahora, vamos a calcular en Matlab como sería la representación de la función [math] u(x) [/math] para distintos dominios:
clear all
close all
% Definimos la discretización de x1 y x2
num_puntos=200;
x1s=linspace(-10^10,10^10,num_puntos); % Esta discretización la hacemos para estudiar el comportamiento en el infinito
x2s=linspace(-10^10,10^10,num_puntos);
y1s=linspace(-10,10,num_puntos); % Esta discretización la hacemos para estudiar con más precisión el comportamiento en un dominio más cercano al (x1,x2)=(0,0)
y2s=linspace(-10,10,num_puntos);
% Definimos la discretización para r y theta
erres=linspace(0,1,300);
thetas=linspace(0,2*pi,300);
% Inicializamos las matrices en las que almacenaremos las soluciones para ambas discretizaciones
matrizx=zeros(num_puntos,num_puntos);
matrizy=zeros(num_puntos,num_puntos);
for w=1:num_puntos
x1=x1s(w);
y1=y1s(w);
for z=1:num_puntos
x2=x2s(z);
y2=y2s(z);
integralx_r=[];
integraly_r=[];
for i=1:length(thetas)
theta=thetas(i);
intx_r=trapz(erres,erres.*log(sqrt((x1-erres.*cos(theta)).^2+(x2-erres.*sin(theta)).^2))); % Calculamos la integral respecto de r
inty_r=trapz(erres,erres.*log(sqrt((y1-erres.*cos(theta)).^2+(y2-erres.*sin(theta)).^2))); % Calculamos la integral respecto de r
integralx_r=[integralx_r,intx_r];
integraly_r=[integraly_r,inty_r];
end
matrizx(w,z)=-1/(2*pi)*trapz(thetas,integralx_r); % Añadimos la solución a la matriz haciendo la integral respecto de theta
matrizy(w,z)=-1/(2*pi)*trapz(thetas,integraly_r); % Añadimos la solución a la matriz haciendo la integral respecto de theta
end
end
% Representación gráfica
figure(1)
subplot(1,2,1) % Representación para estudiar el comportamiento en el infinito
surf(x1s,x2s,matrizx,'EdgeColor','none')
xlabel('x_1')
ylabel('x_2')
title('x_1')
title("(x_1, x_2) \in ["+num2str(x1s(1))+","+num2str(x1s(end))+"] \times ["+num2str(x2s(1))+","+num2str(x2s(end))+"]")
subplot(1,2,2)
surf(y1s,y2s,matrizy,'EdgeColor','none') % Representación para estudiar el comportamiento en un dominio más cercano al (x1,x2)=(0,0)
xlabel('x_1')
ylabel('x_2')
title("(x_1, x_2) \in ["+num2str(y1s(1))+","+num2str(y1s(end))+"] \times ["+num2str(y2s(1))+","+num2str(y2s(end))+"]")
sgtitle('Representación de la solución en distintos dominios')
En el código, se discretizan las variables [math]x_{1}, x_{2},y_{1}, y_{2}[/math] para tomar dos conjuntos de puntos, uno espaciado para estudiar el comportamiento en [math] \mathbb{R}^2 [/math] y otro más denso para estudiarlo cerca del origen. Se hace el cambio a polares de la función a integrar y se calcula la solución numérica de esta utilizando el método del trapecio, iterando sobre los puntos discretizados y almacenando los resultados en matrices. Finalmente, se visualizan las soluciones en dos subgráficas utilizando la función “surf”.
En la primera imagen, la función resulta adecuada para estudiar el comportamiento de la solución en valores de [math] x [/math] muy alejados del origen. Esto se podría interpretar como su comportamiento en el “infinito”. Sin embargo, esta gráfica es inadecuada para estudiar el comportamiento con respecto al origen ya que no se representa bien el comportamiento alrededor de este. La inmensa anchura de la región hace que las discretizaciones en las que se toman los valores estén muy espaciadas. Además, la alta pendiente presente en la vecindad del origen requiere una mayor densidad de puntos para capturar su comportamiento. En la segunda imagen se toma una región más limitada cerca del origen, la cual si que facilita un estudio más preciso de su comportamiento.
Por otro lado, hay que destacar que el potencial logarítmico no ‘desaparece’ en el infinito, Su comportamiento asintótico, como indica la observación [math]3.35[/math] del Teorema [math] 3.33 [/math] del libro Partial Differential Equations in Action. Springer. Sandro Salsa ,viene dado por:
Donde definimos [math] M [/math] como:
que aplicado a nuestro problema inicialmente planteado con [math] f [/math] la función característica de la bola de radio [math] 1 [/math], tendríamos que:
Buscamos ahora comparar el comportamiento asintótico teórico con el obtenido numéricamente. En primer lugar, vamos a hacer una representación gráfica de ambas en el dominio [math] (x_1,x_2) \in [-100,100] \times [-100,100] [/math] y para ello, hemos elaborado el siguiente código en Matlab:
clear all
close all
% Definimos la función que describe su comportamiento asintótico
limite=@(x1,x2) -1/2*log(sqrt(x1.^2+x2.^2));
% Definimos la discretización de x1 y x2
num_puntos=400;
x1s=linspace(-1,1,num_puntos);
x2s=linspace(-1,1,num_puntos);
% Definimos la discretización para r y theta
erres=linspace(0,1,1000);
thetas=linspace(0,2*pi,1000);
% Inicializamos las matrices en las que almacenaremos la solución
matriz=zeros(length(x1s),length(x2s));
for w=1:num_puntos
x1=x1s(w);
for z=1:num_puntos
x2=x2s(z);
integral_r=[];
for i=1:length(thetas)
theta=thetas(i);
int_r=trapz(erres,erres.*log(sqrt((x1-erres.*cos(theta)).^2+(x2-erres.*sin(theta)).^2))); % Calculamos la integral respecto de r
integral_r=[integral_r,int_r];
end
integral_theta=trapz(thetas,integral_r); % Calculamos la integral respecto de theta
matriz(w,z)=-1/(2*pi)*integral_theta; % Añadimos la solución a la matriz
end
end
% Representación gráfica
figure(1)
subplot(1,2,1)
surf(x1s,x2s,matriz,'EdgeColor','none') % Representación de la solución
xlabel('x_1')
ylabel('x_2')
title("Representación de la solución")
[X1_meshgrid,X2_meshgrid]=meshgrid(x1s,x2s);
subplot(1,2,2)
surf(x1s,x2s,limite(X1_meshgrid,X2_meshgrid),'EdgeColor','none') % Representación de la función que describe el comportamiento asintótico
xlabel('x_1')
ylabel('x_2')
title("Representación de la función que describe el comportamiento asintótico")
sgtitle("Comparación de la solucion con la función que describe el comportamiento asintótico para (x_1, x_2) \in ["+num2str(x1s(1))+","+num2str(x1s(end))+"] \times ["+num2str(x2s(1))+","+num2str(x2s(end))+"]")Primero de todo, en el código se define la función asintótica. A continuación, se discretizan las variables [math](x_{1},x_{2})[/math] para su representación y se establece otra discretización para las variables [math]r[/math] y [math]\theta[/math]. Seguidamente, mediante un bucle, se calculan las integrales respecto de [math]r[/math] y respecto de [math]\theta[/math] utilizando el método del trapecio. Estos resultados se almacenan en una matriz y, por último, se representan gráficamente la solución numérica y el comportamiento asintótico descrito teóricamente.
En las gráficas obtenidas no se pueden apreciar bien las diferencias entre ambas funciones. Para poder apreciar la magnitud de esta diferencia, hemos ejecutado el código anterior pero con un dominio más cercano al punto [math] (x_1,x_2)=(0,0) [/math]. En particular, para el dominio [math] (x_1,x_2) \in [-1,1] \times [-1,1] [/math] las gráficas obtenidas son las siguientes:
De esta manera, al restringir el dominio de representación, si se aprecian con claridad las diferencias entre ambas. En concreto, se observa una singularidad en el punto [math] (x_1,x_2)=(0,0) [/math] para la función del comportamiento asintótico [math] u(x) = -\frac{1}{2} \cdot log|x| [/math], mientras que la solución calculada de manera numérica no presenta dicha singularidad. Esto último contradice lo citado anteriormente ya que mencionamos que ambas funciones eran singulares en [math] x=0 [/math]. Sin embargo, esta contradicción es consecuencia de realizar los cálculos de manera numérica.
Para poder estudiar el comportamiento asintótico de la solución calculada por el método del trapecio, hemos tenido en cuenta que ésta es simétrica con respecto al eje [math] x_1[/math] y, por lo tanto, basta con estudiar lo que sucede para el semieje positivo. Para ello, hemos elaborado otro código en Matlab:
clear all
close all
% Definimos la función que describe su comportamiento asintótico
limite=@(x1,x2) -1/2*log(sqrt(x1.^2+x2.^2));
% Definimos la discretización de x1 y x2
num_puntos=200;
x1s=[0.001];
for h=1:num_puntos-1
x1s=[x1s,x1s(end)*1.1]; % Hacemos una discretización no equidistribuida para x1s
end
%x1s=linspace(0.0001,10^10,num_puntos)
x2s=linspace(-10,10,num_puntos-1);
x2s=[x2s,0]; % Añadimos el valor x2=0, pues es donde más diferencia se va a producir entre ambas funciones
x2s=sort(x2s);
% Definimos la discretización para r y theta
erres=linspace(0,1,300);
thetas=linspace(0,2*pi,300);
% Inicializamos las matrices en las que almacenaremos la solución
matriz=zeros(length(x1s),length(x2s));
% Crear un vídeo
pelicula=VideoWriter('Ej5','MPEG-4'); % creo el video
pelicula.FrameRate=10; % controla velocidad
open(pelicula); % abrimos el vídeo
error=[];
for w=1:num_puntos
x1=x1s(w);
figura=figure(1);
for z=1:num_puntos
x2=x2s(z);
integral_r=[];
for i=1:length(thetas)
theta=thetas(i);
int_r=trapz(erres,erres.*log(sqrt((x1-erres.*cos(theta)).^2+(x2-erres.*sin(theta)).^2))); % Calculamos la integral respecto de r
integral_r=[integral_r,int_r];
end
integral_theta=trapz(thetas,integral_r); % Calculamos la integral respecto de theta
matriz(w,z)=-1/(2*pi)*integral_theta;
end
error=[error,max(abs(matriz(w,:)-limite(x1,x2s)))]; % Calculamos el máximo del error
plot(x2s,limite(x1,x2s),'r','LineWidth',1.5)
hold on
plot(x2s,matriz(w,:),'b--','LineWidth',1.5)
hold off
ylim([-25,5])
xlabel('Valores de x_2')
title("Representación de ambas funciones para x_1="+num2str(x1))
legend('Función que describe el comportamiento asintótico','Solución usando método del trapecio','Location','southeast')
% Añadir el frame al vídeo
imagen=getframe(figura);
writeVideo(pelicula,imagen); % inserto la imagen
end
% Representación gráfica del error con respecto a los valores de x1
figure(2)
sgtitle('Representación del error máximo para distintos valores de x_1')
semilogx(x1s,error,'bo-','LineWidth',1.5)
xlim([10^(-3),10^5])
xlabel('Valores de x_1 en escala logarítmica')
ylabel('error')
% Cerrar el vídeo
close(pelicula);
Este código representa la solución calculada usando el método del trapecio y la función [math] u(x) = \frac{-1}{2} log|x| [/math] para [math] x_2 \in [-10,10] [/math] y valores fijos de [math] x_1 [/math] entre [math] 0.001 [/math] y [math] 1.72 \cdot 10^5 [/math].
Con respecto al código, en primer lugar definimos la función [math] u(x) [/math], a la que hemos denominado limite . A continuación, hemos discretizado las variables [math] x_1 [/math] y [math] x_2 [/math], así como las variables [math] r [/math] y [math] \theta [/math], recordando que para resolver las integrales correspondientes para calcular la solución hemos hecho uso de coordenadas polares. Cabe destacar que puesto que para valores de [math]x_{1}[/math] cercanos a [math] 0 [/math] la solución sufre cambios mucho más notorios, hemos llevado a cabo una discretización de esta variable para que tenga muchos puntos cercanos a dicho valor y a medida que aumenta el mismo, estén más espaciados. Además, hemos optado por el método del trapecio para calcular numéricamente dichas integrales, la primera con respecto a la variable [math] r [/math] y la segunda con respecto a la variable [math] \theta [/math]. Finalmente, hemos graficado ambas funciones para cada uno de los valores de [math] x_1 [/math] y hemos añadido cada una de esas imagenes a un vídeo, que es el siguiente:
Al comienzo de este vídeo, cuando el valor de [math] x_1 [/math] es próximo a [math] 0 [/math] , se puede apreciar como ambas funciones sí presentan diferencias en un entorno del punto [math] x_2=0 [/math] (para cada valor fijado de [math] x_1 [/math]). En concreto, la solución obtenida al integrar de manera numérica no presenta la singularidad que sí presenta la función [math] u(x) [/math]. Sin embargo, el objetivo de este vídeo es comprobar que la función [math] u(x) [/math] describe el comportamiento asintótico de la solución y, efectivamente, a medida que crece [math]x_1 [/math], la diferencia entre estas funciones disminuye. En particular, a partir del valor [math] x_1=0.86 [/math], ambas funciones se solapan y decrecen hacia [math] - \infty [/math] a la misma velocidad.
Para comprobar este hecho recién explicado, el código también devuelve la siguiente gráfica:
Es imprescindible mencionar que debido a la discretización llevada a cabo para la variable [math] x_1 [/math], hemos considerado más oportuno poner el eje [math]x[/math] de esta gráfica en escala logarítmica, pues de la otra manera la mayoría de los puntos se almacenan en un intervalo muy pequeño y es muy complicado extraer conclusiones.
Con respecto a esta gráfica, es muy llamativo como el error inicial, para [math] x_1=0.001 [/math] es realmente elevado, tomando el valor de [math] 3.2 [/math]. Aunque este valor pueda parecer muy alto, realmente concuerda con los resultados analizados en el vídeo, pues está directamente relacionado con la singularidad que se produce en la función [math] u(x) [/math]. Además, esta última gráfica confirma que a partir de [math] x_1=0.86 [/math], la diferencia entre ambas funciones es realmente pequeña y, por lo tanto, concluimos que efectivamente [math] u(x) [/math] describe el comportamiento asintótico de la solución.
5 Aplicaciones
Hasta el momento, se ha realizado un estudio exhaustivo de la ecuación de Laplace y la ecuación de Poisson, abordando su comportamiento, sus soluciones y sus aproximaciones tanto en [math] \mathbb{R}^{2}[/math] como en [math]\mathbb{R}^{3}[/math]. No obstante, utilizando estas ecuaciones, los científicos e ingenieros pueden modelar y resolver una variedad de problemas prácticos con precisión y eficacia.
- Electrostática y Magnetostática:
En el contexto de la electrostática, la ecuación de Laplace y la ecuación de Poisson se utilizan para modelar y resolver problemas relacionados con el campo y el potencial eléctricos en sistemas sin cargas y con cargas distribuidas, respectivamente. Estas ecuaciones son fundamentales para comprender el comportamiento de circuitos eléctricos, condensadores, líneas de transmisión y muchos otros dispositivos eléctricos.
- Transferencia de Calor:
En problemas de transferencia de calor, la ecuación de Laplace y la ecuación de Poisson se aplican para modelar la distribución de temperatura en sistemas estacionarios y transitorios. Estas ecuaciones son utilizadas en el diseño de sistemas de refrigeración, análisis de aislamiento térmico, y optimización de procesos de fabricación que involucran calor.
- Mecánica de Fluidos y Aerodinámica:
En la mecánica de fluidos, la ecuación de Laplace y la ecuación de Poisson se utilizan para modelar y resolver problemas relacionados con el potencial de flujo. Estas ecuaciones son esenciales en el diseño de cuerpos sumergidos y sistemas de tuberías. En aerodinámica, estas ecuaciones se aplican para calcular el potencial de velocidad y la distribución de presión alrededor de objetos en movimiento, como aviones, automóviles y cuerpos aerodinámicos.
6 Referencias
- Wikipedia. Regla del trapecio
- [Partial differential equations in action from modelling to theory. Sandro Salsa]
- de Laplace | Ecuación. Usos y Ejemplos
