First, turn off the annoying spelling warnings and then load some packages that may be useful:

Off[General::spell]
Off[General::spell1]
<< Graphics`Arrow`
<< Graphics`PlotField`

The Problem: Social Mobility

Consider a town of 30,000 families which have been divided by the local chamber of commerce for planning purposes into three economic brackets (or classes): lower, middle, and upper. A builder wants to know the future populations of these brackets so they can decide what types (read: cost) of houses to build. The city also wants this information to help decide other issues like taxes and social services.

The past year the flow between these populations was the following:
20% of the lower move into the middle.
10% of the middle move back to the lower.
10% of the middle move to the upper.
15% of the upper move down to the middle.
5% of the lower move directly into the upper.
4% of the upper move down to the lower.

Last year there were 12,000 lower, 10,000 middle, and 8,000 upper income families.

The Model

We will model this situation with differential equations. We know the flow (rate of change over time) between the three brackets, so we can construct a model based on the first derivatives of the populations in each bracket.

Formulate the model

Let's use the following notation:
x: number of families in upper income bracket,
y: number of families in the lower income bracket,
z: number of families in the lower income bracket,
t: time measured in years.

A family is considered our basic unit, so there is no inclusion in our model of families which split apart or merge.

The inital number of families in each bracket has been provided:

x0   = 8000;
y0  =  10000;
z0  =  12000;

A Diagram will help us construct the differential equations which model the system. I am constructing this diagram in Mathematica, but you could draw it by hand or with another graphics program if you wish.

ptx  =  {0, 0};
pty  =  {-2, 2};
ptz  =  {2, 2};
cx  =  Graphics[Circle[ptx, 1]];
cy  =  Graphics[Circle[pty, 1]];
cz  =  Graphics[Circle[ptz, 1]];
tx  =  Graphics[Text["x", ptx]];
ty  =  Graphics[Text["y", pty]];
tz  =  Graphics[Text["z", ptz]];
axy  =  
        Graphics[
          Arrow[{-0.5, 0}, {-2.5, 
              1.5}]];
axz  =  
        Graphics[
          Arrow[{+0.5, 0}, {+2.5, 
              1.5}]];
azy  =   Graphics[
          Arrow[{+2, 2.5}, {-2, 2.5}]];
ayx  =  
        Graphics[
          Arrow[{-1.75, 1.75}, {-0.25, 
              0.25}]];
azx  =  
        Graphics[
          Arrow[{+1.75, 1.75}, {0.25, 
              0.25}]];
ayz  =   Graphics[
          Arrow[{-1.75, 2}, {+1.75, 
              2}]];
fxy  =  
        Graphics[
          Text["15%", {-1.5, 0.5}]];
fyx  =  Graphics[Text["10%", {-0.8, 1.1}]]; 
    fxz  =  Graphics[Text["4%", {1.5, 0.5}]]; 
    fzx  =  Graphics[Text["5%", {0.8, 1.0}]]; 
    fzy  =  Graphics[Text["20%", {0, 2.6}]]; 
    fyz  =  Graphics[Text["10%", {0, 2.1}]];
Show[cx, cy, cz, tx, ty, tz, axy, axz, azy, ayx, azx, ayz, fxy, fyx, 
      fxz, fzx, fzy, fyz, Frame -> True]

Show[cx, cy, cz, tx, ty, tz, axy, axz, azy, ayx, azx, ayz, fxy, fyx, 
      fxz, fzx, fzy, fyz]

From this diagram it is easy to construct the differential equations we want to solve. First, let's define all the flows from one population to another:

flowxy  =  15/100;
flowyx  =  10/100;
flowxz  =  4/100;
flowzx  =  5/100;
flowzy  =  20/100;
flowyz  =  10/100;

deqx  =  
      D[x[t], t]  ==  
-(flowxy  +  
                  flowxz)*x[t] + flowyx*y[t]   +  
          flowzx*z[t]
deqy  =  
      D[y[t], t]  ==  
flowxy*
            x[t]  - (flowyz + flowyx)*y[t]   +  
          flowzy*z[t]
deqz  =  
      D[z[t], t]  ==  
