Ecuaciones diferenciales. Transformada de Laplace. Funciones analíticas de una matriz cuadrada.

Francisco Javier Pérez González
Departamento de Análisis Matemático
Universidad de Granada

Comandos básicos para resolver ecuaciones diferenciales

El comando DSolve[ ]

Este comando es uno de los más potentes de Mathematica y puede resolver de forma simbólica una gran cantidad de EDOs. Su sintaxis es

DSolve[listadeEDOs(y)CondicionesIniciales, listadeFunciones, VariableIndependiente]

e intenta resolver simbólicamente las EDO(s) con condiciones iniciales dadas en listadeEDOs(y)CondicionesIniciales para las funciones en listadeFunciones. La variable independiente es VariableIndependiente. Para el caso de una sola EDO sin condiciones iniciales y una sola función incógnita, los argumentos primero y segundo pueden ir sin llaves de lista. Cuando este comando tiene éxito, su salida es una lista de listas de reglas de sustitución (como Solve[ ]). Es frecuente que en la solución intervengan funciones especiales como, por ejemplo, la función seno integral, funciones de Bessel, funciones de Airy y otras (todas ellas documentadas en la ayuda de Mathematica). Además, dependiendo de cómo representemos las funciones, la respuesta dada por Mathematica puede presentar algunos cambios formales. Veamos algunos ejemplos.

