Diferencia entre revisiones de «Ecuación de Ondas»

De MateWiki
Saltar a: navegación, buscar
(→‎Solución de la ecuación de ondas de dimensión 2 mediante el producto de convolución usando la solución fundamental)
(→‎Solución de la ecuación de ondas de dimensión 2 mediante el producto de convolución usando la solución fundamental)
Línea 559: Línea 559:
 
[[Archivo:Trapecio 2.png|700px|thumb|center|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>]]
 
[[Archivo:Trapecio 2.png|700px|thumb|center|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>]]
  
 +
 +
[[Archivo: Trapecio 5.png|700px|thumb|center|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>]]
  
 
[[Archivo: Trapecio 3.png|700px|thumb|center|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>]]
 
[[Archivo: Trapecio 3.png|700px|thumb|center|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>]]
−
 
−
[[Archivo: Integral 1.png|700px|thumb|center|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>]]
 
−
 
−
[[Archivo: Trapecio 5.png|700px|thumb|center|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>]]
 
  
 
[[Archivo: Integral 2.png|700px|thumb|center|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>]]
 
[[Archivo: Integral 2.png|700px|thumb|center|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>]]
Línea 571: Línea 569:
  
 
[[Archivo: Integral 4.png|700px|thumb|center|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>]]
 
[[Archivo: Integral 4.png|700px|thumb|center|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>]]
 +
[[Archivo: Integral 1.png|700px|thumb|center|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>]]
  
 
=== Código ===
 
=== Código ===

Revisión del 22:58 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

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 Ecuación de ondas en dimensión 1

En esta sección, se buscarán soluciones para la ecuación de ondas en dimensión 1 y se estudiarán dichas soluciones.

2.1 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.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.1.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.1.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.1.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]

2.2 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.

2.3 Particularización del sistema de la ecuación de ondas a unos datos iniciales específicos

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.

2.3.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]

2.3.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 desde el centro de la cuerda a un extremo y volver (rebotando simétricamente) 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.

2.3.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


2.4 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 . [/math]

[math] 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án las integrales 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, es periódica de periodo 2 y la velocidad de propagación es 1. Esto último de que la velocidad de propagación es 1 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.

Aquí se observa más claramente que el tiempo que tarda en llegar la onda de un extremo a otro es de una unidad de tiempo. Como la longitud de la cuerda es de una unidad de espacio, la velocidad de propagación es necesariamente 1.

Volviendo a la primera figura, es posible observar que la onda parece que se mueve en el sentido positivo de \( x \) a medida que transcurre el tiempo, de ahí que viaje en un sólo sentido. Además, en el momento en el que llega a la frontera en \( x=1 \), la onda rebota simétricamente cambiando el signo de la tensión. Estando en el sentido negativo de su posición vertical comienza a moverse en este caso hacia el sentido negativo de \( x \) a medida que transcurre el tiempo. De nuevo, cuando llega a la frontera en \( x=1 \) vuelve a rebotar simétricamente cambiando el signo de la tensión. Esto sigue haciéndolo una y otra vez. Además, en el momento exacto en el que llega a la frontera la solución parece anularse por un instante a lo largo de todo el intervalo de espacio \( [0,1] \).

2.4.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


2.5 Condiciones de frontera de Neumann

Finalmente, cambiaremos las condiciones de frontera para que sean de tipo Neumann. Con esta idea, se tiene el siguiente sistema:

[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), obteniendo la solución:

[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 . [/math]
[math] B_n = \frac{2}{n\pi} \int_0^1 200(x - 1/2)e^{-100(x - 1/2)^2} \cos(n\pi x) \, dx. [/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án las integrales 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] \).

De nuevo, se cumple alguna cosa que se cumplía anteriormente. La solución sigue siendo periódica de periodo 2 y la velocidad de propagación es 1 ya que este es el tiempo que tarda en moverse la onda de un extremo a otro de la cuerda. Sin embargo, la amplitud máxima ha dejado de ser 1 y lo más importante, los extremos de la cuerda ya no se encuentran fijos. Podemos ver que ahora los valores de \( u(0,t) \) y \( u(1,t) \) no son cero para todo tiempo \( t \). En este caso, estos valores varían, aunque se puede observar que las derivadas \( u_x(0,t) \) y \(u_x(1,t) \) sí que se mantienen en cero en todo instante \( t \) tal y como indican las condiciones Neumann. Esto se puede ver ya que la pendiente de la recta tangente a \( u(0,t_0) \) contenida en el plano \( t=t_0 \) siempre es 0 (esto se puede ver a ojo) independientemente del valor \( t_0 \) y de igual forma con \( u(1,t_0) \).

