Deprecated: The each() function is deprecated. This message will be suppressed on further calls in /home/mat/public_html/w/includes/diff/DairikiDiff.php on line 434

Revisión del 21:03 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]

Nótese que al ser el seno y el coseno funciones periódicas de periodo \( 2\pi \), la función solución será periódica en tiempo y espacio de periodo \( 2 \) independientemente de las condiciones iniciales. La periodicidad en espacio no es de gran interés pues la cuerda sólo se encuentra en el intervalo \( [0,1] \), sin embargo, la periodicidad en tiempo se observará en las siguientes secciones.

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

La solución final con estas condiciones iniciales es:

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

Una vez hallada, se hará una representación de dicha solución. Para graficar la solución se tomarán los primeros 50 términos de la serie y se aproximará la integral haciendo uso del método del trapecio con una discretización de 300 puntos. Además, cabe destacar que se graficará para el intervalo de tiempo \( t \in [0,2] \). De forma que tomando \( x\in[0,1] \) para observar la cuerda por completo se tiene la siguiente gráfica:

En la figura superior se muestra la gráfica de la solución aproximada del sistema dado por la ecuación de ondas y sus condiciones iniciales y frontera mostradas anteriormente. Esta solución se representa tomando 50 términos de la serie de Fourier y 300 puntos en la discretización del método del trapecio. Además, se representa para \( t \in [0,2] \) y \( x\in[0,1] \).

Con esta imagen, se puede intuir que la solución es periódica en tiempo de periodo 2 tal y como se mencionó en la anterior sección. Por lo tanto, se ha decidido dibujarla para un intervalo de tiempo mayor \( t \in [0,4] \) para observar que esto en efecto es cierto:

En la figura superior se muestra la gráfica de la solución aproximada del sistema dado por la ecuación de ondas y sus condiciones iniciales y frontera mostradas anteriormente. Esta solución se representa tomando 50 términos de la serie de Fourier y 300 puntos en la discretización del método del trapecio. Además, se representa para \( t \in [0,4] \) y \( x\in[0,1] \).

Tal y como se puede observar, esto es cierto, es periódica en tiempo de periodo 2. Además, se puede observar que la amplitud máxima es 1, la cual se alcanza en los valores enteros de tiempo. También es posible ver que el tiempo que tarda la onda en llegar de un extremo a otro es una unidad de tiempo. Como además estamos considerando que la cuerda mide una unidad de espacio, la velocidad de propagación es uno, tal y como se especifica al plantear el sistema. Finalmente, cabe destacar que tal y como se ha impuesto en el sistema, en los extremos \( x=0 \) y \( x=1 \) el valor de la solución es 0, ya que la cuerda está fija en ellos.

En el sentido más físico, la gráfica muestra como la onda se mueve del centro a los extremos y una vez que llega a ellos, como estos se encuentran fijos, la onda rebota de forma simétrica cambiando el signo de la tensión.

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

En esta sección consideraremos unas condiciones iniciales correspondientes a una onda que viaja en un solo sentido. Para ello, basta tomar como datos iniciales los correspondientes a una onda cuya forma en un punto dado \( x \) y tiempo \( t \) está determinada por una función \( f \) tal que \( u(x,t)=f(x-t) \). Por lo tanto, en este caso tomaremos \( u_0(x) = f(x) \) y \( u_1(x) = -f'(x) \), donde se ha escogido \( f(x) = e^{-100(x-1/2)^2} \). El sistema resultante 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) = 200(x - 1/2)e^{-100(x - 1/2)^2}, \end{aligned} \right. [/math]

Una vez planteado el sistema, nos disponemos a calcular la solución de la misma forma que en las anteriores secciones.

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[ 2 \int_0^1 e^{-100(x-1/2)^2} \sin(n\pi x) \, dx \cos(n\pi t) + \frac{2}{n\pi} \int_0^1 200(x - 1/2)e^{-100(x - 1/2)^2} \sin(n\pi x) \, dx \sin(n\pi t) \right] \sin(n\pi x) [/math]

Una vez hallada, se hará una representación de dicha solución. Para graficar la solución se tomarán los primeros 50 términos de la serie y se aproximará la integral haciendo uso del método del trapecio con una discretización de 300 puntos. Además, cabe destacar que se graficará para el intervalo de tiempo \( t \in [0,4] \). De forma que tomando \( x\in[0,1] \) para observar la cuerda por completo se tiene la siguiente gráfica:

En la figura superior se muestra la gráfica de la solución aproximada del sistema dado por la ecuación de ondas y sus condiciones iniciales y frontera mostradas anteriormente. Esta solución se representa tomando 50 términos de la serie de Fourier y 300 puntos en la discretización del método del trapecio. Además, se representa para \( t \in [0,4] \) y \( x\in[0,1] \).

De nuevo, podemos ver que se cumple todo lo que se cumplía en la sección anterior. Las condiciones frontera se cumplen, la amplitud máxima es de 1 y se alcanza en tiempos enteros y la velocidad de propagación es 1. Esto último se puede observar de forma más clara si se observa la gráfica de perfil:

En la figura superior se muestra la gráfica de la solución aproximada del sistema dado por la ecuación de ondas y sus condiciones iniciales y frontera mostradas anteriormente. Esta solución se representa tomando 50 términos de la serie de Fourier y 300 puntos en la discretización del método del trapecio. Además, se representa para \( t \in [0,4] \) y \( x\in[0,1] \). En este caso, se representa de perfil con el fin de que sea posible observar correctamente la velocidad de propagación de la onda.

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 280: / Línea 280: @@
 Una vez hallada, se hará una representación de dicha solución. Para graficar la solución se tomarán los primeros 50 términos de la serie y se aproximará la integral haciendo uso del método del trapecio con una discretización de 300 puntos. Además, cabe destacar que se graficará para el intervalo de tiempo \( t \in [0,4] \). De forma que tomando \( x\in[0,1] \) para observar la cuerda por completo se tiene la siguiente gráfica:
 [[Archivo: Apar 4.png|700px|thumb|center|En la figura superior se muestra la gráfica de la solución aproximada del sistema dado por la ecuación de ondas y sus condiciones iniciales y frontera mostradas anteriormente. Esta solución se representa tomando 50 términos de la serie de Fourier y 300 puntos en la discretización del método del trapecio. Además, se representa para \( t \in [0,4] \) y \( x\in[0,1] \).]]
+De nuevo, podemos ver que se cumple todo lo que se cumplía en la sección anterior. Las condiciones frontera se cumplen, la amplitud máxima es de 1 y se alcanza en tiempos enteros y la velocidad de propagación es 1. Esto último se puede observar de forma más clara si se observa la gráfica de perfil:
-Ahora vemos como el perfil viaja a velocidad 1:
 [[Archivo: Apar 4b.png|700px|thumb|center|En la figura superior se muestra la gráfica de la solución aproximada del sistema dado por la ecuación de ondas y sus condiciones iniciales y frontera mostradas anteriormente. Esta solución se representa tomando 50 términos de la serie de Fourier y 300 puntos en la discretización del método del trapecio. Además, se representa para \( t \in [0,4] \) y \( x\in[0,1] \). En este caso, se representa de perfil con el fin de que sea posible observar correctamente la velocidad de propagación de la onda.]]
 === Código ===

Diferencia entre revisiones de «Ecuación de Ondas»