4
$\begingroup$

I know this code is wrong, but I am not sure why. I have a system of coupled differential equations. The model is given by the following equations:

$\frac{dS_0} {dt} = (1 - \mu_0) \rho_0 S_0 \frac{S_\gamma}{S_\rho} - \gamma_0 S_0$

$\frac{dS_1} {dt} = \mu_0 \rho_0 S_0 \frac{S_\gamma}{S_\rho} - \gamma_1 S_1 + (1 - \mu_1) \rho_1 S_1 \frac{S_\gamma}{S_\rho}$

$\vdots$

$\frac{dS_n}{dt} = \mu_{n-1} \rho_{n-1} S_{n-1} \frac{S_\gamma}{S_\rho} + (1-\mu_n) \rho_n S_n \frac{S_\gamma}{S_\rho} - \gamma_n S_n$.

Caveat: $\frac{S_\gamma}{S_\rho}$ is a constant term that scales cell birth and cell death rates over all life-stages so the population size remains constant over time. I thought it should be simple to solve the system with Mathematica because I have solved the simplest cases of the ODE system by hand (the 2-stage and 3-stage models) and have an idea that if I put it into matrix form correctly, it will have a triangular structure. Here's my code:

(*Change the recursion limit on Wolfram Mathematica on the Cloud*)
(* Define the number of stages *)
n = 2;
(* Define the state vector *)
S = Array[s, n];
(* Define the coefficient vectors *)
\[Mu] = Array[\[Mu], n];
\[Rho] = Array[\[Rho], n];
\[Gamma] = Array[\[Gamma], n];
(* Define the coefficient matrices *)
A = DiagonalMatrix[-\[Gamma]];
B = DiagonalMatrix[(1 - \[Mu]) \[Rho]];
(* Define the derivative of the state vector *)
dSdt = D[S[t], t];
(* Construct the derivative in terms of matrices A and B *)
ode = dSdt == A.S + B.S \[Gamma]/S \[Rho];
(* Define the initial condition vector *)
S0 = {1, 0}; (* First stage = 1, others = 0 *)
(* Solve the ODE system with the initial condition *)
solution = DSolve[{ode, S[0] == S0}, S[t], t]

If there are n=2 stages, there should be two coupled ODEs; if there are three stages, there should be three coupled ODEs. Is this an inefficient way to input everything to make a matrix ODE system? It seems the way I have formed the matrix ODE system here is either incorrect or missing something because when I run DSolve, it returns it unevaluated. The error message output is: "Function S appears with no arguments." I have searched here a few times to find the best way to formulate similar problems but haven't found anything that jumped out at me as a good way to correct my code.

$\endgroup$
3
  • $\begingroup$ you can not type \[Mu] = Array[\[Mu], n] because you will get infinite recursion !Mathematica graphics $\endgroup$
    – Nasser
    Feb 26 at 2:52
  • $\begingroup$ Ah OK! That explains the recursion wall! Thank you! $\endgroup$ Feb 26 at 2:56
  • $\begingroup$ I changed those 3 lines to the format of "Table[[Mu]i, {i, n}];" which could be more sensible... $\endgroup$ Feb 26 at 3:03

1 Answer 1

1
$\begingroup$

Starting all over again as correct equations posted.

Since Mathematica index starts from 1 and not zero, it is easier to start your equations from 1 and not 0 as you showed.

Let $\beta=\frac{S_{\gamma}}{S_{\rho}}$ a constant. Then

\begin{align*} S_{1}^{\prime} & =\left( 1-\mu_{1}\right) \rho_{1}S_{1}\beta-\gamma _{1}S_{1}\\ S_{2}^{\prime} & =\mu_{1}\rho_{1}S_{1}\beta-\gamma_{2}S_{2}+\left( 1-\mu _{2}\right) \rho_{2}S_{2}\beta\\ S_{3}^{\prime} & =\mu_{2}\rho_{2}S_{2}\beta-\gamma_{3}S_{3}+\left( 1-\mu _{3}\right) \rho_{3}S_{3}\beta\\ & \vdots\\ S_{n}^{\prime} & =\mu_{n-1}\rho_{n-1}S_{n-1}\beta-\gamma_{n}S_{n}+\left( 1-\mu_{n}\right) \rho_{n}S_{n}\beta \end{align*}

In Matrix form the above becomes