En el sentido físico, podemos observar cómo la onda se dirige a un extremo de la cuerda y en el momento que llega hace que aumente el valor de la solución en dicho extremo considerablemente. Una vez que llega a su tope, decrece de nuevo enviando la onda de vuelta en la dirección del otro extremo y continúa así de forma periódica.

2.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, 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


3 Ecuación de ondas en más de una dimensión

En esta sección se estudiará la ecuación de ondas en más de una dimensión, representando las soluciones fundamentales para distintas dimensiones y observando un ejemplo en dimensión 2 en el cual se obtiene la solución mediante la solución fundamental y la convolución.

3.1 Representación de la solución fundamental en distintas dimensiones

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]\left \{ \begin{array}{ll} u_{tt}-c^2\Delta u=0, \quad x \in \mathbb{R}^n, t\gt0, \\ u(x,0)=0, u_t(x,0)=\delta (x), \quad x \in \mathbb{R}^n, \end{array} \right. [/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.

3.1.1 Solución Fundamental en 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]

Si la representamos en una gráfica en función de x y de t, tomando \( x \in [-1,1] \) y \( t \in [0,1] \) se obtiene lo siguiente:

En la figura superior se muestra la gráfica de la solución fundamental de la ecuación de ondas de dimensión 1 en función de \( x \) y \( t \) para \( x \in [-1,1] \) y \( t \in [0,1] \).

3.1.2 Solución Fundamental en 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\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 a la hora de representarla se representará la regularización:

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

Sin embargo, para representarla en función de \( x_1, x_2 \) y \( t \) donde \( x_1, x_2 \) son las dos coordenadas de \( x \) se necesitarían 4 dimensiones. Por lo tanto, a la hora de representarla es importante notar que esta función (tanto la solución sin regularizar como la regularizada) es radial. Esto permite que sea representada tan sólo en función del radio \( r = |x| \) y el tiempo \( t \). Si hacemos esto para \( r\in [0,1] \) y \( t \in [0,1] \) tenemos lo siguiente:

En la figura superior se muestra la gráfica de la solución fundamental de la ecuación de ondas de dimensión 2 en función de \( r \) y \( t \) para \( r \in [0,1] \) y \( t \in [0,1] \).

3.1.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]

Para representarla, se sustituirá la delta de Dirac por su aproximación:

[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\) ). De nuevo, esta función es radial, lo cual permite que sea representada en función del radio \( r= |x| \) y el tiempo \( t \). Si hacemos esto para \( r\in[0,1] \) y \( t \in [0,1] \) se tiene lo siguiente:


En la figura superior se muestra la gráfica de la solución fundamental de la ecuación de ondas de dimensión 3 en función de \( r \) y \( t \) para \( r \in [0,1] \) y \( t \in [0,1] \).

3.1.4 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

3.1.5 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

3.1.6 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


3.2 Solución de la ecuación de ondas de dimensión 2 mediante el producto de convolución usando la solución fundamental

La solución del problema en dimensión 2:

[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_{\mathbb R^2} K_2(x - y, t) h(y) \, dy. [/math]

Dado que esta solución es radial (la solución fundamental es radial y h también y la solución es la convolución de dos funciones radiales), se facilita su representación. Sin embargo, en este caso, en vez de representarla en función del radio \( r=|x| \) y el tiempo \( t \) hemos decidido representarla en función de las coordenadas de \(x \) fijando un tiempo \(t \). Para ello, dado que sabemos que es radial, se ha decidido calcular la integral en coordenadas polares y tan sólo obtener los valores de la solución para los puntos de la forma \( x=(r,0) \) variando \( r \in [0,1]\). De esta forma los valores de la solución en el resto de puntos \(x \) en la bola unidad se obtienen girando la solución en los puntos de la forma \( x=(r,0) \).

La integral en coordenadas polares teniendo en cuenta que \( y=(\rho \cos(\beta), \rho \sin(\beta)) \) tiene la siguiente forma (usamos la regularización de la solución fundamental para la representación):

[math] u(r,t)=\int_{0}^{2 \pi} \int_{0}^{1/2} \frac{\chi_{(0,ct)} (r^2+\rho^2-2 r \rho cos(\beta))}{\epsilon + 2 \pi c \sqrt{c^2 t^2 – (r^2+\rho^2-2 r \rho cos(\beta))}} d\rho d\beta [/math]

Se han calculado varios valores de la solución variando \( r\in [0,1] \) y fijando los tiempos \(t=0,0.5,1,2 \). Con dichos valores, sabiendo que la solución es radial se han obtenido el resto de valores para varios puntos de la circunferencia unidad. Finalmente, se ha representado la solución en 4 gráficas variando el tiempo a fijar. Estas vienen representadas a continuación:

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


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

3.2.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

3.2.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
Obtenido de «https://mat.caminos.upm.es/w/index.php?title=Ecuación_de_Ondas&oldid=73199»