Download the Maple worksheet to supplement the paper "Centers of the United States" (The College Mathematics Journal, 36 (November 2005) No. 5, 366-373)

An HTML version is included below.

Maple worksheet to supplement the paper

"The center of the United States and other applications of calculus to geography"

to appear in the College Mathematics Journal

I. Procedures and data

 > with(plots): with(plottools): with(CurveFitting):

Warning, the name changecoords has been redefined

Warning, the assigned name arrow now has a global binding

 > ############# # Constants # ############# milesperbox:=30.3:    #number of miles for each graph-paper unit NSslope:=76/55:     #the slope of a north-south line on the graph paper

 > ########################################################### # The procedure PFit                                      # #                                                         # # Procedure call: PFit([[x1,y1],...,[xm,ym]],n);          # # Input: 1.) any set of points [[x1,y1],...,[xm,ym]]      # #        2.) a positive integer n                         # # Output: 1.) A plot of the m points and and the graph of # #         a polynomial p(x) of degree n that fits the     # #         data.                                           # #         2.) The polynomial p(x)                         # ########################################################### PFit := proc (t,n)        local dataplot, curveplot, p, i, a, j, k, minx, maxx, miny, maxy;        p(x):=a[0];        for i from 0 to n do                 p(x):=p(x)+a[i]*x^i        end do:        f(x):=LeastSquares(t, x, curve=p(x));        dataplot:=plot(t,style=point, symbol=point, color=blue):        minx:=t[1,1];        maxx:=t[1,1];        miny:=t[1,2];        maxy:=t[1,2];        for j from 1 to nops(t) do                minx:=min(minx,t[j,1],minx);                maxx:=max(maxx,t[j,1],maxx);                miny:=min(minx,t[j,2],miny);                maxy:=max(maxy,t[j,2],maxy);        end do;        curveplot:=plot(f(x),x=minx..maxx,y=miny..maxy):        print(f(x));        display(dataplot,curveplot);  end proc:

 > ####################################################### # Points on the northern and souther border of the US # ####################################################### N:=[[0,61.5],[1,63],[2,66],[3,70],[5,73.5],[7,76],[10,79],[13,81.5],[16.5,86.5],[19,84],[20,85],[25,81],[30,77],[35,73],[40,70],[45,66],[50,63],[55,60.5],[60,56],[65,53],[67,49],[67.5,45],[68,39],[70,39],[71.5,39],[73,39.5],[75,40],[78,40.5],[82,41],[85,40.5],[88,42.5],[90,44.5],[92,43],[92.5,39]]: S:=[[0,61.5],[2,55],[2,51],[5,48],[7,45],[10,41.5],[13,39],[17,37],[18,33],[19,29],[22,28.5],[23.5,25],[23.5,21],[25.5,16.5],[28,19.5],[31,20.5],[35,21],[38,20],[41,18],[43,17],[45,18],[48,17],[51,15],[54,14],[55,11],[55,6],[56,2],[57,0],[58.5,3],[59,8],[59,11],[61,16],[64,17.5],[67,18.5],[71,20],[74,24.5],[77,27.5],[79,31],[83,32],[85.5,33],[86,36],[89,37.5],[92.5,39]]:

 > ################################################ # The locations and populations of the 75 most # # populated metropolitan areas                 # ################################################ population:=[[80,32,21199865],[2,56,16373645],[58,42,9157540],[75,30,7608070],[2,66,7039362],[77,32,6188463],[85,35,5819100],[65,40,5456428],[34,30,5221801],[33,22,4669571],[57,22,4112198],[59,3,3876380],[18,83,3554760],[11,47,3251876],[53,53,2968806],[67,37,2945831],[3,51,2813833],[50,37,2603607],[29,52,2581506],[55,8,2395997],[69,34,2358695],[13,78,2265223],[61,34,1979202],[5,66,1796857],[45,40,1776062],[59,45,1689572],[58,9,1644561],[59,37,1607486],[27,23,1592383],[74,25,1569541],[11,56,1563282],[64,35,1540157],[64,23,1499293],[43,18,1337726],[21,60,1333914],[67,25,1251509],[30,25,1249763],[55,29,1231311],[84,33,1188613],[70,23,1187941],[82,34,1183110],[73,39,1170111],[49,30,1135614],[59,5,1131184],[59,13,1100491],[75,39,1098201],[62,43,1088514],[36,36,1083346],[57,32,1025598],[72,27,996512],[61,23,962441],[41.5,35.5,950558],[5,62,922516],[52,23,921106],[81,36,875583],[11,44,843746],[40,36,803235],[76.5,37.5,732117],[45,47,716998],[21,44,712738],[59,26,687249],[18,37,679622],[4,59,661645],[77,33,637958],[74,33,629401],[77,35,624776],[64,39,618203],[42,20,602894],[68,36,594746],[83,35,591932],[55,6.5,589959],[44,31,583845],[24.5,18.5,569463],[4,65,563598],[64,17,549033]]:

