{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 34343 "# Investigate the minimum surface areas obtained\n# when the graph of a parametric curv e is revolved\n# about lines parallel to lines y = mx + beta.\n\n# fil ename = ParametricRevslant.mws\n\n# Notes:\n# The run times are quite \+ long when the full\n# minimizations are done. If you only want\n# to s ee a few surfaces, you can bypass the\n# full minimization near the en d of the\n# worksheet. You can also reduce the number\n# (Nvalm) of sl opes used in the first pass\n# calculations. (I typically use somethin g\n# like Nvalm = 3.) As things stand, you can\n# use any one of 18 pr oblems (my favorite\n# is number 18, the crooked egg). If you wish\n# \+ to define your own equation, search on\n# 'Problem :=' and add your fu nction at the\n# at the end of the block. If you have trouble\n# repro ducing any of the plots in the manuscript,\n# I have problem specific \+ versions of this\n# worksheet in which I may have changed the\n# amoun t of detail to get satisfactory plots.\n# Please feel free to make any changes you like\n# and use the worksheet as you see fit. If you\n# f ind any blunders or better ways to do things,\n# I'll get a kick out o f seeing your changes. \n\nrestart:\n\nSMG2b := proc(mminS,mmaxS)\n# T his is the function procedure for the Golden\n# Search minimization of Smin. SMG2b is called\n# by the main procedure.\n local tol2,tol2_d ummy,zmax2,zmin2,Smin2,kminS2:\n global valx2,valf2,iters2,Use_Adapt mw2,\n From_Where:\n From_Where := 3:\n (tol2,tol2_dummy) := Set _Tols(From_Where):\n lprint(\"SMG2b is calling GOLD2 ...\");\n # U se Adaptmw to do the integration.\n (valf2,valx2,iters2) := GOLD2(SM G3b,mminS,mmaxS,tol2):\n lprint(\"Back in SMG2b after call to GOLD2 \+ ...\");\n lprint(\"Slope m yielding minimum Smin = \", valx2);\n l print(\"Minimum surface area S(valx2) for this m = \", Smin);\n retu rn(valf2,valx2):\nend proc:\n\nSMG3b := proc(mval)\n# This is the func tion procedure for the Golden 2\n# Search minimization based on the pi voting\n# information. SMG3b is called by GOLD2.\n local tol,theta,c ost,sint,atan,wofx,zofx,\n wpofx,zpofx,Pie,betaSc,SminSc,k,m,\n Fr om_Where,abserrA,relerrA,SvalA,errestA,\n flagA,nfeA:\n global kmi nminS,fA2b:\n m := evalf(mval):\n (theta,cost,sint,atan,wofx,zofx, wpofx,zpofx,Pie):=tcs(m):\n # Define the surface area integrand for \+ Adaptmw.\n fA2b := 'fA2b':\n fA2b := x -> 2*Pie*evalf(\n abs(zof x(x)-k)*sqrt((wpofx(x))^2+(zpofx(x))^2)):\n # Define the error toler ances.\n From_Where := 6:\n (abserrA,relerrA) := Set_Tols(From_Whe re):\n flagA := -1:\n # Approximate the integral of the function f A2b(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(\"App roximate 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 fA2b evals in A daptmw\";\n fi;\n fi;\n SvalA := -SvalA:\n lprint(\"In Smg3 b. m, beta, SvalA = \", m, beta, SvalA):\n return(SvalA):\nend proc: \n\nSMG := proc(zmin,zmax)\n# This is the function procedure for the G olden\n# Search minimization of kmin. SMG is called by\n# the main pro cedure and proc SK.\n local tol,tol_dummy,From_Where:\n global val x,valf,iters:\n # lprint(\"In SMG. zmin, zmax = \", zmin, zmax);\n \+ From_Where := 1:\n (tol,tol_dummy) := Set_Tols(From_Where):\n lpr int(\"SMG is calling GOLD ...\");\n (valf,valx,iters) := GOLD(SA,zmi n,zmax,tol):\n # lprint(\"valx = \", valx):\n lprint(\"Minimum sur face area S(valx) = \", abs(valf)):\n lprint(\"Number of Golden Sear ch iterations = \", iters):\n return(valf,valx):\nend proc:\n\nSA := proc(k)\n# Calculate the surface area and multiply by -1.\n# This is \+ the function procedure for the Golden\n# Search minimization if Adaptm w is used to do\n# the integrations. SA is called by SMG. m is a\n# gl obal variable defined in the main procedure.\n local flagA,abserrA,r elerrA,answer,errestA,nfeA,\n x,Sval,SvalA,theta,cost,sint,atan,wofx ,zofx,\n wpofx,zpofx,Pie,From_Where:\n global fA:\n # lprint(\"I n 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 the integral of the function fA(x).\n (SvalA,erre stA,flagA,nfeA) :=\n Adaptmw(fA,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(\"Integrand evaluations = \", nfeA);\n # Check if Adapt mw was successful.\n if (flagA = 1) then\n error \"Insuff icient queue storage in Adaptmw\";\n fi;\n if (flagA = 2) th en\n error \"Too many fA evals in Adaptmw\";\n fi;\n fi ;\n Sval := -SvalA:\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 local Kval,atan,theta,cost,sint,Pie,\n Smin,kminS,wo fx,zofx,wpofx,zpofx:\n (theta,cost,sint,atan,wofx,zofx,wpofx,zpofx,P ie):=tcs(m):\n # Do the Golden Section minimization for\n # this v alue of m.\n # Note:\n # SK assumes that zmin and zmaz were calcul ated\n # before callin SK.\n lprint(\"In SK. m, zmin, zmax = \", m , zmin, zmax);\n Smin := 'Smin':\n (Smin,kminS) := SMG(zmin,zmax): \n Smin := -Smin:\n lprint(\"In SK. m, kminS, S(kminS) = \", m, km inS, Smin);\n Kval := kminS:\n return(Kval):\nend proc:\n\nSMG2 := proc(mminS,mmaxS)\n# This is the function procedure for the Golden\n# Search minimization of Smin. 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 (t ol2,tol2_dummy) := Set_Tols(From_Where):\n lprint(\"SMG2 is calling \+ GOLD2 ...\");\n # Use Adaptmw to do the integration.\n (valf2,valx 2,iters2) := GOLD2(SMG3,mminS,mmaxS,tol2):\n lprint(\"Back in SMG2 a fter 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(\"Number of Golden 2 Search iterations = \",\n iters2);\n return(valf2,valx2):\nend proc:\n\nSMG3 := pro c(m)\n# This is the function procedure for the Golden 2\n# Search mini mization. SMG3 is called by GOLD2.\n local tol,theta,cost,sint,atan, wofx,zofx,\n wpofx,zpofx,Pie,zmin,zmax,Nzm,zminm,zmaxm:\n global v alx3,valf3,iters3,Use_Adaptmw,mg2,\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 := 'N zm':\n Nzm := 21:\n (zminm,zmaxm) := Estimate_Range(Nzm,zofx,a,b): \n tol := 'tol':\n tol := evalf(1/10^4):\n Use_Adaptmw := 1:\n \+ lprint(\"SMG3 is calling GOLD ...\");\n (valf3,valx3,iters3) :=\n \+ GOLD(SA2,zminm,zmaxm,tol):\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 slope = \", m, valx3); \n lprint(\"Minimum surface area S(valx3) for this m = \", abs(valf3 ));\n lprint(\"Number of Golden Search iterations for this m = \", i ters3);\n return(valf3):\nend proc:\n\nSA2 := 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. 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,errestA,nfeA,\n x,Sval,SvalA,thet a,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 integrand for Adaptmw.\n fA2 := 'fA2' :\n fA2 := x -> 2*Pie*evalf(\n abs(zofx(x)-k)*sqrt((wpofx(x))^2+(z pofx(x))^2)):\n # Define the error tolerances.\n From_Where := 6: \n (abserrA,relerrA) := Set_Tols(From_Where):\n flagA := -1:\n # Approximate the integral of the function fA2(x).\n (SvalA,errestA,f lagA,nfeA) :=\n Adaptmw(fA2,a,b,abserrA,relerrA):\n if (flagA <> - 1) then\n # Print the results.\n # lprint(\"Status flagA \+ = \", flagA);\n # lprint(\"Approximate integral = \", Sval A);\n # lprint(\"Approximate error = \", errestA);\n # l print(\"Integrand evaluations = \", nfeA);\n # Check if Adaptmw w as successful.\n if (flagA = 1) then\n error \"Insufficie nt queue storage in Adaptmw\";\n fi;\n if (flagA = 2) then\n error \"Too many fA evals in Adaptmw\";\n fi;\n fi;\n \+ Sval := -SvalA:\n return(Sval):\nend proc:\n\nEstimate_Range := pr oc(n,g,a,b)\n# Estimate the range 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 deltax := evalf((b-a)/(n-1)):\n i := \+ 'i': \n for i from 1 to n do\n xpt := 'xpt': gptx := 'gpt x':\n xpt := evalf(a+(i-1)*deltax):\n gptx := evalf(g( xpt)):\n gmin := min(gptx,gmin):\n gmax := max(gptx,gm ax):\n od:\n i := 'i':\n elif (n = 1) then\n gpta := \+ 'gpta': gptb := 'gptb': \n gpta := evalf(g(a)):\n gptb := ev alf(g(b)):\n gmin := min(gpta,gptb):\n gmax := max(gpta,gptb ): \n else\n error \"Estimate_Range requires more than points. \";\n fi:\n gmin := gmin - deltax:\n gmax := gmax + deltax:\n \+ return(gmin,gmax):\nend proc:\n\nSet_Tols := proc(From_Where)\n# Set A daptmw integration error tolerances and\n# Golden search error toleran ces.\n local abstol,reltol:\n if (From_Where < 1 or From_Where > 1 1) then\n error \"Illegal tolerance 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_W here = 2) then\n # Outer Golden Search to minimize Smin:\n a bstol := 1e-2:\n reltol := 1e-2: # not used\n elif (From_Where \+ = 3) then\n # Surface area 3d plot:\n abstol := 1e-5:\n \+ reltol := 1e-5:\n else\n # Other Adaptmw integrations:\n \+ abstol := 1e-3:#5:\n reltol := 1e-3:#5:\n fi:\n # lprint(\"In Set_Tols. From_Where, abstol, reltol = \",\n # From_Where, absto l, reltol);\n return(abstol,reltol):\nend proc:\n\ntcs := proc(m)\n# Define the necessary transformation of\n# axes information.\n local cost,sint,atan,wofx,zofx,wpofx,\n zpofx,theta,Pie:\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 := t -> cost*xoft(t)+sint*(yoft(t)-beta):\n zofx := t -> -sint*xoft(t )+cost*(yoft(t)-beta):\n wpofx := t -> cost*xpoft(t)+sint*ypoft(t) :\n zpofx := t -> -sint*xpoft(t)+cost*ypoft(t):\n return(theta,co st,sint,atan,wofx,zofx,wpofx,zpofx,Pie):\nend proc:\n\n#______________ _____________________________________________\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 adaptive quadrature scheme based on\n# Gauss-Kronrod (3,7) formulas. \n# \n# Usage:\n# (answer,errest,flag,nfe) := adaptm w(f,a,b,abserr,relerr):\n# Input parameters:\n# f = name of th e function defining f(x).\n# a, b = end points of integration int erval.\n# abserr = absolute error tolerance desired.\n# relerr = relative error tolerance desired.\n# Output parameters:\n# answer \+ = computed estimate of the integral.\n# errest = estimate of the ab solute error in answer.\n# flag = 0 for normal return.\n# \+ = 1 insufficient storage in queue (120).\n# = 2 too ma ny function evaluations (3577).\n#\n# Adaptmw was adapted from the Mat lab script Adapt.m\n# from the book Fundamentals of Numerical Computin g\n# by Shampine, Allen, and Preuss.\n#\n# IMPORTANT NOTE:\n# If you e xerience difficulties with Adaptmw or have\n# questions about it, plea se do not bug the authors\n# of Numerical Computing. Instead, contact \+ Skip\n# Thompson by email at thompson@radford.edu.\n# \nAdaptmw := pro c(f,aval,bval,abserrval,relerrval)\n#\nlocal maxq, maxfe, eps, length, top, bottom,\nanswer, errest, alpha, beta, h, tol, el, er,\nql, qr, q , e, nfe, flag, Quadmw, Addmw, queue,\na, b, abserr, relerr, A, x, Alp ha, Beta, Q, E,\nQkonrod, H, f1, f2, f3, f4, f5, f6, f7, center:\n#\n# Unload the arguments.\n#\na := aval:\nb := bval:\nrelerr := relerrval :\nabserr := abserrval:\n#\n# Set up the necessary arrays.\n#\nmaxq := 120:\n# queue := zeros(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.10465622602646 72:\nA[3] := 0.4013974147759622:\nA[4] := 0.4509165386584744:\nx[1] := 0.7745966692414834:\nx[2] := 0.9604912687080202:\nx[3] := 0.434243749 3468026:\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 Adaptmw.\";\nfi:\nif abserr < 0 then\n error \"abs err <= 0 is not allowed in Adaptmw.\";\nfi:\n#\n# Initialization.\n#\n length := 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 suff iciently\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 := evalf(f(center - H * x[3])):\nf7 := evalf(f(center + H \+ * x[3])):\nQ := evalf(H * (5 * (f2 + f3) + 8 * f1) / 9):\nQkonrod := e valf(H * (A[1] * (f2 + 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 := e valf(e):\n# \nif evalf(abs(errest)) >\n evalf(max(abserr,relerr*abs( answer))) then\n # Addmw(maxq,q,e,a,b,length,bottom):\n queue[bott om,1] := q:\n queue[bottom,2] := e:\n queue[bottom,3] := a:\n qu eue[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 retu rn,\n# else subdivide the top entry.\n#\nwhile length > 0 do\n #\n \+ # Remove the top entries from the queue.\n #\n q := evalf(queue[t op,1]):\n e := evalf(queue[top,2]):\n alpha := evalf(queue[top,3]) :\n beta := evalf(queue[top,4]):\n length := length - 1:\n if to p < 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,alp ha,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(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 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(cen ter)):\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 := 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 qr := Q:\n er := E:\n nf e := nfe + 7:\n #\n # Update answer and the error estimate.\n # \n answer := evalf(answer + ((ql + qr) - q)):\n errest := evalf(er rest + ((el + er) - e)):\n #\n # Test for failures.\n #\n if l ength >= maxq - 1 then\n flag := 1:\n return(answer,errest,f lag,nfe):\n fi:\n if nfe >= maxfe then\n flag := 2:\n re turn(answer,errest,flag,nfe):\n fi:\n #\n # Test for convergence .\n #\n tol := evalf(max(abserr,relerr*abs(answer))):\n if evalf (abs(errest)) <= tol then\n flag := 0:\n return(answer,erres t,flag,nfe):\n fi:\n #\n # Add new subintervals to queue if erro rs\n # are too big.\n #\n tol := evalf(tol * h / (b - a)):\n i f evalf(abs(el)) > tol then\n # Addmw(maxq,ql,el,alpha,alpha+h,l ength,bottom):\n queue[bottom,1] := ql:\n queue[bottom,2] \+ := el:\n queue[bottom,3] := alpha:\n queue[bottom,4] := al pha + h:\n length := 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)) > t ol then\n # Addmw(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 lengt h := length + 1:\n bottom := bottom:\n if bottom < maxq th en\n bottom := bottom + 1:\n else\n bottom := \+ 1:\n fi:\n fi:\nod: # End of while loop\nflag := 0:\nreturn(an swer,errest,flag,nfe):\nend proc:\n\nGOLD := proc(f::procedure,a::nume ric,b::numeric,T::numeric)\n#\n# To perform a the Golden Search maximi zation for a\n# unimodal function do the following.\n#\n# 1) Provide \+ a procedure to evaluate f(x).\n# 2) Type (fmin,xmin,iters) := GOLD(f, a,b,tolerance):\n# for specific values of a, b, and the toleran ce.\n# 3) The last output provided is the midpoint of the\n# f inal interval and the value of f(x) at that\n# point.\n#\n# Not e: To minimize a unimodal function, maximize -f.\n#\n# This Golden Sea rch Algorithm was written originally by\n# Dr. William P. Fox, Depa rtment of Mathematics,\n# Francis Marion University, Florence, S C 29501\n# Dr. Margie Witherspoon, Department of Computer Science, \n# Francis Marion University, Florence, SC 29501\n#\n# This pro cedure performs the Golden Section Search algorithm\n# to find the max imum of a unimodal function, f(x), over an\n# interval, a < x < b. The program calculates the number of\n# iterations required to insure the final interval is within\n# the user-specified tolerance. This is fou nd by solving for\n# the smallest value of k that makes this inequalit y true:\n# (b-a)(0.618^k)< tolerance. The Golden section concept\n# in volves placing two experiments between [a,b] using the\n# Golden secti on ratios. One experiment is placed at position,\n# a + 0.382(b-a), an d the other at position, a + 0.618(b-a).\n# The function, to be maximi zed, is evaluated at these two\n# points and the functional values are compared. We want to\n# keep the larger functional value and its corr esponding\n# opposite endpoint. At the end of the required iterations, \n# the final interval is the answer. At times when the final\n# answe r must be a single point and not an interval, the\n# convention of sel ecting the midpoint is provided.\n#\n# The Golden Search ratios con6 a nd con3 are defined below\n# in the following manner. At each iteratio n two new tests\n# points are added. They are chosen such that the fol lowing\n# two ratios are equal:\n#\n# Length of current interval Le ngth of larger fraction\n# -------------------------- = ------------- ------------ \n# Length of larger fraction Length of smaller fract ion\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 the 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 = co n3 ~ 0.382.\n#\n local x1,x2, N, val;\n global con6, con3;\n N := 0;\n #\n # Define the Golden Search ratios.\n con6 := evalf((sqrt(5)-1)/2) ;\n con3 := evalf(1-con6);\n # \n # Determine the first two endpoints \+ x1 and x2.\n x1 := a+con3*(b-a);\n x2 := a+con6*(b-a);\n # \n # Determ ine the number of iterations to be performed.\n N := ceil((ln(T/(b-a)) /ln(con6)));\n # lprint(\"In GOLD, no. of iterates =\", N);\n # Perfo rm the maximization.\n iterate(f,a,b,N,x1,x2);\n #\n # The value f(mid point) is found. This midpoint is the\n # midpoint of the interval fou nd in the Nth iteration.\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 global mdpt,fkeep,xkeep,pos;\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)+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 interval a midpoint is found. After finding \n \+ # the midpoint, an error check is performed to ensure the\n # midpo int of the interval within the specified tolerance\n # is the max v alue for the given function.\n # \n if (j = N) then\n m dpt := (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(b n(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 := (an(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# Gauss-Kronrod (3,7) formulas. \n# \n# Usage:\n# ( answer,errest,flag,nfe) := Adaptmw2(f,a,b,abserr,relerr):\n# Input par ameters:\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 tolerance desired.\n # Output parameters:\n# answer = computed estimate of the integral. \n# errest = estimate of the absolute error in answer.\n# flag \+ = 0 for normal return.\n# = 1 insufficient storage in que ue (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 Numerical 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 Nu merical Computing. Instead, contact Skip\n# Thompson by email at thomp son@radford.edu.\n# \nAdaptmw2 := proc(f,aval,bval,abserrval,relerrval )\n#\nlocal maxq, maxfe, eps, length, top, bottom,\nanswer, errest, al pha, 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, f 1, f2, f3, f4, f5, f6, f7, center:\n#\n# Unload the arguments.\n#\na : = aval:\nb := bval:\nrelerr := relerrval:\nabserr := abserrval:\n#\n# \+ Set up the necessary arrays.\n#\nmaxq := 120:\n# queue := zeros(maxq,4 ):\nqueue := array(1..maxq,1..4):\nmaxfe := 3577:\n# a = [0.2684880898 683334, 0.1046562260264672, ...\n# 0.4013974147759622, 0.45091653 86584744]:\n# x = [0.7745966692414834, 0.9604912687080202,\n# 0.4 342437493468026]:\nA := array(1..4):\nx := array(1..3):\nA[1] := 0.268 4880898683334:\nA[2] := 0.1046562260264672:\nA[3] := 0.401397414775962 2:\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 Adaptmw2.\"; \nfi:\nif abserr < 0 then\n error \"abserr <= 0 is not allowed in Ad aptmw2.\";\nfi:\n#\n# Initialization.\n#\nlength := 0:\ntop := 1:\nbot tom := 1:\nnfe := 0:\n#\n# Form an initial approximation answer to the \n# integral over [a,b]. If it is not sufficiently\n# accurate, initia lize the queue and begin the\n# main loop.\n# q,e) := Quadmw2(f,a,b): \n#\nAlpha := a:\nBeta := b:\nH := evalf(0.5 * (Beta - Alpha)):\ncente r := evalf(Alpha + H):\nf1 := evalf(f(center)):\nf2 := evalf(f(center \+ - H * x[1])):\nf3 := evalf(f(center + H * x[1])):\nf4 := evalf(f(cente r - H * x[2])):\nf5 := evalf(f(center + H * x[2])):\nf6 := evalf(f(cen ter - H * x[3])):\nf7 := evalf(f(center + H * x[3])):\nQ := evalf(H * \+ (5 * (f2 + f3) + 8 * f1) / 9):\nQkonrod2 := evalf(H * (A[1] * (f2 + f3 ) +\n A[4] * f1 + A[2] * (f4 + f5)\n + A[3] * (f6 + f7))):\nE := evalf(Qkonrod2 - Q):\nq := Q:\ne := E:\nnfe := nfe + 7:\nanswer := evalf(q):\nerrest := evalf(e):\n# \nif evalf(a bs(errest)) >\n evalf(max(abserr,relerr*abs(answer))) then\n # Add mw2(maxq,q,e,a,b,length,bottom):\n queue[bottom,1] := q:\n queue[b ottom,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:\nf i:\n#\n# Main loop: if queue is empty then return,\n# else subdivide t he top entry.\n#\nwhile length > 0 do\n #\n # Remove the top entri es from the queue.\n #\n q := evalf(queue[top,1]):\n e := evalf( queue[top,2]):\n alpha := evalf(queue[top,3]):\n beta := evalf(que ue[top,4]):\n length := length - 1:\n if top < maxq then\n to p := top + 1:\n else\n top := 1:\n fi:\n h := evalf(0.5*(be ta-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(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 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 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(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 rod2 := 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 := eva lf(answer + ((ql + qr) - q)):\n errest := evalf(errest + ((el + er) \+ - e)):\n #\n # Test for failures.\n #\n if length >= maxq - 1 \+ then\n flag := 1:\n return(answer,errest,flag,nfe):\n fi: \n if nfe >= maxfe then\n flag := 2:\n return(answer,erres t,flag,nfe):\n fi:\n #\n # Test for convergence.\n #\n tol : = evalf(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)) \+ > tol then\n # Addmw2(maxq,ql,el,alpha,alpha+h,length,bottom):\n queue[bottom,1] := ql:\n queue[bottom,2] := el:\n q ueue[bottom,3] := alpha:\n queue[bottom,4] := alpha + h:\n \+ length := length + 1:\n bottom := bottom:\n if bottom < m axq then\n bottom := bottom + 1:\n else\n bott om := 1:\n fi:\n fi:\n if evalf(abs(er)) > tol then\n \+ # Addmw2(maxq,qr,er,alpha+h,beta,length,bottom):\n queue[bottom, 1] := qr:\n queue[bottom,2] := er:\n queue[bottom,3] := al pha + h:\n queue[bottom,4] := beta:\n length := length + 1 :\n bottom := bottom:\n if bottom < maxq then\n b ottom := bottom + 1:\n else\n bottom := 1:\n fi: \n fi:\nod: # End of while loop\nflag := 0:\nreturn(answer,errest,fl ag,nfe):\nend proc:\n\n# Copy of GOLD\nGOLD2 := proc(f::procedure,a::n umeric,b::numeric,T::numeric)\n#\n# To perform a the Golden Search max imization for a\n# unimodal function do the following.\n#\n# 1) Provi de a procedure to evaluate f(x).\n# 2) Type (fmin,xmin,iters) := GOLD 2(f,a,b,tolerance):\n# for specific values of a, b, and the tol erance.\n# 3) The last output provided is the midpoint of the\n# \+ final interval and the value of f(x) at that\n# point.\n#\n# Note: To minimize a unimodal function, maximize -f.\n#\n# This Golden Search Algorithm was written originally by\n# Dr. William P. Fox, \+ Department of Mathematics,\n# Francis Marion University, Florenc e, SC 29501\n# Dr. Margie Witherspoon, Department of Computer Scien ce,\n# Francis Marion University, Florence, SC 29501\n#\n# This \+ procedure performs the Golden Section Search algorithm\n# to find the \+ maximum of a unimodal function, f(x), over an\n# interval, a < x < b. \+ The program calculates the number of\n# iterations required to insure \+ the final interval is within\n# the user-specified tolerance. This is \+ found by solving for\n# the smallest value of k that makes this inequa lity true:\n# (b-a)(0.618^k)< tolerance. The Golden section concept\n# involves placing two experiments between [a,b] using the\n# Golden se ction ratios. One experiment is placed at position,\n# a + 0.382(b-a), and the other at position, a + 0.618(b-a).\n# The function, to be max imized, is evaluated at these two\n# points and the functional values \+ are compared. We want to\n# keep the larger functional value and its c orresponding\n# opposite endpoint. At the end of the required iteratio ns,\n# the final interval is the answer. At times when the final\n# an swer must be a single point and not an interval, the\n# convention of \+ selecting the midpoint is provided.\n#\n# The Golden Search ratios con 62 and con32 are defined below\n# in the following manner. At each ite ration two new tests\n# points are added. They are chosen such that th e following\n# two ratios are equal:\n#\n# Length of current interval Length of larger fraction\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 to the equation r^2+r-1=0. The\n# positive solut ion 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 the Golden Search ratios.\n con62 := evalf((sqrt (5)-1)/2);\n con32 := evalf(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 number 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 maximization.\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 the Nth iteration.\n val := f(mdpt2);\n return(v al,mdpt2,N);\n end proc:\n \n iterate2 := proc(f::procedure,a::numeri c,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 endpoint is assigned to x1n(1) and \n # the second en dpoint 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 following 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 specified.\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)+con62*(b n(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 fou nd. After finding \n # the midpoint, an error check is performed to ensure the\n # midpoint of the interval within the specified toler ance\n # is the max value for the given function.\n # \n i f (j = N) then\n mdpt2 := (an(j) + bn(j))/2;\n if (f(an(j) ) > f(bn(j))) and (f(an(j)) > f(mdpt2)) then \n fkeep2 := f(a n(j));\n xkeep2 := 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 xke ep2 := 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" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 3733 "# Be gin main procedure.\nwith(plots):\nwith(Student[Calculus1]):\nwith(plo ttools):\n\n# Define the numerical precision.\nDigits := 16:\n\n# Defi ne the number of grid points for the\n# plots of the surfaces with min imum area.\nngridw := 'ngridw': ngridz := 'ngridz':\nngridw := 51: ngr idz := 51:\n\nt := 't':\nProblem := 1;\nif (Problem = 1) then\n a := -0.5;\n b := 1;\n xoft := t -> t;\n yoft := t -> t^2;\nelif (Pr oblem = 2) then\n a := 0;\n #b := 2*evalf(Pi);\n b := evalf(Pi); \n # b := evalf(Pi);\n # b := evalf(Pi)/2;\n xoft := t -> sin(2* t);\n yoft := t -> 2*sin(t);\nelif (Problem = 3) then\n a := -5;\n b := 0.1;\n # b := 0.2;\n xoft := t -> cos(t+sin(3*t));\n yof t := t -> sin(t+cos(4*t));\nelif (Problem = 4) then\n a := 0;\n b \+ := 5;\n xoft := t -> t+cos(2*t);\n yoft := t -> t-sin(4*t);\nelif \+ (Problem = 5) then\n # Barbed wire:\n a := -5;\n b := 5;\n xof t := t -> t+cos(2*t);\n yoft := t -> t-sin(4*t);\nelif (Problem = 6) then\n # Pose:\n # Wei-Chi's Example 10\n # Be sure to change o uter minimization\n # interval below.\n a := 0;\n # b := evalf(P i/2); # Case a: 1.5 leafs\n # b := evalf(2*Pi/3); # Case b: 2 le afs\n # b := evalf(Pi); # Case c. 3 leafs\n # b := evalf(Pi/ 3); # Case d: 1 leaf\n b := .4679647278306630; # Case e: 0.5 l eaf\n xoft := t -> sin(3*t)*cos(t);\n yoft := t -> sin(3*t)*sin(t) ;\nelif (Problem = 7) then\n # Wei-Chi's Example 11\n a := 0;\n \+ pie := evalf(Pi);\n b := pie;\n foft := t -> sin(3*t)*cos(t);\n \+ goft := t -> sin(3*t)*sin(t);\n xoft := t -> foft(t-pie/6);\n yoft := t -> goft(t+pie/6);\nelif (Problem = 8) then\n a := -5;\n b := 0.1;\n # b := 0.2;\n xoft := t -> cos(t+sin(5*t));\n yoft := t \+ -> sin(t+cos(5*t));\nelif (Problem = 9) then\n a := -2;\n b := 2; \n xoft := t -> t^2-2;\n yoft := t -> t^3-t;\nelif (Problem = 10) \+ then\n # Folium:\n a := -0.5;\n b := 10;\n xoft := t -> 3*t/(1 +t^3);\n yoft := t -> 3*t^2/(1+t^3);\nelif (Problem = 11) then\n # Spiral:\n a := 0;\n b := 10;\n pie := evalf(Pi);\n xoft := t \+ -> t*cos(t)/(100*pie);\n yoft := t -> t*sin(t)/(100*pie);\nelif (Pro blem = 12) then\n # Spiral:\n a := 0;\n b := 20;\n pie := eval f(Pi);\n xoft := t -> t*cos(t)/(100*pie);\n yoft := t -> t*sin(t)/ (100*pie);\nelif (Problem = 13) then\n # Spiral:\n a := 0;\n b : = 30;\n pie := evalf(Pi);\n xoft := t -> t*cos(t)/(100*pie);\n y oft := t -> t*sin(t)/(100*pie);\nelif (Problem = 14) then\n # Rose: \n a := 0;\n b := 2*evalf(Pi);\n # b := evalf(Pi);\n acon := 1 ; bcon := 2; kcon := 4;\n xoft := t -> (acon+bcon*cos(kcon*t))*cos(t );\n yoft := t -> (acon+bcon*cos(kcon*t))*sin(t);\nelif (Problem = 1 5) then\n # Rose:\n a := 0;\n #b := 2*evalf(Pi);\n b := evalf( Pi);\n acon := 2; bcon := 1; kcon := 4;\n xoft := t -> (acon+bcon* cos(kcon*t))*cos(t);\n yoft := t -> (acon+bcon*cos(kcon*t))*sin(t); \nelif (Problem = 16) then\n # Ellipse:\n a := 0;\n b := 2*evalf (Pi);\n acon := 2;\n bcon := 5;\n xoft := t -> acon*cos(t);\n \+ yoft := t -> acon*sin(t);\nelif (Problem = 17) then\n # Upper right \+ branch of hyperbola:\n a := 0;\n b := 2;\n acon := 1;\n bcon : = 1;\n xoft := t -> acon*(exp(t)+exp(-t))/2;\n yoft := t -> bcon*( exp(t)-exp(-t))/2;\nelif (Problem = 18) then\n # Crooked egg:\n a \+ := 0;\n b := evalf(Pi);#3;\n xoft := t -> ((sin(t))^3 + (cos(t))^3 ) * cos(t);\n yoft := t -> ((sin(t))^3 + (cos(t))^3) * sin(t);\nelse \n # Supply other function information here.\n error \"No problem \+ number was specified.\";\nfi;\na := evalf(a):\nb := evalf(b):\n\nprint level := 1:\n\n# Define x'(t) and y'(t).\nxpoft := t -> D(xoft)(t);\ny poft := t -> D(yoft)(t);\nxpoft(t);\nypoft(t);\n\n# Plot (x(t),y(t)). \nplot([xoft(t),yoft(t),t=a..b],title=`(x(t),y(t))`,\n scaling=con strained,thickness=3);\n# Plot (x'(t),y'(t)).\n" }}}{EXCHG {PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 672 "# Extimate the ranges of xoft and yoft(t).\nn xy := 'nxy':\nxtmin := 'xtmin': xtmax := 'xtmax':\nytmin := 'ytmin': y tmax := 'ytmax':\ndeltaxy := 'deltaxy':\nnxy := 101:\nxtmin := 1e20: \nxtmax := -1e20:\nytmin := 1e20:\nytmax := -1e20:\ndeltaxy := evalf( (b-a)/(nxy-1)):\ni := 'i':\nfor i from 1 to nxy do\n txy := 'txy':\n txy := evalf(a+(i-1)*deltaxy):\n xtmin := min(xtmin,xoft(txy)):\n xtmax := max(xtmax,xoft(txy)):\n ytmin := min(ytmin,yoft(txy)):\n ytmax := max(ytmax,yoft(txy)):\nod:\ni := 'i':\nxtmin := xtmin-delt axy:\nxtmax := xtmax+deltaxy:\nytmin := ytmin-deltaxy:\nytmax := ytmax +deltaxy:\nlprint(\"xtmin, xtmax = \", xtmin, xtmax);\nlprint(\"ytmin, ytmax = \", ytmin, ytmax);\n" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 5309 "# Start Nvalm loop.\n# Define the number of and values of slop es m\n# for which to determine the minimum surface\n# areas.\nprintlev el := 0:\nNvalm := 'Nvalm': Mbegin := 'Mbegin':\nMQuit := 'MQuit': MDe lta := 'MDelta':\nNvalm := 21:#3;\nim := 'im':\nMBegin := -5:\nMQuit : = 5:\nif (Nvalm > 1) then\n MDelta := evalf((MQuit-MBegin)/(Nvalm-1) ):\nfi;\n\nfor im from 1 to Nvalm do\n beta := 0:\n # Define the s lope 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(MBe gin+(im-1)*MDelta):\n fi:\n printlevel := 1:\n lprint(\"Start of Nvalm loop calculations for m = \", m);\n printlevel := 0:\n (the ta,cost,sint,atan,wofx,zofx,wpofx,zpofx,Pie):=tcs(m):\n # Find lower and upper bounds zmin and zmax\n # for the parametric function z(x) .\n Ns := 'Ns': zmin := 'zmin': zmax := 'zmax':\n deltax := 'delta x':\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 the Golden Sectio n minimization for\n # this value of m.\n Smin := 'Smin':\n (Smi n,kminS) := SMG(zmin,zmax):\n Smin := -Smin:\n betaS := evalf(kmin S*sqrt(m^2+1)+beta):\n PlotSk := \"PlotSk\":\n PlotSk := 0: \n i f (PlotSk = 1) then\n # Plot the surface area function S(k).\n \+ printlevel := 2:\n lprint(\"Plot of -SA(k) follows ...\");\n \+ lprint(\" (It should be unimodal.)\");\n lprint(\"Please be patient. This plot takes a while.\");\n plot(-SA,zmin..zmax,titl e=`S(k)`);\n fi;\n # Save the surface area data in order to genera te\n # an animated plot of the lines below.\n mdataS[im] := m:\n \+ kdataS[im] := kminS:\n betaSdata[im] := betaS:\n Smindata[im] := \+ Smin:\n lprint(\"End Golden Search calculations for m = \", m);\n \+ t := 't':\n\n # Plot the axis of revolution.\n # plot(mdataS[im]*t +betaSdata[im],\n # t=a..b,color=[maroon,navy],\n # scal ing=constrained,thickness=3);\n printlevel := 0:\n pline :=\n pl ot(mdataS[im]*t+betaSdata[im],\n t=xtmin..xtmax,y=ytmin..ytmax, \n color=[maroon,navy],\n scaling=constrained,thickness= 3):\n # Plot (x(t),y(t)).\n pcurve :=\n plot([xoft(t),yoft(t),t= a..b],scaling=constrained,\n thickness=3):\n # Display both t he axis of revolution and the curve.\n printlevel := 1:\n display( pcurve,pline);\n\n # Plot the surface of revolution.\n printlevel \+ := 0:\n Both_Parm :=\n plot3d([wofx(t),(zofx(t)-kminS)*cos(th),(zo fx(t)-kminS)*sin(th)],\n t=a..b,th=0..2*Pi,axes=framed,\n \+ labels=[\"w\",\"z\",\"\"],scaling=CONSTRAINED,\n numpoi nts=ngridw*ngridz,\n glossiness=1,lightmodel=light2):\n pri ntlevel := 1:\n lprint(\"Full solid ...\");\n plots[display](Both_ Parm);\n printlevel := 0:\n SFirstPassFrame[im] := Both_Parm:\n \+ printlevel := 0:\n Half_Parm :=\n plot3d([wofx(t),(zofx(t)-kminS)* cos(th),(zofx(t)-kminS)*sin(th)],\n t=a..b,th=0..Pi,axes=fram ed,\n labels=[\"w\",\"z\",\"\"],scaling=CONSTRAINED,\n \+ numpoints=ngridw*ngridz,\n glossiness=1,lightmodel=light2) :\n printlevel := 1:\n lprint(\"Half solid ...\");\n plots[displ ay](Half_Parm);\n printlevel := 0:\n SFirstPassHalfFrame[im] := Ha lf_Parm:\n printlevel := 0:\n Quarter_Parm :=\n plot3d([wofx(t), (zofx(t)-kminS)*cos(th),(zofx(t)-kminS)*sin(th)],\n t=a..b,th =0..Pi/2,axes=framed,\n labels=[\"w\",\"z\",\"\"],scaling=CON STRAINED,\n numpoints=ngridw*ngridz,\n glossiness=1 ,lightmodel=light2):\n printlevel := 1:\n lprint(\"Quarter solid . ..\");\n plots[display](Quarter_Parm);\n\n DoThis := 'DoThis':\n \+ DoThis := 0:\n if (DoThis = 1) then\n printlevel := 2:\n ndevel op := 'ndevelop': thinc := 'thinc':\n ndevelop := 6;\n thinc := ev alf(2*Pi/ndevelop);\n idevelop := 'idevelop': thb := 'thb':\n for \+ idevelop from 1 to ndevelop do\n thb := evalf(idevelop*thinc);\n \+ lprint(\"idevelop, thb := \", idevelop, thb);\n plot3d([wofx (t),(zofx(t)-kminS)*cos(th),(zofx(t)-kminS)*sin(th)],\n t= a..b,th=0..thb,axes=framed,\n labels=[\"w\",\"z\",\"\"],sc aling=CONSTRAINED,\n numpoints=ngridw*ngridz,\n \+ glossiness=1,lightmodel=light2);\n # printlevel := 1:\n # \+ plots[display](Quarter_Parm);\n od;\n idevelop := 'idevelop':\n \+ printlevel := 1:\n ndevelop := 'ndevelop': thinc := 'thinc':\n nde velop := 6;\n thinc := evalf(2*Pi/ndevelop);\n idevelop := 'idevel op': thb := 'thb':\n for idevelop from 1 to ndevelop do\n thb : = evalf(idevelop*thinc);\n lprint(\"idevelop, thb := \", idevelop , thb);\n pdevelop[idevelop] :=\n plot3d([wofx(t),(zofx(t)-k minS)*cos(th),(zofx(t)-kminS)*sin(th)],\n t=a..b,th=0..thb ,axes=framed,\n labels=[\"w\",\"z\",\"\"],scaling=CONSTRAI NED,\n numpoints=ngridw*ngridz,\n glossiness=1 ,lightmodel=light2):\n # printlevel := 1:\n # plots[display] (Quarter_Parm);\n od;\n idevelop := 'idevelop':\n printlevel := \+ 1:\n # Animate the sections of the surface of revolution.\n lprint (\"To animate the first pass minimum surfaces,\");\n lprint(\" cli ck on the plot and choose\");\n lprint(\" Animation/Play/Slow ... \");\n i := 'i':\n display([seq(pdevelop[i],i=1..ndevelop)],\n i nsequence=true,glossiness=1,lightmodel=light2);\n fi;\nod;\n# End of Nvalm loop.\n" }{TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 965 "printlevel := 2:\n# Animate the minimum surface lines.\nt := 't ':\ni := 'i':\n#for i from 1 to Nvalm do\n# l := t -> mdataS[i]*t+be taSdata[i];\n# plot([t,l(t)],\n# t=a..b,scaling=constrained,thickn ess=3);\n#od;\nt := 't':\nprintlevel := 0:\npcurve :=\nplot([xoft(t),y oft(t),t=a..b],thickness=3,\n title=`(x(t),y(t)) and y=mx+beta_S`, \n scaling=constrained):\nprintlevel := 2:\nt := 't':\nfor i from \+ 1 to Nvalm do\n l := t -> mdataS[i]*t+betaSdata[i]:\n # l(t);\n \+ # pline := plot([t,l(t)],\n # pline := plot(l(t),t=a..b,color=[maroo n,navy],\n pline :=\n plot(l(t),t=xtmin..xtmax,y=ytmin..ytmax,\n \+ color=[maroon,navy],\n scaling=constrained,thickness=3): \n printlevel := 2:\n FrameS[i] := display(pcurve,pline):\nod:\ni \+ := 'i':\nprintlevel := 2:\nlprint(\"To animate the surface area lines \");\nlprint(\" y = mx + betaS, click on the\");\nlprint(\" plot and choose Animation/Play/Slow ...\");\ni := 'i':\ndisplay([seq(Frame S[i],i=1..Nvalm)],insequence=true);" }{TEXT -1 0 "" }}{PARA 13 "" 0 " " {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 304 "printlev el := 1:\n# Animate the minimum surfaces of revolution.\nlprint(\"To a nimate the first pass minimum surfaces,\");\nlprint(\" click on the \+ plot and choose\");\nlprint(\" Animation/Play/Slow ...\");\ni := 'i' :\ndisplay([seq(SFirstPassFrame[i],i=1..Nvalm)],\ninsequence=true,glos siness=1,lightmodel=light2);\n" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 296 "# Animate the minimum half surfaces of revolution.\nlprint(\" To animate the first pass minimum surfaces,\");\nlprint(\" click on \+ the plot and choose\");\nlprint(\" Animation/Play/Slow ...\");\ni := 'i':\ndisplay([seq(SFirstPassHalfFrame[i],i=1..Nvalm)],\ninsequence=t rue,glossiness=1,lightmodel=light2);\n" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 236 "printlevel := 2:\n# Print a first pass summary.\ni : = 'i':\nprintf(\" m betaS kminS Smin\");\nfor i f rom 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 1390 "DoThis := 'DoThis':\nDoThis := 0: \nif (DoThis = 1) then\n# Plot Smin as a function of m.\nprintlevel := 2:\nNs := 'Ns':\nNs := 51:\nmfirst := -5:\nmlast := 5:\ndeltam := eva lf((mlast-mfirst)/(Ns-1)):\nNz := 'Nz':\nNz := 21:\ndeltax := 'deltax' :\ndeltax := evalf((b-a)/(Ns-1)):\ni := 'i':\nLeastSminEstimate := 1e1 0:\nfor i from 1 to Ns do\n m := 'm':\n m := evalf(mfirst+(i-1)*de ltam):\n (theta,cost,sint,atan,wofx,zofx,wpofx,zpofx,Pie):=tcs(m):\n printlevel := 1:\n (zmin,zmax) := Estimate_Range(Nz,zofx,a,b):\n \+ # Perform a Golden Search for this value of m.\n kminS := SK(m):\n # Smin := S(kminS):\n Smin := -SA(kminS):\n betaS := evalf(kmin S*sqrt(m^2+1)+beta):\n # Compare betaS to the value determined by cp ivot.\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-be taS);\n mdataS[i] := m:\n kdataS[i] := kminS:\n betaminS[i] := e valf(kminS*sqrt(m^2+1)+beta):\n Sdata[i] := Smin:\n LeastSminEstim ate := min(LeastSminEstimate,Smin):\nod;\nprintlevel := 0:\nlprint(\"B allpark Least Smin = \", LeastSminEstimate);\ni := 'i':\nkdata := [[md ataS[i],kdataS[i]] $i=1..Ns]:\ni := 'i':\nsdata := [[mdataS[i],Sdata[i ]] $i=1..Ns]:\ni := 'i': m := 'm':\nprintlevel := 2:\nplot(sdata,m=mda taS[1]..mdataS[Ns],style=point,\nsymbol=circle,color=blue,title=`Smin \+ vs m`);\nfi;\n" }{TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 1418 "# Perform a Golden Search on Smin(m) to determine\n# the minim um value of Smin.\nprintlevel := 1:\nlprint(\"Please be patient. This \+ minimization takes\"):\nlprint(\" a long time.\"):\n# Refine the int erval on which the outer Golden\n# Search minimization will be perform ed.\nmminS := 'mminS': mmaxS := 'mmaxS':\nif (Problem = 1) then\n mm inS := 0.5: mmaxS := 1.5:\nelif (Problem = 6) then\n # mminS := 1.0 : mmaxS := 2.0: # Case a: 1.5 leafs\n # mminS := 1.5: mmaxS := 2.5 : # Case b: 2 leafs\n # mminS := -0.5: mmaxS := 0.5: # Case c: 3 leafs\n # mminS := 0: mmaxS := 1: # Case d: 1 leaf\n \+ mminS := 0: mmaxS := 1: # Case e: 0.5 leaf\nelif (Problem = 11) \+ then\n mminS := -1.5: mmaxS := 1.5:\nelif (Problem = 12) then\n \+ mminS := -1.5: mmaxS := 1.5:\nelif (Problem = 14) then\n mminS : = -0.5: mmaxS := 0.5:\nelif (Problem = 15) then\n mminS := -0.5: m maxS := 0.5:\nelif (Problem = 16) then\n mminS := -0.5: mmaxS := 0 .5:\nelif (Problem = 17) then\n mminS := 0.5: mmaxS := 1.5:\nelif (Problem = 18) then\n mminS := -1.5: mmaxS := -0.5:\nelse\n mmi nS := -5: mmaxS := 5:\nfi;\n# Do the Golden Section minimization. \nSminmin := 'Sminmin': mminminS := 'mminminS':\nbminminS := 'bminminS ':\n(Sminmin,mminminS) := SMG2(mminS,mmaxS):\nSminmin := -Sminmin;\nbm inminS := evalf(kminminS*sqrt(mminminS^2+1)+beta):\nlprint(\"Smallest \+ Smin, slope, and intercept found = \",\nSminmin, mminminS, bminminS); \n " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 772 "# Plot the minimu m surface and half surface of revolution.\nlprint(\"Parametric plot of surface ...\");\nprintlevel := 0:\n(theta,cost,sint,atan,wofx,zofx,wp ofx,zpofx,Pie)\n:= tcs(mminminS):\nBoth_Parm :=\nplot3d([wofx(t),zofx( t)*cos(th),zofx(t)*sin(th)],\n t=a..b,th=0..2*Pi,axes=framed,\n \+ labels=[\"w\",\"z\",\"\"],scaling=CONSTRAINED,\n numpoints =ngridw*ngridz,\n glossiness=1,lightmodel=light2):\nprintlevel : = 1:\nplots[display](Both_Parm);\nlprint(\"Parametric plot of surface \+ ...\");\nprintlevel := 0:\nHalf_Parm :=\nplot3d([wofx(t),zofx(t)*cos(t h),zofx(t)*sin(th)],\n t=a..b,th=0..Pi,axes=framed,\n labe ls=[\"w\",\"z\",\"\"],scaling=CONSTRAINED,\n numpoints=ngridw*ng ridz,\n glossiness=1,lightmodel=light2):\nprintlevel := 1:\nplot s[display](Half_Parm);\n" }}}}{MARK "11" 0 }{VIEWOPTS 1 1 0 1 1 1803 1 1 1 1 }{PAGENUMBERS 0 1 2 33 1 1 }