# Sistema de ecuaciones diferenciales

Un sistema de ecuaciones diferenciales es un conjunto de varias ecuaciones diferenciales con varias funciones incógnitas y un conjunto de condiciones de contorno. Una solución del mismo es un conjunto de funciones diferenciables que satisfacen todas y cada una de las ecuaciones del sistema. Según el tipo de ecuaciones diferenciales puede tenerse un sistema de ecuaciones diferenciales ordinarias o un sistema de ecuaciones en derivadas parciales.

## Sistemas de ecuaciones diferenciales ordinarias

En un sistema de ecuaciones diferenciales ordinarias de cualquier orden, puede ser reducido a un sistema equivalente de primer orden, si se introducen nuevas variables y ecuaciones. Por esa razón en este artículo sólo se consideran sistemas de ecuaciones de primer orden. Un sistema de ecuaciones diferenciales ordinarias de primer orden escrito en forma explícita es un sistema de ecuaciones de la forma:
   ${\begin{cases}{\cfrac {dx_{1}}{dt}}=F_{1}(x_{1},x_{2},\ldots ,x_{n};t)\\{\cfrac {dx_{2}}{dt}}=F_{2}(x_{1},x_{2},\ldots ,x_{n};t)\\\ldots \\{\cfrac {dx_{n}}{dt}}=F_{n}(x_{1},x_{2},\ldots ,x_{n};t)\end{cases}}$

### Reducción a un sistema de primer orden

Dado un sistema de ecuaciones diferenciales de orden n con m ecuaciones:
   $F_{i}\left(x_{j},{\frac {dx_{j}}{dt}},\ldots ,{\frac {d^{n}x_{j}}{dt^{n}}};t\right)=0\qquad {\mbox{con}}\ i,j\in \{1,2,\ldots ,m\}$
Existe un sistema equivalente de primer orden con a lo sumo (n+1) x m ecuaciones. Para ver esto consideremos un sistema en el que intervienen m funciones incógnitas x$_i$ y sus n derivadas, e introduzcamos un nuevo conjunto de variables y$_{i,k}$ definidas de la siguiente manera:
   $y_{i,k}(t):={\frac {d^{k}x_{i}(t)}{dt^{k}}}$
El sistema de primer orden equivalente en las variables y$_{i,k}$ resulta ser:
   ${\begin{cases}y_{i,k+1}={\cfrac {dy_{i,k}}{dt}}&k\in \{0,2,\ldots ,n-1\}\\F_{i}\left(y_{j,0},y_{j,1},\ldots ,y_{j,n};t\right)=0&i,j\in \{1,2,\ldots ,m\}\end{cases}}$
Como ejemplo de reducción de un sistema de ecuaciones diferenciales podemos considerar las ecuaciones de movimiento de la mecánica newtoniana de una partícula que es un sistema de segundo orden con tres ecuaciones:
   ${\begin{cases}m{\cfrac {d^{2}x}{dt^{2}}}=F_{x}\left(x,y,z;{\cfrac {dx}{dt}},{\cfrac {dy}{dt}},{\cfrac {dz}{dt}};t\right)\\m{\cfrac {d^{2}y}{dt^{2}}}=F_{y}\left(x,y,z;{\cfrac {dx}{dt}},{\cfrac {dy}{dt}},{\cfrac {dz}{dt}};t\right)\\m{\cfrac {d^{2}z}{dt^{2}}}=F_{z}\left(x,y,z;{\cfrac {dx}{dt}},{\cfrac {dy}{dt}},{\cfrac {dz}{dt}};t\right)\end{cases}}$
Si se introducen tres funciones incógnita nuevas que representan la velocidad, el sistema anterior se puede reducir a un sistema de primer orden y seis ecuaciones:
   ${\begin{cases}{\cfrac {dx}{dt}}=v_{x}(t),&m{\cfrac {dv_{x}}{dt}}=F_{x}\left(x,y,z;v_{x},v_{y},v_{z};t\right)\\{\cfrac {dy}{dt}}=v_{y}(t),&m{\cfrac {dv_{y}}{dt}}=F_{y}\left(x,y,z;v_{x},v_{y},v_{z};t\right)\\{\cfrac {dz}{dt}}=v_{z}(t),&m{\cfrac {dv_{z}}{dt}}=F_{z}\left(x,y,z;v_{x},v_{y},v_{z};t\right)\end{cases}}$

