1. ## 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?

2. 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.

3. 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?!

4. 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.

5. 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
$
X' = JF(X^*)\cdot\left(X - X^*\right)
$

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.

6. 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