flowxz*x[t]  + 
          flowyz*y[t]   -  (flowzx + flowzy)*z[t]

Solve the Model

This is a constant coefficient linear system of differential equations. The solution depends on the eigenvalues of the coefficient matrix. Mathematica can get the solution directly:

Clear[x, y, z]

solution  =  
        DSolve[{deqx, deqy, deqz, x[0] == x0, y[0] == y0, 
            z[0] == z0}, {x[t], y[t], z[t]}, t];

Let's assign the solution values so we can examine the solution graphically:

x[t_]  =  x[t]  /.  Part[solution, 1, 1]
y[t_]  =  y[t]  /.  Part[solution, 1, 2]
z[t_]  =  z[t]  /.  Part[solution, 1, 3]

As a simple and obvious check, let's verify that the solution satisfies the initial conditions:

{x[0], y[0], z[0]}

It does! This also reminds us that the number of families in the high economic class starts at 8000, in the middle economic class 10000, and in the low economic class 12000. We need to know that to interpret the graph below:

A plot of the solution is:

Plot[{x[t], y[t], z[t]}, {t, 0, 20}, 
      PlotStyle -> {{RGBColor[1, 0, 0]}, {RGBColor[0, 1, 0]}, {RGBColor[
              0, 0, 1]}}]

KEY:
red: high economic class
green: middle economic class
blue: low economic class
time: years.

Answer the Question

According to our model, the populations in the three income brackets we have chosen will level out after about fifteen or twenty years. By this time the model is probably going to be invalid, however. We see that the conditions will result in the lower economic class being reduced. The quantity we are interested in is the stable limit of the population, which is

Limit[{x[t], y[t], z[t]}, t -> Infinity]

N[%]

Remember, we really have a discrete system here, but we are modeling it with a continuous model. Therefore, any time we report results we should convert back to the discrete case.

Round[%]

After a time (from our graph it looks like 15 to 20 years) our model predicts that the number of families in the income classes will be:
High income: 9137
Middle income: 13858
Low income: 7005

The ratios of the economic classes will be:

N[Limit[{1, y[t]/x[t], z[t]/x[t]}, t -> Infinity]]

Sensitivity Analysis on Initial Conditions

We want to know how sensitive our solution is to the initial conditions. In other words, is the solution we found stable or unstable? Will we approach this solution for other values of the initial conditions?

If the eigenvalues of the constant coefficient matrix  are real and negative, then we will have a stable solution. Let's check:

MatrixForm[
      coefficientmatrix = {{-(flowxy  +  flowxz), flowyx, 
            flowzx}, {flowxy, -(flowyz + flowyx), flowzy}, {flowxz, 
            flowyz, - (flowzx + flowzy)}}]

Eigenvalues[coefficientmatrix]

N[%]

The eigenvalue of zero represents a constant solution. This means the our solutions will stay in the same ratio to each other, but the actual values may vary. Let's look at a solution where we change the initial conditions.

\(\(\(Clear[x, y, z, deqx, deqy, deqz, flowxy, flowyx, flowxz, flowzx, 
      flowzy, flowyz]\)\(
\)
    \)\)
flowxy  =  15/100;
flowyx  =  10/100;
flowxz  =  4/100;
flowzx  =  5/100;
flowzy  =  20/100;
flowyz  =  10/100;

    
deqx  =  
        D[x[t], t]  ==  
-(flowxy  +  
                    flowxz)*x[t] + flowyx*y[t]   +  
            flowzx*z[t];
deqy  =  
        D[y[t], t]  ==  
flowxy*
              x[t]  - (flowyz + flowyx)*y[t]   +  
            flowzy*z[t];
deqz  =  
        D[z[t], t]  ==  
flowxz*x[t]  + 
            flowyz*y[t]   -  (flowzx + flowzy)*z[t];
\(\(\(solution\  = \ 
        DSolve[{deqx, deqy, deqz, x[0] == x0 + 1000, 
            y[0] == y0 - 800, z[0] == z0 + 200}, {x[t], y[t], 
            z[t]}, t];\)\(
\)
    \)\)
