Day 23

Differential Equations

U.S. Population Growth problem (p. 245)

Definition of the initial value problem

Given the data of U.S. population from 1790 to 1990 (by 10), we explore various options for a model.

data
       = {3.9,  5.3,  7.24,  9.6,  12.86,  17,  23.2,  31.4,  38.5,  
        50.2,  63,  76.2,  92.2,  106,  123.2,  132.2,  151.3,  
        179.3,  203.3,  226.5,  248.7}years = Table[10  n + 1790, {n, 0, 20}]\(popsdata = Transpose[{years, data}]; 

    p1 = ListPlot[popsdata, PlotStyle -> {PointSize[.02]}, 
        DisplayFunction -> Identity]; \)

Let's start with linear regression

Needs["Statistics`LinearRegression`"]h[x_] = Fit[popsdata, {1, x}, x]\(p2 = Plot[h[x], {x, 1790, 2000}, DisplayFunction -> Identity]; 

    Show[p1, p2, DisplayFunction -> $DisplayFunction]; \)

\(residuals = 
      Table[h[1790 + 10 \((n - 1)\)] - data[\([n]\)], {n, 1, 21}]; 

    modeldiffs = Table[h[1790 + 10 \((n - 1)\)] - Mean[data], {n, 1, 21}]; 

    ListPlot[residuals, PlotStyle -> {PointSize[.02]}]; 

    SSRes = Length[residuals]*Mean[residuals^2]\)SSReg = Length[modeldiffs]*Mean[modeldiffs^2]SSReg/(SSRes + SSReg)

Looks pretty crummy, but we gather that we're missing a quadratic term in the model. Let's have a look at the regression statistics:

Regress[popsdata, {1, x}, x]

If you didn't pay any attention, and use some graphics, you might accidently accept that model! Remember to look at residuals, for any clues about what might be missing, or even how it might be improved.

Let's look for a model of the exponential form: for that reason, we can takes logs of the data and use linear regression:

Logpopsdata = Transpose[{years, Log[data]}]h[x_] = Fit[Logpopsdata, {1, x}, x]\(p2 = 
      Plot[Exp[h[x]], {x, 1790, 2000}, DisplayFunction -> Identity]; 

    Show[p1, p2, DisplayFunction -> $DisplayFunction]; \)

Regress[Logpopsdata, {1, x}, x]

\(residuals = 
      Table[Exp[h[1790 + 10 \((n - 1)\)]] - data[\([n]\)], {n, 1, 21}]; 

    modeldiffs = 
      Table[Exp[h[1790 + 10 \((n - 1)\)]] - Mean[data], {n, 1, 21}]; 

    ListPlot[residuals, PlotStyle -> {PointSize[.02]}]; 

    SSRes = Length[residuals]*Mean[residuals^2]\)SSReg = Length[modeldiffs]*Mean[modeldiffs^2]SSReg/(SSRes + SSReg)

An improvement, but still not good: There's something systematic in those residuals! The trouble is that the rate is not really constant across the last 200 years. Have a look:

\(ratedata = 
      Table[\((data[\([n + 1]\)] - data[\([n - 1]\)])\)/
          \((20*data[\([n]\)])\), {n, 2, 20}]; 

    rateyears = Table[10  n + 1790, {n, 1, 19}]; 

    ratedata = Transpose[{rateyears, ratedata}]\)p3 = ListPlot[ratedata, PlotStyle -> {PointSize[.02]}];

Hmmm...maybe a linear model would suffice. Let's model the rate, using simple linear regression:

r[x_] = Fit[ratedata, {1, x}, x]\(p4 = Plot[r[x], {x, 1790, 2000}, DisplayFunction -> Identity]; 

    Show[p3, p4, DisplayFunction -> $DisplayFunction]; \)

This changes our differential equation considerably: we now have

intr[x_] = Integrate[r[x], x]alpha = data[[1]]/Exp[intr[1790]]pop[x_] = alpha*Exp[intr[x]]\(p5 = Plot[pop[x], {x, 1790, 2000}, DisplayFunction -> Identity]; 

    Show[p1, p5, DisplayFunction -> $DisplayFunction]; \)

Check the residuals, and compute R^2:

