{VERSION 6 0 "IBM INTEL NT" "6.0" } {USTYLETAB {CSTYLE "Maple Input" -1 0 "Courier" 0 1 255 0 0 1 0 1 0 0 1 0 0 0 0 1 }{PSTYLE "Normal" -1 0 1 {CSTYLE "" -1 -1 "Times" 1 12 0 0 0 1 2 2 2 2 2 2 1 1 1 1 }1 1 0 0 0 0 1 0 1 0 2 2 0 1 }{PSTYLE "Maple Plot" -1 13 1 {CSTYLE "" -1 -1 "Times" 1 12 0 0 0 1 2 2 2 2 2 2 1 1 1 1 }3 1 0 0 0 0 1 0 1 0 2 2 0 1 }} {SECT 0 {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 52605 "# Investigate the minimum surface areas obtained\n# when the graph of a function f(x) i s revolved\n# about lines parallel to lines y = mx + beta.\n\n# filena me = Revslant.mws\n\nrestart:\n\n# Notes:\n# The run times are quite l ong when the full\n# minimizations are done. If you only want\n# to se e a few surfaces, you can bypass the\n# full minimization near the end of the\n# worksheet. You can also reduce the number\n# (Nvalm) of slo pes used in the first pass\n# calculations. (I typically use something \n# like Nvalm = 3.) As things stand, you can\n# use any one of 44 pro blems. If you wish\n# to define your own equation, search on\n# 'Probl em :=' and add your function at the\n# at the end of the block. If you have trouble\n# reproducing any of the plots in the manuscript,\n# I \+ have problem specific versions of this\n# worksheet in which I may hav e changed the\n# amount of detail to get satisfactory plots.\n# Please feel free to make any changes you like\n# and use the worksheet as yo u see fit. If you\n# find any blunders or better ways to do things,\n# I'll get a kick out of seeing your changes.\n\nSMG2b := proc(mminS,mm axS)\n# This is the function procedure for the Golden\n# Search minimi zation of Smin. SMG2b is called\n# by the main procedure.\n local to l2,tol2_dummy,zmax2,zmin2,Smin2,kminS2:\n global valx2,valf2,iters2, Use_Adaptmw2,\n From_Where:\n From_Where := 3:\n (tol2,tol2_dumm y) := Set_Tols(From_Where):\n lprint(\"SMG2b is calling GOLD2 ...\") ;\n # Use Adaptmw to do the integration.\n (valf2,valx2,iters2) := GOLD2(SMG3b,mminS,mmaxS,tol2):\n lprint(\"Back in SMG2b after call \+ to GOLD2 ...\");\n lprint(\"Slope m yielding minimum Smin = \", valx 2);\n lprint(\"Minimum surface area S(valx2) for this m = \",\n \+ abs(valf2));\n # lprint(\"Number of Golden 2 Search iterations = \",\n # iters2);\n return(valf2,valx2):\nend proc:\n\nSM G3b := proc(mval)\n# This is the function procedure for the Golden 2\n # Search minimization based on the pivoting\n# information. SMG3b is c alled by GOLD2.\n local tol,theta,cost,sint,atan,wofx,zofx,\n wpof x,zpofx,Pie,betaSc,SminSc,k,m,\n From_Where,abserrA,relerrA,SvalA,er restA,\n flagA,nfeA:\n global kminSc,kminminS,fA2b:\n m := evalf (mval):\n (theta,cost,sint,atan,wofx,zofx,wpofx,zpofx,Pie):=tcs(m): \n betaSc := 'betaSc': kminSc := 'kminSc':\n SminSc = 'SminSc':\n \+ betaSc := evalf(fcpivot - cpivot*m):\n kminSc := evalf((betaSc-bet a)/sqrt(m^2+1)):\n k := kminSc:\n # Save kmin so the main procedur e can get to it.\n kminminS := kminSc:\n # lprint(\"In SMG3b. m, b etaSc, kminSc = \",\n # m, betaSc, kminSc);\n # Define the surf ace area integrand for Adaptmw.\n fA2b := 'fA2b':\n fA2b := x -> 2 *Pie*evalf(\n abs(zofx(x)-k)*sqrt((wpofx(x))^2+(zpofx(x))^2)):\n \+ # Define the error tolerances.\n From_Where := 6:\n (abserrA,rel errA) := Set_Tols(From_Where):\n flagA := -1:\n # Approximate the \+ integral of the function fA2b(x).\n (SvalA,errestA,flagA,nfeA) :=\n \+ Adaptmw(fA2b,a,b,abserrA,relerrA):\n if (flagA <> -1) then\n \+ # Print the results.\n # lprint(\"Status flagA = \", \+ flagA);\n # lprint(\"Approximate integral = \", SvalA);\n # lprint(\"Approximate error = \", errestA);\n # lprint(\"Inte grand evaluations = \", nfeA);\n # Check if Adaptmw was successfu l.\n if (flagA = 1) then\n error \"Insufficient queue sto rage in Adaptmw\";\n fi;\n if (flagA = 2) then\n err or \"Too many fA2b evals in Adaptmw\";\n fi;\n fi;\n SminSc : = -SvalA:\n lprint(\"In Smg3b. m, betaSc, kminSc, SminSc, beta = \", \n m, betaSc, kminSc, SminSc, beta):\n return(SminSc):\nend pro c:\n\nreflectf := proc(a,b,f,m,betarfl)\n# Reflect the graph of the fu nction f(x) through\n# the line y = mx + beta on the interval [a,b].\n local theta,cost,sint,atan,wofx,zofx,wpofx,\n zpofx,Pie,p1,p2,p3, p4,p5,invw,invz,midxa,midxb,\n EndPerpa,EndPerpb,l,fp,x:\n # Retur ns:\n # p1 = plot of line l(x) = mx + betarfl (black)\n # p2 = plot of f(x) (blue)\n # p3 = plot of reflection of f(x) thru l(x) (red)\n # p4 = perpendicular line thru (a,l(a)) (green)\n # p 5 = perpendicular line thru (b,l(b)) (green) \n # printlevel := 0:\n x := 'x':\n p1 := 'p1': p2 := 'p2': p3 := 'p3': p4 := 'p4':\n p 5 := 'p5': EndPerpa := 'EndPerpa':\n EndPerpb := 'EndPerpb': l := 'l ':\n invw := 'invw': invz := 'invz':\n midxa := 'midxa': midxb := \+ 'midxb':\n l := x -> m*x + betarfl:\n fp := x -> D(f)(x):\n (the ta,cost,sint,atan,wofx,zofx,wpofx,zpofx,Pie)\n := tcs2(m,betarfl):\n # Reflect zofx(x) about the w-axis: z -> -z.\n invw := x -> cost* wofx(x) - sint*(-zofx(x)):\n invz := x -> sint*wofx(x) + cost*(-zofx (x)) + betarfl:\n # midxa := evalf((a+invw(a))/2):\n # midxb := ev alf((b+invw(b))/2):\n # p1 := plot([x,l(x),x=midxa..midxb],\n p1 : = plot([x,l(x),x=a..b],\n color=black,scaling=constrained):\n \+ p2 := plot([x,f(x),x=a..b],\n color=blue,scaling=constrained ):\n p3 := plot([invw(x),invz(x),x=a..b],\n color=red,scalin g=constrained,\n title=`Reflection of f(x) thru y=mx+beta`):\n if (abs(invw(a)-a) > 1e-5 and abs(invw(b)-b) > 1e-5) then\n # \+ Plot the endpoint perpendiculars if they\n # aren't vertical.\n \+ EndPerpa := x -> f(a)+(x-a)*(invz(a)-f(a))/(invw(a)-a):\n p4 \+ := plot([x,EndPerpa(x),x=a..invw(a)],\n color=green,scaling =constrained):\n EndPerpb := x -> f(b)+(x-b)*(invz(b)-f(b))/(invw (b)-b):\n p5 := plot([x,EndPerpb(x),x=b..invw(b)],\n c olor=green,scaling=constrained):\n else\n p4 := p1: p5 := p1:\n fi;\n return(p1,p2,p3,p4,p5):\nend proc:\n\nSMG := proc(zmin,zmax )\n# This is the function procedure for the Golden\n# Search minimizat ion of kmin. SMG is called by\n# the main procedure and proc SK.\n l ocal tol,tol_dummy,From_Where:\n global valx,valf,iters,Use_Adaptmw: \n lprint(\"In SMG. zmin, zmax = \", zmin, zmax);\n From_Where := \+ 1:\n (tol,tol_dummy) := Set_Tols(From_Where):\n lprint(\"SMG is ca lling GOLD ...\");\n Use_Adaptmw := 1:\n if (Use_Adaptmw = 1) then \n # Use Adaptmw to do the integrations.\n (valf,valx,iters) := GOLD(SA,zmin,zmax,tol):\n else\n error \"This feature is no t active.\";\n # Use Maple to do the integrations.\n # (valf ,valx,iters) := GOLD(SM,zmin,zmax,tol):\n fi:\n lprint(\"valx = \" , valx):\n lprint(\"Minimum surface area S(valx) = \",\n abs(va lf)):\n lprint(\"Number of Golden Search iterations = \",\n ite rs):\n return(valf,valx):\nend proc:\n\nSA := proc(k)\n# Calculate t he surface area and multiply by -1.\n# This is the function procedure \+ for the Golden\n# Search minimization if Adaptmw is used to do\n# the \+ integrations. SA is called by SMG. m is a\n# global variable defined i n the main procedure.\n local flagA,abserrA,relerrA,answer,errestA,n feA,\n x,Sval,SvalA,theta,cost,sint,atan,wofx,zofx,\n wpofx,zpofx, Pie,From_Where:\n global fA:\n # lprint(\"In SA. k = \", k);\n ( theta,cost,sint,atan,wofx,zofx,wpofx,zpofx,Pie):=tcs(m):\n # Define \+ the surface area integrand for Adaptmw.\n fA := 'fA':\n fA := x -> 2*Pie*evalf(\n abs(zofx(x)-k)*sqrt((wpofx(x))^2+(zpofx(x))^2)): \n # Define the error tolerances.\n From_Where := 5:\n (abserrA, relerrA) := Set_Tols(From_Where):\n flagA := -1:\n # Approximate t he integral of the function fA(x).\n (SvalA,errestA,flagA,nfeA) :=\n Adaptmw(fA,a,b,abserrA,relerrA):\n if (flagA <> -1) then\n \+ # Print the results.\n # lprint(\"Status flagA = \", f lagA);\n # lprint(\"Approximate integral = \", SvalA);\n # \+ lprint(\"Approximate error = \", errestA);\n # lprint(\"Integ rand evaluations = \", nfeA);\n # Check if Adaptmw was successful .\n if (flagA = 1) then\n error \"Insufficient queue stor age in Adaptmw\";\n fi;\n if (flagA = 2) then\n erro r \"Too many fA evals in Adaptmw\";\n fi;\n fi;\n Sval := -Sv alA:\n return(Sval):\nend proc:\n\nSK := proc(m)\n# Determine kminS \+ for a given slope m. SK is\n# called only by the main procedure.\n l ocal Kval,atan,theta,cost,sint,Pie,\n Smin,kminS,wofx,zofx,wpofx,zpo fx:\n (theta,cost,sint,atan,wofx,zofx,wpofx,zpofx,Pie):=tcs(m):\n \+ # Do the Golden Section minimization for\n # this value of m.\n # \+ Note:\n # SK assumes that zmin and zmaz were calculated\n # before callin SK.\n lprint(\"In SK. m, zmin, zmax = \", m, zmin, zmax);\n \+ Smin := 'Smin':\n # (Smin,kminS) := -SMG(zmin,zmax):\n # kminS : = valx:\n (Smin,kminS) := SMG(zmin,zmax):\n Smin := -Smin:\n lpr int(\"In SK. m, kminS, S(kminS) = \", m, kminS, Smin);\n Kval := kmi nS:\n return(Kval):\nend proc:\n\nSMG2 := proc(mminS,mmaxS)\n# This \+ is the function procedure for the Golden\n# Search minimization of Smi n. SMG2 is called\n# by the main procedure.\n local tol2,tol2_dummy, zmax2,zmin2,Smin2,kminS2:\n global valx2,valf2,iters2,Use_Adaptmw2, \n From_Where:\n From_Where := 2:\n (tol2,tol2_dummy) := Set_Tol s(From_Where):\n lprint(\"SMG2 is calling GOLD2 ...\");\n Use_Adap tmw2 := 1:\n if (Use_Adaptmw2 = 1) then\n # Use Adaptmw to do t he integration.\n (valf2,valx2,iters2) :=\n GOLD2(SMG3,mm inS,mmaxS,tol2):\n else\n error \"In SMG2. This feature is not \+ active.\";\n # Use Maple to do the integrations.\n # (valf2, valx2,iters2) :=\n # GOLD2(SM2,mminS,mmaxS,tol2):\n fi:\n \+ lprint(\"Back in SMG2 after call to GOLD2 ...\");\n lprint(\"Slope m yielding minimum Smin = \", valx2);\n lprint(\"Minimum surface area S(valx2) for this m = \",\n abs(valf2));\n lprint(\"Numbe r of Golden 2 Search iterations = \",\n iters2);\n return( valf2,valx2):\nend proc:\n\nSMG3 := proc(m)\n# This is the function pr ocedure for the Golden 2\n# Search minimization. SMG3 is called by GOL D2.\n local tol,theta,cost,sint,atan,wofx,zofx,\n wpofx,zpofx,Pie, zmin,zmax,Nzm,zminm,zmaxm:\n global valx3,valf3,iters3,Use_Adaptmw,m g2,\n kminminS:\n (theta,cost,sint,atan,wofx,zofx,wpofx,zpofx,Pie) :=tcs(m):\n mg2 := 'mg2':\n mg2 := evalf(m):\n # Estimate zmin, \+ zmax for this value of m.\n Nzm := 'Nzm':\n Nzm := 21:\n (zminm, zmaxm) := Estimate_Range(Nzm,zofx,a,b):\n tol := 'tol':\n tol := e valf(1/10^4):\n Use_Adaptmw := 1:\n lprint(\"SMG3 is calling GOLD \+ ...\");\n if (Use_Adaptmw = 1) then\n # Use Adaptmw to do the i ntegration.\n (valf3,valx3,iters3) :=\n GOLD(SA2,zminm,zm axm,tol):\n else\n # Use Maple to do the integration.\n er ror \"This feature is not active.\";\n # (valf3,valx3,iters3) := \n # GOLD(SM,zminm,zmaxm,tol):\n fi:\n # Save kmin so the \+ main procedure can get to it.\n kminminS := valx3:\n lprint(\"Back in SMG3 from GOLD ...\");\n lprint(\"Slope m and kmin for this slop e = \", m, valx3);\n lprint(\"Minimum surface area S(valx3) for this m = \",\n abs(valf3));\n lprint(\"Number of Golden Search iter ations for this m = \",\n iters3);\n return(valf3):\nend proc: \n\nSA2 := proc(k)\n# Calculate the surface area and multiply by -1.\n # This is the function procedure for the Golden\n# Search minimization if Adaptmw is used to do\n# the integrations. SA2 is called by GOLD. \+ The\n# global variable mg2 is the current slope\n# (passed by GOLD2). \+ It is defined in proc SMG3. \n local flagA,abserrA,relerrA,answer,er restA,nfeA,\n x,Sval,SvalA,theta,cost,sint,atan,wofx,zofx,\n wpofx ,zpofx,Pie,From_Where:\n global fA2:\n (theta,cost,sint,atan,wofx, zofx,wpofx,zpofx,Pie):=tcs(mg2):\n # Define the surface area integra nd for Adaptmw.\n fA2 := 'fA2':\n fA2 := x -> 2*Pie*evalf(\n \+ abs(zofx(x)-k)*sqrt((wpofx(x))^2+(zpofx(x))^2)):\n # Define the erro r tolerances.\n From_Where := 6:\n (abserrA,relerrA) := Set_Tols(F rom_Where):\n # abserrA := 1e-5:\n # relerrA := 1e-5:\n flagA := -1:\n # Approximate the integral of the function fA2(x).\n (SvalA ,errestA,flagA,nfeA) :=\n Adaptmw(fA2,a,b,abserrA,relerrA):\n i f (flagA <> -1) then\n # Print the results.\n # lprint(\"Sta tus flagA = \", flagA);\n # lprint(\"Approximate integra l = \", SvalA);\n # lprint(\"Approximate error = \", errestA );\n # lprint(\"Integrand evaluations = \", nfeA);\n # Check if Adaptmw was successful.\n if (flagA = 1) then\n error \"Insufficient queue storage in Adaptmw\";\n fi;\n if (flag A = 2) then\n error \"Too many fA evals in Adaptmw\";\n f i;\n fi;\n Sval := -SvalA:\n return(Sval):\nend proc:\n\nSEval2 \+ := proc(m,beta)\n# SEval2 is used to evaluate the surface\n# area for \+ the 3d plot at the end of the\n# main procedure. SEval2 calculates sur face\n# areas using w-z coordinates. It is called\n# only from the end of the main procedure.\n local atan,theta,cost,sint,wofx,zofx,\n \+ wpofx,zpofx,Pie,flagA,abserrA,relerrA,\n answer,errestA,nfeA,Sval,Sv alA,x,\n From_Where:\n global fA3:\n atan := tan@@(-1):\n the ta := evalf(atan(m)):\n cost := evalf(cos(theta)):\n sint := eva lf(sin(theta)):\n Pie := evalf(Pi):\n wofx := x -> cost*x+(f( x)-beta)*sint:\n zofx := x -> -sint*x+(f(x)-beta)*cost:\n wpofx \+ := x -> cost+fp(x)*sint:\n zpofx := x -> -sint+fp(x)*cost:\n # \+ lprint(\"Calling tcs with m = \", m):\n # (theta,cost,sint,atan,wofx ,zofx,wpofx,zpofx,Pie):=tcs(m):\n # Define the integrand for Adaptmw .\n # lprint(\"(,m,beta,a,b) = \", m, beta, a, b):\n fA3 := 'fA3': \n fA3 := x -> 2*Pie*evalf(\n abs(zofx(x))*sqrt((wpofx(x))^2+(z pofx(x))^2)):\n fA3(x):\n # Define the error tolerances.\n From_ Where := 4:\n (abserrA,relerrA) := Set_Tols(From_Where):\n flagA : = -1:\n # Approximate the integral of the function fA3(x).\n # lpr int(\"Calling Adaptmw ...\"):\n (SvalA,errestA,flagA,nfeA) :=\n \+ Adaptmw(fA3,a,b,abserrA,relerrA):\n # lprint(\"Back from Adaptmw. f lagA = \", flagA):\n if (flagA <> -1) then\n # Print the result s.\n # lprint(\"Status flagA = \", flagA):\n # lpri nt(\"Approximate integral = \", SvalA):\n # lprint(\"Approximate error = \", errestA):\n # lprint(\"Integrand evaluations = \+ \", nfeA):\n # Check if Adaptmw was successful.\n if (flagA \+ = 1) then\n error \"Insufficient queue storage in Adaptmw\";\n fi:\n if (flagA = 2) then\n error \"Too many fA eva ls in Adaptmw\";\n fi: \n fi:\n Sval := SvalA:\n # From_Whe re := -2:\n # How_Many_Digits := Set_Tols_Int(From_Where):\n # Sva l := 2*Pie*evalf(Re(int(\n # abs(zofx(x))*sqrt((wpofx(x))^2+(zpof x(x))^2),\n # x=a..b,digits=How_Many_Digits))):\n # lprint(\"In SEval2. m, beta, Sval = \",\n # m, beta, Sval);\n return(Sval) :\nend proc:\n\n\nEstimate_Range := proc(n,g,a,b)\n# Estimate the rang e of a function g(x)\n# on an interval [a,b] using n equally\n# spaced points.\n local deltax,xval,i,gmin,gmax,\n gptx,gpta,gptb,xpt:\n \+ deltax := 'deltax': xpt := 'xpt':\n gmin := 'gmin': gmax := 'gmax' :\n gmin := 1e10:\n gmax := -1e10:\n if (n > 1) then\n delt ax := evalf((b-a)/(n-1)):\n i := 'i': \n for i from 1 to n d o\n xpt := 'xpt': gptx := 'gptx':\n xpt := evalf(a+(i- 1)*deltax):\n gptx := evalf(g(xpt)):\n gmin := min(gpt x,gmin):\n gmax := max(gptx,gmax):\n od:\n i := 'i': \n elif (n = 1) then\n gpta := 'gpta': gptb := 'gptb': \n \+ gpta := evalf(g(a)):\n gptb := evalf(g(b)):\n gmin := min(gp ta,gptb):\n gmax := max(gpta,gptb): \n else\n error \"Esti mate_Range requires more than points.\";\n fi:\n gmin := gmin - de ltax:\n gmax := gmax + deltax:\n return(gmin,gmax):\nend proc:\n\n Set_Tols := proc(From_Where)\n# Set Adaptmw integration error toleranc es and\n# Golden search error tolerances.\n local abstol,reltol:\n \+ if (From_Where < 1 or From_Where > 11) then\n error \"Illegal to lerance flag in Set_Tols\";\n fi;\n if (From_Where = 1) then\n \+ # Inner Golden Search for Smin:\n abstol := 1e-4:\n reltol := 1e-4: # not used\n elif (From_Where = 2) then\n # Outer Gol den Search to minimize Smin:\n abstol := 1e-2:\n reltol := 1 e-2: # not used\n elif (From_Where = 3) then\n # Surface area 3 d plot:\n abstol := 1e-5:\n reltol := 1e-5:\n else\n \+ # Other Adaptmw integrations:\n abstol := 1e-5:\n reltol := \+ 1e-5:\n fi:\n # lprint(\"In Set_Tols. From_Where, abstol, reltol = \",\n # From_Where, abstol, reltol);\n return(abstol,reltol): \nend proc:\n\nSforpion2 := proc(a,b)\n# Find kmin and Smin for theta \+ = Pi/2.\n local kminforpion2,Sminforpion2,k,From_Where,\n abserrA, relerrA,bbal,ccal,abserral,relerral,\n ff,abserrZ,relerrZ,flag,bb,cb ,residual,nfeal,\n gg,Aal3,Bal3,iflaga3,Queue_Size,answeral3,\n er restal3,iflagal3,Nfeal3,fval,gal,twopie,\n Sforkis0:\n lprint(\"In side proc SFor_pion2\");\n From_Where := 10:\n (abserrZ,relerrZ) : = Set_Tols(From_Where):\n bbal := a:\n ccal := b:\n flag := -1: \n (bb,cb,residual,flag,nfeal) :=\n Zeromw(SfAL,bbal,ccal,abser rZ,relerrZ):\n if (flag <> -1) then\n # lprint(\"Results follow ...\");\n if flag = 0 then\n lprint(\"flag = 0. A zero w as successfully found.\");\n elif flag = 1 then\n lprint( \"flag = 1. Too many function evals in Adaptmw\");\n elif flag = \+ 2 then\n lprint(\"flag = 2. A pole may have been located.\"); \n fi;\n # lprint(\" bbal = \", bbal);\n # lprint(\" cc al = \", ccal);\n # lprint(\" Number of function evaluations = \" , nfeal);\n cb := bb:\n else\n error \" A fatal error occu rred in Zeromw.\";\n fi;\n # A horizontal line z = -x = k correspo nds to a\n # vertical line x = -k.\n kminforpion2 := -cb:\n gg : = x -> abs(-x-kminforpion2)*sqrt(1+(fp(x))^2):\n Aal3 := a:\n Bal3 := b:\n iflaga3 := -1:\n Queue_Size := 0:\n (answeral3,errestal 3,iflagal3,Nfeal3) :=\n Adaptmw(gg,Aal3,Bal3,Abserral,Relerral,Qu eue_Size):\n if (iflagal3 = 1) then\n error \"Insufficient queu e storage in Adaptmw\";\n fi;\n if (iflagal3 = 2) then\n erro r \"Too many gg function evals in Adaptmw\";\n fi;\n if (iflagal3 \+ = -1) then\n error \"flag = -1 return in proc SFor_pion2\";\n f i;\n twopie := evalf(2*Pi):\n Sminforpion2 := twopie*answeral3:\n \n # Also find the surface area when k = 0.\n gg := x -> abs(-x)*s qrt(1+(fp(x))^2):\n Aal3 := a:\n Bal3 := b:\n iflaga3 := -1:\n \+ Queue_Size := 0:\n (answeral3,errestal3,iflagal3,Nfeal3) :=\n \+ Adaptmw(gg,Aal3,Bal3,Abserral,Relerral,Queue_Size):\n if (iflagal3 = 1) then\n error \"Insufficient queue storage in Adaptmw\";\n f i;\n if (iflagal3 = 2) then\n error \"Too many gg function eval s in Adaptmw\";\n fi;\n if (iflagal3 = -1) then\n error \"fla g = -1 return in proc SfAL\";\n fi;\n Sforkis0 := twopie*answeral3 :\n lprint(\"For theta=Pi/2 and k = 0, the surface area is \",\n \+ Sforkis0);\n lprint(\"For theta=Pi/2, kmin, Smin, S(0) = \",\n \+ kminforpion2,Sminforpion2,Sforkis0);\n return(kminforpion2,Sminfor pion2,Sforkis0):\nend proc:\n\nSfAL := proc(k)\n# Evaluate the residua l for Zeromw when it\n# finds the value of k for which the arc\n# leng th of f(x) from a to k is equal to\n# the arc length from k to b. (See proc\n# Sforpion2 below.)\n global Abserral,Relerral,Aal,Bal,iflaga l,\n Queue_Size,answeral,errestal,Nfeal,fval,\n Aal2,Bal2,iflagal2 ,answeral2,errestal2,\n Nfeal2,gal:\n # lprint(\"Inside proc SfAL \");\n gal := x -> sqrt(1+(fp(x))^2):\n Abserral := 1e-7:\n Rele rral := 1e-7:\n Aal := a:\n Bal := k:\n iflagal := -1:\n Queue _Size := 0:\n (answeral,errestal,iflagal,Nfeal) :=\n Adaptmw(ga l,Aal,Bal,Abserral,Relerral,Queue_Size):\n if (iflagal = 1) then\n \+ error \"Insufficient queue storage in Adaptmw\";\n fi;\n if (i flagal = 2) then\n error \"Too many gal function evals in Adaptmw \";\n fi;\n if (iflagal = -1) then\n error \"flag = -1 return in proc SfAL\";\n fi;\n Aal2 := k:\n Bal2 := b:\n iflagal := \+ -1:\n Queue_Size := 0:\n (answeral2,errestal2,iflagal2,Nfeal2) := \n Adaptmw(gal,Aal2,Bal2,Abserral,Relerral,Queue_Size):\n if (i flagal2 = 1) then\n error \"Insufficient queue storage in Adaptmw \";\n fi;\n if (iflagal2 = 2) then\n error \"Too many gal fun ction evals in Adaptmw\";\n fi;\n if (iflagal2 = -1) then\n e rror \"flag = -1 return in proc SfAL\";\n fi;\n fval := answeral - answeral2:\n return(fval):\nend proc:\n\nFindPivotPoint := proc(a,b )\n# Find the point (c,f(c)) for which the arc length\n# of f(x) from \+ a to c is equal to the arc length\n# from c to b.\n local From_Where ,bbal,ccal,abserral,relerral,\n ff,abserrZ,relerrZ,flag,bb,cb,residu al,nfeal,\n cpivot,fcpivot:\n lprint(\"Inside proc FindPivotPoint \");\n From_Where := 10:\n (abserrZ,relerrZ) := Set_Tols(From_Wher e):\n #abserrZ := 1e-5:\n #relerrZ := 1e-5:\n bbal := a:\n cca l := b:\n flag := -1:\n (bb,cb,residual,flag,nfeal) :=\n Zero mw(SfAL,bbal,ccal,abserrZ,relerrZ):\n if (flag <> -1) then\n # \+ lprint(\"Results follow ...\");\n if flag = 0 then\n lpri nt(\"flag = 0. A zero was successfully found.\");\n elif flag = 1 then\n lprint(\"flag = 1. Too many function evals in Adaptmw \");\n elif flag = 2 then\n lprint(\"flag = 2. A pole may have been located.\");\n fi;\n # lprint(\" bbal = \", bbal) ;\n # lprint(\" ccal = \", ccal);\n # lprint(\" Number of fu nction evaluations = \", nfeal);\n cb := bb:\n else\n erro r \" A fatal error occurred in Zeromw.\";\n fi;\n # A horizontal l ine z = -x = k corresponds to a\n # vertical line x = -k.\n cpivot := cb:\n fcpivot := evalf(f(cpivot)):\n lprint(\"cpivot, fcpivot \+ = \", cpivot, fcpivot);\n return(cpivot,fcpivot):\nend proc:\n\ntcs \+ := proc(m)\n# Define the necessary transformation of\n# axes informati on. The function f(x) and\n# its derivative fp(x) are defined at the\n # beginning of the main procedure.\n local cost,sint,atan,wofx,zofx, wpofx,\n zpofx,theta,Pie,My_Infinityp,My_Infinitym:\n My_Infinityp := Float(infinity):\n My_Infinitym := Float(-infinity):\n if (m = My_Infinityp) then\n # Vertical rotation by 90 degrees ...\n \+ atan := tan@@(-1):\n cost := 0:\n sint := 1:\n Pie \+ := evalf(Pi):\n theta := evalf(Pi/2):\n wofx := x -> (f( x)-beta):\n zofx := x -> -x:\n wpofx := x -> fp(x):\n \+ zpofx := x -> 0*x - 1:\n elif (m = My_Infinitym) then\n # V ertical rotation by -90 degrees ...\n atan := tan@@(-1):\n \+ cost := 0:\n sint := -1:\n theta := evalf(atan(m)):\n P ie := evalf(Pi):\n theta := evalf(-Pi/2):\n wofx := x -> -(f(x)-beta):\n zofx := x -> x:\n wpofx := x -> -fp(x) :\n zpofx := x -> 0*x + 1:\n else\n atan := tan@@(-1):\n theta := evalf(atan(m)):\n cost := evalf(cos(theta)):\n \+ sint := evalf(sin(theta)):\n Pie := evalf(Pi):\n wofx \+ := x -> cost*x+(f(x)-beta)*sint:\n zofx := x -> -sint*x+(f(x )-beta)*cost:\n wpofx := x -> cost+fp(x)*sint:\n zpofx := \+ x -> -sint+fp(x)*cost:\n fi:\n return(theta,cost,sint,atan,wofx,z ofx,wpofx,zpofx,Pie):\nend proc:\n\ntcs2 := proc(m,betarfl)\n# Define \+ the necessary transformation of\n# axes information. The function f(x) and\n# its derivative fp(x) are defined at the\n# beginning of the ma in procedure.\n local cost,sint,atan,wofx,zofx,wpofx,\n zpofx,thet a,Pie,My_Infinityp,My_Infinitym:\n My_Infinityp := Float(infinity): \n My_Infinitym := Float(-infinity):\n if (m = My_Infinityp) then \n # Vertical rotation by 90 degrees ...\n atan := tan@@(-1 ):\n cost := 0:\n sint := 1:\n Pie := evalf(Pi):\n \+ theta := evalf(Pi/2):\n wofx := x -> (f(x)-betarfl):\n \+ zofx := x -> -x:\n wpofx := x -> fp(x):\n zpofx := x - > 0*x - 1:\n elif (m = My_Infinitym) then\n # Vertical rotation by -90 degrees ...\n atan := tan@@(-1):\n cost := 0:\n \+ sint := -1:\n theta := evalf(atan(m)):\n Pie := evalf(Pi ):\n theta := evalf(-Pi/2):\n wofx := x -> -(f(x)-betarfl) :\n zofx := x -> x:\n wpofx := x -> -fp(x):\n zpof x := x -> 0*x + 1:\n else\n atan := tan@@(-1):\n theta : = evalf(atan(m)):\n cost := evalf(cos(theta)):\n sint := e valf(sin(theta)):\n Pie := evalf(Pi):\n wofx := x -> co st*x+(f(x)-betarfl)*sint:\n zofx := x -> -sint*x+(f(x)-betarfl) *cost:\n wpofx := x -> cost+fp(x)*sint:\n zpofx := x -> - sint+fp(x)*cost:\n fi:\n return(theta,cost,sint,atan,wofx,zofx,wpo fx,zpofx,Pie):\nend proc:\n\n# Adaptmw := proc(f,a,b,abserr,relerr)\n# \n# Estimate the definite integral of f(x) from a to b\n# using an ad aptive quadrature scheme based on\n# Gauss-Kronrod (3,7) formulas. \n # \n# Usage:\n# (answer,errest,flag,nfe) := adaptmw(f,a,b,abserr,reler r):\n# Input parameters:\n# f = name of the function defining \+ f(x).\n# a, b = end points of integration interval.\n# abserr \+ = absolute error tolerance desired.\n# relerr = relative error tole rance desired.\n# Output parameters:\n# answer = computed estimate \+ of the integral.\n# errest = estimate of the absolute error in answ er.\n# flag = 0 for normal return.\n# = 1 insufficien t storage in queue (120).\n# = 2 too many function evaluati ons (3577).\n#\n# Adaptmw was adapted from the Matlab script Adapt.m\n # from the book Fundamentals of Numerical Computing\n# by Shampine, Al len, and Preuss.\n#\n# IMPORTANT NOTE:\n# If you exerience difficultie s with Adaptmw or have\n# questions about it, please do not bug the au thors\n# of Numerical Computing. Instead, contact Skip\n# Thompson by \+ email at thompson@radford.edu.\n# \nAdaptmw := proc(f,aval,bval,abserr val,relerrval)\n#\nlocal maxq, maxfe, eps, length, top, bottom,\nanswe r, errest, alpha, beta, h, tol, el, er,\nql, qr, q, e, nfe, flag, Quad mw, Addmw, queue,\na, b, abserr, relerr, A, x, Alpha, Beta, Q, E,\nQko nrod, H, f1, f2, f3, f4, f5, f6, f7, center:\n#\n# Unload the argument s.\n#\na := aval:\nb := bval:\nrelerr := relerrval:\nabserr := abserrv al:\n#\n# Set up the necessary arrays.\n#\nmaxq := 120:\n# queue := ze ros(maxq,4):\nqueue := array(1..maxq,1..4):\nmaxfe := 3577:\n# a = [0. 2684880898683334, 0.1046562260264672, ...\n# 0.4013974147759622, \+ 0.4509165386584744]:\n# x = [0.7745966692414834, 0.9604912687080202,\n # 0.4342437493468026]:\nA := array(1..4):\nx := array(1..3):\nA[1 ] := 0.2684880898683334:\nA[2] := 0.1046562260264672:\nA[3] := 0.40139 74147759622:\nA[4] := 0.4509165386584744:\nx[1] := 0.7745966692414834: \nx[2] := 0.9604912687080202:\nx[3] := 0.4342437493468026:\n#\n# Test \+ the input error tolerances.\n#\neps := evalf(1.0 / 10.0^(Digits)):\nif relerr < 10*eps then\n error \"relerr < 10*eps is not allowed in Ad aptmw.\";\nfi:\nif abserr < 0 then\n error \"abserr <= 0 is not allo wed in Adaptmw.\";\nfi:\n#\n# Initialization.\n#\nlength := 0:\ntop := 1:\nbottom := 1:\nnfe := 0:\n#\n# Form an initial approximation answe r to the\n# integral over [a,b]. If it is not sufficiently\n# accurate , initialize the queue and begin the\n# main loop.\n# q,e) := Quadmw(f ,a,b):\n#\nAlpha := a:\nBeta := b:\nH := evalf(0.5 * (Beta - Alpha)): \ncenter := evalf(Alpha + H):\nf1 := evalf(f(center)):\nf2 := evalf(f( center - H * x[1])):\nf3 := evalf(f(center + H * x[1])):\nf4 := evalf( f(center - H * x[2])):\nf5 := evalf(f(center + H * x[2])):\nf6 := eval f(f(center - H * x[3])):\nf7 := evalf(f(center + H * x[3])):\nQ := eva lf(H * (5 * (f2 + f3) + 8 * f1) / 9):\nQkonrod := evalf(H * (A[1] * (f 2 + f3) +\n A[4] * f1 + A[2] * (f4 + f5)\n \+ + A[3] * (f6 + f7))):\nE := evalf(Qkonrod - Q):\nq := Q:\ne := E: \nnfe := nfe + 7:\nanswer := evalf(q):\nerrest := evalf(e):\n# \nif ev alf(abs(errest)) >\n evalf(max(abserr,relerr*abs(answer))) then\n \+ # Addmw(maxq,q,e,a,b,length,bottom):\n queue[bottom,1] := q:\n que ue[bottom,2] := e:\n queue[bottom,3] := a:\n queue[bottom,4] := b: \n length := length + 1:\n bottom := bottom:\n if bottom < maxq \+ then\n bottom := bottom + 1:\n else \n bottom := 1:\n fi :\nfi:\n#\n# Main loop: if queue is empty then return,\n# else subdivi de the top entry.\n#\nwhile length > 0 do\n #\n # Remove the top e ntries from the queue.\n #\n q := evalf(queue[top,1]):\n e := ev alf(queue[top,2]):\n alpha := evalf(queue[top,3]):\n beta := evalf (queue[top,4]):\n length := length - 1:\n if top < maxq then\n \+ top := top + 1:\n else\n top := 1:\n fi:\n h := evalf(0.5 *(beta-alpha)):\n #\n # (ql,el) := Quadmw(f,alpha,alpha+h):\n Al pha := alpha:\n Beta := evalf(alpha + h):\n H := evalf(0.5 * (Beta - Alpha)):\n center := evalf(Alpha + H):\n f1 := evalf(f(center)) :\n f2 := evalf(f(center - H * x[1])):\n f3 := evalf(f(center + H \+ * x[1])):\n f4 := evalf(f(center - H * x[2])):\n f5 := evalf(f(cen ter + H * x[2])):\n f6 := evalf(f(center - H * x[3])):\n f7 := eva lf(f(center + H * x[3])):\n Q := evalf(H * (5 * (f2 + f3) + 8 * f1) \+ / 9): \n Qkonrod := evalf(H * (A[1] * (f2 + f3) +\n \+ A[4] * f1 + A[2] * (f4 + f5)\n + A[3] * (f6 + f7 ))):\n E := evalf(Qkonrod - Q):\n ql := Q:\n el := E:\n nfe := nfe + 7:\n # \n # (qr,er) := Quadmw(f,alpha+h,beta):\n Alpha := alpha + h:\n Beta := beta:\n H := evalf(0.5 * (Beta - Alpha)):\n \+ center := evalf(Alpha + H):\n f1 := evalf(f(center)):\n f2 := ev alf(f(center - H * x[1])):\n f3 := evalf(f(center + H * x[1])):\n \+ f4 := evalf(f(center - H * x[2])):\n f5 := evalf(f(center + H * x[2] )):\n f6 := evalf(f(center - H * x[3])):\n f7 := evalf(f(center + \+ H * x[3])):\n Q := evalf(H * (5 * (f2 + f3) + 8 * f1) / 9):\n Qkon rod := evalf(H * (A[1] * (f2 + f3) +\n A[4] * f1 + \+ A[2] * (f4 + f5)\n + A[3] * (f6 + f7))):\n E := e valf(Qkonrod - Q):\n qr := Q:\n er := E:\n nfe := nfe + 7:\n # \n # Update answer and the error estimate.\n #\n answer := evalf (answer + ((ql + qr) - q)):\n errest := evalf(errest + ((el + er) - \+ e)):\n #\n # Test for failures.\n #\n if length >= maxq - 1 th en\n flag := 1:\n return(answer,errest,flag,nfe):\n fi:\n \+ if nfe >= maxfe then\n flag := 2:\n return(answer,errest,f lag,nfe):\n fi:\n #\n # Test for convergence.\n #\n tol := e valf(max(abserr,relerr*abs(answer))):\n if evalf(abs(errest)) <= tol then\n flag := 0:\n return(answer,errest,flag,nfe):\n fi: \n #\n # Add new subintervals to queue if errors\n # are too big .\n #\n tol := evalf(tol * h / (b - a)):\n if evalf(abs(el)) > t ol then\n # Addmw(maxq,ql,el,alpha,alpha+h,length,bottom):\n \+ queue[bottom,1] := ql:\n queue[bottom,2] := el:\n queue [bottom,3] := alpha:\n queue[bottom,4] := alpha + h:\n len gth := length + 1:\n bottom := bottom:\n if bottom < maxq \+ then\n bottom := bottom + 1:\n else\n bottom : = 1:\n fi:\n fi:\n if evalf(abs(er)) > tol then\n # Ad dmw(maxq,qr,er,alpha+h,beta,length,bottom):\n queue[bottom,1] := qr:\n queue[bottom,2] := er:\n queue[bottom,3] := alpha + h:\n queue[bottom,4] := beta:\n length := length + 1:\n \+ bottom := bottom:\n if bottom < maxq then\n bottom := bottom + 1:\n else\n bottom := 1:\n fi:\n f i:\nod: # End of while loop\nflag := 0:\nreturn(answer,errest,flag,nfe ):\nend proc:\n\nGOLD := proc(f::procedure,a::numeric,b::numeric,T::nu meric)\n#\n# To perform a the Golden Search maximization for a\n# unim odal function do the following.\n#\n# 1) Provide a procedure to evalu ate f(x).\n# 2) Type (fmin,xmin,iters) := GOLD(f,a,b,tolerance):\n# \+ for specific values of a, b, and the tolerance.\n# 3) The last \+ output provided is the midpoint of the\n# final interval and th e value of f(x) at that\n# point.\n#\n# Note: To minimize a uni modal function, maximize -f.\n#\n# This Golden Search Algorithm was wr itten originally by\n# Dr. William P. Fox, Department of Mathematic s,\n# Francis Marion University, Florence, SC 29501\n# Dr. Ma rgie Witherspoon, Department of Computer Science,\n# Francis Mar ion University, Florence, SC 29501\n#\n# This procedure performs the G olden Section Search algorithm\n# to find the maximum of a unimodal fu nction, f(x), over an\n# interval, a < x < b. The program calculates t he number of\n# iterations required to insure the final interval is wi thin\n# the user-specified tolerance. This is found by solving for\n# \+ the smallest value of k that makes this inequality true:\n# (b-a)(0.61 8^k)< tolerance. The Golden section concept\n# involves placing two ex periments between [a,b] using the\n# Golden section ratios. One experi ment is placed at position,\n# a + 0.382(b-a), and the other at positi on, a + 0.618(b-a).\n# The function, to be maximized, is evaluated at \+ these two\n# points and the functional values are compared. We want to \n# keep the larger functional value and its corresponding\n# opposite endpoint. At the end of the required iterations,\n# the final interva l is the answer. At times when the final\n# answer must be a single po int and not an interval, the\n# convention of selecting the midpoint i s provided.\n#\n# The Golden Search ratios con6 and con3 are defined b elow\n# in the following manner. At each iteration two new tests\n# po ints are added. They are chosen such that the following\n# two ratios \+ are equal:\n#\n# Length of current interval Length of larger fracti on\n# -------------------------- = ------------------------- \n# Len gth of larger fraction Length of smaller fraction\n#\n# If [a,b] is the current interval, this gives\n#\n# (b-a) r*(b-a)\n# ---- --- = ---------- \n# r*(b-a) (1-r)*(b-a)\n#\n# This leads to t he equation r^2+r-1 =0. The\n# positive solution of this equation is\n # r = con6 ~ 0.618.\n# Also, 1-r = 1 - con6 = con3 ~ 0.382.\n#\n local x1,x2, N, val;\n global con6, con3;\n N := 0;\n #\n # Define the Gold en Search ratios.\n con6 := evalf((sqrt(5)-1)/2);\n con3 := evalf(1-co n6);\n # \n # Determine the first two endpoints x1 and x2.\n x1 := a+c on3*(b-a);\n x2 := a+con6*(b-a);\n # \n # Determine the number of iter ations to be performed.\n N := ceil((ln(T/(b-a))/ln(con6)));\n # lprin t(\"In GOLD, no. of iterates =\", N);\n # Perform the maximization.\n iterate(f,a,b,N,x1,x2);\n #\n # The value f(midpoint) is found. This \+ midpoint is the\n # midpoint of the interval found in the Nth iteratio n.\n val := f(mdpt);\n return(val,mdpt,N);\n end proc:\n \n iterate : = proc(f::procedure,a::numeric,b::numeric,\n N::posint ,x1::numeric,x2::numeric)\n #\n local x1n,x2n,an,bn,fx1,fx2,j;\n globa l mdpt,fkeep,xkeep,pos;\n #\n # The first endpoint is assigned to x1n( 1) and \n # the second endpoint is assigned to x2n(1).\n x1n(1) := x1; \n x2n(1) := x2;\n # \n # Assign the original endpoints to variables \n # an(1) and bn(1).\n an(1) := a;\n bn(1) := b;\n # \n # The followi ng loops runs from 2 to N (the number\n # of iterations for the length of the last interval\n # to have length less than the tolerance speci fied.\n for j from 2 to N do\n fx1(j-1) := f(x1n(j-1));\n fx2(j- 1) := f(x2n(j-1));\n if fx1(j-1) <= fx2(j-1) then\n an(j) := \+ x1n(j-1);\n bn(j) := bn(j-1);\n x1n(j) := x2n(j-1);\n \+ x2n(j) := an(j)+con6*(bn(j)-an(j));\n else\n an(j) := an(j- 1);\n bn(j) := x2n(j-1);\n x2n(j) := x1n(j-1);\n x1n (j) := an(j)+con3*(bn(j)-an(j));\n fi;\n #\n # For the Nth in terval a midpoint is found. After finding \n # the midpoint, an err or check is performed to ensure the\n # midpoint of the interval wi thin the specified tolerance\n # is the max value for the given fun ction.\n # \n if (j = N) then\n mdpt := (an(j) + bn(j)) /2;\n if (f(an(j)) > f(bn(j))) and (f(an(j)) > f(mdpt)) then \n \+ fkeep := f(an(j));\n xkeep := an(j);\n pos \+ := \"left end\";\n else\n if (f(bn(j)) > f(mdpt)) and ( f(bn(j)) < f(an(j)))\n then\n fkeep := f(bn(i)); \n xkeep := bn(i);\n pos := \"right end\";\n \+ else\n fkeep := f(mdpt);\n xkeep := (a n(j) + bn(j))/2;\n pos := \"midpoint\";\n fi;\n \+ fi;\n fi;\n od;\n end proc:\n\n# Copy of Adaptmw.\n# Adaptmw2 \+ := proc(f,a,b,abserr,relerr)\n# \n# Estimate the definite integral of \+ f(x) from a to b\n# using an adaptive quadrature scheme based on\n# Ga uss-Kronrod (3,7) formulas. \n# \n# Usage:\n# (answer,errest,flag,nfe ) := Adaptmw2(f,a,b,abserr,relerr):\n# Input parameters:\n# f \+ = name of the function defining f(x).\n# a, b = end points of int egration interval.\n# abserr = absolute error tolerance desired.\n# relerr = relative error tolerance desired.\n# Output parameters:\n # answer = computed estimate of the integral.\n# errest = estima te of the absolute error in answer.\n# flag = 0 for normal retur n.\n# = 1 insufficient storage in queue (120).\n# \+ = 2 too many function evaluations (3577).\n#\n# Adaptmw2 was adapted from the Matlab script Adapt.m\n# from the book Fundamentals of Numer ical Computing\n# by Shampine, Allen, and Preuss.\n#\n# IMPORTANT NOTE :\n# If you exerience difficulties with Adaptmw2 or have\n# questions \+ about it, please do not bug the authors\n# of Numerical Computing. Ins tead, contact Skip\n# Thompson by email at thompson@radford.edu.\n# \n Adaptmw2 := proc(f,aval,bval,abserrval,relerrval)\n#\nlocal maxq, maxf e, eps, length, top, bottom,\nanswer, errest, alpha, beta, h, tol, el, er,\nql, qr, q, e, nfe, flag, Quadmw2, Addmw2, queue,\na, b, abserr, \+ relerr, A, x, Alpha, Beta, Q, E,\nQkonrod2, H, f1, f2, f3, f4, f5, f6, f7, center:\n#\n# Unload the arguments.\n#\na := aval:\nb := bval:\nr elerr := relerrval:\nabserr := abserrval:\n#\n# Set up the necessary a rrays.\n#\nmaxq := 120:\n# queue := zeros(maxq,4):\nqueue := array(1.. maxq,1..4):\nmaxfe := 3577:\n# a = [0.2684880898683334, 0.104656226026 4672, ...\n# 0.4013974147759622, 0.4509165386584744]:\n# x = [0.7 745966692414834, 0.9604912687080202,\n# 0.4342437493468026]:\nA : = array(1..4):\nx := array(1..3):\nA[1] := 0.2684880898683334:\nA[2] : = 0.1046562260264672:\nA[3] := 0.4013974147759622:\nA[4] := 0.45091653 86584744:\nx[1] := 0.7745966692414834:\nx[2] := 0.9604912687080202:\nx [3] := 0.4342437493468026:\n#\n# Test the input error tolerances.\n#\n eps := evalf(1.0 / 10.0^(Digits)):\nif relerr < 10*eps then\n error \+ \"relerr < 10*eps is not allowed in Adaptmw2.\";\nfi:\nif abserr < 0 t hen\n error \"abserr <= 0 is not allowed in Adaptmw2.\";\nfi:\n#\n# \+ Initialization.\n#\nlength := 0:\ntop := 1:\nbottom := 1:\nnfe := 0:\n #\n# Form an initial approximation answer to the\n# integral over [a,b ]. If it is not sufficiently\n# accurate, initialize the queue and beg in the\n# main loop.\n# q,e) := Quadmw2(f,a,b):\n#\nAlpha := a:\nBeta \+ := b:\nH := evalf(0.5 * (Beta - Alpha)):\ncenter := evalf(Alpha + H): \nf1 := evalf(f(center)):\nf2 := evalf(f(center - H * x[1])):\nf3 := e valf(f(center + H * x[1])):\nf4 := evalf(f(center - H * x[2])):\nf5 := evalf(f(center + H * x[2])):\nf6 := evalf(f(center - H * x[3])):\nf7 \+ := evalf(f(center + H * x[3])):\nQ := evalf(H * (5 * (f2 + f3) + 8 * f 1) / 9):\nQkonrod2 := evalf(H * (A[1] * (f2 + f3) +\n \+ A[4] * f1 + A[2] * (f4 + f5)\n + A[3] * (f6 + f7))):\n E := evalf(Qkonrod2 - Q):\nq := Q:\ne := E:\nnfe := nfe + 7:\nanswer : = evalf(q):\nerrest := evalf(e):\n# \nif evalf(abs(errest)) >\n eval f(max(abserr,relerr*abs(answer))) then\n # Addmw2(maxq,q,e,a,b,lengt h,bottom):\n queue[bottom,1] := q:\n queue[bottom,2] := e:\n que ue[bottom,3] := a:\n queue[bottom,4] := b:\n length := length + 1: \n bottom := bottom:\n if bottom < maxq then\n bottom := bott om + 1:\n else \n bottom := 1:\n fi:\nfi:\n#\n# Main loop: if queue is empty then return,\n# else subdivide the top entry.\n#\nwhil e length > 0 do\n #\n # Remove the top entries from the queue.\n \+ #\n q := evalf(queue[top,1]):\n e := evalf(queue[top,2]):\n alp ha := evalf(queue[top,3]):\n beta := evalf(queue[top,4]):\n length := length - 1:\n if top < maxq then\n top := top + 1:\n else \n top := 1:\n fi:\n h := evalf(0.5*(beta-alpha)):\n #\n \+ # (ql,el) := Quadmw2(f,alpha,alpha+h):\n Alpha := alpha:\n Beta := evalf(alpha + h):\n H := evalf(0.5 * (Beta - Alpha)):\n center := evalf(Alpha + H):\n f1 := evalf(f(center)):\n f2 := evalf(f(cente r - H * x[1])):\n f3 := evalf(f(center + H * x[1])):\n f4 := evalf (f(center - H * x[2])):\n f5 := evalf(f(center + H * x[2])):\n f6 \+ := evalf(f(center - H * x[3])):\n f7 := evalf(f(center + H * x[3])): \n Q := evalf(H * (5 * (f2 + f3) + 8 * f1) / 9): \n Qkonrod2 := ev alf(H * (A[1] * (f2 + f3) +\n A[4] * f1 + A[2] * (f 4 + f5)\n + A[3] * (f6 + f7))):\n E := evalf(Qkon rod2 - Q):\n ql := Q:\n el := E:\n nfe := nfe + 7:\n # \n # \+ (qr,er) := Quadmw2(f,alpha+h,beta):\n Alpha := alpha + h:\n Beta : = beta:\n H := evalf(0.5 * (Beta - Alpha)):\n center := evalf(Alph a + H):\n f1 := evalf(f(center)):\n f2 := evalf(f(center - H * x[1 ])):\n f3 := evalf(f(center + H * x[1])):\n f4 := evalf(f(center - H * x[2])):\n f5 := evalf(f(center + H * x[2])):\n f6 := evalf(f( center - H * x[3])):\n f7 := evalf(f(center + H * x[3])):\n Q := e valf(H * (5 * (f2 + f3) + 8 * f1) / 9):\n Qkonrod2 := evalf(H * (A[1 ] * (f2 + f3) +\n A[4] * f1 + A[2] * (f4 + f5)\n \+ + A[3] * (f6 + f7))):\n E := evalf(Qkonrod2 - Q):\n qr := Q:\n er := E:\n nfe := nfe + 7:\n #\n # Update answer and the error estimate.\n #\n answer := evalf(answer + ((ql + qr) - q)):\n errest := evalf(errest + ((el + er) - e)):\n #\n # Tes t for failures.\n #\n if length >= maxq - 1 then\n flag := 1: \n return(answer,errest,flag,nfe):\n fi:\n if nfe >= maxfe th en\n flag := 2:\n return(answer,errest,flag,nfe):\n fi:\n \+ #\n # Test for convergence.\n #\n tol := evalf(max(abserr,rele rr*abs(answer))):\n if evalf(abs(errest)) <= tol then\n flag := 0:\n return(answer,errest,flag,nfe):\n fi:\n #\n # Add new subintervals to queue if errors\n # are too big.\n #\n tol := e valf(tol * h / (b - a)):\n if evalf(abs(el)) > tol then\n # Ad dmw2(maxq,ql,el,alpha,alpha+h,length,bottom):\n queue[bottom,1] \+ := ql:\n queue[bottom,2] := el:\n queue[bottom,3] := alpha :\n queue[bottom,4] := alpha + h:\n length := length + 1: \n bottom := bottom:\n if bottom < maxq then\n bo ttom := bottom + 1:\n else\n bottom := 1:\n fi:\n fi:\n if evalf(abs(er)) > tol then\n # Addmw2(maxq,qr,er,al pha+h,beta,length,bottom):\n queue[bottom,1] := qr:\n queu e[bottom,2] := er:\n queue[bottom,3] := alpha + h:\n queue [bottom,4] := beta:\n length := length + 1:\n bottom := bo ttom:\n if bottom < maxq then\n bottom := bottom + 1:\n else\n bottom := 1:\n fi:\n fi:\nod: # End of \+ while loop\nflag := 0:\nreturn(answer,errest,flag,nfe):\nend proc:\n\n # Copy of GOLD\nGOLD2 := proc(f::procedure,a::numeric,b::numeric,T::nu meric)\n#\n# To perform a the Golden Search maximization for a\n# unim odal function do the following.\n#\n# 1) Provide a procedure to evalu ate f(x).\n# 2) Type (fmin,xmin,iters) := GOLD2(f,a,b,tolerance):\n# \+ for specific values of a, b, and the tolerance.\n# 3) The last output provided is the midpoint of the\n# final interval and t he value of f(x) at that\n# point.\n#\n# Note: To minimize a un imodal function, maximize -f.\n#\n# This Golden Search Algorithm was w ritten originally by\n# Dr. William P. Fox, Department of Mathemati cs,\n# Francis Marion University, Florence, SC 29501\n# Dr. M argie Witherspoon, Department of Computer Science,\n# Francis Ma rion University, Florence, SC 29501\n#\n# This procedure performs the \+ Golden Section Search algorithm\n# to find the maximum of a unimodal f unction, f(x), over an\n# interval, a < x < b. The program calculates \+ the number of\n# iterations required to insure the final interval is w ithin\n# the user-specified tolerance. This is found by solving for\n# the smallest value of k that makes this inequality true:\n# (b-a)(0.6 18^k)< tolerance. The Golden section concept\n# involves placing two e xperiments between [a,b] using the\n# Golden section ratios. One exper iment is placed at position,\n# a + 0.382(b-a), and the other at posit ion, a + 0.618(b-a).\n# The function, to be maximized, is evaluated at these two\n# points and the functional values are compared. We want t o\n# keep the larger functional value and its corresponding\n# opposit e endpoint. At the end of the required iterations,\n# the final interv al is the answer. At times when the final\n# answer must be a single p oint and not an interval, the\n# convention of selecting the midpoint \+ is provided.\n#\n# The Golden Search ratios con62 and con32 are define d below\n# in the following manner. At each iteration two new tests\n# points are added. They are chosen such that the following\n# two rati os are equal:\n#\n# Length of current interval Length of larger fra ction\n# -------------------------- = ------------------------- \n# \+ Length of larger fraction Length of smaller fraction\n#\n# If [a,b] is the current interval, this gives\n#\n# (b-a) r*(b-a)\n# - ------ = ---------- \n# r*(b-a) (1-r)*(b-a)\n#\n# This leads t o the equation r^2+r-1=0. The\n# positive solution of this equation is \n# r = con62 ~ 0.618.\n# Also, 1-r = 1 - con62 = con32 ~ 0.382.\n#\n \+ local x1,x2, N, val;\n global con62, con32;\n N := 0;\n #\n # Define t he Golden Search ratios.\n con62 := evalf((sqrt(5)-1)/2);\n con32 := e valf(1-con62);\n # \n # Determine the first two endpoints x1 and x2.\n x1 := a+con32*(b-a);\n x2 := a+con62*(b-a);\n # \n # Determine the nu mber of iterations to be performed.\n N := ceil((ln(T/(b-a))/ln(con62) ));\n # lprint(\"In GOLD2, no. of iterates =\", N);\n # Perform the m aximization.\n iterate2(f,a,b,N,x1,x2);\n #\n # The value f(midpoint) \+ is found. This midpoint is the\n # midpoint of the interval found in t he Nth iteration.\n val := f(mdpt2);\n return(val,mdpt2,N);\n end proc :\n \n iterate2 := proc(f::procedure,a::numeric,b::numeric,\n \+ N::posint,x1::numeric,x2::numeric)\n #\n local x1n,x2n,an,bn, fx1,fx2,j;\n global mdpt2,fkeep2,xkeep2,pos2;\n #\n # The first endpoi nt is assigned to x1n(1) and \n # the second endpoint is assigned to x 2n(1).\n x1n(1) := x1;\n x2n(1) := x2;\n # \n # Assign the original e ndpoints to variables\n # an(1) and bn(1).\n an(1) := a;\n bn(1) := b; \n # \n # The following loops runs from 2 to N (the number\n # of iter ations for the length of the last interval\n # to have length less tha n the tolerance specified.\n for j from 2 to N do\n fx1(j-1) := f(x 1n(j-1));\n fx2(j-1) := f(x2n(j-1));\n if fx1(j-1) <= fx2(j-1) t hen\n an(j) := x1n(j-1);\n bn(j) := bn(j-1);\n x1n(j ) := x2n(j-1);\n x2n(j) := an(j)+con62*(bn(j)-an(j));\n else \n an(j) := an(j-1);\n bn(j) := x2n(j-1);\n x2n(j) : = x1n(j-1);\n x1n(j) := an(j)+con32*(bn(j)-an(j));\n fi;\n \+ #\n # For the Nth interval a midpoint is found. After finding \n \+ # the midpoint, an error check is performed to ensure the\n # mid point of the interval within the specified tolerance\n # is the max value for the given function.\n # \n if (j = N) then\n \+ mdpt2 := (an(j) + bn(j))/2;\n if (f(an(j)) > f(bn(j))) and (f(a n(j)) > f(mdpt2)) then \n fkeep2 := f(an(j));\n xkee p2 := an(j);\n pos2 := \"left end\";\n else\n \+ if (f(bn(j)) > f(mdpt2)) and (f(bn(j)) < f(an(j)))\n then\n \+ fkeep2 := f(bn(i));\n xkeep2 := bn(i);\n \+ pos2 := \"right end\";\n else\n fkeep2 := f (mdpt2);\n xkeep2 := (an(j) + bn(j))/2;\n pos2 := \"midpoint\";\n fi;\n fi;\n fi;\n od;\n end proc :\n\n# START PROCEDURE Zeromw\n# Zeromw := proc(f,b,c,abserr,relerr)\n #\n# Zeromw computes a root of the nonlinear equation f(x) = 0 when\n# f(x) is a continuous real function of a single real variable x.\n# Th e method used is a combination of bisection and the secant rule.\n#\n# Input parameters:\n# f = the name of the function f(x ).\n# b, c = values of x such that f(b)*f(c) <= 0.\n# \+ abserr, relerr = absolute and relative error tolerances.\n# \+ The stopping criterion is:\n# abs(b-c) <= \+ 2*max(abserr,abs(b)*relerr).\n# Output parameters:\n# b, c = s ee flag options.\n# residual = value of final residual f(b).\n# \+ flag = 0 for normal return; f(b)*f(c) < 0 and the\n# \+ stopping criterion is met (or f(b) = 0). b is\n# \+ always selected so that abs(f(b)) <= abs(f(c)).\n# = 1 \+ if too many function evaluations were made; in\n# thi s version 500 are allowed.\n# = 2 if abs(f(b)) is more th an 100 times the larger\n# of the residuals for the i nput b and c. b and c\n# probably approximate a pole \+ rather than a root.\n# Zeromw was adapted from the Matlab script Zero. m from the book\n# Fundamentals of Numerical Computing by Shampine, Al len, and Preuss.\n#\n# IMPORTANT NOTE:\n# If you exerience difficultie s with Zeromw or have questions about\n# it, please do not bug the aut hors of Numerical Computing. Instead,\n# contact Skip Thompson by emai l at thompson@radford.edu.\n#\n# Example usage:\n# f := x -> x - co s(x);\n# b := 0.2;\n# c := 0.9;\n# abserr := 1e-5;\n# rele rr := 1e-5;\n# (b,c,residual,flag,nfe) := Zeromw(f,b,c,abserr,reler r);\n#\n\nZeromw := proc(f,bval,cval,abserr,relerr)\n# \n# Declare loc al variables.\n#\nlocal digits, eps, maxfe, iter, Tols_Are_OK, Endpoin t_Zero, count, \\\n width, Bracketed, initial_residual, cmb, tol, p, q, bisect, \\\n acmb, a, fa, b, fb, c, fc, nfe, residual, fla g:\n#\n# Unload the bracketing interval endpoints.\n#\nb := bval:\nc : = cval:\n#\n# Approximate unit roundoff.\n#\ndigits := Digits:\neps := 1 / 10^(digits-1):\n#\n# Define the maximum number of iterations:\n# \nmaxfe := 500:\niter := 0:\nnfe := 0:\n#\n# Test the input tolerances .\n#\nTols_Are_OK := 1:\nif (relerr < 10*eps) then\n Tols_Are_OK := \+ 0:\n residual := -1:\n error \"relerr < 10*eps is not allowed in Z ero.mw\";\nfi:\nif (abserr < 0) then\n residual := -1:\n Tols_Are_ OK := 0:\n error \"abserr < 0 is not allowed in Zero.mw.\";\nfi:\n# \n# Check if f is 0 at one of the endpoints; and\n# perform initializa tion.\n#\nEndpoint_Zero := 0:\ncount := 0:\nwidth := evalf(abs(b - c)) :\na := c:\nfa := evalf(f(a)):\nnfe := 1:\nif (fa = 0) then\n b := a :\n residual := 0:\n flag := 0:\n Endpoint_Zero := 1:\nfi:\nfb : = evalf (f(b)):\nnfe := 2:\nif (fb = 0) then\n residual := 0:\n fl ag := 0:\n Endpoint_Zero := 1:\nfi:\n#\n# Check if [a,b] brackets a \+ zero.\n#\nBracketed := 1:\n# if ((fa > 0) = (fb > 0)) then\nif ((fa > \+ 0 and fb > 0) or (fa <= 0 and fb <= 0)) then\n Bracketed := 0:\n e rror \"Initial interval [b,c] does not bracket a root.\";\nfi:\nif ((T ols_Are_OK = 1) and (Endpoint_Zero = 0) and (Bracketed = 1))\nthen\n \+ initial_residual := max(abs(fa),abs(fb)):\n fc := fa:\n #\n # P erform the iteration.\n #\n for iter from 1 to maxfe do\n if \+ (abs(fc) < abs(fb)) then\n #\n # Interchange b and c s o that abs(f(b)) <= abs(f(c)).\n #\n a := b:\n \+ fa := fb:\n b := c:\n fb := fc:\n c := a:\n \+ fc := fa:\n fi:\n cmb := evalf(0.5*(c-b)):\n acm b := abs(cmb):\n tol := max(abserr, abs(b)*relerr):\n #\n \+ # Test the stopping criterion and function count.\n #\n i f (acmb <= tol) then\n residual := fb:\n if (abs(resid ual) > 100*initial_residual) then\n flag := 2:\n \+ break:\n else\n flag := 0:\n break:\n \+ fi:\n fi:\n #\n # Calculate new iterate implicit ly as b+p/q where we arrange\n # p >= 0. The implicit form is use d to prevent overflow.\n #\n p := evalf((b-a)*fb):\n q \+ := fa-fb:\n if (p < 0) then\n p := -p:\n q := -q: \n fi:\n #\n # Update a; check if reduction in the size of bracketing\n # interval is being reduced at a reasonable rate . If not,\n # bisect until it is.\n #\n a := b:\n \+ fa := fb:\n count := count + 1:\n bisect := 0:\n if (co unt >= 4) then\n if (8*acmb >= width) then\n bisect := 1:\n else\n count := 0:\n width := a cmb:\n fi:\n fi:\n #\n # Test for too small a c hange.\n #\n if (p <= abs(q)*tol) then\n #\n \+ # If the change is too small, increment by tol.\n #\n \+ b := evalf(b + tol * sign(cmb)):\n else\n #\n # \+ The root ought to be between b and (c+b)/2.\n #\n if ( p < cmb*q) then\n # Use secant rule.\n b := eval f(b + p/q):\n else\n # Use bisection.\n \+ bisect := 1:\n fi:\n fi:\n if (bisect = 1) then\n \+ b := c - cmb:\n fi:\n #\n # The computation for n ew iterate b has been completed.\n #\n fb := evalf(f(b)):\n \+ nfe := nfe + 1:\n if (fb = 0) then\n c := b:\n \+ residual := fb:\n flag := 0:\n break:\n fi:\n \+ # if ((fb > 0) = (fc > 0)) then\n if ((fb > 0 and fc > 0) or (fb <= 0 and fc <= 0)) then\n c := a:\n fc := fa:\n \+ fi:\n od:\n if (nfe = maxfe) then\n flag := 1:\n fi:\nf i:\nreturn(b,c,residual,flag,nfe)\nend proc:\n" }}}{EXCHG {PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 6719 "# Begin main procedure.\nwith(plots):\nwith( Student[Calculus1]):\nwith(plottools):\n\n# Note:\n# To minimize surfa ce area Smin for a\n# given slope m,\n# Main -> SMG -> GOLD -> SA\n # To minimize Smin for all slopes m,\n# Main -> SMG2 -> GOLD2 -> SM G3 -> SA2\n\n# Define the numerical precision.\nDigits := 16:\n\n# Def ine the number of grid points for the\n# plots of the surfaces with mi nimum area.\nngridw := 'ngridw': ngridz := 'ngridz':\nngridw := 51: ng ridz := 51:\n \n# Note fmin and fmax are used only for\n# graphical pu rposes when the minimum\n# lines are animated. Any bounds will do.\n\n Problem := 15;\nif (Problem = 1) then\n a := 0;\n b := 1;\n f := x -> x^2;\n fmin := 0:\n fmax := 1:\nelif (Problem = 2) then\n \+ a := 0;\n pie := evalf(Pi):\n b := 2*pie;\n f := x -> x + sin(x) ;\n fmin := 0:\n fmax := b+1:\nelif (Problem = 3) then\n a := 0. 1;\n b := 1;\n f := x -> sqrt(x);\n fmin := sqrt(a):\n fmax := sqrt(b):\nelif (Problem = 4) then\n a := 0;\n b := 1;\n f := x \+ -> cos(x);\n fmin := cos(1):\n fmax := 1:\nelif (Problem = 5) then \n a := 1;\n b := 2;\n f := x -> x^3-3*x^2+2*x+3;\n fmin := -5 :\n fmax := 5:\nelif (Problem = 6) then\n # Extremely long run tim e ...\n a := 0;\n b := 4;\n pie := evalf(Pi):\n f := x -> exp( -x)*cos(2*pie*x);\n fmin := -1:\n fmax := 1;\nelif (Problem = 7) t hen\n a := 0;\n b := 8;\n f := x -> cos(x);\n fmin := -1:\n \+ fmax := 1:\nelif (Problem = 8) then\n a := -1;\n b := 3;\n f := \+ x -> x*(x-1)*(x-2);\n fmin := f(1+sqrt(3)/3):\n fmax := f(1-sqrt(3 )/3):\nelif (Problem = 9) then\n a := -1; #0;\n b := 3; #2;\n f \+ := x -> -x*(x-1)*(x-2);\n fmin := f(1-sqrt(3)/3):\n fmax := f(1+sq rt(3)/3):\nelif (Problem = 10) then\n a := 0;\n b := 3;\n f := x -> (x-1)^3;\n fmin := -1:\n fmax := 8:\nelif (Problem = 11) then \n a := 0;\n b := 3;\n f := x -> (x-1)^2;\n fmin := 0:\n fma x := 4:\nelif (Problem = 12) then\n # Serpentine\n a := -5;\n b \+ := 5;\n f := x -> (4*x)/(x^2+2);\n fmin := -sqrt(2):\n fmax := s qrt(2):\nelif (Problem = 13) then\n # Pear\n a := 0;\n b := 1.9; # 2;\n f := x -> sqrt(x^3*(2-x)/9);\n fmin := 0:\n fmax := f(3/ 2):\nelif (Problem = 14) then\n a := 0;\n pie := evalf(Pi):\n b \+ := 2*pie;\n f := x -> x - sin(x);\n fmin := 0:\n fmax := f(2*pie ):\nelif (Problem = 15) then\n a := 0;\n pie := evalf(Pi):\n b : = 2*pie;\n f := x -> x*sin(x);\n fmax := 2:\n fmin := -5:\nelif \+ (Problem = 16) then\n # For [-1,4], Smin is symmetric about m=0.\n \+ # a := -1;\n a := 0;\n b := 4;\n f := x -> x*(x-1)*(x-2)*(x-3); \n fmin := -5:\n fmax := 25:\nelif (Problem = 17) then\n a := 1. 1;\n b := 2.9;\n f := x -> sqrt(1-(x-2)^2);\n fmin := min(f(a),f (b)):\n fmax := f(2):\nelif (Problem = 18) then\n a := 0;\n b := 2*evalf(Pi);\n f := x -> 2*sin(x) - cos(2*x);\n fmin := -3/2:\n \+ fmax := 3:\nelif (Problem = 19) then\n a := 0;\n b := evalf(Pi); \n f := x -> sin(x+sin(2*x));\n fmin := 0:\n fmax := 1:\nelif (P roblem = 20) then\n a := 0.01; # 0;\n b := evalf(Pi);\n f := x - > Si(x);\n fmin := -1:\n fmax := 2:\nelif (Problem = 21) then\n \+ a := 0;\n b := 4;\n f := x -> x^3 - 12*x;\n fmin := -16:\n fma x := max(f(a),f(b)):\nelif (Problem = 22) then\n a := 0;\n b := 4; \n f := x -> x^6 - x^5;\n fmin := f(5/6):\n fmax := max(f(a),f(b ),f(0)):\nelif (Problem = 23) then\n a := -1;\n b := 4;\n f := x -> 3*x^4 - 16*x^3 + 18*x^2;\n fmin := f(3):\n fmax := max(f(a),f( b),f(1)):\nelif (Problem = 24) then\n a := -1;\n b := 4.1;\n f : = x -> x^4 - 4*x^3;\n fmin := f(3):\n fmax := max(f(a),f(b),f(0)): \nelif (Problem = 25) then\n a := 3/2+0.1;\n b := 5/2-0.1;\n asi n := sin@@(-1);\n f := x -> asin(2*x-4)/2;\n fmin := -evalf(Pi/4): \n fmax := evalf(Pi/4):\nelif (Problem = 26) then\n # cosh x\n c val := 1;\n a := -cval;\n b := cval;\n f := x -> (exp(x/cval)+ex p(-x/cval))/2;\n fmin := f(0):\n fmax := max(f(a),f(b)):\nelif (Pr oblem = 27) then\n # cosh x\n cval := 1;\n a := -2*cval;\n b : = cval;\n f := x -> (exp(x/cval)+exp(-x/cval))/2;\n fmin := f(0): \n fmax := max(f(a),f(b)):\nelif (Problem = 28) then\n # sinh x\n \+ cval := 1;\n a := -2*cval;\n b := cval;\n f := x -> (exp(x/cva l)-exp(-x/cval))/2;\n fmin := f(a):\n fmax := f(b):\nelif (Problem = 29) then\n # tanh x\n cval := 1;\n a := -2*cval;\n b := cva l;\n f := x -> ((exp(x/cval)-exp(-x/cval))/2)\n / ((exp(x/ cval)+exp(-x/cval))/2);\n fmin := f(a):\n fmax := f(b):\nelif (Pro blem = 30) then\n # sech x\n cval := 1;\n a := -2*cval;\n b := cval;\n f := x -> 1 / ((exp(x/cval)+exp(-x/cval))/2);\n fmin := m in(f(a),f(b)):\n fmax := f(0):\nelif (Problem = 31) then\n a := 1. 1;\n b := 5;\n f := x -> x + sqrt(x^2-1);\n fmin := 0:\n fmax \+ := max(f(a),f(b)):\nelif (Problem = 32) then\n a := 1.1;\n b := 5; \n f := x -> x - sqrt(x^2-1);\n fmin := 0:\n fmax := max(f(a),f( b)):\nelif (Problem = 33) then\n a := 0.1;\n b := 10;\n f := x - > (x + 1/x)/2;\n fmin := 1:\n fmax := max(f(a),f(b)):\nelif (Probl em = 34) then\n # Tractrix\n alpha := 4;\n a := 0.01;\n b := 3 .9;\n f := x -> alpha*ln((alpha+sqrt(alpha^2-x^2))/x)\n \+ - sqrt(alpha^2-x^2);\n fmin := 0:\n fmax := f(a):\nelif (Problem = 35) then\n a := 0;\n pie := evalf(Pi);\n b := 3*pie;\n f := x -> x*cos(x);\n fmin := f(b);\n fmax := f(2*pie):\nelif (Problem = 36) then\n a := 0;\n pie := evalf(Pi);\n b := 3*pie;\n f := x -> x*cos(2*x);\n fmin := -f(b);\n fmax := f(b):\nelif (Problem = \+ 37) then\n a := 0;\n b := 2;\n f := x -> x/(x^2+1);\n fmin := \+ f(-1);\n fmax := f(1):\nelif (Problem = 38) then\n a := -5;\n b \+ := 10;\n f := x -> 10*x/(x^2+1) + 3;\n fmin := f(-1);\n fmax := \+ f(1):\nelif (Problem = 39) then\n a := 1;\n pie := evalf(Pi);\n \+ b := 3*pie;\n f := x -> sqrt(x) * cos(2*x);\n fmax := sqrt(b);\n \+ fmin := -sqrt(b):\nelif (Problem = 40) then\n a := 0;\n pie := ev alf(Pi);\n b := 3*pie;\n f := x -> (1+sin(x)) * cos(2*x);\n fmax := 2;\n fmin := -2:\nelif (Problem = 41) then\n a := 0;\n pie : = evalf(Pi);\n b := 3*pie;\n f := x -> exp(sin(2*x)-x/5);\n fmax := exp(1);\n fmin := 1/fmax;\nelif (Problem = 42) then\n a := 0; \n pie := evalf(Pi);\n b := 3*pie;\n f := x -> x * exp(sin(2*x)) ;\n fmax := 3*b;\n fmin := 0;\nelif (Problem = 43) then\n a := 0 ;\n pie := evalf(Pi);\n b := 3*pie;\n f := x -> cos(x) * exp(sin (2*x));\n fmax := 3;\n fmin := -3;\nelif (Problem = 44) then\n a := -1;\n b := 1;\n f := x -> x^2;\n fmin := 0:\n fmax := 1:\n else\n # Supply other function information here.\n error \"No prob lem number was specified.\";\nfi;\n\na := evalf(a):\nb := evalf(b):\nf min := evalf(fmin):\nfmax := evalf(fmax):\n# Define the vertical plot \+ box for animated lines.\nFmin := fmin-1:\nFmax := fmax+1:\n\n# Define \+ f'(x).\nfp := x -> D(f)(x):\nfp(x);\n\n# Plot f(x).\nplot(f(x),x=a..b, title=`y = f(x)`);\n" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 3284 "D oThis := 'DoThis':\nDoThis := 1:\nif (DoThis = 1 and Problem = 44) the n\nlprint(\"Revolve right side of f(x) about the line x = -1/2 ...\"); \nSurfaceOfRevolution(f(x),x=0..b,axis=vertical,\n distancefromaxis= -1/2,output=plot);#\n #title=`Revolve right side of f(x) about the l ine x = -1/2`);\nSurfaceOfRevolution(f(x),x=0..b,axis=vertical,\n di stancefromaxis=-1/2,output=integral,\n title=`Revolve right side of \+ f(x) about the line x = -1/2`);\nSurfaceOfRevolution(f(x),x=0..b,axis= vertical,\n distancefromaxis=-1/2,\n title=`Revolve f(x) about the line x = -1/2`);\n\nlprint(\"Revolve left side of f(x) about the line x = -1/2 ...\");\nSurfaceOfRevolution(f(x),x=a..0,axis=vertical,\n \+ distancefromaxis=-1/2,output=plot);#,\n #title=`Revolve left side of f(x) about the line x = -1/2`);\nSurfaceOfRevolution(f(x),x=a..0,axis =vertical,\n distancefromaxis=-1/2,output=integral,\n title=`Revol ve left side of f(x) about the line x = -1/2`);\nSurfaceOfRevolution(f (x),x=a..0,axis=vertical,\n distancefromaxis=-1/2,\n title=`Revolv e left side of f(x) about the line x = -1/2`);\n\nprintlevel := 0:\nRi ghtSide :=\nSurfaceOfRevolution(f(x),x=0..b,axis=vertical,\n distanc efromaxis=-1/2,output=plot,\n title=`Revolve y=x^2, x=0..1 about the line x = -1/2`):\nprintlevel := 0:\nLeftSide :=\nSurfaceOfRevolution( f(x),x=a..0,axis=vertical,\n distancefromaxis=-1/2,output=plot,\n \+ title=`Revolve y=x^2, x=-1..0 about the line x = -1/2`):\nprintlevel : = 1:\ndisplay(\{RightSide,LeftSide\},\n title=`Revolve y=x^2, x=-1.. 1 about the line x = -1/2`);\n\nlprint(\"Revolve y=x^2, x=-1..1 about \+ the line x = -1/2 ...\");\nSurfaceOfRevolution(f(x),x=a..b,axis=vertic al,\n distancefromaxis=-1/2,output=plot,\n title=`Revolve y=x^2, x =-1..1 about the line x = -1/2`);\nSurfaceOfRevolution(f(x),x=a..b,ax is=vertical,\n distancefromaxis=-1/2,output=integral,\n title=`Rev olve y=x^2, x=-1..1 about the line x = -1/2`);\nSurfaceOfRevolution(f (x),x=a..b,axis=vertical,\n distancefromaxis=-1/2,\n title=`Revolv e y=x^2, x=-1..1 about the line x = -1/2`);\n\nprintlevel := 0:\n#disp lay(line([0,0], [3,4], color=red, linestyle=dash));\nkaxis :=\nline([- 1/2,0], [-1/2,1], color=black,thickness=3):\nfsq :=\nplot(f(x),x=a..b, title=`f(x) = x^2`,color=blue,thickness=3):\nrfl1 :=\nplot(f(-1-x),x=- 2..-1,title=``,color=red,thickness=3):\nrfl2 :=\nplot(f(-1-x),x=a..0,t itle=``,color=red,thickness=3):\nprintlevel := 1:\ndisplay(\{kaxis,fsq ,rfl1,rfl2\});\ndisplay(\{kaxis,fsq\});\n\n# No titles...\nSurfaceOfRe volution(f(x),x=0..b,axis=vertical,\n distancefromaxis=-1/2,output=p lot,\n title=` `);\nSurfaceOfRevolution(f(x),x=a..0,axis=vertical,\n distancefromaxis=-1/2,output=plot,\n title=` `);\nSurfaceOfRevolu tion(f(x),x=a..b,axis=vertical,\n distancefromaxis=-1/2,output=plot, \n title=` `);\n\nDoThis := 'DoThis':\nDoThis := 0:\nif (DoThis = 1) then\nlprint(\"Revolve y=x^2, x=-1..1 about the line x = -1/2 ...\"); \nVolumeOfRevolution(f(x),x=a..b,axis=vertical,\n distancefromaxis=- 1/2,output=plot,\n title=`Revolve y=x^2, x=-1..1 about the line x = -1/2`);\nVolumeOfRevolution(f(x),x=a..b,axis=vertical,\n distancefr omaxis=-1/2,output=integral,\n title=`Revolve y=x^2, x=-1..1 about \+ the line x = -1/2`);\nVolumeOfRevolution(f(x),x=a..b,axis=vertical,\n \+ distancefromaxis=-1/2,\n title=`Revolve y=x^2, x=-1..1 about the l ine x = -1/2`);\nfi;\nfi;\n" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 560 "DoThis := 'DoThis':\nDoThis := 0:\nif (DoThis = 1) then\n# Find t he value of c for which the arc length of\n# f from a to c is equal to the length from c to b.\n(cpivot,fcpivot) := FindPivotPoint(a,b):\nfi ;\n\n# If Use_Pivot = 1, the outer Golden Search will be\n# performed \+ using the pivot point to determine betaS\n# rather than a full blown i nner minimization. The\n# results must be interpreted carefully howeve r.\n# This can be done before doing a full optimization\n# since it gi ves an idea where the actual minimum\n# occurs and is much, much faste r.\nUse_Pivot := 0;\n" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 2601 " DoThis := 'DoThis':\nDoThis := 0:\nif (DoThis = 1) then\n Sofk := k \+ -> int((zofx(x)-k)*sqrt(1+(fp(x))^2),x=a..c) -\n int((z ofx(x)-k)*sqrt(1+(fp(x))^2),x=c..b);\n D(Sofk)(k);\n Sofc := c -> \+ int((zofx(x)-k)*sqrt(1+(fp(x))^2),x=a..c) -\n int((zofx (x)-k)*sqrt(1+(fp(x))^2),x=c..b);\n D(Sofc)(c);\n Sofck := (c,k) - > int((zofx(x)-k)*sqrt(1+(fp(x))^2),x=a..c) -\n in t((zofx(x)-k)*sqrt(1+(fp(x))^2),x=c..b);\n lprint(\"Partial wrt c .. .\");\n D[1](Sofck);\n lprint(\"Partial wrt k ...\");\n D[2](Sof ck);\n Sofc1c2k := (c1,c2,k) -> \n int((zofx(x)-k)*sqrt(1+(fp(x ))^2),x=a ..c1) -\n int((zofx(x)-k)*sqrt(1+(fp(x))^2),x=c1..c2) \+ +\n int((zofx(x)-k)*sqrt(1+(fp(x))^2),x=c2..b );\n lprint(\"Par tial wrt c1 ...\");\n D[1](Sofc1c2k);\n lprint(\"Partial wrt c2 .. .\");\n D[2](Sofc1c2k);\n lprint(\"Partial wrt k ...\");\n D[3]( Sofc1c2k);\n lprint(\"________________________________________\");\n \n Sofmb := (m,beta) -> (1/sqrt(m^2+1))*\n int((f(x)-m*x-beta)*sqr t(1+(fp(x))^2),x=a..b);\n lprint(\"First artial wrt m ...\");\n D[ 1](Sofmb);\n lprint(\"First partial wrt beta ...\");\n D[2](Sofmb) ;\n lprint(\"Second partial wrt m ...\");\n D[1,1](Sofmb);\n lpr int(\"Second partial wrt beta ...\");\n D[2,2](Sofmb);\n lprint(\" ________________________________________\");\n\n Sofmb := (m,beta) - >\n (1/sqrt(m^2+1))*\n int((f(x)-m*x-beta)*sqrt(1+(fp(x))^2),x=a.. c1)-\n (1/sqrt(m^2+1))*\n int((f(x)-m*x-beta)*sqrt(1+(fp(x))^2),x= c1..b);\n lprint(\"First partial wrt m ...\");\n SA1 := D[1](Sofmb );\n lprint(\"Second partial wrt m ...\");\n SA11 := D[1,1](Sofmb) ;\n lprint(\"First partial wrt beta ...\");\n SA2 := D[2](Sofmb); \n lprint(\"Second partial wrt beta ...\");\n SA22 := D[2,2](Sofmb );\n lprint(\"Mixed partial ...\");\n SA12 := D[1,2](Sofmb);\n l print(\"________________________________________\");\n #HSA := (m,be ta) -> SA11(m,beta)*SA22(m,beta)-(SA12(m,beta))^2;\n HSA := (m,beta) -> - (SA12(m,beta))^2;\n HSA(m,beta);\n lprint(\"SA22 is 0. If SA 2 = 0 so is SA12; so the\");\n lprint(\" Hessian will be 0 at any \+ critical point.\");\n\n #Vofmb := (m,beta) -> (1/(sqrt(m^2+1))^3)*\n #int( (f(x)-m*x-beta)^2*(1+m*fp(x)),x=a..c) -\n #int( (f(x)-m*x-b eta)^2*(1+m*fp(x)),x=c..b);\n #lprint(\"First partial wrt m ...\"); \n #D[1](Vofmb);\n #simplify(%);\n #lprint(\"First partial wrt b eta ...\");\n #D[2](Vofmb);\n #simplify(%);\n #lprint(\"Second p artial wrt m ...\");\n #D[1,1](Vofmb);\n #simplify(%);\n #lprint (\"Second partial wrt beta ...\");\n #D[2,2](Vofmb);\n #simplify(% );\n #lprint(\"________________________________________\");\n\nfi;\n " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 502 "# g := (c1,c2) -> int( -sqrt(1+fp(x)^2),x=a ..c1) -\n# int(-sqrt(1+fp(x)^2),x= c1..c2) +\n# int(-sqrt(1+fp(x)^2),x=c2..b );\n# plot(g( c1,c2),c1=a..b,c2=a..b,grid=[11,11]);\n# g := (c1,c2) ->\n# - ArcLen gth(f(x),x= a..c1)\n# + ArcLength(f(x),x=c1..c2)\n# - ArcLength(f( x),x=c2..b );\n# plot(g(c1,c2),c1=a..b,c2=a..b,grid=[5,5];\n# DoThis : = 'DoThis':\n# DoThis := 0:\n# if (DoThis = 1) then\n# eqn := ArcLe ngth(f(x), x=a..c) = ArcLength(f(x), x=c..b);\n# fsolve(eqn,c);\n# \+ fi;\n" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 5075 "# Start Nvalm lo op.\n\n# Define the number of and values of slopes m\n# for which to d etermine the minimum surface\n# areas.\nprintlevel := 0:\nNvalm := 'Nv alm': Mbegin := 'Mbegin':\nMQuit := 'MQuit': MDelta := 'MDelta':\nNval m := 21;\nim := 'im':\nMBegin := -5:\nMQuit := 5:\nif (Nvalm > 1) then \n MDelta := evalf((MQuit-MBegin)/(Nvalm-1)):\nfi;\nprintlevel := 1: \nfor im from 1 to Nvalm do\n beta := 0:\n # Define the slope for \+ this trip thru the first\n # pass Nvalm loop.\n m := 'm':\n if ( Nvalm = 1) then\n m := 0:\n else\n m := evalf(MBegin+(im-1 )*MDelta):\n fi;\n lprint(\"______________________________________ ________\");\n lprint(\"Start of Nvalm loop calculations for m = \", m);\n if (im = 1) then\n # Plot f(x). \n printlevel := 0: \n p := plot(f(x),x=a..b,y=Fmin..Fmax,\n scaling=c onstrained,color=red,\n title=`f(x)`):\n printleve l := 2:\n display(p);\n fi;\n printlevel := 1:\n (theta,cos t,sint,atan,wofx,zofx,wpofx,zpofx,Pie):=tcs(m):\n\n printlevel := 2: \n lprint(\"Study these plots.\");\n if (m = 0) then\n plot( \{f(x),0*x+m\},x=a..b,\n title=`y = f(x) and y = m`);\n e lse\n plot(\{f(x),0*x+m,0*x-1/m\},x=a..b,\n title=`y = f(x), y = m, and y = -1/m`);\n fi;\n plot([wofx(x),zofx(x),x=a..b ],\n color=blue,scaling=constrained,title=`z(w)`);\n plot([x, zofx(x),x=a..b],\n color=red,scaling=constrained,title=`z(x)`); \n plot([x,wofx(x),x=a..b],\n color=green,scaling=constrained ,title=`w(x)`);\n plot([x,wpofx(x),x=a..b],\n color=red,scali ng=constrained,title=`w'(x)`);\n plot([x,zpofx(x),x=a..b],\n \+ color=red,scaling=constrained,title=`z'(x)`);\n\n printlevel := 0:\n # Find lower and upper bounds zmin and zmax\n # for the parametri c function z(x).\n Ns := 'Ns': zmin := 'zmin': zmax := 'zmax':\n d eltax := 'deltax':\n Ns := 21:\n deltax := evalf((b-a)/(Ns-1)):\n \+ (zmin,zmax) := Estimate_Range(Ns,zofx,a,b):\n printlevel := 1:\n \+ lprint(\"Start Golden Search calculations for m = \", m);\n # Do th e Golden Section minimization for\n # this value of m.\n Smin := ' Smin':\n # (Smin,kminS) := -SMG(zmin,zmax):\n # kminS := valx:\n \+ (Smin,kminS) := SMG(zmin,zmax):\n Smin := -Smin:\n betaS := evalf (kminS*sqrt(m^2+1)+beta):\n\n # Compare betaS to the value determine d by cpivot.\n lprint(\"m, kminS, Smin, betaS = \",\n m, kminS, Smin, betaS);\n betaSc := 'betaSc':\n betaSc := fcpivot - cpivot* m:\n lprint(\"betaSc, betaS, diff. = \",\n betaSc, betaS, betaS c-betaS);\n\n # Plot the parametric equations for f(x)\n # and the line z = k.\n printlevel := 0:\n p1 := plot([wofx(x),zofx(x),x=a. .b],color=blue,\n scaling=constrained,title=`w(x) and z(x )`):\n p2 := plot([wofx(x),kminS+0*x,x=a..b],color=blue,\n \+ title=`w(x) and line z=kminS`,\n scaling=constrained) :\n printlevel := 1:\n display(\{p1,p2\},scaling=constrained,\n \+ title=`z(w) and kminS`);\n\n PlotSk := \"PlotSk\":\n PlotS k := 0: \n if (PlotSk = 1) then\n # Plot the surface area funct ion S(k).\n printlevel := 2:\n lprint(\"Plot of -SA(k) follo ws ...\");\n lprint(\" (It should be unimodal.)\");\n lpri nt(\"Please be patient. This plot takes a while.\");\n plot(-SA,z min..zmax,title=`S(k)`);\n fi;\n\n printlevel := 1:\n # Save the surface area data in order to generate\n # an animated plot of the \+ lines below.\n mdataS[im] := m:\n kdataS[im] := kminS:\n betaSda ta[im] := betaS:\n Smindata[im] := Smin:\n lprint(\"End Golden Sea rch calculations for m = \", m);\n\n # Plot the minimum surface area line and the\n # function f(x).\n plot(\{m*x+betaS,f(x)\},x=a..b, color=[red,blue],\n scaling=constrained,title=`y=mx+betaS and f (x)`);\n\n # Plot the solid with minimum surface area for\n # this slope.\n lprint(\"Parametric plot of surface ...\");\n printlevel := 0:\n Both_Parm :=\n plot3d([wofx(x),(zofx(x)-kminS)*cos(t),(zo fx(x)-kminS)*sin(t)],\n x=a..b,t=0..2*Pi,axes=framed,\n \+ labels=[\"w\",\"z\",\"\"],scaling=CONSTRAINED,\n numpoint s=ngridw*ngridz,\n glossiness=1,lightmodel=light2):\n print level := 1:\n plots[display](Both_Parm);\n printlevel := 0:\n SF irstPassFrame[im] := Both_Parm:\n\n Theta0 := 'Theta0': Theta1 := 'T heta1':\n Theta0 := 0; Theta1 := Pi; \n printlevel := 0:\n Half_ Parm :=\n plot3d([wofx(x),zofx(x)*cos(t),zofx(x)*sin(t)],\n \+ x=a..b,t=Theta0..Theta1,axes=framed,\n labels=[\"w\",\"z\", \"\"],scaling=CONSTRAINED,\n numpoints=ngridw*ngridz,\n \+ glossiness=1,lightmodel=light2):\n printlevel := 1:\n plots[di splay](Half_Parm);\n\n printlevel := 0:\n # Plot the line of revol ution, f(x), the reflection\n # of f(x), and the endpoint perpendicu lars. Several\n # views are plotted. Take your pick.\n p1 := 'p1': p2 := 'p2': p3 := 'p3':\n p4 := 'p4': p5 := 'p5':\n (p1,p2,p3,p4, p5) := reflectf(a,b,f,m,betaS):\n printlevel := 1:\n plots[display ](p1,p2,p3);\n plots[display](p1,p2,p3,p4,p5);\nod; # End of Nvalm l oop\nlprint(\"Finished Nvalm loop.\");\nlprint(\"_____________________ __________________\");" }}{PARA 13 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 807 "printlevel := 2:\n# Animate the mi nimum surface lines.\nShowMoreOfLines := 'ShowMoreOfLines':\nShowMoreO fLines := 0:\n# Plot the lines y = mx+betaS that yield the\n# minimum \+ surface areas.\ni := 'i':\n# for i from 1 to Nvalm do\n# plot(\{mda taS[i]*x+betaSdata[i],f(x)\},\n# x=a..b,y=Fmin..Fmax,color=[mar oon,navy],\n# scaling=constrained,thickness=3);\n# od;\n# i := \+ 'i':\nprintlevel := 0:\nfor i from 1 to Nvalm do\n plot(\{mdataS[i]* x+betaSdata[i],f(x)\},\n x=a..b,y=Fmin..Fmax,color=[maroon,navy ],\n scaling=constrained,thickness=3):\n FrameS[i] := %:\nod: \ni := 'i':\nprintlevel := 2:\nlprint(\"To animate the surface area li nes\");\nlprint(\" y = mx + betaS, click on the\");\nlprint(\" p lot and choose Animation/Play/Slow ...\");\ni := 'i':\ndisplay([seq(Fr ameS[i],i=1..Nvalm)],insequence=true);\n" }{TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 312 "printlevel := 1:\n# Animate the mi nimum surfaces of revolution.\nlprint(\"To animate the first pass mini mum surfaces,\");\nlprint(\" click on the plot and choose\");\nlprin t(\" Animation/Play/Slow ...\");\ni := 'i':\ndisplay([seq(SFirstPass Frame[i],i=1..Nvalm)],\n insequence=true,glossiness=1,lightmode l=light2);\n" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 239 "printlevel := 2:\n# Print a first pass summary.\ni := 'i':\nprintf(\" m \+ betaS kminS Smin\");\nfor i from 1 to Nvalm do\n printf (\"% 10.3f % 10.3f % 10.3f % 10.3f \\n\",\n mdataS[i],betaSdata[i ],kdataS[i],Smindata[i]);\nod;\n" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 2599 "DoThis := 'DoThis':\nDoThis := 0:\nif (DoThis = 1) \+ then\n# Minimization of vertical rotations.\n# As theta -> Pi/2 or -Pi /2, Smin doesn't approach\n# approach the surface area obtained by rev olving\n# f(x) about the y-axis. It approaches the area\n# obtained wh en f(x) is revolved about x = -kmin.\n# Find kmin and Smin as m -> +- \+ infinity.\n# Additionally, verify the results obtained by\n# revolving f(x) about the x-axis and about the\n# y-axis agree with the correspo nding Maple values.\nplot(SfAL,a..b,title=`SfAL(k)`);\nlprint(\"Callin g Sforpion2 ...\");\n(kminSforpion2,Sminforpion2,Sforkis0) := Sforpion 2(a,b):\nlprint(\"kminSforpion2,Sminforpion2, Sforkis0 = \",\n kminS forpion2,Sminforpion2,Sforkis0);\n\nlprint(\"Revolve f(x) about the y- axis ...\");\nSurfaceOfRevolution(f(x),x=a..b,axis=vertical,\n dista ncefromaxis=0,output=plot,\n title=`Revolve f(x) about the y-axis`); \nSurfaceOfRevolution(f(x),x=a..b,axis=vertical,\n distancefromaxis= 0,output=integral,\n title=`Revolve f(x) about the y-axis`);\nSabout yaxis :=\nSurfaceOfRevolution(f(x),x=a..b,axis=vertical,\n distancef romaxis=0,\n title=`Revolve f(x) about the y-axis`);\nSaboutyaxis := evalf(\nint(2*Pie*abs(x)*sqrt(1+(fp(x))^2),x=a...b,digits=3));\nlprin t(\"Check that the vertical rotation matches\");\nlprint(\" the resu lt of revolving about the y-axis.\");\nlprint(\"Saboutyaxis, Sforkis0 \+ difference = \",\n Saboutyaxis, Sforkis0, Saboutyaxis-Sforkis0);\n\n lprint(\"Plot f, its reflection, and vertical axis ...\");\nprintlevel := 0:\nk := 'k':\nk := -kminSforpion2;\nkaxis :=line([k,fmin],[k,fmax ],color=black,thickness=3):\nfplot := plot(f(x),x=a..b,title=``,\n \+ color=blue,thickness=3):\nrflplot := plot(f(2*k-xbar),xbar=2 *k-b..2*k-a,\n title=``,color=red,thickness=3):\nprintlev el := 1:\ndisplay(\{kaxis,fplot,rflplot\});\ndisplay(\{kaxis,fplot\}); \n\nDoThis := 'DoThis':\nDoThis := 0:\nif (DoThis = 1) then\nlprint(\" Revolve f(x) about the x-axis ...\");\nSurfaceOfRevolution(f(x),x=a..b ,axis=horizontal,\n distancefromaxis=0,output=plot,\n title=`Revol ve f(x) about the x-axis`);\nSaboutxaxis :=\nSurfaceOfRevolution(f(x), x=a..b,axis=horizontal,\n distancefromaxis=0,\n title=`Revolve f(x ) about the x-axis`);\nSurfaceOfRevolution(f(x),x=a..b,axis=horizontal ,\n distancefromaxis=0,output=integral,\n title=`Revolve f(x) abou t the x-axis`);\nm := 'm': SA0 := 'SA0':\nm := 0:\n(theta,cost,sint,at an,wofx,zofx,wpofx,zpofx,Pie):=tcs(m):\nSA0 := -SA(0):\nlprint(\"Check that the horizontal rotation matches\");\nlprint(\" the result of r evolving about the x-axis.\");\nlprint(\"Saboutxaxis, SA0, difference \+ = \",\n Saboutxaxis, SA0, Saboutxaxis-SA0);\nfi;\nfi;\n" }{TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 3034 "Plot_Smin := Plot_ Smin:\nPlot_Smin := 0:\nif (Plot_Smin = 1 and Use_Pivot = 1) then\n \+ # Plot approximate Smin as a function of m.\n printlevel := 2:\n N s := 'Ns':\n Ns := 51:\n mfirst := -5:\n mlast := 5:\n deltam \+ := evalf((mlast-mfirst)/(Ns-1)):\n Nz := 'Nz':\n Nz := 21:\n del tax := 'deltax':\n deltax := evalf((b-a)/(Ns-1)):\n i := 'i':\n \+ LeastSminEstimate := 1e10:\n for i from 1 to Ns do\n m := 'm': \n m := evalf(mfirst+(i-1)*deltam):\n (theta,cost,sint,atan, wofx,zofx,wpofx,zpofx,Pie):=tcs(m):\n printlevel := 1:\n bet aS := 'betaS': kminS := 'kminS':\n betaS := fcpivot - cpivot*m:\n kminS := (betaS-beta)/sqrt(m^2+1):\n Smin := -SA(kminS):\n \+ lprint(\"i, m, betaS, kminS, Smin = \",\n i, m, betaS, km inS, Smin):\n mdataS[i] := m:\n kdataS[i] := kminS:\n b etaminS[i] := betaS:\n Sdata[i] := Smin:\n LeastSminEstimate := min(LeastSminEstimate,Smin):\n od;\n printlevel := 0:\n lpri nt(\"Ballpark Least Smin = \", LeastSminEstimate);\n i := 'i':\n k data := [[mdataS[i],kdataS[i]] $i=1..Ns]:\n i := 'i':\n sdata := [ [mdataS[i],Sdata[i]] $i=1..Ns]:\n i := 'i': m := 'm':\n printlevel := 2:\n # plot(kdata,m=mdataS[1]..mdataS[Ns],style=point,\n # \+ symbol=circle,color=blue,title=`kmin vs m`);\n plot(sdata,m=mdataS [1]..mdataS[Ns],style=point,\n symbol=circle,color=blue,title=` Smin vs m`);\nfi;\n\nif (Plot_Smin = 1 and Use_Pivot <> 1) then\n # \+ Plot the actual Smin as a function of m.\n printlevel := 2:\n Ns : = 'Ns':\n Ns := 51:\n mfirst := -5:\n mlast := 5:\n deltam := \+ evalf((mlast-mfirst)/(Ns-1)):\n Nz := 'Nz':\n Nz := 21:\n deltax := 'deltax':\n deltax := evalf((b-a)/(Ns-1)):\n i := 'i':\n Lea stSminEstimate := 1e10:\n for i from 1 to Ns do\n m := 'm':\n \+ m := evalf(mfirst+(i-1)*deltam):\n (theta,cost,sint,atan,wofx ,zofx,wpofx,zpofx,Pie):=tcs(m):\n printlevel := 1:\n (zmin,z max) := Estimate_Range(Nz,zofx,a,b):\n\n # Perform a Golden Searc h for this value of m.\n kminS := SK(m):\n # Smin := S(kminS ):\n Smin := -SA(kminS):\n betaS := evalf(kminS*sqrt(m^2+1)+ beta):\n # Compare betaS to the value determined by cpivot.\n \+ lprint(\"i, m, betaS, kminS, Smin = \",\n i, m, betaS, kminS , Smin):\n betaSc := 'betaSc':\n betaSc := fcpivot - cpivot* m;\n lprint(\"betaSc, betaS, diff. = \",\n betaSc, betaS, betaSc-betaS);\n mdataS[i] := m:\n kdataS[i] := kminS:\n \+ betaminS[i] := evalf(kminS*sqrt(m^2+1)+beta):\n Sdata[i] := Sm in:\n LeastSminEstimate := min(LeastSminEstimate,Smin):\n od;\n printlevel := 0:\n lprint(\"Ballpark Least Smin = \", LeastSminEs timate);\n i := 'i':\n kdata := [[mdataS[i],kdataS[i]] $i=1..Ns]: \n i := 'i':\n sdata := [[mdataS[i],Sdata[i]] $i=1..Ns]:\n i := \+ 'i': m := 'm':\n printlevel := 2:\n # plot(kdata,m=mdataS[1]..mdat aS[Ns],style=point,\n # symbol=circle,color=blue,title=`kmin vs m`);\n plot(sdata,m=mdataS[1]..mdataS[Ns],style=point,\n sym bol=circle,color=blue,title=`Smin vs m`);\nfi;" }}{PARA 13 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 3775 "# Perform a Golden Search on Smin(m) to determine\n# the minimum value of Smin. \nprintlevel := 1:\nlprint(\"Please be patient. This minimization take s\"):\nlprint(\" a long time unless Use_Pivot = 1.\"):\n# Refine the interval on which the outer Golden\n# Search minimization will be per formed.\nmminS := 'mminS': mmaxS := 'mmaxS':\nMinForVertical := 'MinFo rVertical':\nMinForVertical := 0:\nif (Problem = 1) then\n mminS := \+ 0.5: mmaxS := 1.5:\nelif (Problem = 2) then\n mminS := 0.5: mmaxS := 1.5:\nelif (Problem = 3) then\n mminS := 0.5: mmaxS := 1.5:\nelif ( Problem = 4) then\n mminS := -1.0: mmaxS := 0.0:\nelif (Problem = 5) then\n mminS := -0.5: mmaxS := 0.5:\nelif (Problem = 6) then\n mm inS := -0.5: mmaxS := 0.5:\nelif (Problem = 7) then\n mminS := -0.5: mmaxS := 0.5:\nelif (Problem = 8) then\n mminS := 2.0: mmaxS := 3.0 :\nelif (Problem = 9) then\n mminS := -3.0: mmaxS := -2.0:\nelif (Pr oblem = 10) then\n mminS := 5.0: mmaxS := 5.0:\n MinForVertical := 1:\nelif (Problem = 11) then\n mminS := 1.5: mmaxS := 2.5:\nelif (P roblem = 12) then\n mminS := 0.0: mmaxS := 1.0:\nelif (Problem = 13) then\n mminS := 0.0: mmaxS := 1.0:\nelif (Problem = 14) then\n mm inS := 1.0: mmaxS := 2.0:\nelif (Problem = 15) then\n mminS := -2.5: mmaxS := -1.5:\nelif (Problem = 16) then\n mminS := 5.0: mmaxS := 5 .0:\n MinForVertical := 1:\nelif (Problem = 17) then\n mminS := -0 .5: mmaxS := 0.5:\nelif (Problem = 18) then\n mminS := -1.5: mmaxS : = 0.0:\nelif (Problem = 19) then\n mminS := -0.5: mmaxS := 0.5:\neli f (Problem = 20) then\n mminS := 0.1: mmaxS := 1.0:\nelif (Problem = 21) then\n mminS := 5.0: mmaxS := 5.0:\n MinForVertical := 1:\nel if (Problem = 22) then\n mminS := 5.0: mmaxS := 5.0:\n MinForVerti cal := 1:\nelif (Problem = 23) then\n mminS := 5.0: mmaxS := 5.0:\n \+ MinForVertical := 1:\nelif (Problem = 24) then\n mminS := 5.0: mma xS := 5.0:\n MinForVertical := 1:\nelif (Problem = 25) then\n mmin S := 0.5: mmaxS := 1.5:\nelif (Problem = 26) then\n mminS := -0.5: m maxS := 0.5:\nelif (Problem = 27) then\n mminS := -1.5: mmaxS := -0. 5:\nelif (Problem = 28) then\n mminS := 1.0: mmaxS := 2.0:\nelif (Pr oblem = 29) then\n mminS := 0.0: mmaxS := 1.0:\nelif (Problem = 30) \+ then\n mminS := 0.0: mmaxS := 1.0:\nelif (Problem = 31) then\n mmi nS := 1.5: mmaxS := 2.5:\nelif (Problem = 32) then\n mminS := -0.5: \+ mmaxS := 0.5:\nelif (Problem = 33) then\n mminS := 0.0: mmaxS := 1.0 :\nelif (Problem = 34) then\n mminS := 5.0: mmaxS := 5.0:\n MinFor Vertical := 1:\nelif (Problem = 35) then\n mminS := -4.5: mmaxS := - 3.5:\nelif (Problem = 36) then\n mminS := 5.0: mmaxS := 5.0:\n Min ForVertical := 1:\nelif (Problem = 37) then\n mminS := -0.5: mmaxS : = 0.5:\nelif (Problem = 38) then\n mminS := 0.0: mmaxS := 1.0:\nelif (Problem = 39) then\n mminS := -0.5: mmaxS := 0.5:\nelif (Problem = 40) then\n mminS := -0.5: mmaxS := 0.5:\nelif (Problem = 41) then\n mminS := -0.5: mmaxS := 0.5:\nelif (Problem = 42) then\n mminS := 5.0: mmaxS := 5.0:\n MinForVertical := 1:\nelif (Problem = 43) then \n mminS := -0.5: mmaxS := 0.5:\nelse\n mminS := -5: mmaxS := 5:\n fi;\n\nif (MinForVertical <> 1) then\n # Do the Golden Section minim ization for this\n # value of m.\n Sminmin := 'Sminmin': mminminS \+ := 'mminminS':\n bminminS := 'bminminS':\n # (Sminmin,mminminS) := -SMG2(mminS,mmaxS):\n # (Sminmin,mminminS) := -SMG2b(mminS,mmaxS): \n # mminminS := valx2:\n if (Use_Pivot = 1) then\n # Approxi mate minimization of Smin.\n (Sminmin,mminminS) := SMG2b(mminS,mm axS):\n else\n # Full minimization of Smin.\n (Sminmin,mmi nminS) := SMG2(mminS,mmaxS):\n fi;\n Sminmin := -Sminmin;\n bmin minS := evalf(kminminS*sqrt(mminminS^2+1)+beta):\n lprint(\"Smallest Smin, slope, and intercept found = \",\n Sminmin, mminminS, bmin minS);\n lprint(\"kminminS = \", kminminS);\nfi;\n " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 1277 "# Plot the axis of revolution, f( x), the reflection\n# of f(x) through the axis, and the endpoint\n# pe rpendiculars from the axis to f(x).\nprintlevel := 0:\n(p1,p2,p3,p4,p5 ) := reflectf(a,b,f,mminminS,bminminS):\nprintlevel := 1:\nplots[displ ay](p1,p2,p3);\nplots[display](p1,p2,p3,p4,p5);\n\nlprint(\"Parametric plot of surface ...\");\nprintlevel := 0:\n(theta,cost,sint,atan,wofx ,zofx,wpofx,zpofx,Pie)\n:= tcs(mminminS):\nBoth_Parm :=\nplot3d([wofx( x),(zofx(x)-kminminS)*cos(t),(zofx(x)-kminminS)*sin(t)],\n x=a.. b,t=0..2*Pi,axes=framed,\n labels=[\"w\",\"z\",\"\"],scaling=CON STRAINED,\n numpoints=ngridw*ngridz,glossiness=1,\n lightm odel=light2):\nprintlevel := 1:\nplots[display](Both_Parm);\n\n# Plot \+ the minimal axis of revolution and f(x).\nlprint(\"First reflection pl ot ...\");\nplot(\{mminminS*x+bminminS,f(x)\},\n x=a..b,color=[red ,blue],\n scaling=constrained,\n title=`Minimal y=mx+betaS, an d f(x)`);\nlprint(\"Second reflection plot ...\");\nplot(\{mminminS*x+ bminminS,f(x)\},\n x=a..b,y=Fmin..Fmax,color=[red,blue],\n sca ling=constrained,\n title=`Minimal y=mx+betaS and f(x)`);\n\n# Plo t z(w) for the solid with minimum surface area.\nx := 'x':\nplot([wofx (x),zofx(x),x=a..b],\n color=blue,scaling=constrained,\n title =`Parametric equations (w(x),z(x))`);\n" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 460 "Surface_Area_Plot_wz := Surface_Area_Plot_wz:\nSurfa ce_Area_Plot_wz := 0:\nif (Surface_Area_Plot_wz = 1) then\n printlev el := 2:\n lprint(\"Surface area 3d plot using (w,z) coordinates ... \");\n lprint(\"Please be patient. This plot takes a long time.\"); \n # Plot the surface area as a function of m and beta.\n m := 'm' : beta := 'beta':\n plot3d(SEval2,-10..10,-5..5,grid=[51,51],\n \+ title=`S(m,beta)`,axes=boxed,\n labels=[\"w\",\"z\",\"S \"]);\nfi;\n" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 402 "Surface_Ar ea_Plot_xy := Surface_Area_Plot_xy:\nSurface_Area_Plot_xy := 0:\nif (S urface_Area_Plot_xy = 1) then\n printlevel := 2:\n lprint(\"Surfac e area 3d plot using (x,y) coordinates ...\");\n lprint(\"Please be \+ patient. This plot takes a long time.\");\n m := 'm': beta := 'beta' :\n plot3d(Smxy,-10..10,-5..5,\n grid=[51,51],title=`S(m,be ta)`,axes=boxed,\n labels=[\"w\",\"z\",\"S\"]);\nfi;\n" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}}{MARK "0 0 0" 0 } {VIEWOPTS 1 1 0 1 1 1803 1 1 1 1 }{PAGENUMBERS 0 1 2 33 1 1 }