# Formulate and solve a system of linear, first-order ODEs in Mathematica

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.

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

## 1 Answer

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}]


Solve

  DSolve[eqs, S, t]


Here it is for n=3

Solve

 DSolve[eqs, S, t]


This is for n=6

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]


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


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


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


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


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


• 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. Feb 26 at 3:59
• 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. Feb 27 at 6:52
• @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. Feb 27 at 8:50
• @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 Sγ is and what is Sρ as these do not show in your code you showed. Feb 28 at 4:57
• Thanks User444 and Nasser for both forcing me to clarify and do some editing! Feb 28 at 5:10