x[t_]  =  x[t]  /.  Part[solution, 1, 1];
y[t_]  =  y[t]  /.  Part[solution, 1, 2];
\(\(\(z[t_]\  = \ z[t]\  /. \ Part[solution, 1, 3];\)\(
\)
    \)\)
Plot[{x[t], y[t], z[t]}, {t, 0, 20}, 
      PlotStyle -> {{RGBColor[1, 0, 0]}, {RGBColor[0, 1, 0]}, {RGBColor[
              0, 0, 1]}}]
N[Limit[{x[t], y[t], z[t]}, 
        t -> Infinity]]
N[Limit[{1, y[t]/x[t], z[t]/x[t]}, t -> Infinity]]

From the above we see the final limiting values have changed but the ratio of the values remained the same. Our solution is stable.

Sensitivity Analysis on the flows

We can perform a sensitivity analysis on a great many of the variables in the model. Let's choose to see how the system is sensititve to the flow from the high economic class to the low economic class.  This means we leave the variable flowxz undetermined:

flowxz

Clear[flowxz, x, y, z, deqx, deqy, deqz]

dxdt  =  -(flowxy  +  flowxz)*x + flowyx*y   +  
        flowzx*z
dydt = 
      flowxy*x  - (flowyz + flowyx)*y   +  
        flowzy*z
dzdt = 
      flowxz*x  + flowyz*y   -  (flowzx + flowzy)*z

Stability is defined in terms of the partial derivatives.

stabilitymatrix  =  {{D[dxdt, x], D[dxdt, y], 
          D[dxdt, z]}, 
{D[dydt, x], D[dydt, y], 
          D[dydt, z]}, 
{D[dzdt, x], D[dzdt, y], 
          D[dzdt, z]}}

Eigenvalues[stabilitymatrix]

If the eigenvalues are distinct and negative real, then the solution is stable. This means that we must have the quantity under the square root being positive.

Plot[-26 - 
        50 flowxz + (Sqrt[11 - 600 flowxz + 2500 (flowxz^2)]), {flowxz, 0, 
        1}]

The solution is not robust for all values of flowzy! (Our eigenvalues will not be real and negative for all flowzy). But maybe we are too hasty.  Let's do a little more work:

Plot[
      Re[-26 - 
          50 flowxz + (Sqrt[11 - 600 flowxz + 2500 (flowxz^2)])], {flowxz, 0, 
        1}]

Since the Real part of this expression is always <0, we will have decay.

Plot[
      Im[-26 - 
          50 flowxz + (Sqrt[11 - 600 flowxz + 2500 (flowxz^2)])], {flowxz, 0, 
        1}]

Solve[11 - 600 flowxz + 2500 (flowxz^2) == 0, flowxz]

For flowxz between 1/50 and 11/50, the eigenvalues will be complex, and we will get an oscillatory solution.  However, since the real part is negative, it will relax to an equilibrium solution after a while. So although our solution behaves in a fundamentally different way, the result is still very similar to what we had before. A similar situation occurs for flowyz, and for we investigate the solution in greater detail for  that case.

flowyz:

Clear[flowyz, x, y, z, deqx, deqy, deqz]
flowxy  =  15/100;
flowyx  =  10/100;
flowxz  =  4/100;
flowzx  =  5/100;
flowzy  =  20/100;

    
dxdt  =  -(flowxy  +  flowxz)*x + flowyx*y   +  
          flowzx*z;
dydt = 
        flowxy*x  - (flowyz + flowyx)*y   +  flowzy*z;
dzdt = 
        flowxz*x  + 
          flowyz*y   -  (flowzx + flowzy)*z;

    
stabilitymatrix  =  {{D[dxdt, x], D[dxdt, y], 
            D[dxdt, z]}, 
{D[dydt, x], D[dydt, y], 
            D[dydt, z]}, 
{D[dzdt, x], D[dzdt, y], 
            D[dzdt, z]}};
Eigenvalues[stabilitymatrix]

Plot[-27 - 50 flowyz + 
        2 (Sqrt[-4 + 75 flowyz + 625 (flowyz^2)]), {flowyz, 0, 
        1}]

