Ecuación de Laplace y de Poisson (Raúl, Sofía, Jaime)
| Trabajo realizado por estudiantes | |
|---|---|
| Título | Ecuación del calor |
| Asignatura | EDP |
| Curso | 2023-24 |
| Autores | Raúl Ortega
Sofía Gómez Jaime Sáenz de Miera |
| Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura | |
En este artículo vamos a estudiar la ecuación de Laplace y la ecuación de Poisson. Para ello vamos a calcular soluciones de distintas formas usando cosas como la fórmula de Poisson o el potencial logarítmico, graficándolas y modificando algunos de los parámetros.
Contenido
1 Ecuación de Laplace
En este apartado vamos a estudiar la exactitud que tienen dos formas diferentes de resolver la ecuación de Laplace. Para hacer esto tomamos el siguiente problema:
[math] \begin{cases} \Delta u=0,& x\in B_1\\ u=g,& x\in\partial B_1 \end{cases} [/math]
Con [math] B_1\subset\mathbb{R}^2 [/math] la bola de radio 1 centrada en el (0,0).
1.1 Fórmula de Poisson
Una de las formas de resolver este tipo de problemas es mediante la fórmula de Poisson, que es: [math] u=\frac{R^2-|x|^2}{2\pi R}\int_{\partial B_R}\frac{g(\sigma)}{|x-\sigma|^2} d\sigma, [/math] donde R es el radio de la bola y [math] |x| [/math] la distancia al centro de la bola.
Para ver cómo funciona esta fórmula vamos a resolver el problema tomando [math] g(\theta)=max \{0,1-\frac{2}{\pi}|\theta-\pi|\}[/math]. Como esta función, que es la que representa las condiciones frontera, está en coordenadas polares vamos a expresar la fórmula de Poisson en polares.
Una vez tenemos la fórmula expresada de esta manera, podemos aplicarla a nuestro problema tomando R=1 y sustituyendo g por nuestra función. Teniendo esto en cuenta vamos a resolver el problema realizando la integral mediante la regla del trapecio y a representar la solución.
% Calcular la solución de la ecuación de Laplace en la bola de radio
% unidad usando la fórmula de Poisson. Con condición frontera la función g
g = @(theta) max(0,1-2/pi*abs(theta-pi)); %Definimos la función g(theta)
N = 3*10^2; %Número de particiones en theta y r
h = 2*pi/N;
theta = 0:h:2*pi; %Partición en theta
r = 0:N^(-1):0.9; %Partición en r (no llegamos hasta el borde de la bola)
w = ones(N+1,1);
w(1) = 1/2; w(N+1) = 1/2;
u = zeros(length(r)+1,length(theta));
for j = 1:length(theta)
for i = 1:length(r)
f = g(theta)*(1-r(i).^2)./(1+r(i).^2-2*r(i)*cos(theta-theta(j))); %Función a integrar
u(i,j) = h*w'*f'/(2*pi); %Fórmula de Poison usando la fórmula del trapecio para integrar.
end
u(length(r)+1,j) = g(theta(j)); %Añadimos la condicón frontera directamente
end
r(length(r)+1) = 1;
[Theta, R] = meshgrid(theta, r); % Crear una malla de theta y r
% Convertir coordenadas polares a cartesianas
X = R.*cos(Theta);
Y = R.*sin(Theta);
% Graficar la superficie en coordenadas cartesianas
colormap jet
surf(X, Y, u, 'EdgeColor', 'none');
xlabel('x');
ylabel('y');
zlabel('u(x, y)');
title('Solución de la ecuación de Laplace');
En la imagen podemos observar que la fórmula de Poisson nos da problemas en la frontera, esto se debe al carácter singular de la integral cerca del borde. Para solucionar esto vamos a tomar directamente el valor de g en el borde de la bola. Al hacer esto obtenemos la siguiente gráfica.
Al comparar ambas gráficas podemos observar que efectivamente sustituir g en el borde hace que evitemos la singularidad de la integral en la frontera.
Como hemos visto con este ejemplo la fórmula de Poisson tiene problemas para aproximar la solución cerca de la frontera. Por lo que para estudiar esto, vamos a tomar [math] g(x,y)=xy [/math] como condición frontera y [math] u(x,y)=xy [/math] como solución del problema para calcular los diferentes errores que genera la fórmula de Poisson.
Para hacerlo vamos a tomar tanto u como g en coordenadas polares, ya que al estudiar el problema en una bola esto nos permite hacernos una mejor idea de lo que ocurre.
[math] \begin{cases} U(r,\theta)=r^2cos(\theta)sin(\theta)\\ U(1,\theta)=G(\theta)=cos(\theta)sin(\theta) \end{cases} [/math]
Primero vamos a tomar un punto "lejos" de la frontera [math] (r,\theta)=(0.9,\frac{\pi}{4}) [/math] y a ver cómo evoluciona el error, es decir, la diferencia entre el valor real de u y el valor de la aproximación por la fórmula de Poisson en ese punto, en función de n. Siendo n el número que determina el número de intervalos en [math] 10^n[/math].
Por el ejemplo anterior, cabe esperar que la fórmula en este punto aproxime bien a la función ya que no es un punto cercano a la frontera, y por lo tanto el error vaya disminuyendo según aumenta el número de intervalos.
% Calculamos el error que incorpora la aproximación de la integral mediante
% el método del trapecio en función del tamaño de la discretización
g = @(theta) sin(theta).*cos(theta); %Definimos la función g(theta)
N_max = 8; %Discretización máxima: 10^N_max
err = zeros(N_max,1);
for m = 1:N_max
N = 10^m; %Tamaño de la discretización
h = 2*pi/N;
theta = 0:h:2*pi; %Partición en theta
w = ones(N+1,1);
w(1) = 1/2; w(N+1) = 1/2;
f = g(theta)*(1-0.9.^2)./(1+0.9.^2-2*0.9*cos(theta-pi/4)); %Función a integrar
u = h*w'*f'/(2*pi); %Fórmula de Poisson usando la fórmula del trapecio en el punto (0.9,pi/4)
% Valor de la solución exacta (x*y=R^2*cos(theta)*sin(theta))
% en el punto (0.9, pi/4)
exacta = 0.9^2*cos(pi/4)*sin(pi/4);
% Calculamos el error como la diferencia en valor absoluto entre
% la solución exacta y la calculada con la fórmula de Poisson
err(m) = abs(exacta-u);
end
plot(1:N_max,log10(err),color='b',LineWidth=1.5) %Dibujamos la gráfica del error en escala logarítmica
xlabel('n');
ylabel('log_{10}(error(10^n))');
Como hemos comentado antes, al observar la gráfica vemos que el error va disminuyendo según aumenta n. Esto tiene sentido ya que cuantos más intervalos hay en la discretización más exacta va a ser la aproximación. También podemos observar que se estabiliza cuando llega aproximadamente a [math] n=3 [/math], esto se debe a que con [math] 10^3 [/math] intervalos el método del trapecio ya aproxima lo suficientemente bien a la función original, lo que hace que aunque aumente el número de intervalos la aproximación no tenga mucho margen de mejora. Con esto hemos visto que en un punto que no esté cerca del borde, la fórmula de Poisson va a funcionar bien.
También podemos usar la fórmula del error máximo del método del trapecio para obtener una cota para el error y ver que los errores calculados numéricamente se encuentran por debajo.
Donde [math] [b-a] [/math] es el intervalo de integración, [math]f[/math] la función a integrar, [math]n[/math] el número de particiones de la discretización y [math] \xi [/math] un cierto valor entre [math] [b-a] [/math]. Para nuetro problema, [math][a,b]=[0,2\pi][/math] y [math]f(\theta)=\frac{(1-r^2) sen(\theta)cos(\theta)}{2\pi(1+r^2-2rcos(\alpha - \theta))}[/math].
Derivando [math]f[/math] dos veces y maximizando su valor en [math][0,2\pi][/math], obtenemos la siguiente cota:
La dibujamos en función de [math]n[/math] junto con el error calculado numéricamente:
% Usando el vector err calculado en el anterior código, representamos la misma gráfica junto
% con la gráfica de la cota estimada del error
%Definimos la cota del error estimada mediante la fórmula del trapecio
cota = @(n) abs(-pi^2*(r^6-r^5+6*r^3+r^2-r+2)./(3*(10.^(2*n))*(1-r)^5));
%Dibujamos las gráficas del error en escala logarítmica
hold on
plot(1:N_max,log10(err(:,1)),color='b',LineWidth=1.5)
plot(1:N_max,log10(cota(1:N_max)),color='r',LineWidth=1.5)
legend('Error absoluto', 'Cota del error'); % Leyenda
hold off
Comprobamos que el error de la aproximación queda por debajo de esta cota (que no es muy ajustada).
Ahora vamos a ver cómo se comporta el error según nos acercamos a la frontera. Para hacer esto vamos a tomar el punto [math] (r,\theta)=(1-10^{-n},\frac{\pi}{4}) [/math] y a ver cómo varía el error (la diferencia entre el valor real de u y el valor de la aproximación por la fórmula de Poisson) cuando aumentamos la n, es decir, cuando nos vamos acercando a la frontera. Como hemos comentado anteriormente cabe esperar que al aumentar la n el error va a aumentar, ya que como hemos visto en el primer ejemplo cerca de la frontera la fórmula de Poisson da problemas.
% Calculamos el error que incorpora la aproximación de la integral según nos acercamos al
% borde de la región.
g = @(theta) sin(theta).*cos(theta); %Definimos la función g(theta)
N_max = 20;
err = zeros(N_max);
for n = 1:N_max
N = 10^2; %Tamaño de la discretización
h = 2*pi/N;
theta = 0:h:2*pi; %Partición en theta
w = ones(N+1,1);
w(1) = 1/2; w(N+1) = 1/2;
f=g(theta)*(1-(1-10^(-n)).^2)./(1+(1-10^(-n)).^2-2*(1-10^(-n))*cos(theta-pi/4)); %Función a integrar
u=h*w'*f'/(2*pi); %Fórmula de Poisson usando la fórmula del trapecio en el punto (1-10^-n,pi/4)
err(n)=abs((1-10^(-n))^2*cos(pi/4)*sin(pi/4)-u);
end
plot(1:N_max,err,color='b',LineWidth=1.5) %Dibujamos la gráfica del error
xlabel('n');
ylabel('Error');
En la gráfica podemos ver cómo efectivamente para valores de n pequeños, es decir, puntos que están relativamente lejos de la frontera, el error es pequeño al igual que en el caso anterior. Por el contrario, según aumenta la n, es decir, según nos acercamos a la frontera el error aumenta, lo cual tiene sentido ya que como hemos comentado antes la fórmula de Poisson da problemas en la frontera.
Además, podemos observar que el error se estabiliza aproximadamente en [math] n=4 [/math], esto se debe a que para este valor de n el punto ya esta suficientemente cerca de la frontera como para representar todo el error que se da en esta y por ello aunque nos acerquemos más a la frontera (aumente n) el error se mantiene estable en [math] 0.5[/math].
1.2 Serie de Fourier
La otra forma que vamos a ver para resolver la ecuación de Laplace es mediante series de Fourier. Para hacer esto vamos a tomar nuestro problema en coordenadas polares:
[math] \begin{cases} \Delta U=U_{rr}+\frac{1}{r}U_r+\frac{1}{r^2}U_{\theta\theta}=0,& r\in(0,1),\theta\in[0,2\pi]\\ U(1,\theta)=G(\theta)=sin(\theta)cos(\theta),& \theta\in[0,2\pi] \end{cases} [/math]
y a aplicar el método de separación de variables. Al hacer esto obtenemos que: [math] U(r,\theta)=\frac{a_0}{2}+\sum_{k=1}^\infty a_kr^kcos(k\theta)+b_kr^ksin(k\theta).[/math] Para que esto sea solución falta hacer que cumpla la condición frontera, para ello vamos a tomar la serie de Fourier de [math] G(\theta) [/math] en [math] L^2([-\pi,\pi]) [/math], que es [math] G(\theta)=\frac{\alpha_0}{2}+\sum_{k=1}^\infty \alpha_kcos(k\theta)+\beta_ksin(k\theta)[/math] y a igualar ambas series cuando r=1. Al hacer esto vemos que para que ambas series sean iguale los coeficientes tienen que ser iguales, por lo que vamos a calcular los coeficientes de [math] G(\theta) [/math] y a ver su equivalencia con los de [math] U [/math].
Al realizar estas integrales obtenemos que [math] a_0=0,a_k=0\forall k [/math] y todos los [math] b_k [/math] son 0 excepto [math] b_2=\frac{1}{2} [/math]. Por lo que la solución nos queda de la siguiente manera:
[math]U(r,\theta)=\frac{r^2}{2}sin(2\theta) [/math]
Ahora vamos a dibujar la solución, cuya gráfica va a ser exacta ya que la serie tiene un único término.
% Dibujamos la solución calculada por serie de Fourier
u = @(r,theta) r.^(2).*sin(2.*theta)/2; %Solución de la ecuación
N = 3*10^2; %Número de particiones en theta y r
h = 2*pi/N;
theta = 0:h:2*pi; %Partición en theta
r = 0:N^(-1):1; %Partición en r
for j = 1:N+1
for i = 1:N+1
U(i,j) = u(r(i),theta(j));
end
end
[Theta, R] = meshgrid(theta, r); % Crear una malla de theta y r
% Convertir coordenadas polares a cartesianas
X = R.*cos(Theta);
Y = R.*sin(Theta);
% Graficar la superficie en coordenadas cartesianas
colormap jet
surf(X, Y, U, 'EdgeColor', 'none');
xlabel('x');
ylabel('y');
zlabel('u(x, y)');
title('Solución exacta calculada por serie de Fourier');
Al haber resuelto el problema por ambos métodos, observamos que este último nos da la solución de forma exacta, por lo que en este caso va a ser más efectivo que el anterior. Esto no tendría porque ser así en todos los casos ya que la serie de Fourier para otra función G distinta puede tener infinitos términos, lo que haría que la aproximación no sea exacta, pero aun así sería más exacta que con el método anterior en los puntos cercanos a la forntera.
1.3 Desigualdad de Harnack
En esta sección vamos a usar la desigualdad de Harnack para obtener la región que contiene a todas las soluciones armónicas de la ecuación de Laplace que valen lo mismo en el origen que la solución que hemos calculado.
La desigualdad de Harnack establece que dada una función [math] u\geq 0 [/math] armónica en [math] \Omega \subset \mathbb{R}^n[/math] y [math] B_R(z) \subset \Omega [/math] con [math] z \in \Omega [/math], entonces:
Para poder aplicarla a nuestro problema, vamos a tomar v=U+M, siendo [math] U(r,\theta)=\frac{r^2}{2R^2}sin(2\theta) [/math] la solución al problema en [math] B_R [/math], y M el mínimo de [math] G(\theta)=R^2sin(\theta)cos(\theta) [/math] en [math] B_R [/math]. Al hacer esto v va a ser una función armónica positiva a la que le podemos aplicar la desigualdad de Harnack.
Primero vamos a aplicar la desigualdad variando el radio de la bola y a observar como varia la región que determina en función de esto. Para ello vamos a representar las cotas en escala logarítmica, lo cual nos permite observar mejor su comportamiento.
% Representamos las cotas que nos da la desigualdad de Harnack
% para nuestro ejemplo.
u = @(r,theta) r.^2*sin(theta).*cos(theta); %Definimos la función u(r,theta)
g=@(theta,r) r^2*sin(theta).*cos(theta); %Definimos la función g
N = 3*10^2; %Número de particiones en theta y r
R = [1,2,10]; %Radios de las bolas
n = 2; %Dimensión del espacio
h = 2*pi/N;
theta = 0:h:2*pi; %Partición en theta
C=["#0000FF","#FF0000","#00FF00"];
for i=1:3
r = 0:N^(-1):R(i); %Partición en r
gR=g(theta,R(i));
M1=min(gR); %Calculamos el mínimo de g en la frontera de la bola de radio R
v1=@(r,theta) u(r,theta)-M1;
A=@(r)R(i).^(n-2).*(R(i)-r)./(R(i)+r).^(n-1)*v1(0,0);
B=@(r)R(i).^(n-2).*(R(i)+r)./(R(i)-r).^(n-1)*v1(0,0);
%Graficamos las cotas para los diferentes radios
subplot(1,3,i)
plot(r,log10(A(r)),Color=C(i),LineWidth=1.5)
hold on
plot(r,log10(B(r)),Color=C(i),LineWidth=1.5)
axis([0 10,-3 6])
xlabel('r');
ylabel('log(Cotas)');
title(['Cotas para r= ' num2str(R(i))]);
end
Lo más destacable en las 3 gráficas es que en todas ellas ambas restricciones se van al infinito cunado llegan al radio correspondiente a la bola en la que están calculadas. Esto tiene sentido ya que podemos observar que el denominador de la cota superior se anula cuando r=R y por lo tanto se va a infinito y su logaritmo también, y la cota inferior es 0 cuando r=R y por tanto su logaritmo se va a menos infinito. También podemos observar, que debido a esto, el pico que forman ambas restricciones es más estrecho cuando aumentamos el radio de la bola.
Ahora vamos a ver cómo varían las restricciones si variamos la dimensión.
Al observar estas gráficas podemos observar que independientemente del radio cuando aumentamos la dimensión aumenta la región que definen ambas restricciones. Además, en n=3 el comportamiento sigue siendo el mismo que el comentado para n=2.
2 Ecuación de Poisson
Ahora vamos a buscar una solución de la ecuación de Poisson en [math] \mathbb{R}^2 [/math]:
[math] \Delta u=f,\quad x\in\mathbb{R}^2. [/math]
Para ello vamos a utilizar la solución fundamental del laplaciano, que en [math] \mathbb{R}^2 [/math] se define como [math] \Phi(x)=-\frac{1}{2\pi}log|x| [/math] y cumple que [math] \Delta \Phi(x)=-\delta_2(x) [/math]. Donde [math]\delta_2(x) [/math] es la delta de Dirac en dimensión 2.
Si ahora consideramos el potencial logarítmico de la función f, que es la la convolución de la solución fundamental con f, entonces podríamos expresar la solución como: [math] u(x)=\int_{\mathbb{R}^2}\Phi(x-y)f(y)dy [/math] y formalmente tenemos que:
[math] \Delta u(x)=\int_{\mathbb{R}^2}\Delta_x\Phi(x-y)f(y)dy=- \int_{\mathbb{R}^2}\delta_2f(y)dy=-f(x),[/math]
pero esto no es correcto ya que el laplaciano de la solución fundamental no es integrable. A pesar de no ser correcto nos da la intuición de que el potencial logarítmico podría ser una posible solución del problema.
Ahora para solucionar esto vamos a tomar f como la función característica de la bola de radio 1, que como es una función con soporte compacto podemos asegurar que una solución para la ecuación de Poisson es el potencial logarítmico. (Página 150 Partial Differential Equations in Action. S. Salsa, G. Verzini).
Por tanto, la solución vendrá dada por la siguiente expresión:
Haciendo el cambio de variable [math] \textbf{y}= ( r cos (\theta), r sin( \theta)) [/math] obtenemos:
Representamos la solución integrando mediante el método del trapecio.
N1=300; N2=300; %Número de puntos
a=0; b=1; c=0; d=2*pi; %Extremos de los intervalos
h1=(b-a)/N1; h2=(d-c)/N2;
r=a:h1:b; theta=c:h2:d; %Particiones de r y theta
[rr,tt]=meshgrid(r,theta);
w1=ones(N1+1,1);
w1(1)=1/2; w1(N1+1)=1/2;
w2=ones(N2+1,1);
w2(1)=1/2; w2(N2+1)=1/2;
x11=-1:2/100:1; % Partición de x_1
x21=-1:2/100:1; %Partición de x_2
for j=1:length(x11)
for i=1:length(x21)
f=-log(sqrt((x11(j)-rr.*cos(tt)).^2+(x21(i)-rr.*sin(tt)).^2)).*rr; %Potencial logaritmico
u1(i,j)=h1*h2*w2'*f*w1/(2*pi); %Solución del método del trapecio en 2 variable
end
end
%Graficamos la solución
colormap jet
surf(x11,x21, u1, 'EdgeColor', 'none');
xlabel('x_1');
ylabel('x_2');
zlabel('u(x_1, x_2)');
title('Solución con el potencial logarítmico');
El comportamiento asintótico del potencial logarítmico es:
Donde se define M como:
En nuestro problema, donde [math] f [/math] es la función característica de la bola unidad, tenemos:
%Representación de la solución asintótica
M=pi; %Area del circulo de radio 1
f=@(x1,x2) -M/(2*pi).*log(sqrt(x1.^2+x2.^2)); % Función que describe el comportamiento asintótico
N=200; % Tamaño de la partición
a=-100; b=100;
h=(b-a)/N;
x11=a:h:b; %Partición de x_1
x21=a:h:b; %Partición de x_2
u=zeros(length(x11),length(x21));
for j=1:length(x11)
for i=1:length(x21)
u(i,j)=f(x11(i),x21(j)); %Evaluamos f en la malla de puntos
end
end
%Graficamos el resultado
colormap jet
surf(x11, x21, u, 'EdgeColor', 'none');
xlabel('x_1');
ylabel('x_2');
zlabel('u(x_1, x_2)');
title('Solución cuando x tienda a \infty');
Si lo comparamos con el comportamiento asintótico de la solución calculada observamos que ambas gráficas se asemejan bastante.
