Sample Tank Problem

Setup

Consider the following three tank setup, where the volume of water in each tank is 1 kL (kilolitre).
Fresh water flows into Tank 1 at a rate of k1 kL/min, and water is dumped out of Tank 2 at a rate of k1 kL/min.  Water circulates through the system with the rates as indicated (these are designed so the volume of water in each tank remains constant).   (This drawing was made by selecting “Graphics” in the menu, then choosing “New Graphic,” and then also selecting “Drawing Tools” from under the “Graphics” menu.

We may model this situation as a system of linear differential equations.
Let x1, x2, x3 be the masses of a solute (e.g. salt or chlorine) present in Tanks 1,2,3 (respectively) after t seconds.
Assume the initial solute masses  m1, m2, m3 in each tank.

Our system of linear differential equations may be written as
({{{, ({{x1' =, -k2 x1     + (k2-k1) x3}, {x2' = , k2 x1  - k1 x2 - (k2-k1)x3}, {x3' =,        (k2-k1)x2  - (k2-k1)x3,}})}})
and this can be rewritten as the matrix differential equation X’ = AX, X(0) = X0, where A is defined below.  The solution is given by X = e^(At) X0.

Minimal analysis.  (A good start)

Here are solutions for a few specific cases.

{k1, k2} = {.01, .02};      (*Flow rates*)
{m1, m2, m3} = {50, 20, 50};     (*Initial solutes*)

a = {{-k2, 0, k2-k1},{k2, -k1, -(k2-k1)},{0,k2-k1,-(k2-k1)}};
X0 = {{m1},{m2},{m3}};
MatrixExp[a t].X0 // MatrixForm

{k1, k2} = {.03, .02};      (*Flow rates*)
{m1, m2, m3} = {40, 30, 50};     (*Initial solutes*)

a = {{-k2, 0, k2-k1},{k2, -k1, -(k2-k1)},{0,k2-k1,-(k2-k1)}};
X0 = {{m1},{m2},{m3}};
Chop[FullSimplify[MatrixExp[a t].X0]] // MatrixForm   (*The solution,as a function of t.  Chop replaces numbers that are close to 0 by 0. *)

While this sort of analysis give correct numbers, it is difficult to get an understanding of the situation (for several reasons).  However, getting these specific solutions are often an important first step in understanding a more general problem.

A much better analysis.

Clear[a,k1,k2,m1,m2,m3];  (*Now, k1 is just k1, as opposed to .03.*)
a[k1_,k2_] := {{-k2, 0, k2-k1},{k2, -k1, -(k2-k1)},{0,k2-k1,-(k2-k1)}};
X0 = {{m1},{m2},{m3}};
MatrixExp[a[k1,k2] t].X0 // MatrixForm

Our first attempt at a general solution did not yield anything very helpful.  Let’s investigate the characteristic polynomial to see what our solutions will look like.

CharacteristicPolynomial[a[k1,k2],z]

Hmm... this is tough too.  What about looking at the roots for specific values within the range that would probably be relevant to the problem?

Manipulate[
Roots[CharacteristicPolynomial[a[k1,k2],z]==0, z], 
{k1, .001, .5}, {k2, .001, .5}]

Upon inspection, it appears we will have two complex eigenvalues with negative real part.  This means we will have solutions of the form e^(-** t) (c1 cos(** t) + c2 sin(** t)).  The third eigenvalue appears to be real; it is positive if k1 > k2, and it is negative if k2 > k1.  But, we implicitly assumed that k2 > k1, since k2-k1 is the flow rate from Tank 2 to 3, and from Tank 3 to 1.  Therefore, the third eigenvalue will be negative, giving us a solution of the form e^(-** t).  

Hence, we expect for our final solution to exhibit some sort of oscillation while converging towards the long-term solution of 0 solvent in each tank.

Here are graphs of the solution curves for various parameters with the specific m1,m2,m3=(50,20,30).  Notice that we should disregard the cases where k1>k2.

soln[k1_,k2_,t_] :=MatrixExp[ a[k1,k2] t].{{50}, {20}, {30}};

Manipulate[
Plot[ {soln[k1,k2,t][[1,1]],soln[k1,k2,t][[2,1]] ,soln[k1,k2,t][[3,1]]} , {t, 0, 500}, PlotRange -> {0, 50}, PlotLegends -> {"x1", "x2", "x3"}],
 {k1, .001, .05}, {k2, .001, .05}]
(* [[i,j]] is being used to obtain the (i,j) entry of the matrix.*)

Here’s a slightly nicer implementation, and we include the ability to modify initial quantities.

table[k1_,k2_,t_, m1_,m2_,m3_] :=Table[(MatrixExp[ a[k1,k2] t].Transpose[{{m1, m2, m3}}])[[i,1]], {i, 1, 3}];

Manipulate[
Plot[ Evaluate[table[k1,k2,t,m1,m2,m3] ], {t, 0, 500}, PlotRange -> {0, 50}, PlotLegends -> {"x1", "x2", "x3"}],
 {k1, .001, .05}, {k2, .001, .05}, {m1, 0, 50}, {m2, 0, 50}, {m3, 0 ,50}]