### Existencia y unicidad de la solución

Un sistema de ecuaciones diferenciales general tiene la forma
   (*)${\dot {X}}(t)=A(t)X(t)+f(t),\qquad X(t_{0})=X_{0},\ t_{0}\in [t_{1},t_{2}]$,
donde:
    ${f}(t)\in \mathbb {R} ^{n}$ es una función vectorial.
    $A(t)\in {\mathcal {L}}(\mathbb {R} ^{n},\mathbb {R} ^{n})$ es una función matricial.
Aplicando el Teorema de Picard-Lindelöf a los sistemas de ecuaciones diferenciales ordinarias de la forma (*), se establece la existencia y unicidad de la solución de dicho sistema en las que tanto la matriz $A(t)$ como la función $f(t)$ sean continuas en un intervalo compacto $[t_{1},t_{2}]$.
La demostración se basa en el hecho de que la función $F(t,X)=A(t)X+f(t)$ es Lipschitz respecto a $X$, para todo $X\in \mathbb {R} ^{n}$ si $t\in [t_{1},t_{2}]$. A partir del resultado del teorema se puede concluir que la sucesión de Picard
   $X_{k+1}=X_{0}+\int _{t_{0}}^{t}(A(s)X_{k}(s)+f(s))\ ds$
converge a la solución.

## Sistemas lineales de ecuaciones diferenciales ordinarias

