Diferencia entre revisiones de «Ecuaciones de Laplace y de Poisson»

De MateWiki
Saltar a: navegación, buscar
(→‎Errores de la fórmula de Poisson)
(→‎Ecuación de Laplace)
 
(No se muestran 11 ediciones intermedias de otro usuario)
Línea 1: Línea 1:
 +
{{ TrabajoED | Ecuación de Laplace y de Poisson. Grupo ALA | [[:Categoría:EDP|EDP]]|[[:Categoría:EDP23/24|2023-24]] | Lucía Amores, Aitana Guill y Andrea Navarro}}
 
== Introducción ==
 
== Introducción ==
 
En este documento, nos centraremos en dos ecuaciones que tienen un amplio uso en diversos ámbitos como electrostática, mecánica de fluidos, física teórica y magnetostática: la ecuación de Laplace y la ecuación de Poisson.  
 
En este documento, nos centraremos en dos ecuaciones que tienen un amplio uso en diversos ámbitos como electrostática, mecánica de fluidos, física teórica y magnetostática: la ecuación de Laplace y la ecuación de Poisson.  
 
Ambas ecuaciones las estudiaremos en el plano y las veremos aplicadas en problemas concretos. Veremos las limitaciones que presenta la fórumla de Poisson así como diferentes métodos analíticos para aproximar soluciones, y raíz de esto, analizaremos errores de aproximación. Por último, estudiaremos el comportamiento de soluciones tanto en un dominio dado, utilizando la desigualdad de Harnack, como asintóticamente en infinito.
 
Ambas ecuaciones las estudiaremos en el plano y las veremos aplicadas en problemas concretos. Veremos las limitaciones que presenta la fórumla de Poisson así como diferentes métodos analíticos para aproximar soluciones, y raíz de esto, analizaremos errores de aproximación. Por último, estudiaremos el comportamiento de soluciones tanto en un dominio dado, utilizando la desigualdad de Harnack, como asintóticamente en infinito.
−
 
−
CREO QUE AQUÍ ME FALTAN PUNTOS Y COMAS
 
  
 
== Ecuación de Laplace ==
 
== Ecuación de Laplace ==
−
La ecuación de Laplace, cuyo nombre nombre honra al distinguido físico-matemático Pierre-Simon Laplace, es una ecuación en derivadas parciales de tipo elíptico. REFERENCIA?
+
La ecuación de Laplace, cuyo nombre nombre honra al distinguido físico-matemático Pierre-Simon Laplace, es una ecuación en derivadas parciales de tipo elíptico.  
 
Construyamos un problema de derivadas parciales a partir de la ecuación de Laplace. La ecuación es,
 
Construyamos un problema de derivadas parciales a partir de la ecuación de Laplace. La ecuación es,
 
<center><math> \Delta u = 0 </math></center>
 
<center><math> \Delta u = 0 </math></center>
Línea 172: Línea 171:
 
A continuación, atendiendo a la fórmula del error "teórica" dada por el método del trapecio, la cual sigue la siguiente expresión:
 
A continuación, atendiendo a la fórmula del error "teórica" dada por el método del trapecio, la cual sigue la siguiente expresión:
 
<center><math>Error:=\left|-\cfrac{f''(\xi)(b-a)^3}{12n^2}\right|, </math></center>.
 
<center><math>Error:=\left|-\cfrac{f''(\xi)(b-a)^3}{12n^2}\right|, </math></center>.
 +
 +
siendo <math>\xi</math> un valor dentro del intervalo <math>[a,b]</math>, y n el valor que genera la discretización para la aplicación del método del trapecio. En nuestro caso, <math>[a,b]=[0,2\pi]</math> y <math>f</math> la función a integrar en la fórmula de Poisson.
 +
Ahora bien, el código de MatLab presentado al final de la sección implementa un cálculo del error máximo, tomando para cada n el error máximo en la discretización y acumulándolo en la lista <math>error\_max</math> de tal forma que en la gráfica se presenta una representación del error máximo teórico de la fórmula del trapecio.
 +
 +
[[Archivo: comparacion_errores_ejercicio2.png|400px|miniaturadeimagen|center|Comparación errores con su máximo]]
 +
 +
Como se puede apreciar de manera bastante clara en la gráfica, a medida que aumentan los puntos de la discretización, la cota va disminuyendo. Aún así, se puede apreciar que el error calculado para un punto alejado de la frontera se sigue manteniendo por debajo de dicha cota.
  
 
Estudiemos ahora el error cometido en puntos "cercanos" a la frontera. En este caso vamos utilizar una única discretización de 100 puntos de la fórmula del trapecio y puntos que cada vez se acercan más a la frontera. Partiendo del punto <math>x = (r,\theta)=(0.9,\pi/4)</math> vamos a ir avanzando hacia la frontera con puntos de la forma <math>(r,\theta)=(1-10^{-n}, \pi/4)</math>. La gráfica entonces resulta ser:
 
Estudiemos ahora el error cometido en puntos "cercanos" a la frontera. En este caso vamos utilizar una única discretización de 100 puntos de la fórmula del trapecio y puntos que cada vez se acercan más a la frontera. Partiendo del punto <math>x = (r,\theta)=(0.9,\pi/4)</math> vamos a ir avanzando hacia la frontera con puntos de la forma <math>(r,\theta)=(1-10^{-n}, \pi/4)</math>. La gráfica entonces resulta ser:
Línea 429: Línea 435:
  
 
=== Comportamiento de la solución ===
 
=== Comportamiento de la solución ===
 +
En esta sección, y atendiendo al comportamiento de la solución estudiaremos en qué región se han de encontrar las soluciones armónicas dentro  de la bola <math>B_R</math>, siendo <math>R</math> el radio de la bola, que tengan el mismo valor que nuestra solución en el (0,0).
 +
 +
La '''desigualdad de Harnack''' establece: Sea <math>u</math> armónica y <math>0\leq u</math> en <math>\Omega \in \mathbb{R}^n</math>. Supongamos <math>\overline{B_{R}(z)}\subset \Omega</math>. Entonces:
 +
 +
<center><math> \frac{R^{n-2} (R-r)}{(R+r)^{n-1}} u(z)\leq u(x) \leq \frac{R^{n-2}  (R+r)}{(R-r)^{n-1}}u(z), \forall x \in B_{R}(z).</math></center>
 +
 +
