Diferencia entre revisiones de «ECUACION LOGÍSTICA»

De MateWiki
Saltar a: navegación, buscar
(→‎. Ejercicio 5)
(→‎.Comportamiento asintótico)
Línea 526: Línea 526:
  
  
−
[[Archivo: Foto apar5 b.png|700px|thumb|center| pie de foto]]
+
[[Archivo: Foto apar5 b.png|700px|thumb|center| Gráfica de la solución en función de <math>\rho \in [0,10]</math> superpuesta con la gráfica de la función asintótica <math> -\frac{1}{2} \log\rho</math>.]]
  
  

Revisión del 22:00 19 abr 2024

Trabajo realizado por estudiantes
Título Ecuación de Laplace y de Poisson. Grupo ABMR
Asignatura EDP
Curso 2023-24
Autores Arturo Barrena y Mario Ríos
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura


1 . Introducción

Las ecuaciones de Laplace y Poisson son fundamentales en el campo de la física matemática y la ingeniería, especialmente en el estudio de fenómenos de difusión, electrostática y flujo de calor. Ambas ecuaciones son ecuaciones diferenciales parciales (EDP) que describen el comportamiento de campos escalares en un dominio dado. La ecuación de Laplace representa un caso especial de la ecuación de Poisson, donde la función fuente es cero.

La ecuación de Laplace se expresa matemáticamente como: [math]\Delta u = \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} + \frac{\partial^2 u}{\partial z^2}= 0[/math]

Donde [math]\phi[/math] es el campo escalar y [math]\nabla^2[/math] es el operador Laplaciano.

Por otro lado, la ecuación de Poisson se formula como: [math]\Delta u = \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} + \frac{\partial^2 u}{\partial z^2}= f[/math]

Donde f es una función fuente que puede representar, por ejemplo, densidad de carga, densidad de masa o fuentes térmicas.

2 . Ejercicio 1

En esta sección nos centraremos en el problema:

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

donde [math]B_R\subset \mathbb R^2[/math] es la bola de radio [math]R[/math] centrada en [math](0,0)[/math] y [math]g[/math] es una función cualquiera. En este caso nos enfocaremos en el caso [math]R=1[/math].

En primer lugar, compararemos su solución usando la fórmula de Poisson frente a su solución usando la serie de Fourier, pintando sus gráficas y viendo los errores que tienen al aproximar una solución exacta. Finalmente, usaremos la desigualdad de Harnack.

2.1 . Apart 1

Comenzamos con la solución dada por la fórmula de Poisson. La fórmula de Poisson nos dice que la solución a un problema de la forma anterior es:

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

donde [math]R[/math] es el radio de la bola donde se plantea el problema y [math] ||x|| [/math] la distancia al centro de la bola.

Consideramos en coordenadas polares [math](r,\theta)[/math] la función [math]g(\theta)=\max \left\{0, 1 - \frac{2}{\pi} \left| \theta - \pi \right| \right\} [/math] que impone la condición frontera del problema. La fórmula de Poisson escrita de otra forma, dice que la solución del problema escrita en coordenadas polares es:

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

Por lo tanto, sustituyendo nuestra función [math]g[/math] en coordenadas polares en esta fórmula con [math]R=1[/math], podemos hallar la solución mediante la fórmula de Poisson. Para dibujarla, aproximaremos la integral mediante el método del trapecio.

Antes de dibujarla, destacaremos el problema que surge con esta fórmula a la hora de aproximar la solución en puntos de la frontera de la forma [math](R,\theta)[/math] con [math]\theta \in [0,2\pi)[/math]. Es fácil ver que para los puntos de esta forma, cuando el coseno de la integral se hace 1 (para ciertos valores de [math]s[/math] entre [math]0[/math] y [math]2\pi[/math]), el denominador de dentro de la integral se anula. Por lo tanto, la integral diverge. Esto provoca que si intentamos dibujar la gráfica sin tomar esto en cuenta, se llega a lo siguiente:

Solución aproximada por la fórmula de Poisson en toda la bola de radio 1 usando una discretización de 300 puntos para la fórmula del trapecio.

Se puede observar claramente la irregularidad que se tiene en la frontera. Para corregir esto, simplemente basta con usar la fórmula de Poisson para estimar [math] U(r,\theta)[/math] en puntos [math](r,\theta)[/math] un poco alejados de la frontera e imponer la condición de la frontera dada por [math] g [/math] para los puntos de la frontera. Si hacemos esto (en nuestro caso estimando con la fórmula de Poisson tan sólo los puntos de radio [math]r\in [0,0.9][/math]), se tiene la siguiente gráfica:

Solución aproximada por la fórmula de Poisson en la bola de radio 0.9 (para no llegar a la frontera) e imponiendo la condición en la frontera de la bola de radio 1 manualmente. Se ha usado una discretización de 300 puntos para la fórmula del trapecio.

En esta gráfica se puede ver cómo se ha corregido un poco la irregularidad en la frontera.

2.1.1 . Código

% Limpiar el espacio de trabajo y cerrar todas las figuras previas
clear all
close all