The solution is not robust for all values of flowzy! (Our eigenvalues will not be real and negative for all flowzy). But maybe we are too hasty.  Let's do a little more work:

Plot[
      Re[-27 - 50 flowyz + 
          2 (Sqrt[-4 + 75 flowyz + 625 (flowyz^2)])], {flowyz, 0, 
        1}]

Since the Real part of this expression is always <0, we will have decay.

Plot[
      Im[-27 - 50 flowyz + 
          2 (Sqrt[-4 + 75 flowyz + 625 (flowyz^2)])], {flowyz, 0, 
        1}]

Notice that there is a region where the eigenvalue is complex--what this means is that for flowyz in this region the solution will be an oscillatory function. However, since the real part is negative, it will relax to an equilibrium solution after a while. So although our solution behaves in a fundamentally different way, the result is still very similar to what we had before.

Solve[-4 + 75 flowyz + 625 (flowyz^2) == 0, 
      flowyz]

For flowyz < 1/25, the soution will oscillate.  Let's solve such a system and see what the solutions look like:

\(\(\(Clear[x, y, z, deqx, deqy, deqz, flowxy, flowyx, flowxz, flowzx, 
      flowzy, flowyz]\)\(
\)
    \)\)
flowxy  =  15/100;
flowyx  =  10/100;
flowxz  =  8/100;
flowzx  =  5/100;
flowzy  =  20/100;
\(\(\(flowyz\  = \ 2/100;\)\(
\)
    \)\)
coefficientmatrix = {{-(flowxy  +  flowxz), flowyx, 
            flowzx}, {flowxy, -(flowyz + flowyx), flowzy}, {flowxz, 
            flowyz, - (flowzx + 
                  flowzy)}};
N[Eigenvalues[coefficientmatrix]]

    
deqx  =  
        D[x[t], t]  ==  
-(flowxy  +  
                    flowxz)*x[t] + flowyx*y[t]   +  
            flowzx*z[t];
deqy  =  
        D[y[t], t]  ==  
flowxy*
              x[t]  - (flowyz + flowyx)*y[t]   +  
            flowzy*z[t];
deqz  =  
        D[z[t], t]  ==  
flowxz*x[t]  + 
            flowyz*y[t]   -  (flowzx + flowzy)*z[t];
\(\(\(solution\  = \ 
        DSolve[{deqx, deqy, deqz, x[0] == x0, y[0] == y0, 
            z[0] == z0}, {x[t], y[t], z[t]}, t];\)\(
\)
    \)\)
x[t_]  =  x[t]  /.  Part[solution, 1, 1];
y[t_]  =  y[t]  /.  Part[solution, 1, 2];
\(\(\(z[t_]\  = \ z[t]\  /. \ Part[solution, 1, 3];\)\(
\)
    \)\)
Plot[{x[t], y[t], z[t]}, {t, 0, 20}, 
      PlotStyle -> {{RGBColor[1, 0, 0]}, {RGBColor[0, 1, 0]}, {RGBColor[
              0, 0, 1]}}]
N[Limit[{x[t], y[t], z[t]}, 
        t -> Infinity]]
N[Limit[{1, y[t]/x[t], z[t]/x[t]}, t -> Infinity]]

The behaviour of the solutions has changed to an oscillation rather than exponential growth or decay to a constant value. But the oscillation is damped, and very slight.  We conclude
that although the nature of the solution has changed (exponential growth/decay to an equilibrium solution to oscillation about an equilibrium solution), the solutions are not radically different (approaching constant values), and so the model is robust.

flowxy:

Clear[flowxy, x, y, z, deqx, deqy, deqz]
flowyx  =  10/100;
flowxz  =  4/100;
flowzx  =  5/100;
flowzy  =  20/100;
flowyz  =  10/100;

    
dxdt  =  -(flowxy  +  flowxz)*x + flowyx*y   +  
          flowzx*z;
dydt = 
        flowxy*x  - (flowyz + flowyx)*y   +  flowzy*z;
dzdt = 
        flowxz*x  + 
          flowyz*y   -  (flowzx + flowzy)*z;

    
