Revisión del 20:17 27 may 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

Contenido

1 Introducción
2 Planteamiento del problema
3 Modelización de los desplazamientos transversales de la cuerda
4 Particularización del sistema de ecuación de ondas a unos datos iniciales dados
5 Condiciones iniciales de una onda que viaja en un único sentido
- 5.1 Código
6 Condiciones de frontera de Neumann
- 6.1 Código
7 Representación solución fundamental de dimensión 3 y regularización para la solución fundamental de dimensión 2
8 Solución fundamental de dimensión 2 mediante el producto de convolución
- 8.1 Código
- 8.2 Código

1 Introducción

En este trabajo se estudiará la ecuación que modeliza una onda, llamada la ecuación de ondas.

En primer lugar, se buscarán soluciones para esta ecuación en el caso de una dimensión. Se especificarán para ello algunos datos iniciales como el intervalo que ocupa la cuerda y su velocidad de propagación, además de sus condiciones frontera de tipo Dirichlet, aunque la solución inicialmente se dejará en función de las condiciones iniciales. A continuación, se elegirán unas condiciones iniciales específicas y se dibujará el aspecto de la solución, comentando algunas cosas interesantes de recalcar. En la siguiente sección, se elegirán unas nuevas condiciones iniciales correspondientes a una onda que viaja en un solo sentido y se repetirá la representación y el análisis de nuevo. Finalmente, se concluirá cambiando las condiciones frontera de tipo Dirichlet por unas de tipo Neumann, para observar qué cambios se producen.

En la segunda parte, se representarán las soluciones fundamentales de la ecuación de ondas para dimensión 1,2 y 3 y se usará la de dimensión 2 para representar la solución de un problema en dimensión 2 haciendo uso de la convolución.

2 Planteamiento del problema

El problema que se va a estudiar inicialmente es el de una cuerda vibrante que ocupa el intervalo [math][0,1][/math] con densidad [math]d[/math] y tensión [math]\tau_0[/math] constante de manera que la velocidad de propagación es [math]c = \tau_0/d = 1[/math]. Además, supondremos que la cuerda está fija en los extremos. Llamaremos [math]u_0(x)[/math] y [math]u_1(x)[/math] su posición e impulso iniciales respectivamente. Vamos a plantear el sistema de ecuaciones que modeliza el comportamiento de los desplazamientos transversales de la cuerda, describiendo cada término.

2.1 Ecuación de ondas

La ecuación de ondas en una dimensión, con una cuerda de densidad \(d\) y tensión \(\tau_0\), donde la velocidad de propagación es \(c = \tau_0/d\), se escribe como:

[math] \frac{\partial^2 u}{\partial t^2} = c \frac{\partial^2 u}{\partial x^2} .[/math]

Dado que en este caso \(c = 1\), la ecuación se simplifica a:

[math] \frac{\partial^2 u}{\partial t^2} = \frac{\partial^2 u}{\partial x^2} .[/math]

2.2 Condiciones de frontera

Dado que la cuerda está fija en los extremos, tenemos:

[math] u(0, t) = 0 .[/math]

[math] u(1, t) = 0 .[/math]

Estas condiciones indican que la posición de la cuerda en los puntos \(x = 0\) y \(x = 1\) siempre es cero, es decir, la cuerda no se mueve en los extremos.

2.3 Condiciones iniciales

Las condiciones iniciales especifican la posición e impulso inicial de la cuerda:

Posición inicial:

[math] u(x, 0) = u_0(x) .[/math]

Esto describe la forma inicial de la cuerda en \(t = 0\).

Velocidad inicial:

[math] \frac{\partial u}{\partial t}(x, 0) = u_1(x) .[/math]

Esto describe la velocidad inicial de cada punto de la cuerda en \(t = 0\).

2.4 Descripción de cada término

\(u(x,t)\): Desplazamiento de la cuerda en la posición \(x\) y tiempo \(t\).
\(\frac{\partial^2 u}{\partial t^2}(x,t)\): Aceleración de la cuerda en la posición \(x\) y tiempo \(t\).
\(\frac{\partial^2 u}{\partial x^2}(x,t)\): Curvatura de la cuerda en la posición \(x\) y tiempo \(t\).
\(u_0(x)\): Impulso inicial de la cuerda en la posición \(x\).
\(u_1(x)\): Velocidad inicial de la cuerda en la posición \(x\).