% Datos
R = 1;                                            % Radio del círculo
g = @(theta) max(0, 1 - 2/pi * abs(theta - pi));  % Condición inicial
numvalores_theta = 400;                            % Número de valores de el ángulo theta
numvalores_r = 400;                                % Número de valores de el radio r
numvalores_trapecio = 400;                         % Número de valores al discretizar para la fórmula del trapecio

% Función a integrar en la fórmula de Poisson
integrando = @(s, r, theta) g(s) ./ (R^2 + r^2 - 2 * R * r * cos(theta - s));

% Generación de puntos de evaluación
theta = linspace(0, 2*pi, numvalores_theta);        % Vector de valores de el ángulo theta
r = linspace(0, R-R/10, numvalores_r);              % Vector de valores de el radio r (sin llegar a la frontera) 
s = linspace(0, 2*pi, numvalores_trapecio);         % Discretización fórmula del trapecio

% Cálculo de la solución en función de r y theta con la fórmula de Poisson (sin calcularla 
% en los valores en los que r=R, la frontera)
U = zeros(length(r), length(theta)); % Inicialización de matriz de valores de la solución
for i = 1:length(r)-1 
    for j = 1:length(theta)
        % Integral numérica para calcular la solución en cada punto r(i),theta(j)
        U(i, j) = (R^2 - r(i)^2) / (2 * pi) * trapz(s, integrando(s, r(i), theta(j)));
    end
end

% Valores de la frontera puestos manualmente
for j = 1:length(theta)
    U(length(r),j)= g(theta(j));
end

% Gráficas
X = r' * cos(theta); % Coordenadas x tras deshacer polares
Y = r' * sin(theta); % Coordenadas y tras deshacer polares
surf(X, Y, U, 'EdgeColor', 'none'); % Gráfica de la solución en función de x e y
xlabel('x'); ylabel('y'); zlabel('u(x,y)'); % Etiquetas de ejes
legend('u(x,y)'); % Leyenda


2.2 . Apartado 2

Tal y como hemos observado, la fórmula de Poisson tiene ciertos problemas al aproximar la solución en la frontera debido al carácter singular de la integral. Estos problemas se pueden estudiar con mayor precisión calculando el error que tiene al aproximar una solución exacta. Consideramos por ejemplo la función en coordenadas cartesianas [math]u(x,y)=xy[/math], la cual es armónica y su valor en la frontera es [math]g(x,y)=xy[/math]. Para usar la fórmula de Poisson, escribimos [math]u[/math] en coordenadas polares:

[math] \begin{cases} U(r,\theta)=r^2\cos(\theta)\sin(\theta), r\in[0,1), \theta\in [0,2\pi)\\ U(1,\theta)=G(\theta)=\cos(\theta)\sin(\theta), \theta\in [0,2\pi) \end{cases} [/math]

Si usamos de nuevo la fórmula de Poisson para calcular la solución aproximada haciendo uso de la fórmula del trapecio, podemos compararla con la solución exacta y obtener el error de aproximación que se tiene.

En primer lugar, vamos a observar el error de la fórmula en puntos "alejados" de la frontera. Para estudiar este error tomamos varias discretizaciones en la fórmula del trapecio y dibujamos los distintos errores en una gráfica. En este caso hemos dibujado en una gráfica la función:

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

donde [math]Error(n)[/math] es el error en un punto "alejado" de la frontera (en nuestro caso [math](r,\theta)=(0.9,\pi/4)[/math] es el que hemos escogido) cuando se toman [math]10^n[/math] puntos en la discretización de la fórmula del trapecio. Tal y como vimos visualmente en las gráficas anteriores, la fórmula en el punto escogido "alejado" de la frontera debería aproximar bien a la solución exacta y por lo tanto el error debería ir disminuyendo según aumenta el número de puntos en la discretización de la fórmula del trapecio (lo que hace a la fórmula del trapecio más exacta). La gráfica de esta función tiene la forma:

Gráfica del error (en escala logarítmica) en el punto [math](r,\theta)=(0.9,\pi/4)[/math] cuando se toman [math]10^n[/math] puntos en la discretización de la fórmula del trapecio, en función de n.

Justo como esperábamos, el error disminuye hasta cierto punto en el cual se estabiliza en el cual tiene un valor de alrededor de [math]10^{-15}[/math]. Con esto concluimos que la fórmula de Poisson aproxima bastante bien los puntos relativamente "alejados" de la frontera.

Usando la fórmula del error máximo del método del trapecio podemos obtener una cota para el error y confirmar que los errores anteriores están por debajo de esta cota. Esta fórmula dice lo siguiente:

[math] Error Máximo=\left|-\frac{f’’(\xi)(b-a)^3}{12n^2} \right| [/math]

siendo [math] \xi [/math] un valor perteneciente al intervalo [math] [b-a] [/math] (intervalo de integración), [math]f[/math] la función a integrar y [math]n[/math] el número de subintervalos tomados en la discretización del método del trapecio. En este caso, tenemos [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].

Usando la fórmula, podemos obtener una cota para el error máximo, maximizando su valor para [math] \xi \in [0,2\pi][/math]:

[math]Error Máximo \leq \left|- \frac{\pi^2(r^6 - r^5 + 6 r^3 + r^2 -r +2) }{3 n^2(1-r)^5} \right|[/math]

Es fácil ver que esta cota es más pequeña a medida que se aumentan el número de subintervalos en la discretización de la fórmula del trapecio (ya que [math]n[/math] se encuentra en el denominador), lo cual indica que en efecto el error decrece al aumentar el número de puntos en la discretización. Una vez obtenida esta cota, veamos que los errores anteriores en efecto están por debajo de esta cota:

Gráfica del error (en escala logarítmica) en el punto [math](r,\theta)=(0.9,\pi/4)[/math] cuando se toman [math]10^n[/math] puntos en la discretización de la fórmula del trapecio, en función de n.

Podemos ver claramente que esto ocurre y la cota del error decrece con el valor de [math]n[/math] tal y como habíamos visto.

Ahora, veamos qué tal funciona la fórmula para puntos "cercanos" a a frontera. Para ello, vamos a estudiar el error en distintos puntos "cercanos" a la frontera con una discretización fija de la fórmula del trapecio y vamos a observar gráficamente cómo se comporta a medida que dichos puntos se acercan a la frontera. En este caso hemos dibujado en una gráfica el error en función de [math]n[/math] donde para cada [math]n[/math], [math]Error(n)[/math] es el error en el punto [math](r,\theta)=(1-10^{-n}, \pi/4)[/math] dado por la fórmula de Poisson con una discretización en la fórmula del trapecio de [math]100[/math] puntos. La gráfica queda de la siguiente forma:

Gráfica del error en el punto [math](r,\theta)=(1-10^{-n}, \pi/4)[/math] en función de n con una discretización fija de 100 puntos.

En la gráfica se puede ver cómo el error aumenta a medida que los puntos se acercan a la frontera hasta llegar a cierto punto en el cual se estabiliza, manteniéndose constante igual a [math]0,5[/math].

2.2.1 .Código

% Limpiar el espacio de trabajo y cerrar todas las figuras previas
clear all
close all

% Datos
R = 1;                                            % Radio del círculo
g = @(theta) R^2*cos(theta).*sin(theta);          % Condición inicial
numdiscretizaciones_trapecio=8;                   % Número de discretizaciones al usar la fórmula del trapecio
r=0.9; theta=pi/4;                                % Punto particular en el que se quiere ver el error
sol=r^2*cos(theta)*sin(theta);                    % Solución exacta en dicho punto particular


% Función a integrar en la fórmula de Poisson
integrando = @(s, r, theta) g(s) ./ (R^2 + r^2 - 2 * R * r * cos(theta - s));

error=zeros(1,numdiscretizaciones_trapecio);
for n=1:numdiscretizaciones_trapecio
    numvalores_trapecio = 10^n;  % Número de valores en cada discretización para la fórmula del trapecio
    % Generación de puntos de evaluación
    s = linspace(0, 2*pi, numvalores_trapecio);         % Discretización para la fórmula del trapecio

    % Cálculo de la solución en el punto particular (r, theta) con la fórmula de Poisson
    u = (R^2 - r^2) / (2 * pi) * trapz(s, integrando(s, r, theta));
    
    error(n)=abs(u-sol);                % Error
end

% Gráficas
plot(1:numdiscretizaciones_trapecio,log10(error)); % Gráfica de los errores en función de n
xlabel('n'); ylabel('Log_{10}(Error)'); % Etiquetas de ejes
legend('Log_{10}(Error)'); % Leyenda

2.2.2 .Código

% Limpiar el espacio de trabajo y cerrar todas las figuras previas
clear all
close all

% Datos
R = 1;                                            % Radio del círculo
g = @(theta) R^2*cos(theta).*sin(theta);          % Condición inicial
numdiscretizaciones_trapecio=7;                   % Número de discretizaciones al usar la fórmula del trapecio
r=0.9; theta=pi/4;                                % Punto particular en el que se quiere ver el error
sol=r^2*cos(theta)*sin(theta);                    % Solución exacta en dicho punto particular


% Función a integrar en la fórmula de Poisson
integrando = @(s, r, theta) g(s) ./ (R^2 + r^2 - 2 * R * r * cos(theta - s));

cota=zeros(1,numdiscretizaciones_trapecio);
error=zeros(1,numdiscretizaciones_trapecio);
for n=1:numdiscretizaciones_trapecio
    numvalores_trapecio = 10^n;  % Número de valores en cada discretización para la fórmula del trapecio
    % Generación de puntos de evaluación
    s = linspace(0, 2*pi, numvalores_trapecio);         % Discretización para la fórmula del trapecio

    % Cálculo de la solución en el punto particular (r, theta) con la fórmula de Poisson
    u = (R^2 - r^2) / (2 * pi) * trapz(s, integrando(s, r, theta));
    
    cota(n)=abs(-pi^2*(r^6-r^5+6*r^3+r^2-r+2)/(3*numvalores_trapecio^2*(1-r)^5));
    error(n)=abs(u-sol);                % Error
