La resolución de un sistema matricial recursivo no es del todo correcta

Question

Tengo un problema con un proyecto personal en el que quiero determinar la dinámica de una población.

Supongamos que la siguiente matriz de transición viene dada por

\begin{equation} T = \begin{pmatrix} 0.9 & 0 & 0 \\ 0.05 & 0.95 & 0 \\ 0.05 & 0.05 & 1 \end{pmatrix} \end{equation} La población se determina normalmente por \begin{equation} P_{n} = TP_{n-1} \end{equation} Sin embargo, en este problema la población en $n$ se determina mediante una ecuación diferente que implica una especie de matriz de entrada extra manual, es decir \begin{equation} P_{n} = TP_{n-1} + S \end{equation} Ahora, este problema será diferente de los problemas regulares de la matriz de transición (por lo que sé). Después de $n$ años, quiero tener una población óptima predeterminada, digamos a $n=f$ (último año). Además, la población inicial es, por supuesto, conocida. Ahora podemos simular este sistema matricial hasta $n=f$ de manera que podamos ver si hemos alcanzado la población óptima, digamos $P_f$ . Sin embargo, como esto también es interesante, quiero resolver este sistema matricial recursivo para $S$ de manera que encuentre el $S$ matriz tal que después de $n=f$ años se obtiene la población óptima. He intentado hacer lo siguiente

\begin{align*} P_{n} &= AP_{n-1} + S\\ &= A(AP_{n-2} + S) + S \\ &= A^2P_{n-2} + AS + S \\ &= A^2(AP_{n-3} + S) + AS + S \\ &= A^3P_{n-3} + A^2S + AS + S \end{align*} o expresado de otra manera \begin{align*} P_{n} = A^n P_0 + \sum_{j=1}^n A^{j-1}S \end{align*} Entonces, simplemente resolviendo para $S$ produciría

\begin{equation} S = \bigg(\sum_{j=1}^n A^{j-1}\bigg)^{-1} \bigg(P_n - A^n P_0\bigg) \end{equation}

Sin embargo, cuando intento simular mi respuesta, no parece funcionar del todo. Supongamos que tengo los siguientes datos

\begin{equation} P_0 = \begin{pmatrix} 85 \\ 110 \\ 0 \end{pmatrix} \end{equation}

\begin{equation} T = \begin{pmatrix} 0.9 & 0 & 0 \\ 0.05 & 0.95 & 0 \\ 0.05 & 0.05 & 1 \end{pmatrix} \end{equation}

\begin{equation} P_f = P_4 = \begin{pmatrix} 80 \\ 130 \\ 0 \end{pmatrix} \end{equation}

sustituyendo esta última en nuestra solución encontrada se obtiene

\begin{equation} S = \begin{pmatrix} 7 \\ 10.8 \\ -1 \end{pmatrix} \end{equation} y ejecutando una simulación se obtendrá \begin{equation} P_4 = \begin{pmatrix} 80.000 \\ 145.564 \\ 3.68 \end{pmatrix} \neq P_f \fin{según la ecuación} Sólo mi la entrada parece correcta. La segunda no. Sin embargo, cuando transpongo la matriz, la primera entrada se vuelve incorrecta mientras que la segunda se vuelve correcta. ¿Cómo es posible? ¿Hay algún error o simplemente mi petición no es posible?

EDIT : Mi problema ha sido resuelto, y se puede encontrar un método numérico en python como respuesta a esta pregunta. Sin embargo, también quería publicar mi código R para las personas que están interesadas.

T_matrix = t(matrix(c(0.9,  0,    0,
                      0.05, 0.95, 0,
                      0.05, 0.05, 1), nrow=3))

P = matrix(c(85,110,0), nrow = 3)
GP = matrix(c(80,130,0), nrow = 3)

years <- 4
T_mat <- NULL
for (j in 1:years){
        T_mat[[j]] <- T_matrix%^%(j-1)
}

S  <- (Reduce("+", T_mat) %>% inv()) %*% (GP - T_matrix%^%years %*% P)

pop <- NULL
for(i in 1: (years+1)){
   if (i == 1){
        pop[[1]] <- P
   }else{
        pop[[i]] <- T_matrix %*% pop[[i-1]] + S  
   }
}

final_population <- pop[[years + 1]] %>% print()
```

Answer 1

Debes estar cometiendo algún tipo de error de cálculo. Obtengo un valor completamente diferente para $S$ que sí produce el valor deseado cuando se ejecuta la iteración.

Aquí está mi script en python:

import numpy as np

T= np.array([[.9,0,0],[.05,.95,0],[.05,.05,1]])
P0 = np.array([85,110,0])
P4 = np.array([80,130,0])
A = sum(np.linalg.matrix_power(T, j) for j in range(4))
A1 = np.linalg.inv(A)
S=A1@(P4-np.linalg.matrix_power(T,4)@P0)
print("S=",S)

R = P0
for n in range(4):
    R= T@R+S
print("PF=",R)

Esto produce la salida:

S= [7.04608898   6.74717393 -10.04326291]
PF= [8.00000000e+01  1.30000000e+02 -5.32907052e-15]

Vemos que PF coincide con el valor dado P4 excepto por el error de redondeo.

Es imposible decir cuál es su error sin ver los detalles de su cálculo. El álgebra matricial es correcta.