Para proceder al estudio planteado, comenzamos calculando el mínimo de la función <math>g(x,y)=xy</math> en la frontera de la bola <math>B_R</math> al que llamaremos <math>M</math>. Para facilitar los cálculos de aquí en adelante trabajaremos en coordenadas polares, por lo que nuestra función <math>g(x,y)=xy</math> pasa a ser <math>g(\theta)=R^2cos(\theta)sin(\theta)</math>, siendo <math>R</math> el radio de la bola. Para el cálculo del mínimo se ha empleado un código de Matlab presentado al final de esta sección.
 +
 +
A continuación, definimos la función <math>v:=u-M </math> , siendo u nuestra solución del problema. Esta función <math>v</math> es positiva y armónica, por lo que podemos aplicar la desigualdad de Harnack para <math>\mathbb{R}</math>, es decir, <math>n=2</math>. Obteniendo lo siguiente:
 +
<center><math> \frac{(R-r)}{(R+r)} v(0,0) \leq v(x,y) \leq \frac{(R+r)}{(R-r)}v(0,0), \forall (x,y) \in B_{R}. </math></center>
 +
 +
Representando gráficamente las cotas planteadas, obtenemos las siguientes gráficas:
 +
 +
[[Archivo: harnack_ejercicio4_n2.png|600px|thumb|center|Desigualdad de Harnack en <math>R^2</math>]]
 +
 +
Las desigualdades de Harnack establecen límites uniformes para todas las soluciones armónicas dentro de una bola, independientemente de su radio. Sin embargo, el valor específico de estas soluciones en un punto dado, como <math>v(0,0)</math>, varía con el tamaño de la bola, lo que sugiere una especie de traslación en la región donde las soluciones armónicas son iguales a <math>v</math> en el punto <math>(0,0)</math>. A medida que aumenta el radio de la bola, tanto las cotas superior e inferior aumentan en amplitud, pero la región cubierta por estas cotas se estrecha. Este fenómeno indica que, aunque la estructura de las cotas permanece constante, el rango de valores en el que pueden variar las soluciones armónicas se reduce a medida que la bola crece en tamaño.
 +
 +
Procediendo a un estudio de la desigualdad en <math>\mathbb{R}^3</math> obtenemos las siguientes desigualdades:
 +
<center><math> \frac{R^3(R-r)}{2(R+r)^{2}} \leq v(x,y,z) \leq \frac{R^3(R+r)}{2(R-r)^2} , \forall (x,y,z) \in B_{R}, </math></center>
 +
 +
Llevando a cabo un procedimiento análogo al anterior obtenemos las siguientes gráficas:
 +
 +
[[Archivo: harnack_ejercicio4_n3.png|600px|thumb|center|Desigualdad de Harnack en <math>R^3</math>]]
 +
 +
Se puede concluir de manera análoga al caso de <math>R^2</math>, pero con una región mayor.
 +
 +
El código de Matlab empleado para la representación de las gráficas y el cálculo de las cotas en <math>R^2</math> es el siguiente:
 +
{{matlab|codigo=
 +
clear
 +
close all
 +
clc
 +
format short
 +
 +
% Definimos los radios de las bolas:
 +
RR=[1 2 10];
 +
for i=1:length(RR)
 +
    R=RR(i);
 +
    % Definimos la dimensión del espacio:
 +
    n=2;
 +
   
 +
    % Definimos la función g definida en la frontera:
 +
    g=@(theta) R.^2.*cos(theta).*sin(theta);
 +
   
 +
    % Definimos la solución exacta u:
 +
    u_exacta=@(r,theta) r.^2.*cos(theta).*sin(theta);
 +
   
 +
    % Calculamos el mínimo de esta función:
 +
    N=100;
 +
    theta=linspace(0,2*pi,N);
 +
    r_con_frontera=linspace(0,R,N+1);
 +
    r=r_con_frontera(1:N);
 +
    M=min(g(theta));
 +
   
 +
   
 +
    % Consideramos la función armónica v=u-M:
 +
    v=@(r,theta) u_exacta(r,theta)-M;
 +
   
 +
    % Definimos las cotas de la desigualdad de Harnack:
 +
    cota_inferior=@(r)(R.^(n-2).*(R-r)./((R+r).^(n-1))).*v(0,0);
 +
    cota_superior=@(r)(R.^(n-2).*(R+r)./((R-r).^(n-1))).*v(0,0);
 +
   
 +
    % Dibujamos las cotas para encontrar la región en la que deben estar las
 +
    % soluciones armónicas de la B_1, que valen lo mismo que u en el (0,0):
 +
    subplot(1,3,i)
 +
    hold on
 +
    plot(r,log10(cota_inferior(r)),"Color",'red')
 +
    plot(r,log10(cota_superior(r)),"Color",'blue')
 +
    xlabel('r')
 +
    legend('Cota inferior','Cota superior','Location','southwest')
 +
    title('Región de soluciones en la Bola de radio',num2str(R))
 +
    grid on
 +
    grid minor
 +
    hold off
 +
end
 +
}}
 +
El código de Matlab empleado para la representación de las gráficas y el cálculo de las cotas en <math>R^3</math> es el siguiente:
 +
{{matlab|codigo=
 +
clear
 +
close all
 +
clc
 +
format long
 +
 +
% Definimos los radios de las bolas:
 +
RR=[1 2 10];
 +
for i=1:length(RR)
 +
    R=RR(i);
 +
    % Definimos la dimensión del espacio:
 +
    n=3;
 +
   
 +
    % Definimos la función g definida en la frontera:
 +
    g=@(theta) R.^2.*cos(theta).*sin(theta);
 +
   
 +
    % Definimos la solución exacta u:
 +
    u_exacta=@(r,theta) r.^2.*cos(theta).*sin(theta);
 +
   
 +
    % Calculamos el mínimo de esta función:
 +
    N=300;
 +
    theta=linspace(0,2*pi,N);
 +
    r_con_frontera=linspace(0,R,N+1);
 +
    r=r_con_frontera(1:N);
 +
    M=min(g(theta));
 +
   
 +
    % Consideramos la función armónica v=u-M:
 +
    v=@(r,theta) u_exacta(r,theta)-M;
 +
   
 +
    % Definimos las cotas de la desigualdad de Harnack:
 +
    cota_inferior=@(r)R.^(n-2).*(R-r)./(R+r).^(n-1).*v(0,0);
 +
    cota_superior=@(r)R.^(n-2).*(R+r)./(R-r).^(n-1).*v(0,0);
 +
   
 +
    % Dibujamos las cotas para encontrar la región en la que deben estar las
 +
    % soluciones armónicas de la B_1, que valen lo mismo que u en el (0,0):
 +
    subplot(1,3,i)
 +
    hold on
 +
    plot(r,log10(cota_inferior(r)),"Color",'red')
 +
    plot(r,log10(cota_superior(r)),"Color",'blue')
 +
    xlabel('r')
 +
    legend('Cota inferior','Cota superior','Location','southwest')
 +
    title('Región de soluciones en la Bola de radio',num2str(R))
 +
    grid on
 +
    grid minor
 +
    hold off
 +
end
 +
}}
  
 
== Ecuación de Poisson ==
 