end

% Gráficas
plot(1:numdiscretizaciones_trapecio,log10(error)); % Gráfica de los errores en función de n
hold on
plot(1:numdiscretizaciones_trapecio,log10(cota));  % Gráfica de las cotas en función de n
xlabel('n'); ylabel('Log_{10}(Error)'); % Etiquetas de ejes
legend('log_{10}(Error)', 'log_{10}(Cota)'); % Leyenda
hold off


2.2.3 .Código

% Limpiar el espacio de trabajo y cerrar todas las figuras previas
clear all
close all

% Datos
R = 1;                                    % Radio del círculo
g = @(theta) R^2*cos(theta).*sin(theta);  % Condición inicial
numvalores_trapecio = 100;                % Número de valores al discretizar para la fórmula del trapecio
numpuntos_cercafrontera=100;              % Número de puntos cercanos a la frontera en los que se verá el error
theta=pi/4;                               % Valor de theta particular en el que se quiere ver el error


% Función a integrar en la fórmula de Poisson
integrando = @(s, r, theta) g(s) ./ (R^2 + r^2 - 2 * R * r * cos(theta - s));

% Generación de puntos de evaluación
s = linspace(0, 2*pi, numvalores_trapecio);         % Discretización fórmula del trapecio

error=zeros(1,numpuntos_cercafrontera);
for n=1:numpuntos_cercafrontera
    r=1-10^(-n);
    % Cálculo de la solución en el punto particular (r, theta) con la fórmula de Poisson
    u = (R^2 - r^2) / (2 * pi) * trapz(s, integrando(s, r, theta));
    
    error(n)=abs(u-r^2*cos(theta)*sin(theta));
end

% Gráficas
plot(1:numpuntos_cercafrontera,error); % Gráfica de los errores en función de n
xlabel('n'); ylabel('Error'); % Etiquetas de ejes
ylim([0,0.7])                   % Intervalo de valores en el eje y
legend('Error'); % Leyenda


3 . Apartado 3

Una vez estudiado el problema que se tiene con la fórmula de Poisson para aproximar la solución cerca de la frontera, veamos qué ocurre con la solución obtenida por la serie de Fourier. Empezaremos considerando que el problema se ha planteado en la bola de radio [math]R[/math] y luego particularizaremos la solución para [math]R=1[/math]. Para ello, transformamos nuestro problema a 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,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]

Una vez transformado, empleamos el método de separación de variables, obteniendo así que:

[math] U(r,\theta)=\frac{a_0}{2}+\sum_{k=1}^\infty a_kr^kcos(k\theta)+b_kr^ksin(k\theta).[/math]

Finalmente, imponemos que se cumpla la condición frontera haciendo el desarrollo en serie de fourier de [math] G(\theta) [/math] en [math] L^2([-\pi,\pi]) [/math]:

[math] U(1,\theta)=G(\theta)=\frac{\alpha_0}{2}+\sum_{k=1}^\infty \alpha_k cos(k\theta)+\beta_k sin(k\theta).[/math]

de forma que:

[math] a_0=\alpha_0=\frac{1}{\pi}\int_{-\pi}^\pi G(\theta)d\theta. [/math]
[math] a_k=\frac{\alpha_k}{R^k}=\frac{\frac{1}{\pi}\int_{-\pi}^\pi G(\theta)cos(k\theta)d\theta}{R^k}. [/math]
[math] b_k=\frac{\beta_k}{R^k}=\frac{\frac{1}{\pi}\int_{-\pi}^\pi G(\theta)sin(k\theta)d\theta}{R^k}. [/math]

Pero los coeficientes de Fourier de [math]G(\theta)=R^2\cos(\theta)\sin(\theta)=\frac{R^2sin(2\theta)}{2}[/math] son claramente todos [math]0[/math] menos [math] \beta_2=\frac{R^2}{2}[/math], por tanto, la solución obtenida por la serie de Fourier es la siguiente:

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

la cual no depende del radio [math]R[/math] (sólo depende de él su dominio de definición). En este caso no hace falta ni siquiera estudiar los errores debido a que la solución obtenida por la serie de Fourier es la propia solución exacta. En efecto, pasando a cartesianas obtenemos que [math]u(x,y)=xy[/math]. Aun así se dejará a continuación una gráfica de dicha función para [math]R=1[/math]:

Gráfica de la función [math]u(x,y)=xy[/math] en la bola de radio 1.

3.1 .Código

clear all
close all
% Datos
R = 1;                                            % Radio del círculo
numvalores_theta = 300;                            % Número de valores de el ángulo theta
numvalores_r = 300;                                % Número de valores de el radio r

% Definición de la función
u=@(x,y) x.*y;

% Generación de puntos de evaluación
theta = linspace(0, 2*pi, numvalores_theta);        % Vector de valores de el ángulo theta
r = linspace(0, R, numvalores_r);                   % Vector de valores de el radio r
X = r' * cos(theta);                                % Coordenadas x tras deshacer polares
Y = r' * sin(theta);                                % Coordenadas y tras deshacer polares