II. Finding the bordering functions

 > USPoints:=plot([N,S],style=point, symbol=circle,symbolsize=3,color=black): NorthVector:=[10/sqrt(1+NSslope^2),10*NSslope/sqrt(1+NSslope^2)]: EWseg:=[10/sqrt(1+NSslope^2),10*NSslope/sqrt(1+NSslope^2)]: NorthArrow:=arrow([80,48],[80+10/sqrt(1+NSslope^2),48+10*NSslope/sqrt(1+NSslope^2)], 1,3,.2,harpoon): BigN:=textplot([80+18/sqrt(1+NSslope^2),48+18*NSslope/sqrt(1+NSslope^2),'N']): display(USPoints,NorthArrow,BigN);

 > ############################################################ # Define each of the components of the piecewise functions # ############################################################ PFit([[0,61.5],[1,63],[2,66],[3,70],[5,73.5],[7,76],[10,79],[13,81.5],[16.5,86.5]],3);

 > n1:=x->60.62419002+3.58920658224591093*x-.256191440600515907*x^2+0.808133348309122111e-2*x^3:

 > PFit([[16.5,86.5],[19,84],[20,85]],2);

 > n2:=x->282.1428572-21.2857142857145370*x+.571428571428577947*x^2:

 > PFit([[20,85],[25,81],[30,77],[35,73],[40,70],[45,66],[50,63],[55,60.5],[60,56]],3);

 > n3:=x->108.0952381-1.43143338143012455*x+0.165440115439259884e-1*x^2-0.117845117844413367e-3*x^3:

 > PFit([[60,56],[65,53],[67,49]],2);

 > n4:=x->-688+122/5*x-1/5*x^2:

 > PFit([[67,49],[67.5,45],[68,39],[70,39],[71.5,39],[73,39.5],[75,40]],5);

 > n5:=x->12160424.39-843086.510295552085*x+23374.1045072612578*x^2-323.926751781483630*x^3+2.24392578836979917*x^4-0.621599603007937270e-2*x^5:

 > PFit([[75,40],[78,40.5],[82,41],[85,40.5],[88,42.5],[90,44.5],[92,43],[92.5,39]],6);

 > n6:=x->-15184078.32+1106488.69106954709*x-33566.2742753724960*x^2+542.580817632036769*x^3-4.92891126535511947*x^4+0.238580683105850295e-1*x^5-0.480732439466748370e-4*x^6:

 > PFit([[0,61.5],[2,55],[2,51],[5,48],[7,45],[10,41.5],[13,39],[17,37],[18,33],[19,29]],3);

 > s1:=x->61.02554932-4.34123363851550702*x+.359216714805441926*x^2-0.113921969478251346e-1*x^3:

 > PFit([[19,29],[22,28.5],[23.5,25],[23.5,21],[25.5,16.5]],2);

 > s2:=x->-124.3016951+15.6042199201759360*x-.396026727655951704*x^2:

 > PFit([[25.5,16.5],[28,19.5],[31,20.5],[35,21],[38,20],[41,18],[43,17],[45,18],[48,17],[51,15],[54,14],[55,11]],4);

 > s3:=x->-406.5776352+43.0266023518660177*x-1.59551079328126422*x^2+0.258881095815548905e-1*x^3-0.156155892831574134e-3*x^4:

 > PFit([[55,11],[55,6],[56,2],[57,0],[58.5,3],[59,8],[59,11]],3);

 > s4:=x->-73771.19130+4023.93289727696038*x-73.0468504146236910*x^2+.441352447487981248*x^3:

 > PFit([[59,11],[61,16],[64,17.5],[67,18.5],[71,20]],4);

 > s5:=x->-67792.36764+4121.29764444619923*x-93.8578653188096439*x^2+.949187208987014584*x^3-0.359625202314874668e-2*x^4:

 > PFit([[71,20],[74,24.5],[77,27.5],[79,31],[83,32],[85.5,33]],2);

 > s6:=x->-391.2799118+9.85541683626756892*x-0.572418918722479436e-1*x^2:

 > PFit([[85.5,33],[86,36],[89,37.5],[92.5,39]],2);

 > s7:=x->-910.5146864+20.6100000000001238*x-.111846689895471155*x^2:

 > ############################################## # Define the northern and southern functions # ############################################## north(x):=piecewise(x<0,0,x>=0 and x<=16.5,n1(x),x>16.5 and x<=20,n2(x),x>20 and x<=60,n3(x),x>60 and x<=67,n4(x),x>67 and x<=75,n5(x),x>75 and x<=92.5,n6(x),x>92.5,39): south(x):=piecewise(x<=0,0,x>0 and x<=19,s1(x),x>19 and x<=25.5,s2(x),x>25.5 and x<=55,s3(x),x>55 and x<=59,s4(x),x>59 and x<=71,s5(x),x>71 and x<=85.5,s6(x),x>85.5 and x<=92.5,s7(x),x>92.5,39):

 > USCurve:=plot([north(x),south(x)],x=0..92.5,y=0..92.5,color=black): display(USPoints,USCurve,NorthArrow,BigN);

III. Area of the United States

 > Area:=evalf(Int(north(x)-south(x),x=0..92.5)): RealArea:=Area*milesperbox*milesperbox;