sol1 = DSolve[y'[x] == 3  y[x], y[x], x] (* una ecuación sin condiciones iniciales y una función incógnita por lo que podemos prescindir de las llaves en los dos primeros argumentos *)

sol2 = 
          DSolve[y'[x] == 3  y[x], y, x]

Mathematica nos dice en ambos casos que la solución de la ED es la función y[x] = C[1] (e^(3  x)) donde C[1]  representa una constante arbitraria. Mathematica introduce tantas constantes como el orden de la ED. Cuando se especifican condiciones iniciales el número de constantes se va reduciendo. Observa que, en ambos casos, la solución es una lista de listas de reglas de sustitución (en este caso una única regla de sustitución porque sólo hay una función incógnita). La diferencia entre una y otra forma de presentar la solución está en su comportamiento cuando se aplican reglas de sustitución. Mientras que sol1 solamente sustituye y[x] por su valor, sol2 puede sustituir valores como y[0] o derivadas y'[x]. Observa la diferencia.

y[x] + y[0] + y'[x] /. sol1 [[1]]  
y[x] + y[0] + y'[x] /. sol2[[1]]

Si sustuimos la solución en la ED original obtenemos True.

y'[x] == 3  y[x]/.sol2[[1]]

Por tanto, si quieres sustituir directamente la solución proporcionada por Mathematica de una ED para hacer cálculos con ella, es preferible la segunda forma a la primera. Salvo este detalle, usar una u otra forma de presentar la solución es cuestión de gusto.

La solución antes obtenida no toma valores numéricos porque depende de C[1] y no es posible representarla gráficamente ni hacer cálculos numéricos con ella. Solamente podemos hacer representaciones gráficas o cálculos numéricos con soluciones particulares de una ED. Debes tener en cuenta que no es posible asignar directamente valores numéricos a las constantes porque son símbolos protegidos, aunque podemos hacerlo con reglas de sustitución apropiadas. En general, la forma usual de obtener una solución particular es especificar condiciones iniciales.

sol=DSolve[{y'[x] == 3  y[x], y[1] == 1}, y[x], 
        x] (* una ecuación con condiciones iniciales y una función incógnita por lo que no podemos prescindir de las llaves en el primer argumento *)

Para hacer cálculos y representaciones gráficas es conveniente definir la solución obtenida en la forma usual en que se definen las funciones en Mathematica.

y[x_] = y[x] /. sol[[1]] 
(* Aunque y[
                x] se repite en los dos lados de la definición, no entramos en un bucle infinito porque Mathematica ejecuta primero el comando de sustitución en el lado de la derecha, con lo cual sustituye y[
                x] por la solución obtenida, y después asigna la función y[x] a esta solución *)

Plot[y[x], {x, -1, 1}, PlotStyle -> Hue[ .5]];

Mathematica puede resolver EDLs y sistemas de EDLs con coeficientes constantes.

ecuacion = y''''[t] - 5  y'''[t] + 10  y''[t] - 
            10  y'[t] + 4  y[t] == 
          Cos[t];
y[t] /. DSolve[ecuacion, y, t][[1]]

En este caso Mathematica ha podido calcular de forma simbólica exacta las raíces del polinomio característico, que proporcionan la solución general de la ED homogénea, y también ha podido calcular una solución particular de la ED completa. La solución obtenida es la suma de ambas.

y''''[t] - 5  y'''[t] + 10  y''[t] - 10  y'[t] + 
            4  y[t] == 0 /. 
        y -> 
          Function[{t}, (E^t) C[
                  3] + (E^(2 t)) C[
                  4] + (E^t) C[2] Cos[
                  t] + (E^t) C[1] Sin[
                  t]]//Simplify

ecuacion /. 
        y -> 
          Function[{t}, 
            (1/20) (Cos[t] - 3 Cos[t] Cos[2 t] - 3 Sin[t] - 
                  Cos[2 t] Sin[t] + Cos[t] Sin[2 t] - 
                  3 Sin[t] Sin[2 t])]//Simplify

Si cambiamos un poquito lo anterior pueden pasar cosas como esta.

y[t] /. DSolve[y''''[t] - 5  y'''[t] + 10  y''[t] - 
              10  y'[t] + 4  y[t] == Exp[t^2], y, 
          t][[1]]

No es posible expresar la solución por medio de funciones elementales y Mathematica usa la función especial Erfi[ ]. La parte homogénea de la solución permanece igual. 
La solución general de una EDL homogénea con coeficientes constantes puede presentarse de formas diferentes dependiendo de que las raíces del polinomio característico sean más o menos fáciles de expresar por medio de radicales (cuando ello sea posible).

DSolve[y''''[t] + 5  y'''[t] + 10  y''[t] - 
            10  y'[t] + 4  y[t] == 0, y, t][[1]]

En este caso Mathematica puede calcular de forma simbólica todas las raíces del polinomio carecterístico porque éste es de grado 4, pero estas raíces son muy complicadas y Mathematica prefiere representarlas en la forma anterior. Si queremos que las muestre de forma explícita hay que pedírlo con el comando ToRadicals[ ].

ToRadicals[%]

Cuando las raíces del polinomio característico no pueden calcularse por medio de radicales, cosa frecuente cuando su grado es mayor o igual que 5, Mathematica las deja indicadas de forma simbólica. Para estos casos es preferible usar el comando NDSolve[ ] que veremos más adelante.

DSolve[y'''''[t] - 5  y''[t] + 4  y'[t] - 
            2  y[t] == 0, y, t][[1]]  

    ToRadicals[%]

También podemos, para el caso de ED lineales con coeficientes constantes, introducir algún valor aproximado en nuestra ED y así Mathematica calcula numéricamente las raíces del polinomio caracterísitco. Observa la diferencia entre trabajar simbólicamente y hacerlo numéricamente.

DSolve[y'''''[t] - 5*y''[t] + 4  y'[t] - 2  y[t] == 
          0.0, y, t][[1]]

Hemos hecho un pequeño cambio: hemos sustituido 0 por 0.0, pero ahora Mathematica proporciona la solución de forma explícita. Naturalmente, esta solución es la que interesa a efectos de cálculo pero, matemáticamente, es una solución aproximada de nuestra ED y Mathematica distingue siempre entre cálculos simbólicos (es decir, matemáticamente exactos) y cálculos aproximados (es decir, "exactos" desde un punto de vista práctico).

Para EDLs con coeficientes no constantes las cosas pueden complicarse. Aquí tienes algunos ejemplos.

DSolve[y'[x] - x y[x] == 1, y[x], x]
DSolve[x y''[x] + y[x] == 0, y[x], x]
DSolve[x^2 y''[x] + y[x] == 0, y[x], x]
DSolve[y''[x] - x y[x] == 0, y[x], x]

Salvo la tercera ecuación cuya solución general puede expresarse por medio de funciones elementales, en los demás casos las soluciones hacen intervenir funciones especiales.

Para resolver sistemas de EDs hay que escribir todas las ecuaciones y las condiciones iniciales si las hay.

DSolve[{x'[t] + y'[t] + y[t] == 
          1, x'[t] - y'[t] - x[t] == 1, x[0] == 0, 
        y[0] == 0}, {x, y}, t]

    {x'[t] + y'[t] + y[t] == 1, x'[t] - y'[t] - x[t] == 
            1, x[0] == 0, y[0] == 0} /. % // 
      Simplify

Clear["@"]

Con Mathematica es muy fácil estudiar la dependencia de la solución de una ED con respecto a las condiciones iniciales. Para ello basta con que definamos la solución de forma apropiada. Observa el siguiente ejemplo.

y[a_, b_, t_] =    
      y[t] /. DSolve[{y''[t] +  y'[t] + 2  y[t] == 0, 
              y[0] == a, y'[0] == b}, y, 
            t][[1]]

Ahora podemos usar esta función para representar cualquier solución particular dando valores a los parámetros "a" y "b".

Do[
      Plot[y[1, b, t], {t, 0, 10}, 
        PlotRange -> {- .7, 1.7}], {b, -2, 2, 1}]

Pero en una ED también podemos considerar variables los coeficientes de la ecuación. Esto es interesante porque, en la práctica, estos coeficientes se conocen con aproximación e interesa saber cómo la solución obtenida depende de ellos. Podemos considerar la situación más general: coeficientes y condiciones iniciales arbitrarias. Aquí tienes definida la solución más general del problema de Cauchy para la EDL de orden 1 completa con coeficientes constantes (pero arbitrarios) y término independiente también arbitrario.

Clear["@"]
solEDL1[a_, b_, c_, func_] :=    
      y[t] /. DSolve[{y'[t] +  c y[t] == func[t], 
              y[a] == b}, y, t][[1]]

Todo lo que hay que hacer ahora es elegir una función para sustituirla por "func".

f[t_] = t*Sin[t];
h[a_, b_, c_, t_] = solEDL1[a, b, c, f]

Si queremos soluciones particulares de esta ED podemos obtenerlas fácilmente.

h[0, 1, 3, t]

Solución de la ED y''[t] + c y'[t] + d y[t] = 0, con condiciones iniciales y[0] = a, y'[0] = b

Esta ED es un modelo para oscilaciones libres con amortiguamiento. Vamos a realizar un programa que proporcione las soluciones de esta ED dependiendo de la naturaleza de las raíces de la ecuación característica. A veces Mathematica presenta las soluciones de forma poco apropiada y hay que darle una ayudita para que las escriba como nosotros queremos.

Forma general de la solución sin considerar la naturaleza de las raíces características

Empezamos calculando la solución general de la ED y'' (t) + c y' (t) + d y(t) = 0. Sabemos que la solución debe ser de la forma \(TraditionalForm\`\(C\_1\) e\^\(\(a\_1\) t\) + \(C\_2\) e\^\(\(a\_2\) t\)\ donde\ C\_1\)y (Subscript[C, 2])  son constantes que dependen de las condiciones iniciales (¡y de los coeficientes de la ED!) y (Subscript[α, 1]) , (Subscript[α, 2]) son las raíces de la ecuación característica (z^2) + c z +  d  = 0. Dichas raíces son \(TraditionalForm\`\(\(-c\)\  ± \ \@\(c\^2 - 4  d\)\)\/2\).

Clear["@"]
sol[t_] = 
      y[t] /. DSolve[{y''[t] + c y'[t] + d y[t] == 0, 
              y[0] == a, y'[0] == b}, y, 
            t][[1]]

Vamos a agrupar las exponenciales y simplificamos.

solgeneral[a_, b_, c_, d_, t_] = 
      Collect[%, E^(x_), Simplify]

Forma de la solución cuando (c^2) - 4 d < 0 (raíces complejas)

Cuando (c^2) -  4 d < 0, se tiene que (Sqrt[(c^2) - 4  d]) = I (Sqrt[ 4  d - (c^2)]), las raíces de la ecuación característica son complejas e interesa obtener la solución anterior expresada por medio de las funciones trigonométricas seno y coseno. Para ello vamos a usar ComplexExpand[ ] que expresa las exponenciales complejas en su forma trigonométrica. Antes debemos quitar las raíces cuadradas  √((c^2) - 4 d) para evitar que este comando las desarrolle. Para ello podemos usar la regla  (Sqrt[(c^2) - 4 d])^(1/
                    2) -> I w que se aplica a las raíces del numerador y la regla  (Sqrt[(c^2) - 4 d])^(-1/
              2) -> -I w^(-1) que se aplica a la raíz del denominador. Naturalmente, nosotros sabemos que w = (Sqrt[4  d - (c^2)])y después de simplificar haremos la sustitución apropiada.

ComplexExpand[
      solgeneral[a, b, c, d, 
          t] /. {((c^2) - 4 d)^(1/2) -> 
            I*w, ((c^2) - 4 d)^(-1/2) -> -I*
              w^(-1)}]

Seguidamente agrupamos los términos en Sin[ ] y simplificamos.

Collect[%, Sin[x_], Simplify]

Finalmente, con la regla w  -> √(-(c^2) + 4 d) restituimos el valor de w.

tmp = % /. w -> (Sqrt[4 d - (c^2)])

Esta expresión ya es mejor. Observa que es de la forma λ Cos[ω] + μ Sin[ω]. Es conveniente escribir estas soluciones de otra forma.

ComplexExpand[
      Im[(μ +  I λ) (Cos[ω*t] +  
              I Sin[ω*t])]]

      
ComplexExpand[
        Im[(Sqrt[(λ^2) + (μ^2)]) 
            Exp[I Arg[μ +  I λ]]*Exp[I ω*t]]] // 
      TrigReduce

Hemos probado así que \(TraditionalForm\`λ\ Cos[ω\ t] + μ\ Sin[ω\ t] = \@\(λ\^2 + μ\^2\)\ Sin[
            t\ ω + Arg[i\ λ + μ]]\). Podemos definir un comando que pase a "forma senoidal" una función de la forma  λ Cos[t ω] + μ Sin[t ω]. Incluyo unas reglas de sustitución para que las exponenciales se saquen fuera de las raíces cuadradas, y uso ComplexExpand[ ] para simplificar los argumentos aunque esto presupone que todos ellos son números reales por lo que este comando debe usarse con cuidado.

formasenoidal[λ_*Cos[ω_] + μ_*
            Sin[ω_]] := 
((Sqrt[Simplify])[(λ^2) + (μ^2)]*
            Sin[ω + 
                ComplexExpand[
                  Arg[Simplify[μ +  
                        I λ]]]]) /. {(x___*E^y_*
                  z___)^(1/
                  2) :> (x*z)^(1/
                    2) E^(y/2), (E^
                  u_)^(1/2) :>  E^(u/
                  2)}

Expresemos en forma senoidal la solución antes obtenida.

formasenoidal[tmp]

Podemos definir esta solución en la forma siguiente.

solraicesimag[a_, b_, c_, d_, t_] = 
        2 (Sqrt[-(((b^2) + a b c + (a^2) d)/((c^2) - 4 d))]) (E^(-((c t)/2))) Sin[
            (1/2) (Sqrt[-(c^2) + 4 d]) t + 
              Arg[I a + ((2 b + a c)/(Sqrt[-(c^2) + 4 d]))]];

Hagamos c = 0 para eliminar el amortiguamiento. Obtenemos así la solución general de la ED y'' (t) + d y(t) = 0, que corresponde a oscilaciones libres sin amortiguamiento, con condiciones iniciales y(0) = a,  y' (0) = b.

solraicesimag[a, b, 0, d, t]

Se trata de un movimiento armónico simple con amplitud (Sqrt[(((b^2) + (a^2) d)/d)]), frecuencia angular (Sqrt[d])(radianes por segundo), periodo 2  π/(Sqrt[d])y fase inicial Arg[i a + (b/(Sqrt[d]))].

Forma de la solución cuando (c^2) - 4 d = 0 (una raíz real doble)

Este caso corresponde a una raíz real doble de la ecuación característica y es muy sencillo.

solraizdoble[a_, b_, c_, t_] = 
      y[t] /. DSolve[{y''[t] + c y'[t] + ((c^2)/4) y[t] == 0, 
              y[0] == a, y'[0] == b}, y, 
            t][[1]]

Forma general de la solución considerando la naturaleza de las raíces características

Podemos reunir los resultados obtenidos en una función que en cada caso proporcione la solución en la forma adecuada. Observa la forma en que escribo las condiciones usando un triple = para que siempre se evalúen.

solucion[a_, b_, c_, d_, t_] := 
      If[Negative[(c^2) - 4  d] === True ,     solraicesimag[a, b, c, d, t], 
        If[(c^2) - 4  d === 0,    
          solraizdoble[a, b, c, t],    
          solgeneral[a, b, c, d, t]]]

Comprobemos este comando en los varios casos posibles.

solucion[2, -3/2, 2, 3, t]

solucion[1, -3/2, 1, 15/4, t]

solucion[1, -1, 2, 2, t]

solucion[1, 3, 2, 1, t]

solucion[1, 3, 3, 1, t]

solucion[1, 1, c, 2, t] (* 
      solución general dependiendo del parámetro c *)

Parece que funciona correctamente.

Las soluciones de la ED y''[t] + c y'[t] +  
          d y[t] = λ Cos[ω t] cuando (c^2) - 4  d < 0

En este caso, se trata de oscilaciones forzadas. La fuerza externa aplicada es f(t) = λ cos(ω t). Calculemos la solución general.

y[t] /. DSolve[{y''[t] + c y'[t] +  
                d y[t] == λ Cos[ω t]}, y, 
          t][[1]]

Agrupemos términos.

Collect[%, {C[1], C[2], Cos[x_], Sin[x_]}, Simplify]

Observa que es de la forma "solución general de la ED homogénea + solución particular de la ED completa". Suponiendo que (c^2) - 4  d < 0, pasemos a forma senoidal las dos partes de la solución.

formasenoidal[
        Collect[(ComplexExpand[(E^(-(1/2) (c + (Sqrt[(c^2) - 4 d])) t)) C[
                        1] + (E^((1/2) (-c + (Sqrt[(c^2) - 4 d])) t)) C[2] /. ((c^2) - 4  d )^(1/2) -> I*w]), {Cos[
                x_], Sin[x_]}, Simplify] /. 
          w -> (Sqrt[4  d - (c^2)])];
formasenoidal[((λ (d - (ω^2)) Cos[t ω])/((d^2) + (c^2) (ω^2) - 2 d (ω^2) + (ω^4))) + ((c λ ω Sin[t ω])/((d^2) + (c^2) (ω^2) - 2 d (ω^2) + (ω^4)))];
FullSimplify[%% + %, λ > 0]

Vemos que el movimiento está formado por la superposición de un movimiento oscilatorio amortiguado que prácticamente cesa al cabo de un tiempo, pues para t suficientemente grande, (e^(-((c t)/2)))es muy pequeño, y de un movimiento movimiento armónico simple que permanece a lo largo del tiempo. El sistema acaba oscilando con la misma frecuencia ω que la fuerza exterior aplicada. Además la amplitud del movimiento depende de dicha frecuencia. Para valores grandes de t dicha amplitud es aproximadamente igual a \(TraditionalForm\`λ\ \@\(1\/\(c\^2\ ω\^2 + \((d - ω\^2)\)\^2\)\)\)y, cuando c es muy pequeño, la amplitud se hace muy grande cuando d = (ω^2), es decir cuando ω = (Sqrt[d]), lo que ocurre cuando la frecuencia de la fuerza exterior aplicada es igual a la frecuencia del sistema cuando oscila libremente sin amortiguamiento. Esto es lo que se conoce como resonancia.

El comando NDSolve[ ]

Este comando es quizás el más importante de los comandos numéricos de Mathematica y puede resolver de forma numérica una gran cantidad de EDOs. Su sintaxis es

NDSolve[listadeEDOsyCondicionesIniciales, listadeFunciones, VariableIndependienteyIntervalo, Opciones]

e intenta resolver numéricamente las EDO(s) con condiciones iniciales dadas, que es preciso especificar siempre, en listadeEDOsyCondicionesIniciales para las funciones en listadeFunciones. La variable independiente es VariableIndependiente.  Cuando este comando tiene éxito, su salida es una lista de funciones que Mathematica representa como InterpolatingFunction[ ].

El resultado de algunos comandos de Mathematica como NDSolve[ ] o Interpolation[listadeDatos], entre otros, se muestra como "InterpolatingFunction[intervalo,listadeDatos]" que representa una función que interpola los datos dados en el intervalo especificado. Veamos un ejemplo.

coseno = 
      Interpolation[Table[{x, Cos[x]}, {x, -4. , 4. ,  .25}]]

Si quieres ver los detalles de esta función (ocultos en el símbolo < >) puedes hacerlo como sigue.

InputForm[coseno]

Puedes evaluar, derivar y representar esta función de la forma usual. Sin olvidar que la función está definida solamente en el intervalo [-4,4]. Si pides un valor de la función fuera de este intervalo Mathematica muestra un mensaje de error e intenta hacer una extrapolación para calcularlo de forma aproximada.

coseno[Pi/2]

coseno[4.5]

coseno'[Pi/2]

Plot[{coseno[x], coseno'[x]}, {x, -Pi, Pi}, 
        PlotStyle -> {Hue[0], Hue[ .65]}];

Veamos cómo trabaja NDSolve[ ].

sol = 
      NDSolve[{x''[t] == -x[t], 
          x[0] == 1, x'[0] == 0}  , 
        x, {t, {-4, 4}}   ]

Para trabajar con esta función conviene representarla de forma más manejable.

x = x /. sol[[1]]

Plot[{x[t], x'[t]}, {t, -Pi, Pi}, 
        PlotStyle -> {Hue[0], Hue[ .65]}];

    Plot[{coseno[t], x[t]}, {t, -Pi, Pi}];

y = 
      y /. NDSolve[{y'''''[t] - 5  y''[t] + 4  y'[t] - 
                  2  y[t] == 0, 
              y[0] ==  .5, y'[0] == 
                1, y''[0] == -2, y'''[
                  0] == -3, y''''[0] == 4}, 
            y, {t, -1, 1}][[1]]

Plot[y[t], {t, -1, 1}];

Clear["@"]

y[t_] = 
      y[t] /. NDSolve[{y'[t] == 3.5  Sin[y[t]] -  .5  t^2, 
              y[0] ==  .5}, y, {t, -3, 3}][[1]]

Plot[y[t], {t, -3, 3}];

Clear["@"]

Naturalmente, podemos usar NDSolve[ ] para sistemas de EDs.

NDSolve[{x'[t] == -y[t] - x[t]^2, y'[t] == 
          2  x[t] - y[t], x[0] == 1, y[0] == 1}, {x, y}, {t, 0, 
        10}]

Podemos expresar las soluciones obtenidas de forma más cómoda.

{x[t_], 
        y[t_]} = {x[t], y[t]}/.%[[1]]

Plot[{x[t], y[t]}, {t, 0, 10}, PlotStyle -> {Hue[0], Hue[ .65]}, 
        AspectRatio -> 1/GoldenRatio];

Clear["@"];
aproxsols = 
        NDSolve[{y'[t] == 
              0.7 y[t] - 0.3 y[t] x[t], x'[t] == -0.44 x[t] + 
                0.08 x[t] y[t], y[0] == 4, x[0] == 1.2}, {x, 
            y}, {t, 0, 50}];
\(\(\({x[t_], y[t_]} = {x[t], y[t]} /. aproxsols[\([1]\)];\)\(
\)
    \) (*\ 
      este\ es\ un\ modelo\ poblacional "depredador(y)-presa(x)"\ *) \)
Plot[{x[t], y[t]}, {t, 0, 50}, 
        PlotStyle -> {{Hue[ .65], Thickness[0.01]}, {Hue[0], 
              Thickness[0.01]}}, 
        AspectRatio -> 1/GoldenRatio];

Clear["@"]

Transformada de Laplace

El comando

LaplaceTransform[función[variable1],variable1,variable2]

calcula la transformada de Laplace de "función" y la expresa como una función de "variable2".  Este comando no indica el campo de convergencia de la transformada de Laplace.

LaplaceTransform[(t^2) Exp[a*t], t, s]
LaplaceTransform[t*Sin[a*t], t, s]

A partir de la versión 4 de Mathematica podemos especificar condiciones (con "Assumptions") para ayudar a Mathematica en los cálculos.

LaplaceTransform[t*Sin[a*t], t, s, 
      Assumptions -> a > 0]

Para saber el campo de convergencia debemos hacer el cálculo directamente.

(∫_0^∞) (t^3) 
        Exp[a*t]*(E^(-s*
              t)) \[DifferentialD]t

Mathematica nos dice que el dominio de convergencia de la integral es el semiplano Re[s] > Re[a]. En las aplicaciones de la transformada de Laplace a la resolución de EDs no hay que preocuparse por el dominio de convergencia.

Mathematica conoce las propiedades de la transformada de Laplace.

\(\(\(LaplaceTransform[c*a[t] + d*b[t], t, s] (*\ linealidad\ *) 

    LaplaceTransform[Integrate[h[u], {u, 0, t}], t, s] (*\ 
      transformada\ de\ una\ integral\ *) 

    LaplaceTransform[\(f'\)[t], t, s] (*\ 
      transformada\ de\ una\ derivada\ primera\ *) 

    LaplaceTransform[\(f''\)[t], t, s] (*\ 
      transformada\ de\ una\ derivada\ segunda\ *) 

    D[LaplaceTransform[f[t], t, s], s] (*\ 
      derivada\ primera\ de\ una\ transformada\ *) 

    D[LaplaceTransform[f[t], t, s], {s, 2}]\)\( (*\ 
      derivada\ segunda\ de\ una\ transformada\ *) \)\)\)

La transformada inversa se calcula con el comando InverseLaplaceTransform[ ].

InverseLaplaceTransform[((s^2)/((s^4) + 1)), s, t]

Podemos comprobar este resultado.

LaplaceTransform[%, t, s] // Simplify

Podemos simplificar bastante el resultado obtenido.

%% // FullSimplify // ComplexExpand

Podemos usar la transformada de Laplace para resolver EDLs.

Clear["@"]
  

      
Solve[(LaplaceTransform[y''[t] + y'[t] + y[t] - Cos[t], 
              t, s] /. {y[0] -> 0, y'[0] -> 0}) == 0, 
      LaplaceTransform[y[t], t, s]]

      
LaplaceTransform[y[t], t, s] /. %[[1]]

      
y[t_] = InverseLaplaceTransform[%, s, t]

      
{y''[t] + y'[t] + y[t] == Cos[t], 
        y[0] == 0, y'[0] == 0} // Simplify

Cálculo de funciones analíticas de una matriz cuadrada

Mathematica dispone de un comando, MatrixExp[matriz], que calcula la exponencial de una matriz.

A = {{2,-1,2,0}, {0,3,-1,0}, {0,1,1,0}, {0,1,-3,5}}; 
MatrixExp[A*t]//MatrixForm

Que yo sepa no hay comandos que permitan calcular directamente otras funciones analíticas de una matriz. Supongamos que queremos calcular f(A t) donde f es una función analítica y A es la matriz anterior. El procedimiento que se sigue para ello es el siguiente.  
1) Se calculan los valores propios de la matriz.

Eigenvalues[A]

Hay un valor propio simple λ = 5 y un valor propio triple λ=2.

2) Se calcula el polinomio interpolador cuyo valor en λ = 5 coincide con f(5 t) y cuyo valor y el de sus dos primeras derivadas en λ = 2 coincide con el calor de f(z t) y el de sus dos primeras derivadas respecto a z en el punto z = 2. Para ello Mathematica dispone de un comando

InterpolatingPolynomial[{{(Subscript[x, 1]),  {(Subscript[f, 1]),  (Subscript[df, 1]),  
          (Subscript[ddf, 1]),  ... }} ,{(Subscript[x, 2]), {(Subscript[f, 2]),  (Subscript[df, 2]), (Subscript[ddf, 2]) ,  ... }},{(Subscript[x, 3]), {(Subscript[f, 3]),  (Subscript[df, 3]), (Subscript[ddf, 3]) ,  ... }}... }, var]

que calcula un polinomio en la variable "var" cuyo valor en (Subscript[x, j]) es igual a (Subscript[f, j]), el valor de su primera derivada en (Subscript[x, j]) es igual a (Subscript[df, j]), el valor de su segunda derivada en (Subscript[x, j]) es igual a (Subscript[ddf, j]) ,  etcétera. Cuando solamente se especifica el valor (Subscript[f, j]) la correspondiente lista debe ser de la forma {(Subscript[x, j]),  (Subscript[f, j])}.

(* 
      En nuestro caso sería InterpolatingPolynomial[{{5, 
              f[5 t]}, {2, {f[2 t],  D[f[z t], z] /. z -> 2,  
                D[f[z t], {z, 2}] /. z -> 2}}}, 
          w] *) 
 (* 
      Supongamos que queremos calcular f[A t] donde f[
            z] = ((Exp[(z^2)] - 1)/(z^2)),  
      f[0] = 1.  Esta función es analítica y su serie de potencias converge en todo el plano complejo *) 
f[
          z_] = ((Exp[(z^2)] - 1)/(z^2));

    InterpolatingPolynomial[{{5, 
            f[5 t]}, {2, {f[2 t],  D[f[z t], z] /. z -> 2,  
              D[f[z t], {z, 2}] /. z -> 2}}}, w];

    p[w_] = Collect[%, w, Simplify]

{p[5], p[2], p'[2], p''[2]} - {f[5 t], f[2 t],  
          D[f[z t], z] /. z -> 2,  
          D[f[z t], {z, 2}] /. z -> 2}//Simplify

3) Se evalúa el polinomio interpolador en la matriz A. Ten en cuenta que Mathematica considera A como una lista, por tanto (A^2) no es el producto matricial de A por A sino el cuadrado de los elementos de la lista A.

(A^2)

Por eso, para calcular p(A)  hay que hacerlo de forma apropiada. Para ello extraemos los coeficientes del polinomio y los multiplicamos por las respectivas potencias de A y sumamos.

\(f[A\ t_] = 
      CoefficientList[p[w], w] . 
          Table[MatrixPower[A, k\  - \ 1], 
		{k, 1, 4}] // 
        Simplify\)

En este caso concreto, podemos calcular de forma más sencilla f[A t] usando el comando MatrixExp[ ]. Teniendo en cuenta cómo está definida la función f y que las operaciones en la variable z se corresponden con operaciones (productos y sumas) en la matriz A, podemos proceder como sigue.

(MatrixExp[A . A (t^2)] - IdentityMatrix[4]) . Inverse[A t] . 
        Inverse[A t] // Simplify

Vamos a hacer este cálculo por medio de transformadas de Laplace. Todo lo que debemos calcular es la exponencial de la matriz A . A (t^2). Sabemos que, dada una matriz M, la exponencial (e^(M t)) viene dada por la transformada de Laplace inversa de la matriz (((s I - M))^(-1)), donde el cáculo de dicha transformada se hace elemento a elemento.

B = 
        Inverse[s*IdentityMatrix[4] - A . A] // 
          Simplify;
U[t_] = 
        Table[InverseLaplaceTransform[B[[j, k]] // Simplify, s, t] // 
            Simplify, {j, 1, 4}, {k, 1, 4}]  ;
(U[(t^2)] - IdentityMatrix[4]) . Inverse[A t] . Inverse[A t] // 
      Simplify