Instrucciones
- Guarde todos los archivos en la carpeta EXAMENES.
- No olvide incluir su nombre y su número de puesto en los archivos de la entrega
- No olvide entregar todo el código necesario para replicar sus resultados:
- matlab: debe entregar tanto las funciones como los scripts.
- python: es suficiente con entregar uno o varios cuadernos jupyter.
- No olvide responder a las preguntas planteadas. Formas aceptadas:
- cuaderno jupyter
- documento office
- Hay 7 preguntas, sólo es necesario responder a 5 de ellas. Cada pregunta se evalúa sobre 2 puntos.
- Pregunte todas sus dudas.
Una matriz¶
Consideramos $A_n$, que es una matriz $n\times n$:
- Tiene 1 en la diagonal principal.
- Tiene -1 en la diagonal debajo de la diagonal principal, y también en la (primera fila, última columna).
- Tiene 0 en el resto de la matriz.
Por ejemplo
- si $n=2$:
- si $n=4$:
Ejercicio 1¶
- Escribe una función que reciba como argumento un número
ny devuelva un arrayAde dimensiones $n\times n$ con la estructura descrita antes. - Describa qué factorizaciones admite esta matriz, y calcúlelas para n=4.
function A=An(n)
A = eye(n) - diag(ones(1,n-1),-1);
A(1,n) = -1;
end
A = An(4)
La matriz es cuadrada, luego podría tener factorizaciones LU o PLU, pero no es simétrica, luego no puede tener factorización de Cholesky ni LDLT.
La matriz no es invertible, luego no está garantizado que tenga factorizaciones LU ni PLU, pero el comando lu la encuentra sin problema:
[P,L,U] = lu(A)
P*L*U
Ejercicio 2¶
- Calcule el determinante y los autovalores de $A_n$ para valores de $n$ desde 4 hasta 20.
- Dibuje en una gráfica el máximo de la parte real de los autovalores de $A_n$ para valores de $n$ desde 4 hasta 20. Repita para el mínimo.
- Decida qué factorizaciones de $A_n$ se pueden utilizar para resolver el sistema : $$ A_n\cdot \mathbf{x} = \mathbf{b} $$ donde $\mathbf{b}$ es un vector con 1 en todas sus entradas.
min_n = 4;
max_n = 20;
% Almacenamos los determinantes, autovalor min y max en arrays distintos
% de tamaño max_n - min_n
dets = zeros(1,max_n - min_n);
max_aval = zeros(1,max_n - min_n);
min_aval = zeros(1,max_n - min_n);
for n = [min_n: max_n]
% Los valores de det, max_aval, min_aval correspondientes a la matriz de tamaño
% n los almcenamos en la posición n - min_ n + 1
% por ejemplo, n=4 va en la primera posición (i=1),
% n=20 va en la última posición (i=20 - min_n + 1),
i = n - min_n + 1;
A = An(n);
dets(i) = det(A);
avals_real = real(eig(A));
max_aval(i) = max(avals_real);
min_aval(i) = min(avals_real);
end
sizes = [min_n: max_n];
hold on;
plot(sizes, dets);
plot(sizes, max_aval);
plot(sizes, min_aval);
hold off;
xlabel('matrix size');
legend('dets', 'max aval', 'min aval')
title('eigenvalues and det of matrix of different size');
- La matriz es singular (la matriz tiene determinante 0), luego ninguna factorización puede ser de ayuda para resolver el sistema de ecuaciones. El sistema no tiene solución, pero aunque la hubiera tenido, las factorizaciones no pueden ayudar a resolver el sistema. En general, en este curso sólo hemos aprendido a resolver sistemas de ecuaciones lineales con solución única. Algunos métodos devuelven una "solución", que es una solución de mínimos cuadrados, y he dado esas soluciones por buenas si usaban alguna factorización.
A = An(4);
b = ones(4,1)
A\b
Un sistema de ODEs¶
Consideramos el siguiente sistema de ecuaciones diferenciales para n funciones desconocidas almacenadas en un vector $\overrightarrow{u}(t)=(u_1(t),\dots,u_{n}(t))$:
$$ \partial_t\overrightarrow{u}(t)= - n\cdot A_n\cdot\overrightarrow{u}(t) $$Por ejemplo, para n=4:
$$\partial_t\left(\begin{array}{c} u_1\\ u_2\\ u_3\\ u_4 \end{array}\right) = -4 \left(\begin{array}{cccc} 1 & 0 & 0 & - 1\\ - 1 & 1 & 0 & 0\\ 0 & - 1 & 1 & 0\\ 0 & 0 & - 1 & 1 \end{array}\right)\cdot \left(\begin{array}{c} u_1\\ u_2\\ u_3\\ u_4 \end{array}\right) $$Ejercicio 3¶
- Consideramos el sistema de EDOs anterior para n=2: $$ \partial_t\overrightarrow{u}(t)= - 2\cdot A_2\cdot\overrightarrow{u}(t) $$ para $$\overrightarrow{u}(t)=\left(\begin{array}{c} u_1 (t)\\ u_2 (t) \end{array}\right)$$ o, equivalentemente $$\begin{array}{ccc} \partial_t u_1 & = & - 2 u_1 + 2 u_2\\ \partial_t u_2 & = & 2 u_1 - 2 u_2 \end{array} $$
con las condiciones adicionales: $$ u_1(0)=0, u_2(0)=1 $$ en el intervalo temporal $[0,3]$.
- Identifique el tipo de problema planteado y proponga un método numérico para resolverlo de forma aproximada.
- Aproxime la solución del problema anterior y muestre en una gráfica la evolución de $u_1$ y $u_2$.
- Es un problema de valor inicial, y se puede resolver con el método de Euler o cualquier método de Runge-Kutta, por ejemplo. Usamos
solve_ivpcon el método adaptativo Runge-Kutta-Fehlberg de orden (4,5).
fun = @(t,u)([-2*u(1) + 2*u(2), 2*u(1) - 2*u(2)]);
y0 = [0,1];
N = 120;
t0 = 0;
tf = 1;
tspan = linspace(t0, tf, N);
[ts, ys] = ode45(fun, tspan, y0);
u1s = ys(:,1);
u2s = ys(:,2);
hold on;
plot(ts, u1s)
plot(ts, u2s)
hold off;
xlabel('t')
ylabel('u')
title('Evolucion de u1 y u2 a lo largo del tiempo')
Ejercicio 4¶
Consideramos la siguiente función periódica de periodo L=1: $$ f(x) = e^{-\left(\frac{x-0.5}{10}\right)^2}\;\text{ si } x \in[0,1] $$ Es decir, para x que no está en el intervalo $[0,1]$, se escribe x como suma de un entero n y un número real r que sí está en el intervalo $[0,1]$, y entonces, por ser periódica,
$$f(x)=f(n+r)=f(r) = e^{-\left(\frac{r-0.5}{10}\right)^2}$$
- Calcule de forma aproximada, por el método que prefiera (pero indique qué método ha usado), los coeficientes de la serie de Fourier de $f$.
- Aségurese de que puede reconstruir la función $f$ de forma aproximada a partir de los coeficientes que acaba de calcular.
- Mi intención era poneros la función que aparece en la gráfica, que es $f(x) = e^{-\left(10(x-0.5)\right)^2}$ en vez de $f(x) = e^{-\left(\frac{x-0.5}{10}\right)^2}$. El ejercicio no cambia por usar una u otra función. Podéis ver una animación de la solución de la ecuación para el otro perfil. Más adelante enlazaré a animaciones de la solución usando distintos métodos numéricos.
- Usamos código del cuaderno "fourier". Lo únioco que hay que cambiar es ajustar el periodo al intervalo [0,1] (observad la llamada
complex_quadrature(integrando,0, L), dondeL=1). Si definimos la función de forma periódica (mirad la funciónf_fullmás abajo), entonces se puede usar el intervalo [-1/2,1/2] y la funcióninit_fn.
L = 1;
init_fn = @(x)(exp(-(((x-.5)./10).^2)));
f_full = @(x)(exp(-(((x- floor(x)-.5)./10).^2)));
xs = linspace(-1,3,200);
us = f_full(xs);
plot(xs,us)
function integral = complex_quadrature(func, a, b)
real_func = @(x)(real(func(x)));
imag_func = @(x)(imag(func(x)));
[real_integral, error_real] = quad(real_func, a, b);
[imag_integral, error_imag] = quad(imag_func, a, b);
integral = real_integral + 1j*imag_integral;
end
f = init_fn;
orden = 6;
ks = [-orden:orden];
Xk = zeros(2*orden+1);
integrando = f;
integral = complex_quadrature(integrando,0,L);
Xk(1) = (1/L)*integral;
for k = ks
integrando = @(t) (f(t).*exp(-1j*(2*pi*k/L)*t));
integral = complex_quadrature(integrando,0,L);
if k<0
Xk(2*orden+2+k) = (1/L)*integral;
else
Xk(1+k) = (1/L)*integral;
end
end
hold on;
Ntimes = 200;
t_eval = linspace(0, L, Ntimes);
f_eval = f(t_eval);
plot(t_eval,f_eval)
%reconstruccion
fourier = zeros(1, Ntimes);
for k=ks
onda = @(t)(exp(1j*(2*pi*k/L)*t));
onda_eval = onda(t_eval);
if k<0
fourier = fourier + Xk(2*orden+2+k)*onda_eval;
else
fourier = fourier + Xk(1+k)*onda_eval;
end
end
legends = cell(1, 2);
legends{1} = 'f';
legends{2} = sprintf('Serie de Fourier hasta grado %d',orden);
max(abs(imag(fourier)))%comprobamos que la parte imaginaria es cero
plot(t_eval, real(fourier), 'g-')
title('Reconstruccion con la serie de Fourier')
legend(legends)
hold off;
Ecuación del transporte¶
Queremos resolver la ecuación del transporte: $$ \partial_{t}u + c \partial_{x}u = 0. $$ Es equivalente a escribir: $$ \frac{\partial u}{\partial t}(t,x) + c \frac{\partial u}{\partial x}(t,x) = 0. $$
El dato inicial es la función $f$ definida antes: $$ u(0,x) = f(x) = e^{-\left(\frac{x-0.5}{10}\right)^2}\;\text{ si } x \in[0,1] $$ y en la variable espacial ponemos condiciones de frontera periódicas: $$ u(t,0) = u(t,1) $$
Método espectral¶
Para cada $t$ fijo, descomponemos la función $u$ en su serie de Fourier, en la variable $x$ (aprovechando que es periódica en la variable $x$ con periodo $L=1$): $$ u(t,x)= \sum_{k=-\infty}^\infty \hat{u}_k(t) e^{i2\pi \frac{k}{L} x} $$
Sustituyendo esta expresión en la ecuación obtenemos:
$$ \partial_{t}\left(\sum_{k=-\infty}^\infty \hat{u}_k(t) e^{i2\pi \frac{k}{L} x}\right) + c\partial_{x}\left(\sum_{k=-\infty}^\infty \hat{u}_k(t) e^{i2\pi \frac{k}{L} x}\right)=0 $$Es decir:
$$ \sum_{k=1}^\infty \left(\partial_{t}\hat{u}_k(t) e^{i2\pi \frac{k}{L} x} +c \hat{u}_k(t) \partial_{x}e^{i2\pi \frac{k}{L} x}\right)= \sum_{k=1}^\infty \left(\partial_{t}\hat{u}_k(t) +c \hat{u}_k(t) i2\pi \frac{k}{L}\right) e^{i2\pi \frac{k}{L} x}= 0 $$y como los coeficientes de la serie de Fourier son únicos, tenemos que, para cada $k$:
$$ \partial_{t}\hat{u}_k(t) + c \left(i2\pi \frac{k}{L}\right)\hat{u}_k(t)=0 \quad (1) $$Estas ecuaciones son fáciles de resolver. La solución de la ecuación (1) es
$$ \hat{u}_k(t) = \hat{u}_k(0) e^{-i2\pi c\frac{k}{L}t} $$y la solución completa es:
$$ u(t,x)= \sum_{k=-\infty}^\infty \hat{u}_k(0) e^{i2\pi \frac{k}{L} (x-ct)} $$donde $\hat{u}_k(0)$ es el coeficiente de Fourier de orden k de la función $f$:
$$ f(x) = u(0,x)= \sum_{k=-\infty}^\infty \hat{u}_k(0) e^{i2\pi \frac{k}{L} x} $$Ejercicio 5¶
- Para la solución $u$ de la ecuación del transporte, dibuje las funciones $u(t,x)$, para $x$ en el intervalo $[0,1]$ y $t$ fijo. Utilice los siguientes valores de t: 0,0.1,0.2,0.3,0.4,0.5.
- Tras el largo preámbulo, y al igual que hicimos para la ecuación de ondas, el único cambio necesario para obtener la solución para un tiempo concreto es reemplazar la línea
onda_eval = onda(x_eval)poronda_eval = onda(x_eval - c*t).
c=1;
% para un tiempo concreto
t = 0.2;
fourier = zeros(1, Ntimes);
for k=ks
onda = @(t)(exp(1j*(2*pi*k/L)*t));
onda_eval = onda(t_eval - c*t);
if k<0
fourier = fourier + Xk(2*orden+2+k)*onda_eval;
else
fourier = fourier + Xk(1+k)*onda_eval;
end
end
legends = cell(1, 1);
legends{1} = sprintf('Solución para t=%s',t);
max(abs(imag(fourier)))%comprobamos que la parte imaginaria es cero
plot(t_eval, real(fourier), 'g-')
title(sprintf('Solución para t=%s',t))
legend(legends)
hold off;
- Opcional, pero mejora mucho si dibujamos la función para varios instantes de tiempo en la misma gráfica:
% para un tiempo concreto
Nt = 6;
times = linspace(0, 0.5, Nt);
legends = cell(1, Nt);
hold on;
for j = [1:Nt]
t = times(j);
fourier = zeros(1, Ntimes);
for k=ks
onda = @(t)(exp(1j*(2*pi*k/L)*t));
onda_eval = onda(t_eval - c*t);
if k<0
fourier = fourier + Xk(2*orden+2+k)*onda_eval;
else
fourier = fourier + Xk(1+k)*onda_eval;
end
end
legends{j} = sprintf('t=%.2f',t);
plot(t_eval, real(fourier), '-')
end
hold off;
title('Solución para varios tiempos')
legend(legends)
hold off;
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-1})\right),\quad x_j=j\cdot h = \frac{j}{n} $$No incluímos $u(x_n)$ porque $u(x_0)=u(x_n)$.
Reemplazamos la ecuación en derivadas parciales:
$$ \partial_{t}u + c \partial_{x}u = 0 $$por la siguiente ecuación diferencial ordinaria
$$ \partial_{t}\overrightarrow{u} + c D\cdot \overrightarrow{u} = 0 $$donde $D$ es un operador de diferencias finitas.
Por ejemplo, si usamos diferencias finitas de dos puntos hacia atrás (backwards 2-point finite differences) para aproximar $\partial_x u$ (aprovechamos que las condiciones de frontera son periódicas: $u(x_0)=u(x_n)$).
$$ D\cdot \overrightarrow{u}=\left( \frac{u(x_{0})-u(x_{n-1})}{h}, \frac{u(x_{1})-u(x_{0})}{h}, \dots, \frac{u(x_{n-1})-u(x_{n-2})}{h}\right) $$y la matriz $D$ es
$$D=\frac{1}{h}\left(\begin{array}{rrrrr} 1 & 0 & 0 & 0 & -1 \\ -1 & 1 & 0 & 0 & 0 \\ 0 & -1 & 1 & 0 & 0 \\ 0 & 0 & -1 & 1 & 0 \\ 0 & 0 & 0 & -1 & 1 \end{array}\right)$$
- Este método numérico se llama upwind.
Si usamos diferencias centradas (centered 3-point finite difference) para aproximar $\partial_x u$: $$ D\cdot \overrightarrow{u}=\left( \frac{u(x_{1})-u(x_{n-1})}{2h}, \frac{u(x_{2})-u(x_{0})}{2h}, \dots, \frac{u(x_{j+1})-u(x_{j-1})}{2h}, \dots\right) $$ y la matriz $D$ es $$D=\frac{1}{2h}\left(\begin{array}{rrrrr} 0 & 1 & 0 & 0 & -1 \\ -1 & 0 & 1 & 0 & 0 \\ 0 & -1 & 0 & 1 & 0 \\ 0 & 0 & -1 & 0 & 1 \\ 1 & 0 & 0 & -1 & 0 \end{array}\right)$$
Las condiciones adicionales vienen dadas por la función f anterior:
$$ u_i(0) = f(x_i) = f(i/n) $$Observaciones:
- para $n=2$, es el mismo sistema del ejercicio 3.
- La matriz $D$ se muy similar a la matriz $A$ del ejercicio 1.
Ejercicio 6¶
- Resuelva el sistema anterior de ODEs para n=100:
con condiciones iniciales
$$ u_i(0) = f(x_i) = f(i/n). $$Puede usar la solución exacta para este problema lineal, usando la exponencial de matrices:
$$ \overrightarrow{u}(t) = \exp\left(-ctD\right)\cdot \overrightarrow{u}(0) $$o usar un método numérico de su elección, pero de orden superior a 1 (no se permite el método de Euler que es el objeto del siguiente apartado).
- Dibuje los vectores $\overrightarrow{u}_j(t)$, con el índice $j$ en el eje de las x (recordemos que $u_j$ aproxima a $u(j/n)$) y la cantidad $u$ en el eje de las $y$. La variable $t$ está fija. Utilice los siguientes valores de t: 0,0.1,0.2,0.3,0.4,0.5.
- De forma similar a lo ocurrido para la ecuación de ondas, reducimos nuestra EDP a un sistema de EDOs, que podemos resolver de forma exacta o usando por ejemplo
ode45.
- Usando integración temporal exacta, la solución de $\partial_{t}\overrightarrow{u} = -c D\cdot \overrightarrow{u}$, es $\partial_{t}\overrightarrow{u} =\exp\left(- c D\right)\cdot\overrightarrow{u}_0$.
L = 1;
Nx = 50;
h = L/Nx;
t0 = 0;
tf = 0.5;
Nt = 6;
ts = linspace(t0, tf, Nt);
xs = linspace(0, L - h, Nx);
u0 = init_fn(xs);
D = (1/h)*An(Nx);
nframes = 6;
hold on;
legends = cell(1, nframes);
ts = linspace(t0, tf, nframes);
for j = [1:Nt]
t = ts(j);
u_t = expm(-t*c*D)*u0';
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')
- Si usamos
ode45, los cambios necesarios respecto al ejercicio 3:
- Cambiamos la función
funque define el sistema de EDOs.- Cambiamos el array de valores iniciales
u0que define el sistema de EDOs.- Dibujamos la solución aproximada de otra forma distinta...
- y observamos que la solución no es muy buena...
L = 1;
Nx = 50;
h = L/Nx;
xs = linspace(0, L - h, Nx);
u_ini = init_fn(xs);
D = (1/h)*An(Nx);
fun = @(t,u)(-c*D*u);
t0 = 0;
tf = 0.5;
ts = linspace(t0, tf, nframes);
[ts, ys] = ode45(fun, ts, u_ini);
nframes = 6;
hold on;
legends = cell(1, nframes);
ts = linspace(t0, tf, nframes);
for j = [1:nframes]
t = ts(j);
u_t = ys(j,:);
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')
- La solución mejora sensiblemente, y se acerca a la solución exacta (del sistema de EDOs), si reducimos la tolerancia (por ejemplo
RelTol) en la llamada aode45.
L = 1;
Nx = 50;
h = L/Nx;
xs = linspace(0, L - h, Nx);
u_ini = init_fn(xs);
D = (1/h)*An(Nx);
fun = @(t,u)(-c*D*u);
t0 = 0;
tf = 0.5;
nframes = 6;
ts = linspace(t0, tf, nframes);
opts = odeset('RelTol',1e-6);
[ts, ys] = ode45(fun, ts, u_ini, opts);
hold on;
legends = cell(1, nframes);
ts = linspace(t0, tf, nframes);
for j = [1:nframes]
t = ts(j);
u_t = ys(j,:);
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')
Ejercicio 7¶
- Resuelva el sistema anterior de ODEs:
utilizando el método de Euler. Utilice los siguientes datos:
- Intervalo temporal $[t_0,t_f]=[0, 0.5]$
- Intervalo espacial $[0, 1]$
- Número de divisiones espaciales:
Nx=100. - Repita para tres valores del número de divisiones temporales
Nt: 25, 50, 100.- Dibuje en cada caso los vectores $\overrightarrow{u}_j(t)$, con el índice $j$ en el eje de las x (recordemos que $u_j$ aproxima a $u(j/n)$) y la cantidad $u$ en el eje de las $y$. La variable $t$ está fija. Utilice los siguientes valores de t: 0,0.1,0.2,0.3,0.4,0.5.
- ¿Cuál es el número de Courant ?
- Usamos un código muy similar al código del cuaderno "ondas". Un detalle que requiere atención es que tenemos que usar un paso pequeño para calcular con el método de Euler, pero al dibujar, no podemos dibujar todos los pasos usados en el cálculo. En el cuaderno de "ondas" podéis encontrar la solución (también podéis encontrar ejercicios similares a todos los ejercicios de este examen resueltos).
L = 1;
Nx = 100;
h = L/Nx;
D = (1/h)*An(Nx);
xs = linspace(0, L - h, Nx);
u_ini = init_fn(xs);
nframes = 6;
t0 = 0;
tf = 0.5;
for Nt = [25,50,200]
dt = (tf - t0)/Nt;
numero_de_Courant = c*dt/h;
u = zeros(Nt+1,Nx);
u(1,1:Nx) = u_ini;
uj = u_ini';
for j =[1:Nt]
uj1 = uj - c*dt*(D*uj);
u(j+1,:) = uj1;
uj = uj1;
end
legends = cell(1, nframes);
nplot = 1;
ts = linspace(t0, tf, Nt);
figure;
hold on;
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(sprintf('Numero de Courant: %.3f', numero_de_Courant))
end
Conclusiones¶
Al discretizar la variable espacial, pero usar un método exacto en la parte temporal, aparece difusión. Si usamos un método adaptativo con una tolerancia baja para la parte temporal, el resultado es casi el mismo.
Si c es negativo, entonces entonces el método es inestable. Si en vez de diferencias hacia atrás usamos diferencias hacia delante, el comportamiento se invierte (difusión para c negativo, inestable para c positivo). Por eso esta elección de métodos se llama upwind , porque es importante elegir la diferencia hacia atrás o hacia delante dependiendo de la dirección de propagación. Hay que ir "a favor del viento".
- Al discretizar tanto la variable espacial como la temporal, usando en ambos casos diferencias de dos puntos, observamos que el método es inestable para número de Courant mayor que 1, y tiene difusión si es menor que 1. Esta difusión es numérica, no es parte del problema sino consecuencia de la elección del método numérico, como se explica 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
- En el directorio con la solución del examen puedes ver varios vídeos mostrando la solución exacta y las soluciones aproximadas con varios métodos.