Un sistema lineal de ecuaciones diferenciales con coeficientes constantes es un sistema de la forma:
   ${\begin{bmatrix}{x_{1}}'(t)\\{x_{2}}'(t)\\\dots \\{x_{n}}'(t)\end{bmatrix}}={\begin{bmatrix}a_{1,1}(t)&a_{1,2}(t)&a_{1,3}(t)&\dots &a_{1,n}(t)\\\dots &\dots &\dots &\dots &\dots \\a_{n,1}(t)&a_{n,2}(t)&a_{n,3}(t)&\dots &a_{n,n}(t)\end{bmatrix}}{\begin{bmatrix}{x_{1}}(t)\\{x_{2}}(t)\\\dots \\{x_{n}}(t)\end{bmatrix}}+{\begin{bmatrix}{b_{1}}(t)\\{b_{2}}(t)\\\dots \\{b_{n}}(t)\end{bmatrix}}$
Que de forma abreviada denotaremos como:
$X^{'}(t)=A(t)\cdot X(t)+{\overrightarrow {b(t)}}$
Cuando el sistema anterior verifica que ${\overrightarrow {b(t)}}=0$ decimos que el sistema anterior es homogéneo, en caso contrario diremos que es no homogéneo.
Llamamos sistema fundamental de ecuaciones al conjunto de $n$ soluciones independientes del sistema homogéneo. (Obsérvese que dichas soluciones son vectores).
La matriz que cuyas columnas representan cada una de las soluciones anteriores recibe el nombre de matriz fundamental $\phi (t)$ y verifica que $\phi (t)^{'}=A\cdot \phi (t)$.
Si dada una matriz fundamental existe un $t$ tal que la matriz resultante es la identidad diremos que dicha matriz es la matriz principal en $t$.
La resolución de estos sistemas es muy diversa. Sin embargo, la resolución de sistemas de ecuaciones diferenciales se simplifica mucho si las ecuaciones son de coeficientes constantes, o lo que es lo mismo, un sistema autónomo.

## Resolución de sistemas de n ecuaciones lineales de primer orden con coeficientes constantes
### El problema homogéneo

En el caso de los sistemas lineales, si la matriz del sistema es de coeficientes constantes podemos resolver el sistema usando el método de los autovalores.
Sea el sistema lineal homogéneo:
$X^{'}(t)=A\cdot X(t)$, con $A\in M_{n}$
Comenzamos encontrando los puntos críticos del sistema, es decir, aquellos valores $X$ donde:
$X^{'}(t)=0\Leftrightarrow A\cdot X(t)=0\Leftrightarrow X(t)\in Ker(A)$
Si $Ker(A)={0}$ entonces solo hay un punto crítico. Pero si$Ker(A)\neq 0$, entonces existe ${\overrightarrow {u}}\neq 0$ tal que $A{\overrightarrow {u}}=A(\beta {\overrightarrow {u}})=0$, $\beta {\overrightarrow {u}}$ es una recta por lo que tenemos rectas de puntos críticos. Estas son las soluciones fijas, no cambian con el paso del tiempo.
Por otra parte, si $\lambda$ es el autovalor de $A$ correspondiente al vector ${\overrightarrow {v}}$ tenemos $A{\overrightarrow {v}}=\lambda {\overrightarrow {v}}$
Cada uno de los $k$ autovalores de la matriz nos dará información para construir una solución $x_{k}(t)$ para el sistema.
La solución anterior depende de la naturaleza del autovalor:
Raíz real única: En este caso la solución para un valor $\lambda$ concreto viene dada directamente por $x_{\lambda }(t)=e^{\lambda t}$.
Raíz real de multiplicidad $j$: En este caso no podemos expresar la solución como antes, ya que si lo hacemos tendríamos información redundante.
En este caso las soluciones para un valor de $r$ se construyen multiplicando la función exponencial por $t$ como sigue: $x_{\lambda ,1}(t)=e^{\lambda t},\ x_{\lambda ,2}(t)=t\cdot e^{\lambda t},\ \dots \ ,\ x_{\lambda ,j}(t)=t^{j-1}\cdot e^{\lambda t}$
La solución es la combinación lineal de las anteriores: $x_{\lambda }(t)=\sum _{i=1}^{j}{\tilde {c}}_{i}\cdot x_{\lambda ,i}(t)$
Raíces complejas: Hay que tener en cuenta que siempre vienen en parejas de números conjugados.
En este caso usamos distinguimos la parte real y la imaginaria de la raíz y posteriormente usamos la fórmula de Euler y añadiendo constantes.
$x_{\lambda }(t)=e^{\lambda t}=e^{[Re(\lambda )+iIm(\lambda )]t}=e^{Re(\lambda )t}\cdot e^{iIm(\lambda )t}=e^{Re(\lambda )t}\ \cdot [c_{1}cos(Im(\lambda ))+c_{2}isen(Im(\lambda ))]$
La solución general al sistema homogéneo será, finalmente una combinación lineal de todas ellas multiplicadas cada una de ellas por su autovector correspondiente:
$X(t)=c_{1}\cdot x_{\lambda _{1}}{\overrightarrow {v_{1}}}+c_{2}\cdot x_{\lambda _{2}}{\overrightarrow {v_{2}}}+...+c_{k}\cdot x_{\lambda _{k}}{\overrightarrow {v_{k}}}$, donde $c_{1},c_{2},...,c_{n}$ son constantes.
La solución anterior puede reescribirse de forma matricial como:
Si se trata de un problema de valor inicial, los valores de $c_{1},c_{2},...,c_{n}$ se determinan a partir del dato inicial.
Ejemplo del problema homogéneo
Consideremos el siguiente problema de valor inicial:
   $X^{'}(t)={\begin{bmatrix}x_{1}^{'}(t)\\x_{2}^{'}(t)\end{bmatrix}}={\begin{bmatrix}1&2\\3&2\end{bmatrix}}{\begin{bmatrix}x_{1}(t)\\x_{2}(t)\end{bmatrix}}={\begin{bmatrix}1&2\\3&2\end{bmatrix}}\cdot X(t)\qquad X(0)={\begin{bmatrix}1\\1\end{bmatrix}}$
Los valores propios de la matriz son $a=4$ y $b=-1$ con los autovectores asociados ${\overrightarrow {v_{a}}}={\begin{bmatrix}2\\3\end{bmatrix}},{\overrightarrow {v_{b}}}={\begin{bmatrix}1\\-1\end{bmatrix}}$;
La matriz fundamental asociada a dicho sistema homogéneo sería: $\phi (t)={\begin{bmatrix}2e^{4t}&e^{-t}\\3e^{4t}&-e^{-t}\end{bmatrix}}$
Y la matriz fundamental principal quedaría: $\phi (t)={\begin{bmatrix}{\tfrac {2e^{5t}+3}{5e^{t}}}&{\tfrac {2e^{5t}-2}{5e^{t}}}\\{\tfrac {3e^{5t}-3}{5e^{t}}}&{\tfrac {3e^{5t}+2}{5e^{t}}}\end{bmatrix}}$
Otro método para resolver este tipo de sistemas es utilizar la exponencial matricial:
Sea $X'=A\cdot X,A\in M_{n}$:
Entonces la solución viene dada por la exponencial de la matriz A; para calcularla, utilizando las iteraciones de Picard, distinguiremos los siguientes casos:
* A es diagonal:
   $\ e^{A}={\begin{bmatrix}e^{\lambda _{1}}&0&0&\dots &0\\0&e^{\lambda _{2}}&0&\dots &0\\\dots &&&\dots \\0&0&0&\dots &e^{\lambda _{k}}\end{bmatrix}}$
La solución general será $\ e^{At}$.
* A no es diagonal, pero diagonaliza:
Si la matriz A diagonaliza puede expresarse de la forma: $A=SDS^{-1}$ donde D es una matriz diagonal y así tendríamos $e^{A}=Se^{D}S^{-1}$.
En este caso la solución general a nuestro problema será $e^{At}=Se^{Dt}S^{-1}$
* A no diagonaliza:
Si la matriz A no diagonaliza en general se puede escribir de la forma $A=SJS^{-1}$donde J es la forma de Jordan de la matriz A (Forma canónica de Jordan) y ahora sólo se necesita calcular la exponencial de un bloque de Jordan para terminar calculando la exponencial de la matriz J como la matriz que tiene por bloques la exponencial de cada uno de los bloques de Jordan.
Primero descompondremos cada bloque de Jordan en dos matrices más simples:
   $\ J_{r}(\lambda )=D_{r}(\lambda )+N_{r}$
donde $\ D_{r}$ es una matriz diagonal $r\cdot r$ que contiene a $\lambda$ en sus elementos diagonales y$\ N_{r}$ una matriz nilpotente $r\cdot r$ de la forma:
$N_{r}={\begin{bmatrix}0&1&0&\dots &0\\0&0&1&\dots &0\\\dots &&&\dots \\0&0&0&\dots &1\\0&0&0&\dots &0\end{bmatrix}}$ Utilizando lo anterior, podemos definir: $e^{N_{r}}={\begin{bmatrix}1&1&1/2&\dots &1/(r-1)!\\0&1&1&\dots &1/(r-2)!\\\dots &&&\dots \\0&0&0&\dots &1\end{bmatrix}}$
Por tanto, $e^{J_{r}(\lambda )}=e^{D_{r}(\lambda )}e^{N_{r}}$
Quedando así:
$e^{A}=Se^{J}S^{-1}$. Por tanto, la solución general al problema será: $e^{At}=Se^{Jt}S^{-1}$

### El problema no homogéneo

Sea el problema no homogéneo:
$X^{'}(t)=A\cdot X(t)+{\overrightarrow {b(t)}}$, con $A\in M_{n}$
Entonces la solución general es de la forma $X_{g}(t)=X_{gh}(t)+X_{p}(t)$, donde $X_{gh}(t)$ es la solución general del problema homogéneo y $X_{p}(t)$ es una solución particular del no homogéneo.
Una vez obtenida una solución general para el problema homogéneo $X_{gh}(t)$, (recordemos que: $X_{gh}(t)=\phi (t)\cdot {\overrightarrow {c}}$, siendo ${\overrightarrow {c}}$ un vector de constantes), para calcular $X_{p}(t)$ utilizamos el método de variación de las constantes, donde buscamos convertir los coeficientes $c_{i}$ del vector ${\overrightarrow {c}}$ en funciones.
Queremos soluciones de la forma $\phi (t)\cdot {\overrightarrow {\alpha (t)}}$, por lo que imponemos que sea solución de nuestro problema.
$[\phi \cdot {\overrightarrow {\alpha }}]'=A\cdot \phi \cdot {\overrightarrow {\alpha }}+{\overrightarrow {b}}$
Como $\phi (t)$ es solución fundamental, $\phi '(t)=A\cdot \phi (t)$ .
Si ahora calculamos dicha derivada por la regla de la cadena y utilizando la propiedad anterior de $\phi (t)$, tenemos que:
$[\phi \cdot {\overrightarrow {\alpha }}]'=\phi '\cdot {\overrightarrow {\alpha }}+\phi \cdot {\overrightarrow {\alpha '}}=A\cdot \phi \cdot {\overrightarrow {\alpha }}+\phi \cdot {\overrightarrow {\alpha '}}$.
Luego volviendo a nuestra ecuación y comparando, nos queda que $\phi (t)\cdot {\overrightarrow {\alpha '(t)}}={\overrightarrow {b(t)}}$. Despejamos:
${\overrightarrow {\alpha '(t)}}=(\phi (t))^{-1}\cdot {\overrightarrow {b(t)}}$. Integramos:
${\overrightarrow {\alpha (t)}}=\int _{t_{0}}^{t}(\phi (s))^{-1}\cdot {\overrightarrow {b(t)}}\,{\text{d}}s$
Por lo que hemos obtenido que la solución particular del problema no homogéneo es$X_{p}(t)=\phi (t)\cdot {\overrightarrow {\alpha (t)}}=\phi (t)\cdot \int _{t_{0}}^{t}(\phi (s))^{-1}\cdot {\overrightarrow {b(t)}}\,{\text{d}}s$.
Finalmente: $X_{g}(t)=\phi (t)\cdot {\overrightarrow {c}}+\phi (t)\cdot \int _{t_{0}}^{t}(\phi (s))^{-1}\cdot {\overrightarrow {b(t)}}\,{\text{d}}s$
(Recordemos que los valores de $c_{1},c_{2},...,c_{n}$ se determinan a partir del dato inicial del PVI).
Ejemplo del problema no homogéneo
Por último, estudiaremos un problema no homogéneo con autovalores imaginarios para ilustrar otro ejemplo.
   $X(t)={\begin{bmatrix}x_{1}^{'}(t)\\x_{2}^{'}(t)\end{bmatrix}}={\begin{bmatrix}2&-5\\1&-2\end{bmatrix}}{\begin{bmatrix}x_{1}(t)\\x_{2}(t)\end{bmatrix}}+{\begin{bmatrix}\csc t\\\sec t\end{bmatrix}}$
Primero hallamos la solución del caso homogéneo:
Los valores propios de la matriz son $\pm {i}$. La matriz fundamental principal asociada a este problema es: $\phi (t)={\begin{bmatrix}2\sin t+\cos(t)&-5\sin(t)\\\sin(t)&\cos(t)-2\sin(t)\end{bmatrix}}$
Ahora resolvemos la parte no homogénea:
Usando las fórmulas vistas, nuestro vector ${\overrightarrow {\alpha (t)}}$ resulta: ${\overrightarrow {\alpha (t)}}={\begin{bmatrix}{\tfrac {2}{3}}t+ln(({\tfrac {\operatorname {sen} t}{\cos t}})^{1/5})\\{\tfrac {-4}{5}}t+{\tfrac {2}{5}}ln(\operatorname {sen} t)\end{bmatrix}}$
Por lo visto anteriormente nuestra solución general del problema no homogéneo tendrá la siguiente forma sustituyendo los datos antes calculados:
   $X_{g}(t)=\phi (t)\cdot {\overrightarrow {c}}+\phi (t)\cdot {\overrightarrow {\alpha (t)}}$