$$ \begin{pmatrix} S_{1}^{\prime}\\ S_{2}^{\prime}\\ S_{3}^{\prime}\\ \vdots\\ S_{n}^{\prime} \end{pmatrix} = \begin{pmatrix} -\gamma_{1}+\left( 1-\mu_{1}\right) \rho_{1}\beta & 0 & 0 & \cdots & \cdots & 0\\ \mu_{1}\rho_{1}\beta & -\gamma_{2}+\left( 1-\mu_{2}\right) \rho_{2}\beta & 0 & \cdots & \cdots & 0\\ 0 & \mu_{2}\rho_{2}\beta & -\gamma_{3}+\left( 1-\mu_{3}\right) \rho_{3}\beta & 0 & \cdots & 0\\ 0 & 0 & \cdots & \cdots & \cdots & 0\\ 0 & 0 & 0 & \cdots & \mu_{n-1}\rho_{n-1}\beta & -\gamma_{n}+\left( 1-\mu _{n}\right) \rho_{n}\beta \end{pmatrix} \begin{pmatrix} S_{1}\\ S_{2}\\ S_{3}\\ \vdots\\ S_{n}% \end{pmatrix} $$

To generate this in Mathematica, a simple loop will do.

n = 2;
S = Array[s, n]
S = Map[#[t] &, S]
dSdt = Map[D[#, t] &, S]
μArray = Array[μ, n]
ρArray = Array[ρ, n]
γArray = Array[γ, n]
eqs = First@Last@Reap@Do[z = (-γArray[[m]] + 
     (1 - μArray[[m]])*ρArray[[m]]*β)*S[[m]]; 
     If[m > 1, z = z + (μArray[[m - 1]]*ρArray[[m - 1]]*β)*S[[m - 1]]]; 
     Sow[dSdt[[m]] == z], {m, 1, n}]

Mathematica graphics

Solve

  DSolve[eqs, S, t]

Mathematica graphics

Here it is for n=3

Mathematica graphics

Solve

 DSolve[eqs, S, t]

Mathematica graphics

This is for n=6

Mathematica graphics

and so on.

old answer

Try this and see if this is what you meant

n = 2;
S = Array[s, n]
S = Map[#[t] &, S]

Mathematica graphics

μArray = Array[μ, n]
ρArray = Array[ρ, n]
γArray = Array[γ, n]

Mathematica graphics

A = DiagonalMatrix[-μArray]
B = DiagonalMatrix[(1 - μArray)  ρArray]
dSdt = Map[D[#, t] &, S]

Mathematica graphics

ode = Thread[dSdt == A . S + B . S  γArray/S  ρArray]

Mathematica graphics

S0 = {1, 0}
ic = MapThread[(#1 == #2) /. t -> 0 &, {S, S0}]

Mathematica graphics

solution = DSolve[{ode, ic}, S, t]

Mathematica graphics

$\endgroup$
6
  • $\begingroup$ This looks like it! And I can seen now that the "map" and "thread" parts are very useful here. I'll have to read up on them in the tech specs. Thank you, Nasser! I don't use this site often enough to upvote your answer, yet. But I will as I soon as my reputation goes up! Rifles through a pad of computations, looking for a particular computation... Yeah these look a lot like what I was expecting but I will compare to my original pen-and-paper work and see if I myself got something similar for the small number of cases that are tractable to calculate by hand. $\endgroup$ Feb 26 at 3:59
  • $\begingroup$ This is very close to what I was thinking of, I just can't figure out the correction so it's actually solving my model. What I've got here, is almost my complete model. In stage 2, here s[2][t], there should be a coefficient term from s[1][t] since they're coupled differential equations. So, I'm at the drawing board trying to figure out the corrections I need to make to the coefficient vectors. $\endgroup$ Feb 27 at 6:52
  • $\begingroup$ @BenHardisty if you look at your ode's. they are not coupled. s1 ode does depend on s2, and s2 ode does not depends on s1. I just used the same equations you had. May be if you can post the ode's in Latex to see what is the problem. if these odes are from a book, you can post screen shot. $\endgroup$
    – Nasser
    Feb 27 at 8:50
  • $\begingroup$ @BenHardisty I can't read your equations easily. Better posting them in your question so they show more clearly. I also do not know what is and what is as these do not show in your code you showed. $\endgroup$
    – Nasser
    Feb 28 at 4:57
  • $\begingroup$ Thanks User444 and Nasser for both forcing me to clarify and do some editing! $\endgroup$ Feb 28 at 5:10

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.