== Ecuación de Poisson ==
Línea 543: Línea 672:
  
 
Tal y como se ve, se observa que el comportamiento es el esperado.
 
Tal y como se ve, se observa que el comportamiento es el esperado.
 +
 +
[[Categoría:Grado en Matemáticas]]
 +
[[Categoría:EDP]]
 +
[[Categoría:EDP23/24]]

Revisión actual del 23:21 19 abr 2024

Trabajo realizado por estudiantes
Título Ecuación de Laplace y de Poisson. Grupo ALA
Asignatura EDP
Curso 2023-24
Autores Lucía Amores, Aitana Guill y Andrea Navarro
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

1 Introducción

En este documento, nos centraremos en dos ecuaciones que tienen un amplio uso en diversos ámbitos como electrostática, mecánica de fluidos, física teórica y magnetostática: la ecuación de Laplace y la ecuación de Poisson. Ambas ecuaciones las estudiaremos en el plano y las veremos aplicadas en problemas concretos. Veremos las limitaciones que presenta la fórumla de Poisson así como diferentes métodos analíticos para aproximar soluciones, y raíz de esto, analizaremos errores de aproximación. Por último, estudiaremos el comportamiento de soluciones tanto en un dominio dado, utilizando la desigualdad de Harnack, como asintóticamente en infinito.

2 Ecuación de Laplace

La ecuación de Laplace, cuyo nombre nombre honra al distinguido físico-matemático Pierre-Simon Laplace, es una ecuación en derivadas parciales de tipo elíptico. Construyamos un problema de derivadas parciales a partir de la ecuación de Laplace. La ecuación es,

[math] \Delta u = 0 [/math]

con [math]u:\mathbb{R}^3 \rightarrow \mathbb{R}[/math] como incógnita. Esta ecuación se define en un dominio [math] \Omega [/math], y en su frontera [math] \partial \Omega [/math] se pueden añadir condiciones de contorno Dirichlet, Neuman o mixtas. En este trabajo nos centraremos en las condiciones de contorno de tipo Dirichlet, las cuales igualan [math] u [/math] a una función específica [math] g [/math]. Con todo esto, llegamos al siguiente problema de derivadas parciales.

[math] \begin{cases} \Delta u = 0, & \text{x} \in \Omega \\ u = g, & \text{x} \in \partial \Omega \end{cases}[/math]

A fin de comprender mejor este problema, podemos examinar el siguiente ejemplo en concreto.

2.1 Ejemplo bola unidad

Sea [math] B_1 ⊂ R^2 [/math] la bola unidad centrada en el origen. Planteamos el problema,

[math] \begin{cases} \Delta u = 0, & \text{x} \in B_1 \\ u = g, & \text{x} \in \partial B_1 \end{cases}[/math]

A lo largo de este documento, en lo que tiene que ver con la ecuación de Laplace, visitaremos este ejemplo varias veces.


Dentro de la resolución del problema planteado existen varios métodos por los que proceder. En este trabajo estudiaremos dos: la fórmula de Poisson y la serie de Fourier.

2.2 Solución por la fórmula de Poisson

Encontrar la solución del problema mediante la fórmula de Poisson viene determinado por el siguiente teorema.