% Gráfica
surf(X, Y, u(X,Y), "FaceColor", "interp",'EdgeColor', 'none')      % Gráfica
xlabel('x'); ylabel('y'); zlabel('u(x,y)')                          % Etiquetas de ejes
legend('Solución exacta');                                          % Leyenda


4 . Apartado 4

En esta sección analizaremos qué dice la desigualdad de Harnack. Usaremos dicha desigualdad para encontrar la región en la que deben estar todas las soluciones armónicas en la [math]B_R[/math] (para [math]R=1,2,10[/math]), que valen lo mismo que u en el punto [math](0,0)[/math].

La desigualdad de Harnack establece lo siguiente: sea [math] u [/math] una función armónica tal que [math] u \geq 0 [/math] en un dominio [math]\Omega \subset \mathbb{R}^n[/math]. Suponemos [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), r=||z-x||.[/math]

Vamos a comenzar hallando el mínimo de la función [math] g(x,y) = xy [/math] en la frontera de la bola [math] \partial B_R [/math], al cual llamaremos [math]M[/math]. Si expresamos la función en la frontera en coordenadas polares tenemos [math] g(\theta)=R^{2}cos(\theta)sen(\theta)[/math], de modo que:

[math]M := \min_{(x, y) \in B_R} \{g(x,y)\} = \min_{\theta \in [0,2\pi)} \{g(\theta)\} = \min_{\theta \in [0,2\pi)} \{R^2 cos(\theta)sen(\theta)\}= \min_{\theta \in [0,2\pi)} \{R^2 \frac{\sin(2\theta)}{2}\} =\frac{R^2}{2} \min_{\theta \in [0,2\pi)} \{\sin(2\theta)\}= - \frac{R^2}{2} .[/math]


Sea [math]u[/math] la solución al problema en [math] B_R [/math] en coordenadas cartesianas. Definimos la función:

[math] v := u-M = u+\frac{R^2}{2} \geq 0.[/math]
.

Como [math]v \geq 0 [/math] podemos aplicarle la desigualdad de Harnack para [math] n=2 [/math], de forma que se tiene:

[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] con [math] r =||(x,y)||.[/math]

para un radio [math]R[/math] genérico. Teniendo en cuenta que [math] u(0,0)=0 [/math] tendremos que [math] v(0,0)= \frac{R^2}{2} [/math]. Sabiendo esto se llega a:

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

Por lo tanto, sabiendo que [math] v = u+\frac{R^2}{2}[/math], tenemos que la región en la que deben estar todas las soluciones armónicas en la [math]B_R[/math], que valen [math]0[/math] en el punto [math](0,0)[/math] es:

[math] \frac{R^2(R-r)}{2(R+r)} - \frac{R^2}{2} = -\frac{R^2r}{R+r}\leq u(x, y) \leq \frac{R^2r}{R-r}=\frac{R^2(R+r)}{2(R-r)}- \frac{R^2}{2}, \forall (x,y) \in B_{R}, [/math] con [math] r =||(x,y)||.[/math]

Como se puede observar, una vez fijado el radio de la bola [math]R[/math] las cotas solo dependen de [math]r[/math]. Sin embargo, a la hora de representar estas cotas visualmente, al ser la cota inferior negativa no es posible usar escala logarítmica. Por lo tanto, para estudiar el comportamiento de las cotas al variar el [math]R[/math] gráficamente, usaremos las cotas de la función [math]v[/math]. Si representamos las funciones cota [math]v_1(r):=\frac{(R-r)}{(R+r)}[/math] y [math]v_2(r):= \frac{R^2(R+r)}{2(R-r)}[/math] en escala logarítmica en función de [math]r\in [0,R][/math] para varios valores de [math]R[/math] (en este caso [math]R=1,2,10[/math]) se tienen las siguientes gráficas:

Gráficas de las funciones cota [math]v_1[/math] y [math]v_2[/math] en escala logarítmica en función del radio [math]r[/math] para distintos valores [math]R=1,2,10[/math] por separado.
Gráficas de las funciones cota [math]v_1[/math] y [math]v_2[/math] en escala logarítmica en función del radio [math]r[/math] para distintos valores [math]R=1,2,10[/math] superpuestas.

En estas gráficas se puede ver que el valor de las cotas superior e inferior va aumentando a medida que aumentamos el radio [math]R[/math], desplazándose ambas en el sentido positivo del eje y. Además, la región que abarcan ambas cotas va haciéndose más pequeña (se va estrechando) a medida que se aumenta el radio [math]R[/math], tal y como se ve en la segunda figura. Esta región, tal y como vimos en clase se va estrechando cada vez más y más, de tal manera que si tomamos el límite cuando [math]R[/math] tiende a infinito ambas cotas serán iguales a [math]v(0,0)=\frac{R}{2}[/math] (teorema de Liouville).

Si pasamos el problema a dimensión 3 (teniendo en cuenta que el mínimo no cambia ya que se tiene la misma función [math]g(x,y,z)=xy[/math]), las cotas serían:

[math]v_1(r):= \frac{R^3(R-r)}{2(R+r)^{2}} \leq v(x,y,z) \leq \frac{R^3(R+r)}{2(R-r)^2} =: v_2(r), \forall (x,y,z) \in B_{R}, [/math] con [math] r =||(x,y,z)||.[/math]

Por lo tanto, sabiendo que [math] v = u+\frac{R^2}{2}[/math], tenemos que la región en la que deben estar todas las soluciones armónicas en la [math]B_R[/math], que valen [math]0[/math] en el punto [math](0,0,0)[/math] es:

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

Representando de nuevo las funciones cota [math]v_1[/math] y [math]v_2[/math] (definidas de forma análoga) en escala logarítmica en función de [math]r\in [0,R][/math] para varios valores de [math]R[/math] (en este caso [math]R=1,2,10[/math]) junto con las gráficas anteriores se llega a lo siguiente:

Gráficas de las funciones cota [math]v_1[/math] y [math]v_2[/math] en dimensión 2 y 3 superpuestas en escala logarítmica en función del radio [math]r[/math], para distintos valores [math]R=1,2,10[/math] por separado.
Gráficas de las funciones cota [math]v_1[/math] y [math]v_2[/math] en dimensión 2 y 3 superpuestas en escala logarítmica en función del radio [math]r[/math], para distintos valores [math]R=1,2,10[/math] superpuestas.

Se tienen las mismas condiciones que para el caso en dimensión 2. El valor de las cotas superior e inferior va aumentando a medida que aumentamos el radio [math]R[/math], desplazándose ambas en el sentido positivo del eje y la región que abarcan las cotas se estrecha a medida que se aumenta [math]R[/math]. Sin embargo, tal y como se puede ver en la segunda figura, la región abarcada es mayor para dimensión 3. Esto se debe a que las nuevas cotas son las mismas que en dimensión 2 pero multiplicadas cada una por un factor (en el caso de [math]v_1[/math] se multiplica por [math]\frac{R}{R+r}\leq 1[/math] haciéndola menor y en el caso de [math]v_2[/math] se multiplica por [math]\frac{R}{R-r} \geq 1[/math] haciéndola mayor) que hace que la región se extienda.

4.1 .Código

clear all
close all

% Datos
VectR = [1,2,10];                                  % Radios de la bola
numvalores_r = 300;                                % Número de valores de el radio r

figure(1)
for i=1:length(VectR)
    R=VectR(i);       % Radio de la bola

    % Definición de las cotas en función de r
    v1_2=@(r) R*(R-r)./(2*(R+r));                       %Cota 1, dimensión 2
    v2_2=@(r) R*(R+r)./(2*(R-r));                       %Cota 2, dimensión 2
    
    % Generación de puntos de evaluación
    r = linspace(0, R, numvalores_r);                   % Vector de valores de el radio r
    
    % Gráfica
    subplot(1,length(VectR),i)
    plot(r,log10(v1_2(r)))                               % Gráfica v1, dimensión 2
    hold on
    plot(r,log10(v2_2(r)))                               % Gráfica v2, dimensión 2

    xlabel('r'); ylabel('log_{10}(v_i(r))');            % Etiquetas de ejes
    legend('Cota inferior, n=2','Cota superior, n=2');  %Leyenda                   
    title("R = "+ num2str(VectR(i)))                    % Título
    hold off
end

figure(2)
for i=1:length(VectR)
    R=VectR(i);       % Radio de la bola

    % Definición de las cotas en función de r
    v1_2=@(r) R*(R-r)./(2*(R+r));                       %Cota 1, dimensión 2
    v2_2=@(r) R*(R+r)./(2*(R-r));                       %Cota 2, dimensión 2

    v1_3=@(r) v1_2(r)*R./(R+r);                            %Cota 1, dimensión 3
    v2_3=@(r) v2_2(r)*R./(R-r);                            %Cota 2, dimensión 3
    
    % Generación de puntos de evaluación
    r = linspace(0, R, numvalores_r);                   % Vector de valores de el radio r
    
    % Gráfica
    subplot(1,length(VectR),i)
    plot(r,log10(v1_2(r)))                               % Gráfica v1, dimensión 2
    hold on
    plot(r,log10(v2_2(r)))                               % Gráfica v2, dimensión 2
    plot(r,log10(v1_3(r)))                               % Gráfica v1, dimensión 3
    plot(r,log10(v2_3(r)))                               % Gráfica v2, dimensión 3

    xlabel('r'); ylabel('log_{10}(v_i(r))');            % Etiquetas de ejes
    legend('Cota inferior, n=2','Cota superior, n=2','Cota inferior, n=3','Cota superior, n=3');   %Leyenda                 
    title("R = "+ num2str(VectR(i)))                    % Título
    hold off
end


4.2 .Código

clear all
close all

% Datos
VectR = [1,2,10];                                  % Radios de la bola
numvalores_r = 300;                                % Número de valores de el radio r

