Los números de Fibonacci se definen así:
$$F_0=0, F_1=1$$$$F_{n} = F_{n-1} + F_{n-2}$$y tenemos $F_2=1,F_3=2,F_4=3,F_5=5,F_6=8\dots$
También se puede definir el vector $$\mathbf{v}_n=\left( \begin{array}{c} F_{n}\\F_{n-1} \end{array} \right)$$ con dos números de Fibonacci consecutivos, y entonces tenemos las ecuaciones $$ \left( \begin{array}{c} F_{n}\\F_{n-1} \end{array} \right)= \left(\begin{array}{cc} 1 & 1\\ 1 & 0 \end{array}\right) \cdot \left( \begin{array}{c} F_{n-1}\\F_{n-2} \end{array} \right) $$
Iterando, podemos calcular el número de Fibonacci $F_{n}$ como $$ \left( \begin{array}{c} F_{n}\\F_{n-1} \end{array} \right)= \left(\begin{array}{cc} 1 & 1\\ 1 & 0 \end{array}\right)^{n-1} \cdot \left( \begin{array}{c} F_{1}\\F_0 \end{array} \right)= \left(\begin{array}{cc} 1 & 1\\ 1 & 0 \end{array}\right)^{n-1} \cdot \left( \begin{array}{c} 1\\0 \end{array} \right) $$
Vamos a comparar varias técnicas para calcular las potencias de la matriz $$ M = \left(\begin{array}{cc} 1 & 1\\ 1 & 0 \end{array}\right) $$
Usaremos en ambos casos funciones que reciben dos argumentos:
M = [1,1; 1,0]
F10 = [1; 0]
#Método 1
# ATENCION: para usar este código desde matlab
# es necesario compiar este código en un archivo de nombre
# potencia1.m y que ese archivo este en el PATH
function P = potencia1(A,n)
%Calcula la potencia A^n de una matriz A a un exponente n
%
%Utiliza una iteración sencilla, multiplicando la identidad por A n veces
%
%Para calcular A^n tiene que hacer n productos de matrices 2x2.
%Cada producto de matrices 2x2 necesita 8 multiplicaciones luego el
%número de multiplicaciones necesarias es 8*n
[d,d] = size(A);
P = eye(d);
for i = [1:n]
P = A*P;
end
end
for i = [1:10]
sprintf('M a la potencia %d',i)
potencia1(M,i)
end
%Los numeros de Fibonacci F_n y F_{n-1}
n = 6
potencia1(M,n-1)*F10
%Método 2
% ATENCION: para usar este código desde matlab
% es necesario compiar este código en un archivo de nombre
% potencia2.m y que ese archivo este en el PATH
function P = potencia2(A,n)
%Calcula la potencia A^n de una matriz A a un exponente n
%
%TAREA: Explica cómo funciona ...
if n==0
[d,d] = size(A);
P = eye(d);
else
resto = mod(n, 2);
cociente = floor(n/2);
An2 = potencia2(A, cociente);
if resto==0
P = An2 * An2;
else
P = An2 * An2 * A;
end
end
end
for i = [1:10]
sprintf('M a la potencia %d',i)
potencia2(M,i)
end
%Los numeros de Fibonacci F_n y F_{n-1}
n = 6
potencia2(M,n-1)*F10
Usa ahora la diagonalización de M. Debajo tienes los comandos necesarios
%Método 3
%TAREA: escribe el cuerpo de la funcion
function P = potencia3(A,n)
%Calcula la potencia A^n de una matriz A a un exponente n
%
% Utiliza la diagonalización A = P*D*Pinv
% y la formula A^n = P*(D^n)*Pinv
end
B = [0,1; 4,0]
[avecs, avals] = eig(B);
P = avecs
D = avals
Pinv = inv(P)
%Recuperamos la matriz B mediante P*D*Pinv (con error numérico)
P*D*Pinv
%Comparamos B@B@B con P*(D*D*D)*Pinv
B*B*B
% de nuevo error numérico
P*(D*D*D)*Pinv
% Para obtener potencias de matrices diagonales podemos usar esta función
function Dn = potencia_diagonal(D,n)
%Devuelve D^n, para D una matriz diagonal y n un número entero
%
%Si D no es diagonal, o n no es entero, devuelve un resultado incorrecto
%Cuando el argumento de diag es una matriz, devuelve los elementos de la diagonal
valores_diagonal = diag(D);
d = length(D);
Dn = zeros(d);
for i = [1:d]
Dn(i,i) = valores_diagonal(i)^n;
end
end
potencia_diagonal(D,3)
Amanece en nuestro crucero con las 1000 personas del pasaje en los camarotes, en el nivel 2. Cada hora, su ubicación puede haber cambiado:
Definimos la cantidad de personas $x_n$ en cubierta, $y_n$ en el nivel 1 y $z_n$ en el nivel 2, transcurridas $n$ horas.
Planteamos las ecuaciones de cantidad de personas en cada nivel como un sistema matricial $$\left(\begin{array}{c} x_{n + 1}\\ y_{n + 1}\\ z_{n + 1} \end{array}\right) = \left(\begin{array}{ccc} 0.4 & 0.3 & 0.4\\ 0.3 & 0.5 & 0.4\\ 0.3 & 0.2 & 0.2 \end{array}\right) \cdot \left(\begin{array}{c} x_n\\ y_n\\ z_n \end{array}\right)$$ Este modelo puede predecir una cantidad de personas en cubierta que no es un entero. No es un problema, porque debe interpretarse como un número medio.
Al igual que antes podemos expresar la cantidad de gente en la hora n directamente en función de la cantidad de gente en la hora inicial $n=0$.
$$\left(\begin{array}{c} x_{n }\\ y_{n}\\ z_{n} \end{array}\right) = \left(\begin{array}{ccc} 0.4 & 0.3 & 0.4\\ 0.3 & 0.5 & 0.4\\ 0.3 & 0.2 & 0.2 \end{array}\right)^n \cdot \left(\begin{array}{c} x_0\\ y_0\\ z_0 \end{array}\right) = \left(\begin{array}{ccc} 0.4 & 0.3 & 0.4\\ 0.3 & 0.5 & 0.4\\ 0.3 & 0.2 & 0.2 \end{array}\right)^n \cdot \left(\begin{array}{c} 0\\ 0\\ 1000 \end{array}\right) $$Usando potencias de matrices, calcula la cantidad esperada de gente en cubierta al cabo de 1,2,3 y 4 horas.
M = [0.4, 0.3, 0.4; 0.3, 0.5, 0.4; 0.3, 0.2, 0.2]
Una matriz con entradas positivas, tal que la suma de las entradas de cada columna es 1, se denomina "matriz estocástica". Las matrices estocásticas aparecen yosal modelizar procesos estocásticos en los que un sistema puede moverse de uno a otro estado (una persona cambia de cubierta, un barco puede estar en espera, cargando, descargado, repostando o en ruta, etc).
Es un resultado conocido que una matriz estocástica tiene el autovalor 1, y el resto de autovalores, sean reales o complejos, tienen módulo menor que 1.
Podemos definir la exponencial de un número $x$ mediante la serie:
$$e^x = \exp(x) = \sum_{n=0}^\infty \frac{1}{n!}x^n$$y podemos definir la exponencial de una matriz $x$ mediante la serie:
$$e^M =\exp(M) = \sum_{n=0}^\infty \frac{1}{n!}M^n$$Como vistéis en clase, es más práctico diagonalizar la matriz para calcular:
$$e^M = \sum_{n=0}^\infty \frac{1}{n!}P\cdot D^n \cdot P^{-1} = P\cdot \left(\begin{array}{ccc} \sum_{n=0}^\infty \frac{(d_1)^n}{n!} & \dots & 0\\ \vdots & \dots & \vdots\\ \vdots & \dots & \sum_{i=0}^\infty \frac{(d_3)^n}{n!} \end{array}\right) \cdot P^{-1} = P\cdot \left(\begin{array}{ccc} \exp(d_1) & \dots & 0\\ \vdots & \dots & \vdots\\ \vdots & \dots & \exp(d_3) \end{array}\right) \cdot P^{-1} $$function eM = exponencial1(M, truncar=10)
%Recibe como argumento una matriz M y devuelve la exponencial de M
%
%Suma varios términos de la serie (1/n!)*M^n, terminando en el valor "truncar"
%que se pasó como argumento. Calcula las potencias M^n sin diagonalizar M.
# TAREA: completa la función
end
function eM = exponencial2(M, truncar=10)
%Recibe como argumento una matriz M y devuelve la exponencial de M
%
%Utiliza la diagonalización A = P*D*Pinv
%y la formula exp(A) = P*exp(D)*Pinv
# TAREA: completa la función
end
M1 = [1,1; 1,0]
M2= [0.4, 0.3, 0.4; 0.3, 0.5, 0.4; 0.3, 0.2, 0.2]
Busca el desarrollo en serie del coseno, y repite el ejercicio anterior para el cálculo del coseno de una matriz.