4
$\begingroup$

I'm writing a module on population dynamics for my elementary differential equations class's MMA lab. I am pulling in the US population data with no trouble:

cd = CountryData["US", {"Population", {1790 , 2000}}]

Then the plot goes smoothly:

rp = DateListPlot[cd, PlotStyle -> RGBColor[1, 0, 0]]

yielding:

![US Population Plot](http://imgur.com/gallery/aS1I1)

Then I solve the Malthusian model using the 1790 census as my initial condition:

malthus = DSolve[{p'[t] == k p[t], p[1790] == 3.93*10^6}, p[t], t]

and find the growth constant using the 2000 census:

Solve[(malthus[[1, 1, 2]] /. t -> 2000) == 281.42*10^6, k]

and plot it too:

malthusplot = Plot[malthus[[1, 1, 2]], {t, 1790, 2000}, PlotStyle -> RGBColor[0, 1, 0]]

Multhusian Plot

However, when I attempt to display these two graphs together the result is a correct display of whichever plot is listed first in the Show command, and the other plot is displayed as either a vertical or horizontal line. My best guess is that time series data does not play well with others. Any ideas?

$\endgroup$
2
$\begingroup$

The units of the time values used by DateListPlot and Plot are different. The former uses seconds, with the initial time fixed in on the first of January 1900,

cd["Times"] // Short

(* {-3471206400, -3155673600, -2997907200, -2840140800, 
    <<122>>, 3092601600, 3124137600, 3155673600} *)

It reinterprets them as years only for the ticks labels. The latter symbol uses years for the time values:

(times = First@ Cases[malthusplot, Line[list_] :> list[[All, 1]], Infinity]) // Short

(* {1790., 1790.06, <<73>>, 1999.93, 2000.} *)

A possible way is then to convert the time values of the output DateListPlot to years, and first call malthusplot in Show so as to use the unit years for both graphics:

rpBis = rp /. Line[list_] :> Line[MapAt[DateValue[#, "Year"] &, list, {All, 1}]];

Show[malthusplot, rpBis]

enter image description here


The choice of the initial time 1st of January 1900 for the seconds also explains why the vertical line in Show[rp, malthusplot] is located at 1900.

The values given by times (see definition above) are interpreted as seconds by DateListPlot and then converted to years for the ticks labels:

DateString /@ MinMax[DateObject /@ times]
(* {"Mon 1 Jan 1900 00:29:50", "Mon 1 Jan 1900 00:33:19"} *)
|improve this answer|||||
$\endgroup$
  • $\begingroup$ Excellent. Thanks, Xavier, for your detailed explanation. $\endgroup$ – ChrisB Nov 11 '16 at 1:39
  • $\begingroup$ @ChrisB Glad this is useful. The converse approach is to modify the other plot with malthusplotBis = malthusplot /. Line[list_] :> Line[MapAt[(# - 1900) 31536000 &, list, {All, 1}]], and then evaluate Show[rp, malthusplotBis]. $\endgroup$ – user31159 Nov 11 '16 at 1:45
  • 1
    $\begingroup$ @ChrisB and @xavier Don't forget that cd has a "Dates" property as well. :) $\endgroup$ – Edmund Nov 11 '16 at 1:57
0
$\begingroup$

The ordinates of the two plots are inconsistent in that cd is measured in seconds, and malthusplot is measured in years. Using the conversion factor, 3.154 10^7 gives,

dd = cd; dd[[2, 2]] = N[dd[[2, 2]]/(3.154 10^7) + 1900];
dp = ListLinePlot[dd, PlotStyle -> Red];
Show[malthusplot, dp]

enter image description here

|improve this answer|||||
$\endgroup$
  • $\begingroup$ Thanks, bbgodfrey. That's a nice efficient fix for the disagreement. $\endgroup$ – ChrisB Nov 11 '16 at 1:41
0
$\begingroup$

You may make use of Quantity and other functions in the Units guide.

Solve Malthus with quantities.

qMalthus = 
 DSolve[{p'[t] == k p[t], p[Quantity[1790, "Years"]] == Quantity[3.93*10^6, "People"]}, 
  p[t], t]

Mathematica graphics

Growth constant with quantities.

qSol = Solve[(p[t] /. qMalthus /. t -> Quantity[2000, "Years"]) == 
         Quantity[281.42*10^6, "People"], k]

Mathematica graphics

The function with quantities is:

p[t] /. First@qMalthus /. First@qSol

Mathematica graphics

When we input t as a "Years" Quantity then its unit will cancel with the per years ( 1 / "Years" ) and leave a "People" Quantity.

A data series compatible with DateListPlot can be obtained by a DateObject for the year and the "People" quantities. Use QuantityMagnitude to extract the year value from the t quantities.

(qRes = Table[{DateObject@QuantityMagnitude@{t}, 
               p[t] /. First@qMalthus /. First@qSol}, 
          {t, Quantity[1790, "Years"], Quantity[2000, "Years"]}]) // Short

Mathematica graphics

Then

qMalthusPlot = DateListPlot[qRes, PlotStyle -> RGBColor[0, 1, 0]]

Mathematica graphics

Which can be given directly to Show with rp as defined in the OP.

Show[rp, qMalthusPlot]

Mathematica graphics

Before all these fancy newfangled computrons including units in calculations was a good partial check that you had gotten the workings of your model correct.

Hope this helps.

|improve this answer|||||
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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