> #sfb -- 13 Nov 2001
> us1:=[[1790,3.9],[1800,5.3],[1810,7.2],[1820,9.6],[1830,12.9],[1840,17.1],[1850,23.1]];
> us2:=[[1850,23.1],[1860,31.4],[1870,38.6],[1880,50.2],[1890,62.9],[1900,76.0],[1910,92.0],[1920,105.7],[1930,122.8],[1940,131.7]];
> us3:=[[1790,3.9],[1800,5.3],[1810,7.2],[1820,9.6],[1830,12.9],[1840,17.1],[1850,23.1],[1860,31.4],[1870,38.6],[1880,50.2],[1890,62.9],[1900,76.0],[1910,92.0],[1920,105.7],[1930,122.8],[1940,131.7]];
> us4:=[[1790,3.9],[1800,5.3],[1810,7.2],[1820,9.6],[1830,12.9],[1840,17.1],[1850,23.1],[1860,31.4],[1870,38.6],[1880,50.2],[1890,62.9],[1900,76.0],[1910,92.0],[1920,105.7],[1930,122.8],[1940,131.7],[1950,150.7],[1960,179.3],[1970,203.3],[1980,226.5],[1980,226.5],[1990,248.7],[2000,281.4]];
> world:=[[1650,0.510],[1700,0.625],[1750,0.710],[1850,1.130],[1900,1.6],[1950,2.565],[1960,3.05],[1970,3.721],[1980,4.476],[1990,5.32]];#
> plot(us1);
> plot(us2);
> plot(us3);
> plot(us4);
> plot(world);
> dpoverp:=proc(f)
> local i, n;
> n:=nops(f);
> return [seq( [f[i][1],(f[i+1][2]-f[i-1][2])/(f[i+1][1]-f[i-1][1])/f[i][2]], i=2..n-1 )];
> end;
> datax:=proc(f)
> local i, n;
> n:=nops(f);
> return [seq(f[i][1],i=1..n)];
> end;
> datay:=proc(f)
> local i, n;
> n:=nops(f);
> return [seq(f[i][2],i=1..n)];
> end;
> datax(us1);datay(us1);
> #
> dpoverp(us1);dpoverp(us2);dpoverp(us3);dpoverp(us4);
> plot(dpoverp(us1));
> with(stats);
> fit[leastsquare[[x,y],y=b]]([datax(dpoverp(us1)),datay(dpoverp(us1))]);
> plot([dpoverp(us1),0.0298],x=1790..1850,y=0..0.035);
> minus1790:=x->x-1790;
> fit[leastsquare[[x,y],y=a*x+b]]([map(minus1790,datax(us1)),map(ln,datay(us1))]);
> exp(1.371974653);
> evalf(dsolve({diff(P(t),t)=0.298*P(t),P(0)=us1[1][2]}));
> plot([us1,3.9*exp(0.0298*(x-1790))],x=1790..1850);
> plot([us4,3.9*exp(0.0298*(x-1790))],x=1790..2000,y=0..250);
> dpoverpvsp:=proc(f)
> local i, n;
> n:=nops(f);
> return [seq( [f[i][2],(f[i+1][2]-f[i-1][2])/(f[i+1][1]-f[i-1][1])/f[i][2]], i=2..n-1 )];
> end;#
> plot(dpoverp(us2));
> plot(dpoverpvsp(us2));
> fit[leastsquare[[x,y],y=a*x+b]]([datax(dpoverpvsp(us2)),datay(dpoverpvsp(us2))]);
> plot([dpoverpvsp(us2),-.1571110018*10^(-3)*x+.3065615312*10^(-1)],x=0..125);
> eqn:=evalf(dsolve({diff(P(t),t)=P(t)*(0.0307-0.000157*P(t)),P(0)=us2[1][2]}));
> plot([us2,30700/(157+1172*exp(-0.0307*(x-1850)))],x=1850..1940);
> plot([us3,30700/(157+1172*exp(-0.0307*(x-1850)))],x=1790..2000,0..250);
> plot([us4,30700/(157+1172*exp(-0.0307*(x-1850)))],x=1790..2000,0..250);
> plot(dpoverp(us3));
> plot(dpoverpvsp(us3));
> plot(dpoverp(us4));
> plot(dpoverpvsp(us4));
> # From  Slicing Pizzas, Racing Turtles and Futher adventures in Applied Mathematics
> # time scale problems
> plot(dpoverp(world));
> plot(dpoverpvsp(world));
> povern:=proc(f)
> local i, n;
> n:=nops(f);
> return [seq( [f[i][1],1/f[i][2]], i=1..n )];
> end;
> #This doesn't work, ask why. :povern(world);
> plot(povern(world));
> minus1650:=x->x-1650;#
> fit[leastsquare[[x,y],y=a*x+b]]([map(minus1650,datax(world)),datay(povern(world))]);
> evalf(1/1.913);evalf((-0.00512/1.913));#
> N0:=0.525;alpha:=0.002675;line:=1/N0-alpha*(x-1650)/N0;
> plot([povern(world),line],x=1650..2000);
> dsolve(diff(P(t),t)=a*P(t)^2);
> plot([world,N0/(1-alpha*(x-1650))],x=1650..2000);
> solve(1-alpha*(x-1650)=0,x);
>