IV. Geographic center of the United States

 > xbar:=evalf(Int(x*(north(x)-south(x)),x=0..92.5)/Int(north(x)-south(x),x=0..92.5)); ybar:=evalf(Int(.5*(north(x)^2-south(x)^2),x=0..92.5)/Int(north(x)-south(x),x=0..92.5));

 > Center:=plot([[xbar,ybar]],style=point, symbol=circle, symbolsize=12,color=black): USCurve:=plot([north(x),south(x)],x=0.1..92.5,y=-2..92.5,color=black,axes=NONE): display(USCurve,Center,NorthArrow,BigN);

 >

V. The geographic median

 > ########################################### # Find the east-west line that splits the # # region in two equal subregions          # ########################################### EWline:=x->(-1/NSslope)*x+74.3: Wint:=solve(north(x)-EWline(x),x): Eint:=solve(south(x)-EWline(x),x): (1/Area)*(evalf(Int(north(x)-EWline(x),x=Wint..Eint))+evalf(Int(north(x)-south(x),x=Eint..92.5)));

 > ############################################# # Find the north-south line that splits the # # region in two equal subregions            # ############################################# NSline:=x->NSslope*x-8.3: Nint:=solve(north(x)-NSline(x),x): Sint:=solve(south(x)-NSline(x),x): (1/Area)*(evalf(Int(north(x)-NSline(x),x=Sint..Nint))+evalf(Int(north(x)-south(x),x=0..Sint)));

 > ######################### # The population median # ######################### xmed:=solve(NSline(x)-EWline(x),x); ymed:=NSline(solve(NSline(x)-EWline(x),x));

 > Lineplot:=plot([NSline(x),EWline(x)],x=0..92.5,color=black,axes=NONE): display(USCurve,Center,Lineplot,NorthArrow,BigN);

VI.  Population center (mean)

 > totalpopulation:=0: xmoment:=0: ymoment:=0: for i from 1 to nops(population) do        totalpopulation:=totalpopulation+population[i,3]:        xmoment:=xmoment+population[i,1]*population[i,3]:        ymoment:=ymoment+population[i,2]*population[i,3]:        j[i]:=[population[i,1],population[i,2]]: end do: PopulationCounted:=totalpopulation; PercentOfUSPop:=evalf(totalpopulation/(281421906-1211537-626932-3808610)); xpop:=evalf(xmoment/totalpopulation); ypop:=evalf(ymoment/totalpopulation); cities:=convert(j,'list'): cityplot:=plot(cities,style=point, symbol=cross,color=black, symbolsize=8,view=[0..90,0..90]): centercity:=plot([[evalf(xmoment/totalpopulation),evalf(ymoment/totalpopulation)]],style=point, symbol=box, symbolsize=12,color=black,axes=NONE): display(USCurve,cityplot,centercity,NorthArrow,BigN);

VII. Population median

 > ############################################# # Find the north-south line that splits the # # population in two                         # ############################################# NSPopline:=x->NSslope*x-48: westpop:=0: for i from 1 to nops(population) do        if (population[i,2]>NSPopline(population[i,1])) then                westpop:=westpop+population[i,3]        end if end do: evalf(westpop/totalpopulation);

 > ########################################### # Find the east-west line that splits the # # population in two                       # ########################################### EWPopline:=x->(-1/NSslope)*x+78.5: northpop:=0: for i from 1 to nops(population) do        if (population[i,2]>EWPopline(population[i,1])) then                northpop:=northpop+population[i,3]        end if end do: evalf(northpop/totalpopulation);

 > PopLineplot:=plot([NSPopline(x),EWPopline(x)],x=0..92.5,color=black): DashLineplot:=plot([EWline(x),NSline(x)],x=0..92.5,color=black,linestyle=DASH): display(USCurve,centercity,PopLineplot,NorthArrow,BigN);

 > ######################### # The population median # ######################### xpopmed:=solve(NSPopline(x)-EWPopline(x),x); ypopmed:=NSPopline(solve(NSPopline(x)-EWPopline(x),x));

VIII. The perimeter of the United States

 > TotalPerimeter:=milesperbox*(evalf(Int(sqrt(1+diff(north(x),x)^2),x=0..92.5))+evalf(Int(sqrt(1+diff(south(x),x)^2),x=0..92.5))); CanadianBorder:=milesperbox*evalf(Int(sqrt(1+diff(north(x),x)^2),x=16.5..92.5)); MexicanBorder:=milesperbox*evalf(Int(sqrt(1+diff(south(x),x)^2),x=2..25.5)); AtlanticCoast:=milesperbox*evalf(Int(sqrt(1+diff(south(x),x)^2),x=57..92.5)); GulfCoast:=milesperbox*evalf(Int(sqrt(1+diff(south(x),x)^2),x=25.5..57)); PacificCoast:=milesperbox*(evalf(Int(sqrt(1+diff(north(x),x)^2),x=0..16.5))+evalf(Int(sqrt(1+diff(south(x),x)^2),x=0..2)));

 >

 >