colores = lines(length(VectR)); %Paleta de colores para las gráficas
figure(1)
for i=1:length(VectR)
    R=VectR(i);       % Radio de la bola

    % Definición de las cotas en función de r
    v1_2=@(r) R*(R-r)./(2*(R+r));                       %Cota 1, dimensión 2
    v2_2=@(r) R*(R+r)./(2*(R-r));                       %Cota 2, dimensión 2
    
    % Generación de puntos de evaluación
    r = linspace(0, R, numvalores_r);                   % Vector de valores de el radio r
    
    % Gráficas

    plot(r,log10(v1_2(r)),'Color', colores(i,:))                     % Gráfica v1, dimensión 2
    hold on
    plot(r,log10(v2_2(r)), 'Color', colores(i,:))                    % Gráfica v2, dimensión 2
end
xlabel('r'); ylabel('log_{10}(v_i(r))');            % Etiquetas de ejes   
legend('Cotas R=1','','Cotas R=2','','Cotas R=10','');  %Leyenda      
hold off
figure(2)
for i=1:length(VectR)
    R=VectR(i);       % Radio de la bola

    % Definición de las cotas en función de r
    v1_2=@(r) R*(R-r)./(2*(R+r));                       %Cota 1, dimensión 2
    v2_2=@(r) R*(R+r)./(2*(R-r));                       %Cota 2, dimensión 2

    v1_3=@(r) v1_2(r)*R./(R+r);                         %Cota 1, dimensión 3
    v2_3=@(r) v2_2(r)*R./(R-r);                         %Cota 2, dimensión 3
    
    % Generación de puntos de evaluación
    r = linspace(0, R, numvalores_r);                   % Vector de valores de el radio r
    
    % Gráfica
    plot(r,log10(v1_2(r)), 'r')                               % Gráfica v1, dimensión 2
    hold on
    plot(r,log10(v2_2(r)), 'r')                               % Gráfica v2, dimensión 2
    plot(r,log10(v1_3(r)), 'b')                               % Gráfica v1, dimensión 3
    plot(r,log10(v2_3(r)), 'b')                               % Gráfica v2, dimensión 3               
end
xlabel('r'); ylabel('log_{10}(v_i(r))');   % Etiquetas de ejes 
legend('Cotas n=2','','Cotas n=3');        %Leyenda    
hold off


5 . Ejercicio 5

Se define la ecuación de Poisson en [math]\mathbb R^n[/math] mediante la siguiente expresión:

[math] \Delta u = f [/math]
donde [math]u:\mathbb{R}^n \rightarrow \mathbb{R}[/math] y [math] f \in C^2(\mathbb{R})[/math].

5.1 Solución fundamental del Laplaciano para dimensión 2.

Se define la solución fundamental del Laplaciano para dimensión 2 como la función [math] \Phi: \mathbb{R}^2 \rightarrow \mathbb{R}[/math] tal que:

[math] \Phi(x) =-\frac{1}{2\pi} log|x|.[/math]

5.2 El potencial logarítmico .

En este apartado buscamos aproximar la ecuación de Poisson [math] \Delta u = f [/math] utilizando la solución fundamental teniendo en cuenta que [math]f[/math] es la función característica de la bola de radio [math]1[/math]. Dicha función tiene la siguiente expresión:

[math]f(x) = \begin{cases} 1, & \text{si } |x| \leq 1 \\ 0, & \text{si } |x| \gt 1 \end{cases}[/math]

Sabemos que la solución de Poisson viene dada por el producto de convolución:

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

Es decir, la solución será:

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

la cual también llamaremos potencial logarítmico. Haciendo el cambio de variable a polares [math] \textbf{y}= ( r \cos (\theta), r \sin( \theta)) [/math] y [math] \textbf{x} = (\rho \cos(\beta), \rho \sin (\beta))[/math] se tiene que la solución es:

[math] u(\rho, \beta) = - \frac{1}{4 \pi} \int _{0}^{2 \pi} \int_{0}^{1} r \log(r^2 + \rho^2 – 2 r \rho \cos(\theta - \beta)) d r d \theta [/math]

Nótese que en la expresión anterior no depende de beta ( se puede ver que su derivada con respecto a beta es 0) de manera que podemos obtener finalmente la siguiente expresión:

[math] u(\rho, \beta) = - \frac{1}{4 \pi} \int _{0}^{2 \pi} \int_{0}^{1} r \log(r^2 + \rho^2 – 2 r \rho \cos(\theta)) d r d \theta [/math]

Si dibujamos esta función para [math]\beta \in [0,2\pi], \rho \in [0,10][/math] se tiene lo siguiente:

Gráfica de la solución para [math]\beta \in [0,2\pi], \rho \in [0,10][/math] en función de x e y.

5.3 .Comportamiento asintótico

Ahora estudiaremos el comportamiento del potencial logarítmico en el infinito. Su comportamiento asintótico es :

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

a medida que [math] |x| \rightarrow +\infty [/math]. Donde definimos [math]M[/math] como:

[math] M := \int_{\mathbb{R}^2} f(y) \, dy [/math]

que aplicado a nuestro problema inicialmente planteado con \( f \) la función característica de la bola de radio 1, tendríamos que:


[math] M := \int_{\mathbb{R}^2} f(y) \, dy = \int_{B_1} 1 \, dy = \pi [/math]

De modo que tenemos que:

[math] \lim_{{|x|\to\infty}} u(x) = -\frac{1}{2} \log|x| [/math]

Es decir, la solución obtenida se comporta en el infinito como la función [math]-\frac{1}{2} \log|x|[/math]. Veamos que esto es cierto.

Para verlo tomaremos valores muy altos de [math] |x|=\rho [/math] y compararemos [math]u(\rho,\beta)[/math] con [math] -\frac{1}{2} \log\rho[/math]. Si dibujamos la gráfica de la solución en función de [math]\rho[/math] junto con la gráfica de [math] -\frac{1}{2} \log\rho[/math] para [math]\rho \in [0,10][/math] observamos lo siguiente:


Gráfica de la solución en función de [math]\rho \in [0,10][/math] superpuesta con la gráfica de la función asintótica [math] -\frac{1}{2} \log\rho[/math].


Tal y como se puede observar, a partir de un radio [math]\rho[/math] entre [math]1[/math] y [math]2[/math], las gráficas se superponen y a partir de ahí se mantienen superpuestas, justo como era de esperar.

5.4 .Código

clear all
close all
% Datos
R = 1;                                             % Radio del círculo
numvalores_theta = 300;                            % Número de valores de la discretización del ángulo theta para la fórmula del trapecio
numvalores_r = 300;                                % Número de valores de la discretización del radio r para la fórmula del trapecio
numvalores_ro = 300;                               % Número de valores del radio ro
numvalores_beta = 300;                             % Número de valores del ángulo beta
Rofinal=10;                                        % Valor hasta el cual llegan los valores del radio ro

% Definición de la función
integrando=@(ro,r,theta)-1/(4*pi)*r.*log(r.^2+ro.^2-2*r.*ro.*cos(theta)) ;

% Generación de puntos de evaluación
theta = linspace(0, 2*pi, numvalores_theta);        % Vector de valores del ángulo theta
r = linspace(0, R-0.000001, numvalores_r);                   % Vector de valores del radio r
ro = linspace(0, Rofinal, numvalores_ro);                 % Vector de valores del radio ro
beta = linspace(0, 2*pi, numvalores_beta);          % Vector de valores del ángulo beta

%Integral
u=zeros(1,length(ro));                                % Vector con los valores de la función u en función de ro
for i=1:length(ro)
    integrando1=zeros(1,length(theta));               % Vector que acumula los valores de la primera integral para distintos theta
    for j=1:length(theta)
        integrando1(j)=trapz(r,integrando(ro(i),r,theta(j))); % Primera integral
    end
    u(i)=trapz(theta,integrando1);  % Segunda integral
end

U=u'*ones(1,length(beta));                          % Matriz con los valores de la función u en función de ro y beta
X = ro' * cos(beta);                                % Coordenadas x tras deshacer polares
Y = ro' * sin(beta);                                % Coordenadas y tras deshacer polares
% Gráfica
surf(X, Y, U, "FaceColor", "interp",'EdgeColor', 'none')      % Gráfica
xlabel('x'); ylabel('y'); zlabel('u(x,y)')                   % Etiquetas de ejes
legend('Solución');                                          % Leyenda


5.5 .Código

clear all
close all
% Datos
R = 1;                                             % Radio del círculo
numvalores_theta = 300;                            % Número de valores de la discretización del ángulo theta para la fórmula del trapecio
numvalores_r = 300;                                % Número de valores de la discretización del radio r para la fórmula del trapecio
numvalores_ro = 300;                               % Número de valores del radio ro
numvalores_beta = 300;                             % Número de valores del ángulo beta
Rofinal=10;                                        % Valor hasta el cual llegan los valores del radio ro

% Definición de la función
integrando=@(ro,r,theta)-1/(4*pi)*r.*log(r.^2+ro.^2-2*r.*ro.*cos(theta)) ;

% Generación de puntos de evaluación
theta = linspace(0, 2*pi, numvalores_theta);        % Vector de valores del ángulo theta
r = linspace(0, R-0.000001, numvalores_r);                   % Vector de valores del radio r
ro = linspace(0, Rofinal, numvalores_ro);                 % Vector de valores del radio ro
beta = linspace(0, 2*pi, numvalores_beta);          % Vector de valores del ángulo beta

%Integral
u=zeros(1,length(ro));                                % Vector con los valores de la función u en función de ro
for i=1:length(ro)
    integrando1=zeros(1,length(theta));               % Vector que acumula los valores de la primera integral para distintos theta
    for j=1:length(theta)
        integrando1(j)=trapz(r,integrando(ro(i),r,theta(j))); % Primera integral
    end
    u(i)=trapz(theta,integrando1);  % Segunda integral
end

% Gráfica
plot(ro,u) %Gráfica solución
hold on
plot(ro,-1/2*log(ro),'--') %Gráfica función asintótica
xlabel('x'); ylabel('y');                   % Etiquetas de ejes
legend('Solución', 'Función asintótica');   % Leyenda
hold off


5.6 Referencias

Obtenido de «https://mat.caminos.upm.es/w/index.php?title=ECUACION_LOGÍSTICA&oldid=71758»