2.2.1 Teorema
La solución [math] u \in C^2(B_R \cup C(\overline{B_R}) [/math]del problema [math] \begin{cases} \Delta u = 0, & \text{x} \in B_R \\ u = g, & \text{x} \in \partial B_R \end{cases}[/math] donde [math]g[/math] es una función continua viene dado por la fórmula de Poisson
[math] u(\vec{x})=\frac{R^2-|\vec{x}|^2}{w_n R}\int_{\partial B_R}\frac{g(\sigma)}{|\vec{x}-\sigma|^2} d\sigma. [/math]


Es importante destacar que el denominador [math]|\vec{x}-\sigma|^2 [/math] se anula en la frontera de la bola llevándonos a una indeterminación y haciendo que la integral diverja. Esto también pasa si expresamos la fórmula en coordenadas polares [math](r, \theta)[/math]. Tomando,

[math]\vec{x} = (x_1, x_2) = (rcos(\theta), rsen(\theta)) [/math]
[math] g(rcos(\theta), rsen(\theta)) = G(\theta)[/math]

llegamos a la fórmula de Poisson,

[math] U(r,\theta)=\frac{R^2-r^2}{2\pi R}\int_{0}^{2\pi}\frac{g(s)}{R^2+r^2-2Rrcos(\theta-s)} ds. [/math]

En este caso, el denominador de la integral también se anula para puntos cercanos a la frontera, concretamente cuando el coseno se hace 1.

Veamos todo esto aplicado a la bola unidad pasado a coordenadas polares,

[math] \begin{cases} U_{rr} + \frac{1}{r} U_r + \frac{1}{r^2} U_{\theta\theta} = 0, & (r, \theta) \in B_1 \\ U(1, \theta) = G(\theta), & (r, \theta) \in \partial B_1 \end{cases}[/math]

Para la resolución de este problema vamos a tomar [math]G(\theta) = máx\{0, 1-\frac{2}{\pi} |\theta - \pi|\}[/math] y utilizaremos la fórmula de Poisson en polares tomando [math]R = 1[/math] y calculando la integral de manera aproximada con la fórmula del trapecio. Como hemos visto antes, la fórmula da problemas para puntos cercanos al borde. Esto lleva a una representación irregular de la frontera si dibujamos la solución sin tener en cuenta este problema.

Solución del problema de Poisson sin aplicar la condición frontera

Es por esto por lo que la condición frontera no es prescindible en el problema. Para estimar [math] U(r,\theta)[/math], utilizaremos la fórmula de Poisson para puntos ligeramente alejados del borde, y en el propio borde utilizaremos la condición frontera [math]U(R, \theta) = G(\theta)[/math].

Solución del problema de Poisson aplicando la condición frontera

En esta nueva gráfica podemos apreciar como la frontera no presenta las irregularidades anteriores. Para conseguir estas gráficas hemos implementado el siguiente código en MatLab.

2.2.1.1 Códigos
% Este código resuelve el problema planteado con la fórmula de Poisson sin aplicar condición frontera:
clear
close all
clc
format long

% Radio de la bola:
R=1;

% Definimos nuestra función g:
g=@(theta)max(0, 1-2/pi*abs(theta-pi));

% Definimos los intervalos r y theta sobre los que vamos a trabajar:
N=200; % Número de divisiones del intervalo
r=linspace(0,R,N); % No podemos llegar hasta la frontera
theta=linspace(0,2*pi,N); % Intervalo theta.

% Calculamos la solución de u con la fórmula de Poisson:
NN=100; % División del intervalo de integración.
s=linspace(0,2*pi,NN); % Intervalo de integración, para el método del trapecio.
u=zeros(N,N); % Inicializamos la matriz que almacena los valores de la solución.
for i=1:length(theta)
    for j=1:length(r) % Calculamos solución a través de la fórmula del trapecio:
        u(i,j)=1/(2*pi).*trapz(s,g(s).*(R^2-r(j).^2)./(R^2+r(j).^2-2.*r(j).*cos(s-theta(i))));
    end
end
% Creamos la malla:
[RR,THETA]=meshgrid(r,theta);
figure
surf(RR.*cos(THETA),RR.*sin(THETA),u,'EdgeColor','flat') % Representación gráfica
xlabel('x')
ylabel('y')
zlabel('u(x,y)')
title('Solución de la ecuación de Laplace')
% Este código resuelve el problema planteado con la fórmula de Poisson aplicando la condición frontera:
clear
close all
clc
format long

% Radio de la bola:
R=1;

% Definimos nuestra función g:
g=@(theta)max(0, 1-2/pi*abs(theta-pi));

% Definimos los intervalos r y theta sobre los que vamos a trabajar:
N=200; % Número de divisiones del intervalo
r=linspace(0,R-N^(-1),N-1); % No podemos llegar hasta la frontera
theta=linspace(0,2*pi,N); % Intervalo theta.

% Calculamos la solución de u con la fórmula de Poisson:
NN=1000; % División del intervalo de integración.
s=linspace(0,2*pi,NN); % Intervalo de integración, para el método del trapecio.
u=zeros(N,N); % Inicializamos la matriz que almacena los valores de la solución.
for i=1:length(theta)
    for j=1:length(r) % Calculamos solución a través de la fórmula del trapecio:
        u(i,j)=1/(2*pi).*trapz(s,g(s).*(R^2-r(j).^2)./(R^2+r(j).^2-2.*r(j).*cos(s-theta(i))));
    end
    % Establecemos condiciones frontera (por presentar problemas en esta):
    u(i,N)=g(theta(i));
end
r(N)=1;
% Creamos la malla:
[RR,THETA]=meshgrid(r,theta);
figure
surf(RR.*cos(THETA),RR.*sin(THETA),u,'EdgeColor','flat') % Representación gráfica
xlabel('x')
ylabel('y')
zlabel('u(x,y)')
title('Solución de la ecuación de Laplace')


2.2.2 Errores de la fórmula de Poisson

Como hemos observado anteriormente, la fórmula de Poisson presenta ciertas dificultades a la hora de aproximar la solución en la frontera, debido a la singularidad inherente de la integral. En este apartado examinaremos estas irregularidades calculando el error de la aproximación frente a la solución exacta.

Volviendo al ejemplo de la bola unidad, esta vez vamos a suponer el problema,

[math] \begin{cases} \Delta u = 0, & \text{x} \in B_1 \\ u = g(x,y) = xy, & \text{x} \in \partial B_1 \end{cases}[/math]

que tiene como solución exacta la función armónica [math]u(x,y) = xy[/math]. Tal y como hemos hecho en el apartado anterior, calculamos la solución aproximada utilizando la fórmula de Poisson, la fórmula del trapecio y pasando a coordenadas polares.

[math] \begin{cases} U_{rr} + \frac{1}{r} U_r + \frac{1}{r^2} U_{\theta\theta} = 0, & (r, \theta) \in B_1 \\ U(1,\theta)=G(\theta)=\cos(\theta)\sin(\theta), & (r, \theta) \in \partial B_1 \end{cases}[/math]

Vamos a distinguir entre dos tipos de puntos, aquellos que estén "alejados" de la frontera y aquellos sean inmediatos a esta. Primero estudiaremos el error en puntos "alejados" de la frontera. Esto lo hacemos empleando distintas discretizaciones con [math] 10^n [/math] puntos en la fórmula del trapecio. Después, calculamos el error para cada discretización en un punto alejado, en este caso lo hacemos evaluando en [math]x = (r,\theta)=(0.9,\pi/4)[/math]. Posteriormente graficamos el error aplicando la fórmula,

[math]f(n):=\log_{10}(Error(10^n)), [/math]
.

Obteniendo así la siguiente gráfica.

Error de un punto alejado a la frontera

Los resultados de esta gráfica cuadran con nuestra intuición, a mayor número de puntos menor es el error que se comete. Además podemos apreciar que el error se acaba estabilizando al rededor de [math]10^{-15}[/math]. Luego podemos concluir que la fórmula de Poisson proporciona una buena aproximación para puntos relativamente alejados de la frontera.

A continuación, atendiendo a la fórmula del error "teórica" dada por el método del trapecio, la cual sigue la siguiente expresión:

[math]Error:=\left|-\cfrac{f''(\xi)(b-a)^3}{12n^2}\right|, [/math]
.

siendo [math]\xi[/math] un valor dentro del intervalo [math][a,b][/math], y n el valor que genera la discretización para la aplicación del método del trapecio. En nuestro caso, [math][a,b]=[0,2\pi][/math] y [math]f[/math] la función a integrar en la fórmula de Poisson. Ahora bien, el código de MatLab presentado al final de la sección implementa un cálculo del error máximo, tomando para cada n el error máximo en la discretización y acumulándolo en la lista [math]error\_max[/math] de tal forma que en la gráfica se presenta una representación del error máximo teórico de la fórmula del trapecio.

Comparación errores con su máximo

Como se puede apreciar de manera bastante clara en la gráfica, a medida que aumentan los puntos de la discretización, la cota va disminuyendo. Aún así, se puede apreciar que el error calculado para un punto alejado de la frontera se sigue manteniendo por debajo de dicha cota.

Estudiemos ahora el error cometido en puntos "cercanos" a la frontera. En este caso vamos utilizar una única discretización de 100 puntos de la fórmula del trapecio y puntos que cada vez se acercan más a la frontera. Partiendo del punto [math]x = (r,\theta)=(0.9,\pi/4)[/math] vamos a ir avanzando hacia la frontera con puntos de la forma [math](r,\theta)=(1-10^{-n}, \pi/4)[/math]. La gráfica entonces resulta ser:

Error de puntos cercanos a la frontera

En este caso el error incrementa según los puntos se acercan a la frontera y se acaba estabilizando en [math]-0,3[/math]. El código utilizado para dibujar estas gráficas es el siguiente. 

% El siguiente código ejecuta todas las gráficas planteadas a lo largo de esta sección
clear
close all
clc
format long

% Definimos el radio de la bola;
R=1;

% Definimos la función g:
g =@(theta) (R.^2.*cos(theta).*sin(theta));

% Definimos la solución exacta:
u_exacta =@(r,theta) (r.^2.*cos(theta).*sin(theta));

% Tomamos el punto r=0.9 y theta=pi/4:
r_valor=0.9;
theta_valor=pi/4;

% Calculamos la solución u en el punto anterior:
n=1:1:8;
error=zeros(1,length(n)); % Inicializamos la matriz que almacena el error
error_log=zeros(1,length(n)); % Inicializamos la matriz que almacena el error logarítmico
u_valor=zeros(1,length(n)); % Inicializamos la matriz que almacena la solución de u para distintas discretizaciones del método del trapecio
for i=1:length(n)
    s=linspace(0,2*pi,10^(n(i))); % Intervalo de integración
    u_valor(i)=1/(2*pi).*trapz(s,g(s).*(R^2-r_valor.^2)./(R^2+r_valor.^2-2.*r_valor.*cos(s-theta_valor))); % Solución a través de la fórmula de Poisson
    
    % Calculamos el error:
    error(i)=abs(u_exacta(r_valor,theta_valor)-u_valor(i));
    % Calculamos el error logarítmico:
    error_log(i)=log10(error(i));
end
% Representamos gráficamente el error logarítmico para el punto r=0.9 y theta=pi/4:
plot(n,error_log,'o-')
xlabel('n')
ylabel('log_{10}(error(10^n))')
grid on 
grid minor
legend('Errores')
title('Gráfica del error logarítmico')


% Observamos la existencia de un error máximo que nos permita confirmar
% que los errores calculados anteriormente se encuentran por debajo siempre:

% Calculamos la segunda derivada (requerida por la fórmula):
syms rr ss thetatheta
% Definimos de forma simbólica la función a derivar:
f=rr^2*cos(ss)*sin(ss)/(1+rr^2-2*rr*cos(thetatheta-ss));
derivada1=diff(f,ss); % Derivamos la función f en función de la variable s
derivada2=diff(derivada1,ss); % Derivamos la primera derivada en función de la variable s
% Calculamos el error máximo del método del trapecio:
nn=1:1:3; % Conforme el número sea más elevado aquí más aumentará el tiempo de compilación del programa.
error_max=zeros(length(nn),1); % Inicializamos a ceros el vector que almacenará los datos del error máximo del trapecio
for i=1:length(nn)
    N=10^i;
    s=linspace(0,2*pi,N);
    error_trapz=zeros(N,1); % Inicializamos a ceros el vector que almacenará los datos del error del trapecio para cada n
    for j=1:N
        s_valor=s(j);
        error_trapz(j)=log10(abs(-(2*pi).^3.*subs(derivada2, [rr,ss,thetatheta] ,[r_valor,s_valor,theta_valor])./(12.*(N.^2))));
    end
    error_max(i)=max(error_trapz);
end
% Representamos gráficamente el error máximo junto con el error calculado
% anteriormente:
figure
hold on
plot(nn,error_log(1:length(nn)),'Color','blue')
plot(nn,error_max,'Color','red')
hold off
legend('Error','Error máximo')
title('Comparación de errores')
grid on
grid minor
xlabel('n')
ylabel('Error')


% Fijamos ahora para N=100 (puntos en la fórmula del trapecio).
% Vamos a calcular el error en los puntos de la forma r=1-10^(-n) y theta=pi/4 y veremos qué ocurre:
N=100;
r_valores=zeros(1,N); % Inicializamos vector que almacena los valores de r.
u_valores_r=zeros(1,N); % Inicializamos el vector que almacena los valores de la solución u para distintos valores de r.
u_exacta_r=zeros(1,N); % Inicializamos el vector que almacena los valores de la solución u exacta para distintos valores de r.
error_r=zeros(1,N); % Inicializamos el vector que almacena los errores.
error_log_r=zeros(1,N); % Inicializamos el vector que almacena los errores logarítmicos.
% theta_valor=pi/4 es el mismo que hemos usado antes, por lo que no lo volvemos a inicializar.
% Definimos el vector de integración de la fórmula del trapecio:
s=linspace(0,2*pi,N);
for j=1:N
    r_valores(j)=1-10^(-j);
    u_valores_r(j)=1/(2*pi).*trapz(s,g(s).*(R^2-r_valores(j).^2)./(R^2+r_valores(j).^2-2.*r_valores(j).*cos(s-theta_valor)));
    u_exacta_r(j)=u_exacta(r_valores(j),theta_valor);
    error_r(j)=abs(u_exacta_r(j)-u_valores_r(j));
    error_log_r(j)=log10(error_r(j));
end

% Representamos gráficamente
figure
plot(1:1:N,error_log_r,'.-')
xlabel('n')
ylabel('log_{10}(error(100))')
grid on 
grid minor
legend('Errores')
title('Errores en los puntos cercanos a la frontera')


2.3 Solución por la serie de Fourier

Una vez estudiada la solución por la fórmula de Poisson, en este apartado calcularemos la solución por serie de Fourier. La obtención de esta solución se basa en una cambio de variable a polares del problema, como el visto previamente, y resolución por separación de variables. Hacemos esto aplicado al problema anterior en una bola de radio [math] R [/math].

[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 sin(\theta)cos(\theta), \theta\in[0,2\pi). \end{cases} [/math]

Lo cual nos conduce a,

[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]

que tiene coeficientes,

[math]\alpha_{0} =\frac{1}{\pi}\int_{-\pi}^{\pi} \cos\Bigl(k \theta \Bigl) dx [/math]
[math]\alpha_{k} =\frac{1}{\pi}\int_{-\pi}^{\pi} cos(\theta) \sin(\theta) \cos\Bigl(k \theta \Bigl) dx [/math]


[math]\beta_{k} =\frac{1}{\pi}\int_{-\pi}^{\pi} cos(\theta) \sin(\theta) \sin\Bigl(k \theta \Bigl) dx [/math]

Sin embargo, todos los coeficientes de Fourier son nulos, salvo por [math] \beta_2=\frac{R^2}{2}[/math], luego la solución por series de Fourier final es,

[math]U(r,\theta)=\frac{r^2}{2}\sin(2\theta)=r^2\cos(\theta)\sin(\theta), r\in [0,R], \theta\in[0,2\pi)[/math]
.

Veamos ahora la gráfica que dibuja esta solución tomando [math] R = 1 [/math].

Representación de la solución
%% Ecuación de Laplace
%% Apartado 3
clear
close all
clc
format long

% Definimos el radio de la bola:
R=1;

% k términos de la serie:
K=[5 10 100];

% Definimos nuestra función g y la función u exacta:
g=@(theta) R.^2.*cos(theta).*sin(theta);
u_exacta=@(r,theta) r.^2.*cos(theta).*sin(theta);

% Definimos los intervalos de theta y r:
N=300; % Número de divisiones de los intervalos theta y r.
Theta=linspace(0,2*pi,N); % Intervalo theta.
r=linspace(0,R,N); % Intervalo r.


% Calculamos los coeficientes de Fourier para distintos k-términos:
for i=1:length(K)
    k=K(i); % k-términos de la serie.
    int=linspace(-pi,pi,k); % Intervalo de integración para los coef de Fourier.
    a0=1/pi.*integral(g,-pi,pi); % Coeficiente de Fourier a_0 (asociado a 1/2).
    a=zeros(k,1); % Inicializamos la matriz de los coeficientes de Fourier asociados al coseno.
    b=zeros(k,1); % Inicializamos la matriz de los coeficientes de Fourier asociados al seno.
    for j=1:k
        f_cos=@(theta) g(theta).*cos(j.*theta); % Definimos la función a integrar asociada al coseno.
        f_sin=@(theta) g(theta).*sin(j.*theta); % Definimos la función a integrar asociada al seno.
        a(j)=1/pi.*integral(f_cos,-pi,pi); % Integramos de manera exacta obteniendo los coeficientes de Fourier del coseno.
        b(j)=1/pi.*integral(f_sin,-pi,pi); % Integramos de manera exacta obteniendo los coeficientes de Fouries del seno.
    end
    % Calculamos la solución u_k asociada al número de términos de la serie:
    suma=@(r,theta)0;
    for j=1:k
        suma=@(r,theta) suma(r,theta) + a(j).*(r/R).^j.*cos(j.*theta)+b(j).*(r/R).^j.*sin(j.*theta);
    end
    eval(['u_' num2str(k) '=@(r,theta) 1/2*a0 + suma(r,theta);']);
end

[RR,THETA]=meshgrid(r,Theta);
figure
for i=1:length(K)
    k=K(i);
    subplot(2,2,i)
    eval(['surf(RR.*cos(THETA),RR.*sin(THETA),u_' num2str(k) '(RR,THETA),"EdgeColor","flat")'])
    title("Solución de u para k=", num2str(k))
end
subplot(2,2,4)
surf(RR.*cos(THETA),RR.*sin(THETA),u_exacta(RR,THETA),"EdgeColor","flat")
title('Solución exacta')



2.3.1 Errores de la solución por serie de Fourier

Al igual que con la solución dada por la fórmula de Poisson, vamos a comparar la solución calculada con la solución exacta [math]u(x,y)=xy[/math] en la bola de radio 1. Vamos a contrastar las soluciones con 5, 10 y 100 términos con la exacta. Calcularemos el error resultante al aproximar la serie de Fourier y lo representamos gráficamente, mostrando el máximo error en función del número de términos utilizados en una escala logarítmica, es decir, pintaremos la función,
[math] \log_{10}(error(n)) [/math]

donde [math] error(n) [/math] es ahora el supremo del error con [math]n [/math]términos. Obetenemos así la siguente gráfica.

Representación de la solución

Podemos apreciar en esta gráfica errores de orden muy pequeño, por lo general, los podemos considerar pequeño. Luego, nos damos cuenta de que, a diferencia de solución por la fórmula de Poisson, esta no nos da problemas en la frontera. Veamos el código utilizado para hacer la gráfica. 

% A continuación se calcularán los errores de la aproximación de Fourier en función de los términos:
KK=0:1:100;
error_log=zeros(length(KK),1); % Inicializamos el vector donde se almacenarán los errores logarítmicos.
for i=1:length(KK)
    k=KK(i);
    int=linspace(-pi,pi,k); % Intervalo de integración para los coef de Fourier.
    a0=1/pi.*integral(g,-pi,pi); % Coeficiente de Fourier a_0 (asociado a 1/2).
    a=zeros(k,1); % Inicializamos la matriz de los coeficientes de Fourier asociados al coseno.
    b=zeros(k,1); % Inicializamos la matriz de los coeficientes de Fourier asociados al seno.
    for j=1:k
        f_cos=@(theta) g(theta).*cos(j.*theta); % Definimos la función a integrar asociada al coseno.
        f_sin=@(theta) g(theta).*sin(j.*theta); % Definimos la función a integrar asociada al seno.
        a(j)=1/pi.*integral(f_cos,-pi,pi); % Integramos de manera exacta obteniendo los coeficientes de Fourier del coseno.
        b(j)=1/pi.*integral(f_sin,-pi,pi); % Integramos de manera exacta obteniendo los coeficientes de Fouries del seno.
    end
    % Calculamos la solución u_k asociada al número de términos de la serie:
    suma=@(r,theta)0;
    for j=1:k
        suma=@(r,theta) suma(r,theta) + a(j).*(r/R).^j.*cos(j.*theta)+b(j).*(r/R).^j.*sin(j.*theta);
    end
    eval(['u_' num2str(i) '=@(r,theta) 1/2*a0 + suma(r,theta);']);
    eval(['u=@(r,theta)u_' num2str(i) '(r,theta) ;']);
    U_exacta=u_exacta(RR,THETA);
    U_fourier=u(RR,THETA);
    valores_max=max(U_exacta-U_fourier);
    sup=max(valores_max,[],'all');
    error_log(i)=log10(sup);
end


% Representamos gráficamente
figure
plot(KK,error_log)
legend('Errores')
title('Errores supremos frente al número de términos')
xlabel('n')
ylabel('log_{10}(error(10^n))')


2.4 Comportamiento de la solución

En esta sección, y atendiendo al comportamiento de la solución estudiaremos en qué región se han de encontrar las soluciones armónicas dentro de la bola [math]B_R[/math], siendo [math]R[/math] el radio de la bola, que tengan el mismo valor que nuestra solución en el (0,0).

La desigualdad de Harnack establece: Sea [math]u[/math] armónica y [math]0\leq u[/math] en [math]\Omega \in \mathbb{R}^n[/math]. Supongamos [math]\overline{B_{R}(z)}\subset \Omega[/math]. Entonces:

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

Para proceder al estudio planteado, comenzamos calculando el mínimo de la función [math]g(x,y)=xy[/math] en la frontera de la bola [math]B_R[/math] al que llamaremos [math]M[/math]. Para facilitar los cálculos de aquí en adelante trabajaremos en coordenadas polares, por lo que nuestra función [math]g(x,y)=xy[/math] pasa a ser [math]g(\theta)=R^2cos(\theta)sin(\theta)[/math], siendo [math]R[/math] el radio de la bola. Para el cálculo del mínimo se ha empleado un código de Matlab presentado al final de esta sección.

A continuación, definimos la función [math]v:=u-M [/math] , siendo u nuestra solución del problema. Esta función [math]v[/math] es positiva y armónica, por lo que podemos aplicar la desigualdad de Harnack para [math]\mathbb{R}[/math], es decir, [math]n=2[/math]. Obteniendo lo siguiente:

[math] \frac{(R-r)}{(R+r)} v(0,0) \leq v(x,y) \leq \frac{(R+r)}{(R-r)}v(0,0), \forall (x,y) \in B_{R}. [/math]

Representando gráficamente las cotas planteadas, obtenemos las siguientes gráficas:

Desigualdad de Harnack en [math]R^2[/math]

Las desigualdades de Harnack establecen límites uniformes para todas las soluciones armónicas dentro de una bola, independientemente de su radio. Sin embargo, el valor específico de estas soluciones en un punto dado, como [math]v(0,0)[/math], varía con el tamaño de la bola, lo que sugiere una especie de traslación en la región donde las soluciones armónicas son iguales a [math]v[/math] en el punto [math](0,0)[/math]. A medida que aumenta el radio de la bola, tanto las cotas superior e inferior aumentan en amplitud, pero la región cubierta por estas cotas se estrecha. Este fenómeno indica que, aunque la estructura de las cotas permanece constante, el rango de valores en el que pueden variar las soluciones armónicas se reduce a medida que la bola crece en tamaño.

Procediendo a un estudio de la desigualdad en [math]\mathbb{R}^3[/math] obtenemos las siguientes desigualdades:

[math] \frac{R^3(R-r)}{2(R+r)^{2}} \leq v(x,y,z) \leq \frac{R^3(R+r)}{2(R-r)^2} , \forall (x,y,z) \in B_{R}, [/math]

Llevando a cabo un procedimiento análogo al anterior obtenemos las siguientes gráficas:

Desigualdad de Harnack en [math]R^3[/math]

Se puede concluir de manera análoga al caso de [math]R^2[/math], pero con una región mayor.

El código de Matlab empleado para la representación de las gráficas y el cálculo de las cotas en [math]R^2[/math] es el siguiente: 

clear
close all
clc
format short

% Definimos los radios de las bolas:
RR=[1 2 10];
for i=1:length(RR)
    R=RR(i);
    % Definimos la dimensión del espacio:
    n=2;
    
    % Definimos la función g definida en la frontera:
    g=@(theta) R.^2.*cos(theta).*sin(theta);
    
    % Definimos la solución exacta u:
    u_exacta=@(r,theta) r.^2.*cos(theta).*sin(theta);
    
    % Calculamos el mínimo de esta función:
    N=100;
    theta=linspace(0,2*pi,N);
    r_con_frontera=linspace(0,R,N+1);
    r=r_con_frontera(1:N);
    M=min(g(theta));
    
    
    % Consideramos la función armónica v=u-M:
    v=@(r,theta) u_exacta(r,theta)-M;
    
    % Definimos las cotas de la desigualdad de Harnack:
    cota_inferior=@(r)(R.^(n-2).*(R-r)./((R+r).^(n-1))).*v(0,0);
    cota_superior=@(r)(R.^(n-2).*(R+r)./((R-r).^(n-1))).*v(0,0);
    
    % Dibujamos las cotas para encontrar la región en la que deben estar las
    % soluciones armónicas de la B_1, que valen lo mismo que u en el (0,0):
    subplot(1,3,i)
    hold on
    plot(r,log10(cota_inferior(r)),"Color",'red')
    plot(r,log10(cota_superior(r)),"Color",'blue')
    xlabel('r')
    legend('Cota inferior','Cota superior','Location','southwest')
    title('Región de soluciones en la Bola de radio',num2str(R))
    grid on
    grid minor
    hold off
end

El código de Matlab empleado para la representación de las gráficas y el cálculo de las cotas en [math]R^3[/math] es el siguiente: 

clear
close all
clc
format long

% Definimos los radios de las bolas:
RR=[1 2 10];
for i=1:length(RR)
    R=RR(i);
    % Definimos la dimensión del espacio:
    n=3;
    
    % Definimos la función g definida en la frontera:
    g=@(theta) R.^2.*cos(theta).*sin(theta);
    
    % Definimos la solución exacta u:
    u_exacta=@(r,theta) r.^2.*cos(theta).*sin(theta);
    
    % Calculamos el mínimo de esta función:
    N=300;
    theta=linspace(0,2*pi,N);
    r_con_frontera=linspace(0,R,N+1);
    r=r_con_frontera(1:N);
    M=min(g(theta));
    
    % Consideramos la función armónica v=u-M:
    v=@(r,theta) u_exacta(r,theta)-M;
    
    % Definimos las cotas de la desigualdad de Harnack:
    cota_inferior=@(r)R.^(n-2).*(R-r)./(R+r).^(n-1).*v(0,0);
    cota_superior=@(r)R.^(n-2).*(R+r)./(R-r).^(n-1).*v(0,0);
    
    % Dibujamos las cotas para encontrar la región en la que deben estar las
    % soluciones armónicas de la B_1, que valen lo mismo que u en el (0,0):
    subplot(1,3,i)
    hold on
    plot(r,log10(cota_inferior(r)),"Color",'red')
    plot(r,log10(cota_superior(r)),"Color",'blue')
    xlabel('r')
    legend('Cota inferior','Cota superior','Location','southwest')
    title('Región de soluciones en la Bola de radio',num2str(R))
    grid on
    grid minor
    hold off
end


3 Ecuación de Poisson

En esta parte del documento se va a proceder a resolver la ecuación de Poisson mediante un nuevo método. En primer lugar, esta ecuación viene dada por [math] \Delta u = f [/math] siendo [math]u:\mathbb{R}^n \rightarrow \mathbb{R}[/math]y [math] f \in C^2(\mathbb{R})[/math]. Para este estudio se particularizará para n=2 y n=3.

Se comenzará definiendo la solución fundamental de esta ecuación, pues se utilizará posteriormente para calcular el potencial newtoniano o logarítimico. Esta viene dada por [math] \phi(x) = -\frac{1}{2\pi} log(|x|) [/math] si n=2 y [math] \phi(x) = \frac{1}{4\pi |x| } [/math] si n=3.

3.1 Potencial newtoniano para [math] \mathbb{R}^3[/math]

Supongamos que [math] f(x) [/math] representa la densidad de una carga contenida en un conjunto compacto dentro del espacio tridimensional [math] \mathbb{R}^3[/math].

Entonces, la expresión [math] \phi(x-y)f(y)dy[/math] denota el potencial en el punto x. De esta manera el potencial total viene dado por la siguiente expresión:

[math] u(\mathbf{x}) = \int_{\mathbb{R}^n} f(\mathbf{y}) \phi(\mathbf{x} - \mathbf{y}) d\mathbf{y} [/math]

Esta fórmula se conoce como potencial Newtoniano de f y se aplica sobre funciones que en el infinito tienden de manera rápida a cero.

A su vez, sea [math]f \in C^2(\mathbb{R}^3)[/math] con soporte compacto. Sea [math]u[/math] el potencial newtoniano de [math]f[/math], definido por el potencial Newtoniano. Entonces, [math]u[/math] es la única solución en [math]\mathbb{R}^3[/math] de [math]\Delta u = -f[/math] que pertenece a [math]C^2(\mathbb{R}^3)[/math] y se anula en el infinito. Es decir, dado:

[math] \begin{cases} \Delta u = -f \\ u(x) \rightarrow 0 & \text{si } |x| \rightarrow \infty \end{cases}[/math]

su solución viene dada por el potencial newtoniano que se presenta anteriormente.

3.2 Potencial logarítmico para [math] \mathbb{R}^2[/math]

Si el problema se presenta en [math] \mathbb{R}^2[/math] el potencial newtoniano se sustituye por el potencial logarítimico:

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

De esta manera, Sea [math]f \in C^2(\mathbb{R}^2)[/math] con soporte compacto, [math]u(\mathbf{x})[/math] es la única solución en [math]\mathbb{R}^2[/math] de:

[math] \begin{cases} \Delta u = -f \\ u(x)=-\frac{M}{2\pi} log |x| + O\left(\frac{1}{|x|} \right) & \text{si } |x| \rightarrow \infty \end{cases}[/math]
donde
[math] M= \int _{\mathbb{R}^2} f(\mathbf y) d\mathbf y[/math]
.

3.3 Ejemplo

A continuación, veremos un ejemplo de lo explicado anteriormente para una correcta comprensión de ello. De esta manera, mediante el potencial logarítmico se aproximará la ecuación de Poisson cuando f sea la función característica de la bola de radio 1 y se estudiará su comportamiento en el infinito.


Para resolver este ejemplo se ha realizado un código en Matlab. La integral se ha simplificado de la siguiente manera:

[math] -\frac{1}{2\pi}\int_{\partial B_1 0} log|\mathbf{x} - \mathbf{y}| d\mathbf{y} = -\frac{\pi}{4}\frac{1}{2\pi}\int_{-1} ^1 \int_{-1} ^1 log|\mathbf{x} - \mathbf{y}| dx dy [/math]
%aSe divide el dominio de cada variable en 200 puntos. Para y1 e y2 se coge
%el intervalo [-1,1] pues cuando este fuera de este intervalo la función vale 0 y por tanto la integral se anulará y para x1 y x2 se coge un intervalo muy grande pues la función tiene dominio todo R2.
i1 = linspace(-1,1,200);
i2 = linspace(-1,1,200);
i3 = linspace(-10^4,10^4,200);
i4 = linspace(-10^4,10^4,200);

% Se define la función f anteriormente obtenida
f = @(y1, y2, x1, x2) -(pi/4)*(1/(2*pi))*log(sqrt((y1-x1).^2 + (y2-x2).^2));

% Generamos una malla de puntos
[A, B] = meshgrid(i1, i2);
[C, D] = meshgrid(i3,i4);

%Creamos una matriz cero donde se irán metiendo los valores que
%próximamente evaluaremos
values = zeros(size(A)); 

% Se crea el siguiente bucle para obtener los valores que necesitamos para
% la representación gráfica
for i = 1:length(i1)
    for j = 1:length(i2)
        fnum = @(x, y) f(x, y, i3(i), i4(j)); % Se crea este bucle para ir evaluando f en los distintos puntos (x1,x2)
        values(i, j) = trapz(i2, trapz(i1, fnum(A, B), 2)); % Evaluar fnum en la malla (A, B) y calcular la integral por el método del trapecio
    end
end

% Se dibujan todos los puntos obtenidos
surf(C,D,values)


La representación obtenida es la siguiente:

Representación de la solución


A continuación se va a comprobar si la solución se comporta de la manera esperada en el infinito. Tal y como se ha especificado anteriormente el comportamiento asintótico sigue la siguiente función:

[math] u(x)=-\frac{M}{2\pi} log |x| + O(\frac{1}{|x|})[/math]

siendo

[math] M= \int _{\mathbb{R}^2} f(\mathbf y) d\mathbf y[/math]


A continuación, se compara esta función con la solución obtenida anteriormente con el fin de observar si se comporta adecuadamente. Se ha creado el siguiente código: 

graf=zeros(length(i1));
xgraf=zeros(length(i1));
for i = 1:length(i1)
        fnum = @(x, y) f(x, y, i3(i), i4(i));
        
        graf(i) = trapz(i2, trapz(i1, fnum(A,B), 2));
        xgraf(i)=sqrt(2*(i3(i)^2));
end


hold on
hol= @(x) -1/2*log(x);
valy = hol(xgraf);
% Crear el gráfico de las líneas y los puntos
hold on
plot(xgraf, valy, '-','DisplayName','Función asintótica'); 
plot(xgraf, graf, '*','DisplayName','Función obtenida');
xlabel('Eje x')
ylabel('Eje y')
title('Gráfico función asintótica y función obtenida')


Se ha obtenido la siguiente gráfica, siendo la verde la función obtenida y la azul la función asintótica:

Representación de la solución

Tal y como se ve, se observa que el comportamiento es el esperado.

Obtenido de «https://mat.caminos.upm.es/w/index.php?title=Ecuaciones_de_Laplace_y_de_Poisson&oldid=71805»