Results 1 to 6 of 6

Math Help - steady state - coupled ODEs

  1. #1
    Member
    Joined
    Jan 2008
    Posts
    114

    steady state - coupled ODEs

    I have a system of 6 coupled ODEs:

     W\frac{dl_E}{dt} = -l_E\rho_F + \psi_Ll_B

     \frac{dl_B}{dt} = \rho_Fl_E - \psi_Ll_B - \chi_Ll_B

     \frac{dl_I}{dt} = \chi_Ll_B - \omega_Ll_I

     \sigma\frac{d\rho_F}{dt} = \gamma_{rr}\sigma\rho_I - \chi_0\rho_F - (ml_E\rho_F - m\psi_Ll_B + \frac{m\chi_Ll_B\rho_F}{1 - \rho_F})

     \sigma\frac{d\rho_I}{dt} = \frac{\gamma_S\sigma}{K + c} + \chi_0f\rho_F + f(1 + \frac{\rho_F}{1 - \rho_F})(m\chi_Ll_B) - \sigma\gamma_{rr}\rho_I

     \frac{dc}{dt} = \gamma(\omega_LRl_I) - \lambda(c - 1)



    And I am trying to find the steady state for each one...I know this involves setting each equation equal to zero, however not sure where to go from there as the equations are coupled so it seems a little complicated. Could anyone please give me a hint as to what method I just approach?

    Thanks in advance!
    Follow Math Help Forum on Facebook and Google+

  2. #2
    Super Member
    Joined
    Aug 2008
    Posts
    903
    Divide through to get just derivatives on left side, set the right side to zero then plug them into Mathematica using Solve or Reduce and solve for the six variables. I'm assuming the right side is all constants and six variables.

    Oh yeah, if Solve or Reduce can't solve it, then I'd assign real values to the constants and then use FindRoot or NSolve to solve it and then if I had to solve it for the general case, I would next start removing the assigned values one by one and replacing them with the general term until Solve or Reduce could solve the semi-general case.
    Last edited by shawsend; December 9th 2009 at 05:22 AM.
    Follow Math Help Forum on Facebook and Google+

  3. #3
    Member
    Joined
    Jan 2008
    Posts
    114
    Thank you for the reply shawsend, I have already solved them in Matlab but I am trying to solve them analytically...is this even possible?!
    Follow Math Help Forum on Facebook and Google+

  4. #4
    Super Member
    Joined
    Aug 2008
    Posts
    903
    Very good then. I'm not familiar with Matlab. However, if you supply the equation set to Mathematica's Solve or Reduce, without specifying particular values for the constants, they MAY return an explicit analytical solution to it. Sometime today, I'll check it ok.
    Follow Math Help Forum on Facebook and Google+

  5. #5
    Senior Member
    Joined
    Mar 2009
    Posts
    378
    The problem with non-linear ODEs is that it is quite difficult to find a solution for the original system. However, if you linearize about an equilibrium you can find a solution that approximates the actual solution about that equilibrium.

    To do this you will rewrite the system as
    <br />