De esta forma se puede escribir el siguiente sistema que recoge toda la información mencionada anteriormente:

[math] \left\{ \begin{aligned} &u_{tt}(x,t) = u_{xx}(x,t) & 0 \lt x \lt 1, t \gt 0, \\ &u(0, t) = u(1, t)=0, & t \gt 0, \\ &u(x, 0) =u_0(x) \\ &u_t(x, 0) = u_1(x), \end{aligned} \right. [/math]

3 Modelización de los desplazamientos transversales de la cuerda

Para resolver la ecuación de ondas utilizando el método de separación de variables y expresar la solución en términos de los coeficientes de Fourier de los datos iniciales, se siguen los siguientes pasos:

1. Plantear la solución mediante separación de variables:

Asumimos una solución de la forma \( u(x,t) = X(x)T(t) \). Sustituyendo en la ecuación de ondas \( \frac{\partial^2 u}{\partial t^2} = \frac{\partial^2 u}{\partial x^2} \) y dividiendo por \( X(x)T(t) \), obtenemos:

[math] \frac{T''(t)}{T(t)} = \frac{X''(x)}{X(x)} = -\lambda [/math]

Esto nos lleva a dos ecuaciones ordinarias:

[math] X''(x) + \lambda X(x) = 0 [/math]

[math] T''(t) + \lambda T(t) = 0 [/math]

2. Resolver las ecuaciones ordinarias:

La solución de \( X(x) \) depende del valor de \( \lambda \). Considerando las condiciones de frontera \( X(0) = 0 \) y \( X(1) = 0 \), obtenemos una familia de soluciones:

[math] X_n(x) = \sin(\sqrt{\lambda_n} x) = \sin(n\pi x) [/math]

con \( \lambda_n = (n\pi)^2 \) y \( n = 1, 2, 3, \ldots \).

Para \( T(t) \), la familia de soluciones obtenida es:

[math] T_n(t) = A_n \cos(n\pi t) + B_n \sin(n\pi t) [/math]

3. Construir la solución general:

La solución general es la suma de todas las soluciones particulares:

[math] u(x,t) = \sum_{n=1}^{\infty} \left[ A_n \cos(n\pi t) + B_n \sin(n\pi t) \right] \sin(n\pi x) [/math]

4. Determinar los coeficientes \( A_n \) y \( B_n \) usando las condiciones iniciales:

Dado \( u(x,0) = u_0(x) \) y \( \frac{\partial u}{\partial t}(x,0) = u_1(x) \), tenemos:

[math] u_0(x) = \sum_{n=1}^{\infty} A_n \sin(n\pi x) [/math]

[math] u_1(x) = \sum_{n=1}^{\infty} n\pi B_n \sin(n\pi x) [/math]

Los coeficientes de Fourier \( A_n \) y \( B_n \) se determinan como:

[math] A_n = 2 \int_0^1 u_0(x) \sin(n\pi x) \, dx [/math]

[math] B_n = \frac{2}{n\pi} \int_0^1 u_1(x) \sin(n\pi x) \, dx [/math]

5. Expresar la solución final:

Sustituimos los coeficientes en la solución general y obtenemos la solución final:

[math] u(x,t) = \sum_{n=1}^{\infty} \left[ \left( 2 \int_0^1 u_0(x) \sin(n\pi x) \, dx \right) \cos(n\pi t) + \left( \frac{2}{n\pi} \int_0^1 u_1(x) \sin(n\pi x) \, dx \right) \sin(n\pi t) \right] \sin(n\pi x) [/math]

4 Particularización del sistema de ecuación de ondas a unos datos iniciales dados

En esta sección, se particulariza el problema de la cuerda vibrante eligiendo unas condiciones iniciales específicas. Para hallar la solución, simplemente se han sustituido las condiciones iniciales en la solución obtenida anteriormente. El sistema correspondiente a dichas condiciones iniciales es el siguiente:

[math] \left\{ \begin{aligned} &u_{tt}(x,t) = u_{xx}(x,t) & 0 \lt x \lt 1, t \gt 0, \\ &u(0, t) = u(1, t) = 0, & t \gt 0, \\ &u(x, 0) = e^{-100(x - 1/2)^2}, \\ &u_t(x, 0) = 0, \end{aligned} \right. [/math]

de forma que la solución tras sustituir los correspondientes datos iniciales se halla a continuación.

4.1 Determinar los coeficientes [math] A_{n} [/math] y [math] B_{n} [/math] usando las condiciones iniciales

Dado \( u_0(x) = e^{-100(x - 1/2)^2} \) y \( u_1(x) = 0 \):

[math] A_n = 2 \int_0^1 e^{-100(x - 1/2)^2} \sin(n\pi x) \, dx [/math]

[math] B_n = 0[/math]

4.2 Expresar la solución final

[math] u(x,t) = \sum_{n=1}^{\infty} \left( 2 \int_0^1 e^{-100(x - 1/2)^2} \sin(n\pi x) \, dx \right) \cos(n\pi t) \sin(n\pi x) [/math]

Para graficar la solución en el intervalo de tiempo \( t \in [0,2] \), tomaremos los primeros 50 términos de la serie.

En la figura superior se muestra la aproximación de la función u, dibujada para [math](x,t)\in[0,1]\times[0,1][/math]. En las dos figuras inferiores se muestra la diferencia entre la aproximación de la solución en [math]x=1/2[/math] y la solución estacionaria en [math]x=1/2[/math] a lo largo del tiempo para [math]t\in[0,1][/math]. En el caso de la figura inferior izquierda se hace esto para [math]\kappa=1[/math], mientras que en el caso de la figura inferior derecha se hace para [math]\kappa=1/2[/math]

4.3 Código

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

% Datos
numterminos_serie = 50;                    % Número de términos para la serie de Fourier
numvalores_t = 400;                        % Número de valores de t
numvalores_x = 400;                        % Número de valores de x
numvalores_trapecio = 300;                 % Número de valores al discretizar para la fórmula del trapecio

% Generación de puntos de evaluación
t = linspace(0, 2, numvalores_t);               % Vector de valores de t
x = linspace(0, 1, numvalores_x);               % Vector de valores de x
s = linspace(0, 1, numvalores_trapecio);        % Vector de valores de s

% Cálculo de la solución con los primeros 50 términos de la serie de Fourier
U = zeros(length(x), length(t));
for n = 1:numterminos_serie
    integrando = @(s) 2*exp(-100*(s-1/2).^2).*sin(n*pi*s);
    U = U + trapz(s,integrando(s))*cos(n*pi*t)'*sin(n*pi*x);
end

%Gráfica
surf(x,t,U, 'EdgeColor', 'none'); % Gráfica de la solución en función de x y t
xlabel('x'); ylabel('t'); zlabel('u(x,t)'); % Etiquetas de ejes
legend('u(x,t)'); % Leyenda

5 Condiciones iniciales de una onda que viaja en un único sentido

Tomar como datos iniciales los correspondientes a una onda que viaja en un solo sentido implica que la forma de la onda en un punto dado \( x \) y tiempo \( t \) está determinada por la función \( f(x - t) \). En otras palabras, la onda se desplaza hacia la derecha a una velocidad constante \( c \), donde \( c \) es la velocidad de propagación de la onda. Esto se logra al modificar apropiadamente las condiciones iniciales \( u_0(x) \) y \( u_1(x) \) para que reflejen esta propiedad de la onda unidireccional. En este caso particular, \( u_0(x) \) representa la forma inicial de la onda y \( u_1(x) \) representa la derivada de \( u_0(x) \) con respecto a \( x \), ajustada negativamente, para reflejar la velocidad y dirección de propagación.

[math] \left\{ \begin{aligned} &u_{tt}(x,t) = u_{xx}(x,t) & 0 \lt x \lt 1, t \gt 0, \\ &u(0, t) = u(1, t) = 0, & t \gt 0, \\ &u(x, 0) = e^{-100(x - 1/2)^2}, \\ &u_t(x, 0) = 200(x - 1/2)e^{-100(x - 1/2)^2}, \end{aligned} \right. [/math]

notése que se toma \( u_0(x) = f(x) \) y \( u_1(x) = -f'(x) \), donde \( f(x) = e^{-100(x-1/2)^2} \).

Calcular los coeficientes de Fourier:

Los coeficientes de Fourier se calculan utilizando las siguientes fórmulas:

[math] A_n = 2 \int_0^1 e^{-100(x-1/2)^2} \sin(n\pi x) \, dx \\ B_n = \frac{2}{n\pi} \int_0^1 200(x - 1/2)e^{-100(x - 1/2)^2} \sin(n\pi x) \, dx [/math]

Expresar la solución final:

Una vez que tenemos los coeficientes de Fourier, la solución se expresa como una serie infinita:

[math] u(x,t) = \sum_{n=1}^{\infty} \left[ A_n \cos(n\pi t) + B_n \sin(n\pi t) \right] \sin(n\pi x) [/math]

Esto representa una suma infinita sobre todos los términos de la serie de Fourier. A continuación se hará la representación gráfica de la serie:

En la figura superior se muestra la aproximación de la función u, dibujada para [math](x,t)\in[0,1]\times[0,1][/math]. En las dos figuras inferiores se muestra la diferencia entre la aproximación de la solución en [math]x=1/2[/math] y la solución estacionaria en [math]x=1/2[/math] a lo largo del tiempo para [math]t\in[0,1][/math]. En el caso de la figura inferior izquierda se hace esto para [math]\kappa=1[/math], mientras que en el caso de la figura inferior derecha se hace para [math]\kappa=1/2[/math]

Ahora vemos como el perfil viaja a velocidad 1:

En la figura superior se muestra la aproximación de la función u, dibujada para [math](x,t)\in[0,1]\times[0,1][/math]. En las dos figuras inferiores se muestra la diferencia entre la aproximación de la solución en [math]x=1/2[/math] y la solución estacionaria en [math]x=1/2[/math] a lo largo del tiempo para [math]t\in[0,1][/math]. En el caso de la figura inferior izquierda se hace esto para [math]\kappa=1[/math], mientras que en el caso de la figura inferior derecha se hace para [math]\kappa=1/2[/math]

5.1 Código

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

% Datos
numterminos_serie = 50;                    % Número de términos para la serie de Fourier
numvalores_t = 400;                        % Número de valores de t
numvalores_x = 400;                        % Número de valores de x
numvalores_trapecio = 300;                 % Número de valores al discretizar para la fórmula del trapecio

% Generación de puntos de evaluación
t = linspace(0, 4, numvalores_t);               % Vector de valores de t
x = linspace(0, 1, numvalores_x);               % Vector de valores de x
s = linspace(0, 1, numvalores_trapecio);        % Vector de valores de s

% Cálculo de la solución con los primeros 50 términos de la serie de Fourier
U = zeros(length(x), length(t));
for n = 1:numterminos_serie
    integrando1 = @(s) 2*exp(-100*(s-1/2).^2).*sin(n*pi*s);
    integrando2 = @(s) 2/(n*pi)*200*(s-1/2).*exp(-100*(s-1/2).^2).*sin(n*pi*s);
    U = U + (trapz(s,integrando1(s))*cos(n*pi*t)' + trapz(s,integrando2(s))*sin(n*pi*t)')*sin(n*pi*x);
end

%Gráfica
surf(x,t,U, 'EdgeColor', 'none'); % Gráfica de la solución en función de x y t
xlabel('x'); ylabel('t'); zlabel('u(x,t)'); % Etiquetas de ejes
legend('u(x,t)'); % Leyenda

6 Condiciones de frontera de Neumann

El nuevo sistema de ecuaciones para la ecuación de ondas con condiciones de frontera de Neumann es:

[math] \left\{ \begin{align*} &u_{tt}(x,t) = u_{xx}(x,t) & \text{para } 0 \lt x \lt 1, t \gt 0, \\ &u_x(0, t) = u_x(1, t) = 0, & \text{para } t \gt 0, \\ &u(x, 0) = e^{-100(x - 1/2)^2}, \\ &u_t(x, 0) = 200(x - 1/2)e^{-100(x - 1/2)^2}, \end{align*} \right. [/math]

Tal y como se hizo en los apartados anteriores se resuelve el sistema por el método de separación de variables (no se desarrollará de nuevo):

[math] u(x,t)= \sum _{n=1}^{\infty} [A_n cos (n \pi t) + B_n sin(n \pi t)] cos(n \pi x)[/math]

donde

[math] A_n = 2 \int_0^1 e^{-100(x-1/2)^2} \cos(n\pi x) \, dx \\ B_n = 2 \int_0^1 200(x - 1/2)e^{-100(x - 1/2)^2} \cos(n\pi x) \, dx [/math]

Resolviendo ambos coeficientes se puede hallar la solución para las condiciones de Newmann y unicamente faltaría por representar la solución final de forma gráfica.

En la figura superior se muestra la aproximación de la función u, dibujada para [math](x,t)\in[0,1]\times[0,1][/math]. En las dos figuras inferiores se muestra la diferencia entre la aproximación de la solución en [math]x=1/2[/math] y la solución estacionaria en [math]x=1/2[/math] a lo largo del tiempo para [math]t\in[0,1][/math]. En el caso de la figura inferior izquierda se hace esto para [math]\kappa=1[/math], mientras que en el caso de la figura inferior derecha se hace para [math]\kappa=1/2[/math]

6.1 Código

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

% Datos
numterminos_serie = 50;                    % Número de términos para la serie de Fourier
numvalores_t = 400;                        % Número de valores de t
numvalores_x = 400;                        % Número de valores de x
numvalores_trapecio = 300;                 % Número de valores al discretizar para la fórmula del trapecio

% Generación de puntos de evaluación
t = linspace(0, 2, numvalores_t);               % Vector de valores de t
x = linspace(0, 1, numvalores_x);               % Vector de valores de x
s = linspace(0, 1, numvalores_trapecio);        % Vector de valores de s

% Cálculo de la solución con los primeros 50 términos de la serie de Fourier
U = zeros(length(x), length(t));
for n = 1:numterminos_serie
    integrando1 = @(s) 2*exp(-100*(s-1/2).^2).*cos(n*pi*s);
    integrando2 = @(s) 2/(n*pi)*200*(s-1/2).*exp(-100*(s-1/2).^2).*cos(n*pi*s);
    U = U + (trapz(s,integrando1(s))*cos(n*pi*t)' + trapz(s,integrando2(s))*sin(n*pi*t)')*cos(n*pi*x);
end

%Gráfica
surf(x,t,U, 'EdgeColor', 'none'); % Gráfica de la solución en función de x y t
xlabel('x'); ylabel('t'); zlabel('u(x,t)'); % Etiquetas de ejes
legend('u(x,t)'); % Leyenda

7 Representación solución fundamental de dimensión 3 y regularización para la solución fundamental de dimensión 2

La solución fundamental es la solución que se obtiene al dar un impulso inicial muy localizado en \( x = 0 \). Matemáticamente resuelve el sistema

[math]u_{tt} - c^2 \Delta u = 0, \quad x \in \mathbb{R}^n, \, t \gt 0, [/math]

con condiciones iniciales [math]u(x,0) = 0, \quad u_t(x,0) = \delta(x), \quad x \in \mathbb{R}^n,[/math] donde [math]\delta(x)[/math] es la delta de Dirac. Formalmente, se define como el siguiente límite [math]\delta(x) = \lim_{r \to 0} \frac{1}{|B(0,r)|} \chi_{B(0,r)}(x) \sim\begin{cases}\infty & \text{si } x = 0, \\0 & \text{si } x \neq 0,\end{cases} [/math] donde [math]\chi_{B(0,r)}(x)[/math] es la función característica de la bola centrada en 0 de radio r y [math]|B(0,r)|[/math] es el volumen de la bola.

7.1 Solución Fundamental dimensión 1

La expresión de la solución fundamental de la ecuación de ondas de dimensión 1 es:

[math]K_1(x,t) = \frac{1}{2c} [H(x + ct) - H(x - ct)], [/math]

donde H(s) es la función de Heaviside,

[math] H(s) = \begin{cases} 0, & \text{si } s \lt 0 \\ 1, & \text{si } s \geq 0 \end{cases} [/math]

En la figura superior se muestra la aproximación de la función u, dibujada para [math](x,t)\in[0,1]\times[0,1][/math]. En las dos figuras inferiores se muestra la diferencia entre la aproximación de la solución en [math]x=1/2[/math] y la solución estacionaria en [math]x=1/2[/math] a lo largo del tiempo para [math]t\in[0,1][/math]. En el caso de la figura inferior izquierda se hace esto para [math]\kappa=1[/math], mientras que en el caso de la figura inferior derecha se hace para [math]\kappa=1/2[/math]

7.2 Solución Fundamental dimensión 2

La expresión de la solución fundamental de la ecuación de ondas de dimensión 2 es:

[math] K_2(x,t) = \frac{1}{2\pi c} \frac{1}{\sqrt{c^2 t^2 - |x|^2}}\chi_{B(0,ct)(x)}, [/math]

donde \(\chi_{B(0,ct)}(x)\) es la función característica de la bola de centro \(0\) y radio \(ct\). Para evitar la singularidad se usará la regularización:

[math] K_{\epsilon, 2}(x, t) = \frac{c^2 t^2 - |x|^2}{(c^2 t^2 - |x|^2) + \epsilon} \chi_{B(0,ct)}(x), \quad \epsilon = 0.01 [/math]

En la figura superior se muestra la aproximación de la función u, dibujada para [math](x,t)\in[0,1]\times[0,1][/math]. En las dos figuras inferiores se muestra la diferencia entre la aproximación de la solución en [math]x=1/2[/math] y la solución estacionaria en [math]x=1/2[/math] a lo largo del tiempo para [math]t\in[0,1][/math]. En el caso de la figura inferior izquierda se hace esto para [math]\kappa=1[/math], mientras que en el caso de la figura inferior derecha se hace para [math]\kappa=1/2[/math]

7.3 Solución Fundamental dimensión 3

La expresión de la solución fundamental de la ecuación de ondas de dimensión 3 es:

[math] K_3(x,t) = \frac{\delta(|x| - ct)}{4\pi c |x|}, \quad t \gt 0 [/math]

se sustituirá la delta de Dirac por su aproximación. Veáse que dicha aproximación es:

[math]\delta(s) \sim \varphi_k(s) = \sqrt{\frac{k}{\pi}} e^{-ks^2}, k \gg 1 [/math]

con k = 1000 (observar que \(\int_{-\infty}^{\infty} \varphi_k(s) \, ds = 1\).

En la figura superior se muestra la aproximación de la función u, dibujada para [math](x,t)\in[0,1]\times[0,1][/math]. En las dos figuras inferiores se muestra la diferencia entre la aproximación de la solución en [math]x=1/2[/math] y la solución estacionaria en [math]x=1/2[/math] a lo largo del tiempo para [math]t\in[0,1][/math]. En el caso de la figura inferior izquierda se hace esto para [math]\kappa=1[/math], mientras que en el caso de la figura inferior derecha se hace para [math]\kappa=1/2[/math]

7.3.1 Código

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

% Datos
c=1;                                       % Velocidad de propagación escogida
numvalores_t = 400;                        % Número de valores de t
numvalores_x = 400;                        % Número de valores de r

% Generación de puntos de evaluación
t = linspace(0, 1, numvalores_t);          % Vector de valores de t
x = linspace(-1, 1, numvalores_x);         % Vector de valores de x

% Crear la cuadrícula de puntos
[X, T] = meshgrid(x, t);

% Definir la función de Heaviside
H = @(s) s >= 0;

% Definir la función K1(x, t)
K1 = @(x, t) 1/(2*c) * ( H(x + c*t) - H(x - c*t) );

%Gráfica
surf(x,t,K1(X,T),'EdgeColor', 'none'); % Gráfica de la solución fundamental para n=1 en función de x y t
shading interp
xlabel('x'); ylabel('t'); zlabel('K_1(x,t)'); % Etiquetas de ejes
legend('K_1(x,t)'); % Leyenda
title('Solución fundamental para n=1') % Título

7.3.2 Código

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

% Datos
c=1;                                       % Velocidad de propagación escogida
epsilon=0.01;                              % Epsilon escogido para evitar la singularidad
numvalores_t = 400;                        % Número de valores de t
numvalores_r = 400;                        % Número de valores de r

% Generación de puntos de evaluación
t = linspace(0, 1, numvalores_t);              % Vector de valores de t
r = linspace(0, 1, numvalores_r);              % Vector de valores de r

% Crear la cuadrícula de puntos
[R, T] = meshgrid(r, t);

% Definir la función K2_epsilon(r, t)
K2_epsilon = @(r, t) (r<c*t)./(epsilon + 2*pi*c*sqrt(c^2*t.^2-r.^2)); 

%Gráfica
surf(r,t,K2_epsilon(R,T),'EdgeColor', 'none'); % Gráfica de la solución fundamental para n=2 en función de r y t
shading interp
xlabel('r'); ylabel('t'); zlabel('K_2^{\epsilon}(r,t)'); % Etiquetas de ejes
legend('K_2^{\epsilon}(r,t)'); % Leyenda
title('Solución fundamental para n=2') % Título

7.3.3 Código

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

% Datos
c=1;                                       % Velocidad de propagación escogida
k=1000;                                    % Valor de k escogido para aproximar la delta de Dirac
numvalores_t = 400;                        % Número de valores de t
numvalores_r = 400;                        % Número de valores de r

% Generación de puntos de evaluación
t = linspace(0, 1, numvalores_t);              % Vector de valores de t
r = linspace(0, 1, numvalores_r);              % Vector de valores de r

% Crear la cuadrícula de puntos
[R, T] = meshgrid(r, t);


% Definir la aproximación de la delta de dirac
Phi_k = @(s) sqrt(k/pi)*exp(-k*s.^2);

% Definir la función K2_epsilon(r, t)
K3_k = @(r, t) Phi_k(r-c*t)./(4*pi*c*r); 

%Gráfica
surf(r,t,K3_k(R,T),'EdgeColor', 'none'); % Gráfica de la solución fundamental para n=3 en función de r y t
shading interp
xlabel('r'); ylabel('t'); zlabel('K_3(r,t)'); % Etiquetas de ejes
legend('K_3(r,t)'); % Leyenda
title('Solución fundamental para n=3') % Título

8 Solución fundamental de dimensión 2 mediante el producto de convolución

La solución del problema en dimensión 2 viene dada a partir del siguiente sistema:

[math] \left \{ \begin{array}{ll} u_{tt}-c^2\Delta u=0, \quad x \in \mathbb{R}^2, t\gt0, \\ u(x,0)=0, u_t(x,0)=h (x) = \chi _{B(0,\frac{1}{2})}(x), \quad x \in \mathbb{R}^2 \end{array} \right. [/math]

viene dada por la convolución

[math] u(x,t) = \int K_2(x - y, t) h(y) \, dy. [/math]

( se debe dibujar las soluciones para t = 0, 0.5, 1, 2.)

En la figura superior se muestra la aproximación de la función u, dibujada para [math](x,t)\in[0,1]\times[0,1][/math]. En las dos figuras inferiores se muestra la diferencia entre la aproximación de la solución en [math]x=1/2[/math] y la solución estacionaria en [math]x=1/2[/math] a lo largo del tiempo para [math]t\in[0,1][/math]. En el caso de la figura inferior izquierda se hace esto para [math]\kappa=1[/math], mientras que en el caso de la figura inferior derecha se hace para [math]\kappa=1/2[/math]

8.1 Código

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

% Datos
c=1;                                     % Velocidad de propagación escogida
epsilon=0.01;                            % Epsilon escogido para evitar la singularidad
numvalores_r = 50;                       % Número de valores del radio r
numvalores_theta = 50;                   % Número de valores del ángulo theta para representar las gráficas
numvalores_trapecio = 600;               % Número de valores al discretizar para la fórmula del trapecio

% Generación de puntos de evaluación
t = [0,0.5,1,2];                               % Vector de valores de t
r = linspace(0, 1, numvalores_r);              % Vector de valores de r
theta = linspace(0,2*pi, numvalores_theta);    % Vector de valores de theta

ro = linspace(0, 1/2, numvalores_trapecio);      % Vector de valores de ro para la fórmula del trapecio
beta = linspace(0,2*pi, numvalores_trapecio);    % Vector de valores de beta para la fórmula del trapecio

% Definir la función K2 (para los puntos de la forma x-y con x=(r,0) e y=(ro*cos(theta),ro*sin(theta)) )
K2 = @(r, t, ro, beta) (sqrt( (r.^2+ro.^2-2*r.*ro.*cos(beta)) )<c*t)./(epsilon + 2*pi*c*sqrt(c^2*t.^2-(r.^2+ro.^2-2*r.*ro.*cos(beta)))) ;

%Integral
U=zeros(length(r),length(t));                            
for i=1:length(r)
    for j=1:length(t)
        integrando1=zeros(1,length(beta));  % Vector que acumula los valores de la primera integral para distintos beta
        for l=1:length(beta)
            integrando1(l)=trapz(ro, K2(r(i),t(j),ro,beta(l)) ); % Primera integral
        end
        % Se almacenan los valores de la integral doble para cada valor de r y t en U
        U(i,j)=trapz(beta,integrando1);  % Segunda integral
    end
end

% Puntos para dibujar las gráficas
X = r' * cos(theta);                                % Coordenadas x tras deshacer polares
Y = r' * sin(theta);                                % Coordenadas y tras deshacer polares

% Gráficas
for j=1:length(t)
    figure(j)
    surf(X, Y, (U(:,j)*ones(1,length(theta)) ), "FaceColor", "interp",'EdgeColor', 'none')      % Gráfica
    xlabel('x'); ylabel('y'); zlabel('u(x,y)')                                                  % Etiquetas de ejes
    legend("Solución en t=" + num2str(t(j)));                                                   % Leyenda
    title("t=" + num2str(t(j)))                                                                 % Título
end

8.2 Código

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

% Datos
c=1;                                     % Velocidad de propagación escogida
epsilon=0.01;                            % Epsilon escogido para evitar la singularidad
numvalores_r = 50;                       % Número de valores del radio r
numvalores_theta = 50;                   % Número de valores del ángulo theta para representar las gráficas

% Generación de puntos de evaluación
t = [0,0.5,1,2];                                % Vector de valores de t
r = linspace(0, 1, numvalores_r);               % Vector de valores de r
theta = linspace(0,2*pi, numvalores_theta);     % Vector de valores de theta

% Definir la función K2 (para los puntos de la forma x-y con x=(r,0) e y=(ro*cos(theta),ro*sin(theta)) )
K2 = @(r, t, ro, theta) (sqrt(r.^2+ro.^2-2*r.*ro.*cos(theta))<c*t)./(epsilon + 2*pi*c*sqrt(c^2*t.^2-(r.^2+ro.^2-2*r.*ro.*cos(theta)))) ;

%Integral
U=zeros(length(r),length(t));                               
for i=1:length(r)
    for j=1:length(t)
        % Se almacenan los valores de la integral para cada valor de r y t en U
        U(i,j)=integral2(@(theta,ro) K2(r(i),t(j),ro,theta), 0,2*pi,  0,1/2);  
    end
end

% Puntos para dibujar las gráficas
X = r' * cos(theta);                                % Coordenadas x tras deshacer polares
Y = r' * sin(theta);                                % Coordenadas y tras deshacer polares

% Gráficas
for j=1:length(t)
    figure(j)
    surf(X, Y, (U(:,j)*ones(1,length(theta)) ), "FaceColor", "interp",'EdgeColor', 'none')      % Gráfica
    xlabel('x'); ylabel('y'); zlabel('u(x,y)')                                                  % Etiquetas de ejes
    legend("Solución en t=" + num2str(t(j)));                                                   % Leyenda
    title("t=" + num2str(t(j)))                                                                 % Título
end

@@ Línea 184: / Línea 184: @@
 </math>
 </center>
-de forma que la solución tras sustituir los correspondientes datos iniciales es la siguiente:
+de forma que la solución tras sustituir los correspondientes datos iniciales se halla a continuación.
 === Determinar los coeficientes <math> A_{n} </math> y <math> B_{n} </math> usando las condiciones iniciales ===

Diferencia entre revisiones de «Ecuación de Ondas»