Ecuación de ondas¶
Vamos a estudiar distintas formas de resolver la ecuación de ondas. Las técnicas que vamos a estudiar sirven para otras ecuaciones, y para otros enfoques, como podéis leer más despacio en estas referencias:
- Leo M. González, Fabricio Macià, Antonio Souto Iglesias: Estabilidad, aproximación numérica y mecánica de fluidos
- J.M. Sanz-Serna: Fourier Techniques in Numerical Methods for Evolutionary Problems
Comenzamos con una cuerda vibrante en una cierta posición inicial. Para simplificar, tomaremos la velocidad inicial igual a cero, y las condiciones de contorno corresponden a los dos extremos fijos. Entre los exámenes y prácticas de otros cursos podrás encontrar más variantes, y si necesitas alguna otra, puedes consultar a los profesores.
Cuando hayas terminado de recorrer el cuaderno, te recomendamos que pruebes a cambiar el perfil inicial.
L = 1;
init_fn = @(x)(exp(-((x-0.5).^2)./0.01));
xs = linspace(0,L,200);
us = init_fn(xs);
plot(xs, us);
xlabel('x')
ylabel('u')
title('initial position of the string (initial speed is zero)')
Serie de senos¶
Comenzamos por recordar la serie de senos de una función f.
Una función $f:[0,L]\rightarrow\mathbb{R}$ tal que $f(0)=f(L)=0$ se puede escribir como suma de una serie de senos $$ f(x)=\sum_{k=1}^\infty \hat{f}_k \sin\left(\frac{\pi k}{l}x\right) $$ donde $$ \hat{f}_k = \frac{2}{L}\int_0^L f(x)\sin\left(\frac{\pi k}{l}x\right)\: d x $$
f = init_fn;
%Primeros coeficientes de la serie de senos de f
for k =[0:6]
k
integral = quad(@(x)(2*f(x).*sin(pi*k*x)), 0,1 )
end
Discrete sine transform¶
Cuando estudiamos la transformada de Fourier, vimos que el algoritmo FFT era una forma muy eficiente de aproximar los coeficientes de la serie de Fourier.
La discrete sine transform es una versión discreta de la serie de senos, análoga a la DFT (Discrete Fourier Transform) para la transformada de Fourier.
En general no se usa en tratamiento de señales, pero es útil precisamente para resolver ecuaciones en derivadas parciales por técnicas espectrales.
%Aproximamos los primeros coeficientes de la serie de senos de f
% mediante la dst
%Compara con los valores obtenidos antes
N = 100;
h = L/(N+1);
xs = linspace(h, L-h, N);
fs = f(xs);
Xk = 2*dst(fs)*(L/(N+1));
Xk(1:6)
% Aproximamos los los primeros coeficientes de la serie de senos de f
% mediante la dst
N = 20;
h = L/(N+1);
xs = linspace(h, L-h, N);
fs = f(xs);
Xk = 2*dst(fs)*(L/(N+1));
%Dibujamos los Coeficientes de la serie de senos
plot(Xk(1:20),'go')
legend('Coeficientes')
title('Coeficientes de la serie de senos')
%Dibujamos la función f
n_eval = 200
x_eval = linspace(0, L, n_eval);
f_eval = f(x_eval);
hold on;
plot(x_eval,f_eval);
%reconstruccion de la función mediante su serie de senos truncada
senos = zeros(1,n_eval);
for k = [1:N]
onda = @(x)(sin((pi*(k)/L)*x));
onda_eval = onda(x_eval);
senos = senos + Xk(k).*onda_eval;
end
plot(x_eval, senos, 'g-')
hold off;
title('Reconstrucción con la serie de senos')
legend('f', sprintf('Serie de senos hasta grado %d',N))
Al sumar sólo algunos de los términos de la serie, tenemos una aproximación a la función $f$:
$$ f(x)\approx\sum_{k=1}^N \hat{f}_k \sin\left(\frac{\pi k}{l}x\right) $$Ejercicio¶
- Dibuja las sumas parciales de la serie de senos, añadiendo cada vez más términos de la serie, para ver cómo mejora la aproximación. Calcula el máximo de la diferencia entre f y su aproximación.
Ejercicio¶
- Aproxima con series de senos otras funciones similares pero que no sean diferenciables: ¿empeora mucho la calidad de la aproximación?
Método espectral¶
Queremos resolver la ecuación de la cuerda vibrante: $$ \partial_{tt}u = \partial_{xx}u $$ (después de absorber las constantes reescalando).
Es equivalente a escribir: $$ \frac{\partial^2u}{\partial t^2}(t,x) = \frac{\partial^2u}{\partial x^2}(t,x). $$
Para cada $t$ fijo, descomponemos la función u en su serie de senos, en la variable $x$: $$ u(t,x)=\sum_{k=1}^\infty \hat{u}_k(t) \sin\left(\frac{\pi k}{l}x\right) $$ Sustituyendo esta expresión en la ecuación obtenemos: $$ \partial_{tt}\left(\sum_{k=1}^\infty \hat{u}_k(t) \sin\left(\frac{\pi k}{l}x\right)\right) =\partial_{xx}\left(\sum_{k=1}^\infty \hat{u}_k(t) \sin\left(\frac{\pi k}{l}x\right)\right) $$ Es decir: $$ \sum_{k=1}^\infty \partial_{tt}\hat{u}_k(t) \sin\left(\frac{\pi k}{l}x\right) =\sum_{k=1}^\infty \hat{u}_k(t) \partial_{xx}\sin\left(\frac{\pi k}{l}x\right) =\sum_{k=1}^\infty \hat{u}_k(t) \left(-(\frac{\pi k}{l})^2\right)\sin\left(\frac{\pi k}{l}x\right) $$ y como los coeficientes de la serie de senos son únicos, tenemos que, para cada $k$: $$ \partial_{tt}\hat{u}_k(t) = -\left(\frac{\pi k}{l}\right)^2\hat{u}_k(t) $$
Estas ecuaciones son fáciles de resolver. La solución de $\partial_{tt}\hat{u}_k(t) = -\left(\frac{\pi k}{l}\right)^2\hat{u}_k(t) $ con $u(0)=u_0$ y $u'(0)=v_0$ es $$ \hat{u}_k(t) = u_0\cos\left(\frac{\pi k}{l}t\right) + v_0\sin\left(\frac{\pi k}{l}t\right) $$
Como dijimos, asumimos que la velocidad inicial es 0. Escribimos la posición inicial de la cuerda $f$ como serie de senos: $$ f(x)=\sum_{k=1}^\infty \hat{f}_k \sin\left(\frac{\pi k}{l}x\right) $$ y deducimos que la solución del problema: $$ \begin{array}{rcl} \partial_{tt}u &=& \partial_{xx}u\\ u(0,x)&=&f(x)\\ \partial_t u(0,x)&=&0\\ u(t,0)&=&0\\ u(t,l)&=&0 \end{array} $$ tiene la forma: $$ u(t,x)=\sum_{k=1}^\infty \hat{f}_k\cos\left(\frac{\pi k}{l}t\right) \sin\left(\frac{\pi k}{l}x\right) $$
Usando la fórmula para el producto de $\sin$ y $\cos$ obtenemos: $$ u(t,x)=\sum_{k=1}^\infty \hat{f}_k\left(\sin\left(\frac{\pi k}{l}t+\frac{\pi k}{l}x\right)+\sin\left(\frac{\pi k}{l}t-\frac{\pi k}{l}x\right)\right) $$ Es decir, la solución es suma de ondas viajeras.
Dibujamos la posición de la cuerda en varios instantes del intervalo $[0,1]$.
t0 = 0;
tf = 1;
N = 100;
h = L/(N+1);
xs = linspace(h, L-h, N);
fs = f(xs);
Xk = 2*dst(fs)*(L/(N+1));
%reconstruccion
nframes = 6;
hold on;
legends = cell(1, nframes);
ts = linspace(t0, 1, nframes);
for j = [1:nframes]
t = ts(j);
fourier = zeros(1, n_eval);
for k = [1: N]
onda = @(x)(sin((pi*k/L)*x));
onda_eval = onda(x_eval);
fourier = fourier + Xk(k)*cos(pi*k*t/L)*onda_eval;
end
plot(x_eval, fourier)
legends{j} = sprintf('t=%.3f',t);
end
hold off;
title('Ondas, solucion exacta (pero truncada)')
xlabel('x')
ylabel('u')
legend(legends)
Puedes generar una animación siguiendo este tutorial
Discretización espacial, solución temporal exacta¶
Comparamos ahora con otro enfoque: Reemplazamos la función inicial $u$ de una variable real por un vector con los valores de $u$ en los puntos de un mallado:
$$ \overrightarrow{u}=\left(u(x_0),\dots,u(x_j),\dots,u(x_n)\right),\quad x_j=x_0+j\cdot h $$y reemplazamos la ecuación en derivadas parciales: $$ \partial_{tt}u = \partial_{xx}u $$ por la siguiente ecuación diferencial ordinaria $$ \partial_{tt}\overrightarrow{u} = D\cdot \overrightarrow{u} $$ donde $D$ es un operador de diferencias finitas (que aprovecha que $u(x_0)=u(x_n)=0$): $$ D\cdot \overrightarrow{u}=\left(0,\frac{-u(x_{0})+2u(x_1)-u(x_{2})}{h^2},\dots,\frac{-u(x_{j-1})+2u(x_j)-u(x_{j+1})}{h^2},\dots,0\right) $$
La matriz $D$ es una vieja conocida: $$D=\left(\begin{array}{rrrrr} 2.0 & -1.0 & 0.0 & 0.0 & 0.0 \\ -1.0 & 2.0 & -1.0 & 0.0 & 0.0 \\ 0.0 & -1.0 & 2.0 & -1.0 & 0.0 \\ 0.0 & 0.0 & -1.0 & 2.0 & -1.0 \\ 0.0 & 0.0 & 0.0 & -1.0 & 2.0 \end{array}\right)$$
diff_fin_matrix = @(n)(diag(2*ones(1,n)) - diag(ones(1,n-1),1) - diag(ones(1,n-1),-1));
D = diff_fin_matrix(5)
Este sistema de EDOs se puede resolver de forma exacta. El primer paso es expresarlo como un sistema de primer orden: $$ \frac{\partial}{\partial t} \left( \begin{array}{c} \overrightarrow{u}\\ \overrightarrow{v} \end{array} \right)= \left( \begin{array}{c|c} \mathbf{0}&\mathbf{Id}\\ \hline D&\mathbf{0} \end{array} \right)\cdot \left( \begin{array}{c} \overrightarrow{u}\\ \overrightarrow{v} \end{array} \right)$$ Llamamos a esta matriz $$ C = \left( \begin{array}{c|c} \mathbf{0}&\mathbf{Id}\\ \hline D&\mathbf{0} \end{array} \right) $$ y la solución es: $$ \left( \begin{array}{c} \overrightarrow{u}(t)\\ \overrightarrow{v}(t) \end{array} \right)= e^{t\cdot C}\cdot \left( \begin{array}{c} \overrightarrow{u}(0)\\ \overrightarrow{v}(0) \end{array} \right) $$
L = 1;
Nx = 50;
h = L/Nx;
t0 = 0;
tf = 1;
xs = linspace(0, L , Nx+1);
xs_int = xs(2:Nx);
u0 = init_fn(xs_int);
v0 = zeros(1,Nx-1);
uv0 = [u0,v0];
C = zeros(2*(Nx-1),2*(Nx-1));
C(Nx:2*(Nx-1),1:Nx-1) = -(1/h^2)*diff_fin_matrix(Nx-1);
C(1:Nx-1,Nx:2*(Nx-1)) = eye(Nx-1);
nframes = 6;
hold on;
legends = cell(1, nframes);
ts = linspace(t0, 1, nframes);
for j = [1:nframes]
t = ts(j);
uv_t = expm(t*C)*uv0';
u_t = zeros(1,Nx+1);
u_t(2:Nx) = uv_t(1:Nx-1);
plot(xs, u_t)
legends{j} = sprintf('t=%.3f',t);
end
hold off;
legend(legends)
xlabel('x')
ylabel('u')
title('Integracion temporal exacta del problema espacial discretizado')
Ecuación de ondas con diferencias finitas¶
Usar la técnica de diferencias finitas siguiente (una aproximación razonable): $$ \begin{array}{lcl} \overrightarrow{u}_{j+1} &=& \overrightarrow{u}_j + h \overrightarrow{v}_j\\ \overrightarrow{v}_{j+1} &=& \overrightarrow{v}_j + h D\cdot \overrightarrow{u}_j \end{array} $$ es equivalente a usar un método de Euler para integrar el sistema de primer orden: $$ \frac{\partial}{\partial t} \left( \begin{array}{c} \overrightarrow{u}\\ \overrightarrow{v} \end{array} \right)= \left( \begin{array}{c|c} \mathbf{0}&\mathbf{Id}\\ \hline D&\mathbf{0} \end{array} \right)\cdot \left( \begin{array}{c} \overrightarrow{u}\\ \overrightarrow{v} \end{array} \right)$$
L = 1;
Nx = 50;
h = L/Nx;
t0 = 0;
tf = 1;
Nt = 700;
k = (tf - t0)/Nt;
xs = linspace(0, L , Nx+1);
xs_int = xs(2:Nx);
u0 = init_fn(xs_int);
v0 = zeros(1,Nx-1);
A = -(1/h^2)*diff_fin_matrix(Nx-1);
u = zeros(Nt+1,Nx+1); %solution to WE
u(1,2:Nx) = u0;
vt = v0;
for j =[1:Nt]
ut_int = u(j,2:Nx);
u_new = ut_int + k*vt;
vt = vt + k*(A*ut_int')';
u(j+1,2:Nx) = u_new;
end
nframes = 6;
hold on;
legends = cell(1, nframes);
nplot = 1;
ts = linspace(t0, tf, Nt);
for t = linspace(1, Nt, nframes)
plot(xs, u(floor(t),:))
legends{nplot} = sprintf('t=%.3f',ts(floor(t)));
nplot = nplot + 1;
end
hold off;
legend(legends)
xlabel('x')
ylabel('u')
title('Integracion temporal con Euler del problema espacial discretizado')
Ejercicio¶
Observa que para Nx=100 divisiones espaciales y Nt = 700 divisiones espaciales, aparecen oscilaciones al final del intervalo temporal. Comprueba que con un intervalo de tiempo un poco mayor o unas pocas divisiones espaciales menos, las oscilaciones crecen de forma exponencial.
Ejercicio (Método adaptativo de Runge-Kutta 4,5)¶
Resuelve el sistema de EDOs anterior usando el método adaptativo de Runge-Kutta (4,5), llamando al método solve_ivp de scipy.integrate.
Diferencia temporal de segundo orden¶
Es más habitual usar una diferencia de segundo orden para aproximar la derivada temporal: $$ \partial_{tt}u(t_{j},x)\approx \frac{-u(t_{j+1},x) + 2u(t_j,x) -u(t_{j-1},x)}{h^2} $$ Al aplicar este método a la ecuación de ondas $$ \partial_{tt}u = c\partial_{xx}u, $$ con un paso espacial $h$ y un paso temporal $k$, la condición de estabilidad es $$ \frac{ck}{h}<1. $$
L = 1;
Nx = 100;
h = L/Nx;
t0 = 0;
tf = 1;
Nt = 200;
k = (tf - t0)/Nt;
c = 1.0;
lamb2 = (c*k/h)^2;
xs = linspace(0, L , Nx+1);
xs_int = xs(2:Nx);
u0 = init_fn(xs_int);
u = zeros(Nt+1,Nx+1);
A = -(1/h^2)*diff_fin_matrix(Nx-1);
u = zeros(Nt+1,Nx+1); %solution to WE
u(1,2:Nx) = u0;
u(2,2:Nx) = (1-lamb2)*u(1,2:Nx) + (lamb2/2)*(u(1,1:Nx-1)+u(1,3:Nx+1));
for j =[2:Nt]
u(j+1,2:Nx) = 2*(1-lamb2)*u(j,2:Nx)-u(j-1,2:Nx)+lamb2*(u(j,1:Nx-1)+u(j,3:Nx+1));
end
nframes = 6;
hold on;
legends = cell(1, nframes);
nplot = 1;
ts = linspace(t0, tf, Nt);
for t = linspace(1, Nt, nframes)
plot(xs, u(floor(t),:))
legends{nplot} = sprintf('t=%.3f',ts(floor(t)));
nplot = nplot + 1;
end
hold off;
legend(legends)
xlabel('x')
ylabel('u')
title('Integracion temporal de segundo orden del problema espacial discretizado')
Ejercicio¶
Prueba a cambiar el intervalo temporal, el número de divisiones espaciales y temporales, de modo que la condición de estabilidad no se cumpla, pero por un margen estrecho (por ejemplo $1<\frac{ck}{h}<2$), y observa si la solución es estable.