stabilitymatrix  =  {{D[dxdt, x], D[dxdt, y], 
            D[dxdt, z]}, 
{D[dydt, x], D[dydt, y], 
            D[dydt, z]}, 
{D[dzdt, x], D[dzdt, y], 
            D[dzdt, z]}};
Eigenvalues[stabilitymatrix]

Plot[-49 - 
        100 flowxy + (Sqrt[561 - 4200 flowxy + 10000 (flowxy^2)]), {flowxy, 
        0, 1}, PlotRange -> {-70, 0}]

The solution is robust for all values of flowxy!

flowyx:

Clear[flowyx, x, y, z, deqx, deqy, deqz]
flowxy  =  15/100;
flowxz  =  4/100;
flowzx  =  5/100;
flowzy  =  20/100;
flowyz  =  10/100;

    
dxdt  =  -(flowxy  +  flowxz)*x + flowyx*y   +  
          flowzx*z;
dydt = 
        flowxy*x  - (flowyz + flowyx)*y   +  flowzy*z;
dzdt = 
        flowxz*x  + 
          flowyz*y   -  (flowzx + flowzy)*z;

    
stabilitymatrix  =  {{D[dxdt, x], D[dxdt, y], 
            D[dxdt, z]}, 
{D[dydt, x], D[dydt, y], 
            D[dydt, z]}, 
{D[dzdt, x], D[dzdt, y], 
            D[dzdt, z]}};
Eigenvalues[stabilitymatrix]

Plot[-27 - 
        50 flowyx + (Sqrt[2]) (Sqrt[17 - 100 flowyx + 1250 (flowyx^2)]), {flowyx, 
        0, 1}]

The solution is robust for all values of flowyx!

flowzx:

Clear[flowzx, x, y, z, deqx, deqy, deqz]
flowxy  =  15/100;
flowyx  =  10/100;
flowxz  =  4/100;
flowzy  =  20/100;
\(\(\(flowyz\  = \ 10/100;\)\(
\)
    \)\)
dxdt  =  -(flowxy  +  flowxz)*x + flowyx*y   +  
          flowzx*z;
dydt = 
        flowxy*x  - (flowyz + flowyx)*y   +  flowzy*z;
\(\(\(dzdt = 
        flowxz*x\  + flowyz*y\ \  - \ \((flowzx + flowzy)\)*z;\)\(
\)
    \)\)
stabilitymatrix  =  {{D[dxdt, x], D[dxdt, y], 
            D[dxdt, z]}, 
{D[dydt, x], D[dydt, y], 
            D[dydt, z]}, 
{D[dzdt, x], D[dzdt, y], 
            D[dzdt, z]}};
Eigenvalues[stabilitymatrix]

Plot[-59 - 
        100 flowzx + (Sqrt[241 - 2200 flowzx + 10000 (flowzx^2)]), {flowzx, 
        0, 1}, PlotRange -> {-80, 0}]

The solution is robust for all values of flowzx!

flowzy:

Clear[flowzy, x, y, z, deqx, deqy, deqz]
flowxy  =  15/100;
flowyx  =  10/100;
flowxz  =  4/100;
flowzx  =  5/100;
flowyz  =  10/100;

    
dxdt  =  -(flowxy  +  flowxz)*x + flowyx*y   +  
          flowzx*z;
dydt = 
        flowxy*x  - (flowyz + flowyx)*y   +  flowzy*z;
dzdt = 
        flowxz*x  + 
          flowyz*y   -  (flowzx + flowzy)*z;

    
stabilitymatrix  =  {{D[dxdt, x], D[dxdt, y], 
            D[dxdt, z]}, 
{D[dydt, x], D[dydt, y], 
            D[dydt, z]}, 
{D[dzdt, x], D[dzdt, y], 
            D[dzdt, z]}};
Eigenvalues[stabilitymatrix]

Plot[-22 - 
        50 flowzy + (Sqrt[79 - 700 flowzy + 2500 (flowzy^2)]), {flowzy, 0, 
        1}, PlotRange -> {-80, 0}]

The solution is robust for all values of flowzy!