Skip to main content
Fixed typo
Source Link
Michael E2
  • 244.7k
  • 18
  • 350
  • 774

Here's a way to automate the solution of a system that, like the OP's, is a linear system whose coefficient matrix happens to have a constant eigenvector (constant in that it is independent of $t$). If $\dot x = Ax$ and $A$ inis an $n \times n$ matrix that has some constant eigenvectors $v_1,v_2,\dots,v_k$, then completing them to a basis $v_1,v_2,\dots,v_k,v_{k+1},\dots,v_n$, we can form the transition matrix $$P = \left( v1 \mid v_2 \mid \cdots \mid v_n \right) \,.$$ Then with $x = Pu$, we get the transformed ODE $\dot u = P^{-1}AP\,u$, which hopefully DSolve can deal with. It should be able to find the first $k$ solutions, but it might fail to find the last $n-k$.

vars = {x1, x2};
ode = {x1'[t] == x2[t], 
   x2'[t] == 2 Cos[t]/(3 + 2 Sin[t]) x1[t] - 
     x2[t] (3 - 2 Cos[t] + 2 Sin[t])/(3 + 2 Sin[t])};
rhs = D[Through[vars[t]], t] /. First@Solve[ode, D[Through[vars[t]], t]];
If[AllTrue[rhs, Internal`LinearQ[#, Through[vars[t]]] &],   (* check linearity *)
 rhsA = CoefficientArrays[rhs, Through[vars[t]]][[2]],
 rhsA = Failure["nonlinear", <||>]
 ]
(*  SparseArray[<<4>>]    < check: was linear > *)

ev = Cases[Eigenvectors[rhsA], v_ /; FreeQ[v, t]]
If[Length[ev] < 1, Print["Bummer. No constant eigenvectors."]]
(*  {{-1, 1}}    < check: got a constant eigenvector > *)

pp = Transpose@Join[ev, NullSpace[ev]];  (* transition matrix *)
aa = Inverse[pp].rhsA.pp // Simplify;    (* transform coefficient matrix *)
dsol = DSolve[{u'[t], v'[t]} == aa.{u[t], v[t]}, {u, v}, t]  (* check DSolve works *)
(*  check: DSolve found a solution
  {{u -> Function[{t}, E^-t C[2] + C[1] (-3 + 2 Cos[t])], 
    v -> Function[{t}, C[1] (3 + 2 Sin[t])]}}
*)

pp.{u[t], v[t]} /. First[dsol] // Simplify;  (* pp.{u, v} yields {x1, x2} *)
xsol = Thread[vars -> (Function @@ {t, #} & /@ %)]
(*
  {x1 -> Function[t, 6 C[1] - E^-t C[2] - 2 C[1] Cos[t] + 2 C[1] Sin[t]], 
   x2 -> Function[t, E^-t C[2] + 2 C[1] Cos[t] + 2 C[1] Sin[t]]}
*)

ode /. xsol // Simplify    (* check solution *)    
(*  {True, True}  *)

Here's a way to automate the solution of a system that, like the OP's, is a linear system whose coefficient matrix happens to have a constant eigenvector (constant in that it is independent of $t$). If $\dot x = Ax$ and $A$ in an $n \times n$ matrix that has some constant eigenvectors $v_1,v_2,\dots,v_k$, then completing them to a basis $v_1,v_2,\dots,v_k,v_{k+1},\dots,v_n$, we can form the transition matrix $$P = \left( v1 \mid v_2 \mid \cdots \mid v_n \right) \,.$$ Then with $x = Pu$, we get the transformed ODE $\dot u = P^{-1}AP\,u$, which hopefully DSolve can deal with. It should be able to find the first $k$ solutions, but it might fail to find the last $n-k$.

vars = {x1, x2};
ode = {x1'[t] == x2[t], 
   x2'[t] == 2 Cos[t]/(3 + 2 Sin[t]) x1[t] - 
     x2[t] (3 - 2 Cos[t] + 2 Sin[t])/(3 + 2 Sin[t])};
rhs = D[Through[vars[t]], t] /. First@Solve[ode, D[Through[vars[t]], t]];
If[AllTrue[rhs, Internal`LinearQ[#, Through[vars[t]]] &],   (* check linearity *)
 rhsA = CoefficientArrays[rhs, Through[vars[t]]][[2]],
 rhsA = Failure["nonlinear", <||>]
 ]
(*  SparseArray[<<4>>]    < check: was linear > *)

ev = Cases[Eigenvectors[rhsA], v_ /; FreeQ[v, t]]
If[Length[ev] < 1, Print["Bummer. No constant eigenvectors."]]
(*  {{-1, 1}}    < check: got a constant eigenvector > *)

pp = Transpose@Join[ev, NullSpace[ev]];
aa = Inverse[pp].rhsA.pp // Simplify;    (* transform coefficient matrix *)
dsol = DSolve[{u'[t], v'[t]} == aa.{u[t], v[t]}, {u, v}, t]  (* check DSolve works *)
(*  check: DSolve found a solution
  {{u -> Function[{t}, E^-t C[2] + C[1] (-3 + 2 Cos[t])], 
    v -> Function[{t}, C[1] (3 + 2 Sin[t])]}}
*)

pp.{u[t], v[t]} /. First[dsol] // Simplify;  (* pp.{u, v} yields {x1, x2} *)
xsol = Thread[vars -> (Function @@ {t, #} & /@ %)]
(*
  {x1 -> Function[t, 6 C[1] - E^-t C[2] - 2 C[1] Cos[t] + 2 C[1] Sin[t]], 
   x2 -> Function[t, E^-t C[2] + 2 C[1] Cos[t] + 2 C[1] Sin[t]]}
*)

ode /. xsol // Simplify    (* check solution *)    
(*  {True, True}  *)

Here's a way to automate the solution of a system that, like the OP's, is a linear system whose coefficient matrix happens to have a constant eigenvector (constant in that it is independent of $t$). If $\dot x = Ax$ and $A$ is an $n \times n$ matrix that has some constant eigenvectors $v_1,v_2,\dots,v_k$, then completing them to a basis $v_1,v_2,\dots,v_k,v_{k+1},\dots,v_n$, we can form the transition matrix $$P = \left( v1 \mid v_2 \mid \cdots \mid v_n \right) \,.$$ Then with $x = Pu$, we get the transformed ODE $\dot u = P^{-1}AP\,u$, which hopefully DSolve can deal with. It should be able to find the first $k$ solutions, but it might fail to find the last $n-k$.

vars = {x1, x2};
ode = {x1'[t] == x2[t], 
   x2'[t] == 2 Cos[t]/(3 + 2 Sin[t]) x1[t] - 
     x2[t] (3 - 2 Cos[t] + 2 Sin[t])/(3 + 2 Sin[t])};
rhs = D[Through[vars[t]], t] /. First@Solve[ode, D[Through[vars[t]], t]];
If[AllTrue[rhs, Internal`LinearQ[#, Through[vars[t]]] &],   (* check linearity *)
 rhsA = CoefficientArrays[rhs, Through[vars[t]]][[2]],
 rhsA = Failure["nonlinear", <||>]
 ]
(*  SparseArray[<<4>>]    < check: was linear > *)

ev = Cases[Eigenvectors[rhsA], v_ /; FreeQ[v, t]]
If[Length[ev] < 1, Print["Bummer. No constant eigenvectors."]]
(*  {{-1, 1}}    < check: got a constant eigenvector > *)

pp = Transpose@Join[ev, NullSpace[ev]];  (* transition matrix *)
aa = Inverse[pp].rhsA.pp // Simplify;    (* transform coefficient matrix *)
dsol = DSolve[{u'[t], v'[t]} == aa.{u[t], v[t]}, {u, v}, t]  (* check DSolve works *)
(*  check: DSolve found a solution
  {{u -> Function[{t}, E^-t C[2] + C[1] (-3 + 2 Cos[t])], 
    v -> Function[{t}, C[1] (3 + 2 Sin[t])]}}
*)

pp.{u[t], v[t]} /. First[dsol] // Simplify;  (* pp.{u, v} yields {x1, x2} *)
xsol = Thread[vars -> (Function @@ {t, #} & /@ %)]
(*
  {x1 -> Function[t, 6 C[1] - E^-t C[2] - 2 C[1] Cos[t] + 2 C[1] Sin[t]], 
   x2 -> Function[t, E^-t C[2] + 2 C[1] Cos[t] + 2 C[1] Sin[t]]}
*)

ode /. xsol // Simplify    (* check solution *)    
(*  {True, True}  *)
Source Link
Michael E2
  • 244.7k
  • 18
  • 350
  • 774

Here's a way to automate the solution of a system that, like the OP's, is a linear system whose coefficient matrix happens to have a constant eigenvector (constant in that it is independent of $t$). If $\dot x = Ax$ and $A$ in an $n \times n$ matrix that has some constant eigenvectors $v_1,v_2,\dots,v_k$, then completing them to a basis $v_1,v_2,\dots,v_k,v_{k+1},\dots,v_n$, we can form the transition matrix $$P = \left( v1 \mid v_2 \mid \cdots \mid v_n \right) \,.$$ Then with $x = Pu$, we get the transformed ODE $\dot u = P^{-1}AP\,u$, which hopefully DSolve can deal with. It should be able to find the first $k$ solutions, but it might fail to find the last $n-k$.

vars = {x1, x2};
ode = {x1'[t] == x2[t], 
   x2'[t] == 2 Cos[t]/(3 + 2 Sin[t]) x1[t] - 
     x2[t] (3 - 2 Cos[t] + 2 Sin[t])/(3 + 2 Sin[t])};
rhs = D[Through[vars[t]], t] /. First@Solve[ode, D[Through[vars[t]], t]];
If[AllTrue[rhs, Internal`LinearQ[#, Through[vars[t]]] &],   (* check linearity *)
 rhsA = CoefficientArrays[rhs, Through[vars[t]]][[2]],
 rhsA = Failure["nonlinear", <||>]
 ]
(*  SparseArray[<<4>>]    < check: was linear > *)

ev = Cases[Eigenvectors[rhsA], v_ /; FreeQ[v, t]]
If[Length[ev] < 1, Print["Bummer. No constant eigenvectors."]]
(*  {{-1, 1}}    < check: got a constant eigenvector > *)

pp = Transpose@Join[ev, NullSpace[ev]];
aa = Inverse[pp].rhsA.pp // Simplify;    (* transform coefficient matrix *)
dsol = DSolve[{u'[t], v'[t]} == aa.{u[t], v[t]}, {u, v}, t]  (* check DSolve works *)
(*  check: DSolve found a solution
  {{u -> Function[{t}, E^-t C[2] + C[1] (-3 + 2 Cos[t])], 
    v -> Function[{t}, C[1] (3 + 2 Sin[t])]}}
*)

pp.{u[t], v[t]} /. First[dsol] // Simplify;  (* pp.{u, v} yields {x1, x2} *)
xsol = Thread[vars -> (Function @@ {t, #} & /@ %)]
(*
  {x1 -> Function[t, 6 C[1] - E^-t C[2] - 2 C[1] Cos[t] + 2 C[1] Sin[t]], 
   x2 -> Function[t, E^-t C[2] + 2 C[1] Cos[t] + 2 C[1] Sin[t]]}
*)

ode /. xsol // Simplify    (* check solution *)    
(*  {True, True}  *)