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)));

>

>