X' = JF(X^*)\cdot\left(X - X^*\right)<br />
    where X = \begin{bmatrix}l_e \\l_b \\l_I \\ \rho_F \\ \rho_I \\ c\end{bmatrix}, F(X) is the RHS, JF is the Jacobian, and X^* is the equilibrium solution.

    EDIT: To find the equilibrium it is probably best to use Newton's Method, since this is non-linear.
    Follow Math Help Forum on Facebook and Google+

  6. #6
    Super Member
    Joined
    Aug 2008
    Posts
    903
    Hi. I solved for the equilibrium points but first I changed some of the variable names since I is for imaginary in Mathematica and also not good to use any variable names starting with capital letters in Mathematica so everywhere a capital letter, I changed it to lower case. Here's the Solve code but since I used the same subscripted names, I have to first change it to raw input form to cut and paste it into this thread so it's hard to read. However if you then cut and paste this code into Mathematica then select it and do a Cell/Convert To/Stanard Form, you can recover all the subscripted formatting.

    Code:
    myeqtpt = First[Solve[{(-Subscript[l, e])*(Subscript[\[Rho], f]/w) + 
           Subscript[\[Psi], l]*(Subscript[l, b]/w) == 0, 
         Subscript[\[Rho], f]*Subscript[l, e] - Subscript[\[Psi], l]*
            Subscript[l, b] - Subscript[\[Chi], l]*Subscript[l, b] == 0, 
         Subscript[\[Chi], l]*Subscript[l, b] - Subscript[\[Omega], l]*
            Subscript[l, i] == 0, 
         Subscript[\[Gamma], rr]*Subscript[\[Rho], i] - Subscript[\[Chi], 0]*
            (Subscript[\[Rho], f]/\[Sigma]) - 
           (m*Subscript[l, e]*Subscript[\[Rho], f] - m*Subscript[\[Psi], l]*
              Subscript[l, b] + (m*Subscript[\[Chi], l]*Subscript[l, b]*
               Subscript[\[Rho], f])/(1 - Subscript[\[Rho], f]))*(1/\[Sigma]) == 
          0, Subscript[\[Gamma], S]/(k + c) + Subscript[\[Chi], 0]*f*
            (Subscript[\[Rho], f]/\[Sigma]) + 
           f*(1 + Subscript[\[Rho], f]/(1 - Subscript[\[Rho], f]))*
            ((m*Subscript[\[Chi], l]*Subscript[l, b])/\[Sigma]) - 
           Subscript[\[Gamma], rr]*Subscript[\[Rho], i] == 0, 
         \[Gamma]*(Subscript[\[Omega], l]*r*Subscript[l, i]) - \[Gamma]*(c - 1) == 
          0}, {c, Subscript[\[Rho], i], Subscript[\[Rho], f], 
         Subscript[l, i], Subscript[l, b], Subscript[l, e]}]]
    Mathematica then returns the following solution:

    l_i\to 0,l_e\to 0,\rho _i\to -\frac{\gamma _S}{(-1+f) (1+k) \gamma _{\text{rr}}},c\to 1,l_b\to 0,\rho _f\to -\frac{\sigma  \gamma _S}{(-1+f) (1+k) \chi _0}

    Which I then checked by back-substituting:

    Code:
    In[240]:=
    (-Subscript[l, e])*(Subscript[\[Rho], f]/w) + 
       Subscript[\[Psi], l]*(Subscript[l, b]/w) /. myeqtpt
    Subscript[\[Rho], f]*Subscript[l, e] - Subscript[\[Psi], l]*
        Subscript[l, b] - Subscript[\[Chi], l]*Subscript[l, b] /. myeqtpt
    Subscript[\[Chi], l]*Subscript[l, b] - Subscript[\[Omega], l]*
        Subscript[l, i] /. myeqtpt
    Subscript[\[Gamma], rr]*Subscript[\[Rho], i] - Subscript[\[Chi], 0]*
        (Subscript[\[Rho], f]/\[Sigma]) - 
       (m*Subscript[l, e]*Subscript[\[Rho], f] - m*Subscript[\[Psi], l]*
          Subscript[l, b] + (m*Subscript[\[Chi], l]*Subscript[l, b]*
           Subscript[\[Rho], f])/(1 - Subscript[\[Rho], f]))*(1/\[Sigma]) /. 
      myeqtpt
    FullSimplify[Subscript[\[Gamma], S]/(k + c) + Subscript[\[Chi], 0]*f*
         (Subscript[\[Rho], f]/\[Sigma]) + 
        f*(1 + Subscript[\[Rho], f]/(1 - Subscript[\[Rho], f]))*
         ((m*Subscript[\[Chi], l]*Subscript[l, b])/\[Sigma]) - 
        Subscript[\[Gamma], rr]*Subscript[\[Rho], i] /. myeqtpt]
    \[Gamma]*(Subscript[\[Omega], l]*r*Subscript[l, i]) - \[Gamma]*(c - 1) /. myeqtpt
    
    Out[240]=
    0
    
    Out[241]=
    0
    
    Out[242]=
    0
    
    Out[243]=
    0
    
    Out[244]=
    0
    
    Out[245]=
    0
    Last edited by shawsend; December 9th 2009 at 12:13 PM.
    Follow Math Help Forum on Facebook and Google+

Similar Math Help Forum Discussions

  1. Steady State Numerical ODEs
    Posted in the Differential Equations Forum
    Replies: 1
    Last Post: September 10th 2010, 03:46 AM
  2. Help with steady-state unemployment
    Posted in the Advanced Applied Math Forum
    Replies: 6
    Last Post: March 31st 2010, 03:13 PM
  3. steady state
    Posted in the Differential Equations Forum
    Replies: 1
    Last Post: March 20th 2010, 12:38 AM
  4. coupled 1st order ODEs, and initial conditions...
    Posted in the Differential Equations Forum
    Replies: 1
    Last Post: January 27th 2010, 03:30 PM
  5. Steady State Vector
    Posted in the Discrete Math Forum
    Replies: 1
    Last Post: December 4th 2009, 05:47 PM

Search Tags


/mathhelpforum @mathhelpforum