\(residuals = 
      Table[pop[1790 + 10 \((n - 1)\)] - data[\([n]\)], {n, 1, 21}]; 

    modeldiffs = Table[pop[1790 + 10 \((n - 1)\)] - Mean[data], {n, 1, 21}]; 
    
ListPlot[residuals, PlotStyle -> {PointSize[.02]}]; 

    SSRes = Length[residuals]*Mean[residuals^2]\)SSReg = Length[modeldiffs]*Mean[modeldiffs^2]SSReg/(SSRes + SSReg)

That was one way to get the quadratic term model: models upon models (first model the rate, then use it to integrate for the population model). In this case, we didn't use any logs. Now let's get the same type of model, only by proceding directly to fit a linear model (including quadratic term) to the log(data):

h[x_] = Fit[Logpopsdata, {1, x, x^2}, x]\(p2 = 
      Plot[Exp[h[x]], {x, 1790, 2000}, DisplayFunction -> Identity]; 

    Show[p1, p2, DisplayFunction -> $DisplayFunction]; \)

Regress[Logpopsdata, {1, x, x^2}, x]

Check the residuals, and compute R^2:

\(residuals = 
      Table[Exp[h[1790 + 10 \((n - 1)\)]] - data[\([n]\)], {n, 1, 21}]; 

    modeldiffs = 
      Table[Exp[h[1790 + 10 \((n - 1)\)]] - Mean[data], {n, 1, 21}]; 

    ListPlot[residuals, PlotStyle -> {PointSize[.02]}]; 

    SSRes = Length[residuals]*Mean[residuals^2]\)SSReg = Length[modeldiffs]*Mean[modeldiffs^2]SSReg/(SSRes + SSReg)

This is just a hair better than the models-upon-models model!

Solving systems of first-order homogeneous linear equations with constant coefficients

Definition of the initial value problem

This is the model shown in Figure 3: Income Class Model. We are solving a problem of the form X'=MX with the initial value X(0) known (and since the equations conserve x+y+z, we will use an initial X(0) whose entries sum to 1). The sum should be 1 when we're done!

M = {{-.25, .1, .04}, {.2, -.2, .15}, {.05, .1, -.19}}

X0 = {1, 0, 0}

Eigenvalues and Eigenvectors

Eigenvalues[M]

To find numeric values of the eigenvalues use the N[] function.

N[%]

eigenlist = Eigensystem[M]

Define the eigenvalues and eigenvector as the appropriate components of eigenlist.  This can be done by entering them directly, but this "automated" approach makes it easy to reuse the code for other systems.

l1 = 
      eigenlist[[1
          ]][[1]]l2 = 
      eigenlist[[1
          ]][[2]]l3 = 
      eigenlist[[1
          ]][[3
        ]]

v1 = 
      eigenlist[[2
          ]][[1]]v2 = 
      eigenlist[[2
          ]][[2]]v3 = 
      eigenlist[[2
          ]][[3
        ]]

Define a vector of arbitrary constants.

c = {c1, c2, c3};

Notice that multiplying a list by a function distributes the function across the list.

(E^(-2 t)) v1

Define the matrix Φ.

Phi[t_] = 
      Transpose[{(E^(l1 t)) v1, (E^(l2 t)) v2, (E^(l3 t)) v3}]

The general solution

Phi[t].c

Solution of the initial value problem

cspec = Inverse[Phi[0]].X0

Phi[t].Simplify[cspec]

soln = Simplify[%]

The functions x, y, and z are defined to be the appropriate components of soln.

x[t_] = soln[[1]]y[t_] = soln[[2]]z[t_] = soln[[3]]

The constant terms above (those terms with E^0.t terms) give the asymptotic behavior, which we see below:

Plot the solution curves

Plot[{x[t], y[t], z[t]}, {t, 0, 100}];

ParametricPlot3D[{x[t], y[t], z[t]}, {t, 0, 20}, 
      ViewPoint -> {2.228, -1.932, 1.659}];

Now let's solve the same problem using Markov chains:

Needs["Statistics`DescriptiveStatistics`"];

MM = Transpose[{{.75, .2, .05}, {.1, .8, .1}, {.04, .15, .81}}]Eigenvalues[MM]v1 = Eigenvectors[MM][[1]](v1/((∑_(i = 1)^3) v1[[i
          ]]))