#if kernelopts(version)='kernelopts(version)' then # ERROR(`This package requires Release 4`) fi; `evalf/INTERVAL`:= proc() map(evalf,INTERVAL(args)) end; gmax := proc (f) options `Maple Advisor Database 1.00 for Maple V Release 4 and 5`, `Copyright (c) 1998 by Robert B. Israel. All rights reserved`; -gmin(-f,args[2..nargs]) end; gmin := proc(f::algebraic, r::(name = realcons .. realcons), m::name, b0::realcons) local x, a, b, fn, fnp, fnpp, fa, fb, bsf, msf, res, fl, n,i, ai, bi, bb, mb, f0; global `gmin/eps1`, `gmin/eps2`, `gmin/eps3`; options `Maple Advisor Database 1.00 for Maple V Release 4 and 5`, `Copyright (c) 1998 by Robert B. Israel. All rights reserved`; Digits:= Digits+5; x := op(1, r); a := evalf(op(1, op(2, r))); b := evalf(op(2, op(2, r))); if a > b then ERROR(`invalid interval `, [a, b]) fi; if has(f,{Heaviside,abs,signum,piecewise,min,max}) then fl:= convert(f,pwlist,x); if has(fl,{Heaviside,abs,signum,piecewise,min,max}) then ERROR(`could not convert to piecewise list`) fi; n:= (nops(fl)+1)/2; bsf:= infinity; msf:= {}; Digits:= Digits-5; for i from 1 to n do if i=1 then ai:= a else ai:= max(a,evalf(fl[2*i-2],Digits+5)) fi; if i=n then bi:= b else bi:= min(b,evalf(fl[2*i],Digits+5)) fi; if ai > bi then next fi; bb:= gmin(fl[2*i-1], x=ai .. bi,'mb',bsf); res:= `gmin/newmin`(mb,bb,bsf,msf); bsf:= res[1]; msf:= res[2]; od; if nargs >= 3 then m:= evalf(msf) fi; RETURN(evalf(bsf)); fi; if has(r, infinity) then bsf:=`gmin/infinite`(f,r,msf); if nargs >= 3 then m:= msf fi; RETURN(bsf) fi; fn := unapply(f, x); fa := traperror(evalf(fn(a))); if fa = lasterror then fa:= evalf(limit(f,x=a,right)) fi; if not type(fa,realcons) then ERROR(`Can not find`,Limit(f,x=a,right)) fi; fb := traperror(evalf(fn(b))); if fb = lasterror then fb:= evalf(limit(f,x=b,left)) fi; if not type(fb,realcons) then ERROR(`Can not find`,Limit(f,x=b,left)) fi; if has([fa,fb],infinity) then `gmin/eps2`:= Float(1,5-Digits) else `gmin/eps2` := Float(1, 5-Digits)*(abs(fa)+abs(fb)+.01)/2 fi; fnp := unapply(diff(f, x), x); fnpp := unapply(diff(f, x $ 2), x); if fa < fb -`gmin/eps2` then bsf:= fa; msf:= {a} elif fa < fb + `gmin/eps2` then bsf:= min(fa,fb); msf:= {a,b} else bsf:= fb; msf:= {b} fi; if (a < 0) and (b > 0) then # check 0 too f0:= evalf(fn(0)); bb:= `gmin/newmin`({0},f0,bsf,msf); bsf:= bb[1]; msf:= bb[2]; fi; if (nargs = 4) and (b0 < bsf) then bsf:= b0; msf:= {} fi; `gmin/eps1`:= Float(1,5-Digits)*(abs(b)+abs(a)); if b = a then `gmin/eps3`:= 1 else `gmin/eps3`:= `gmin/eps2`/(b-a) fi; if bsf = -infinity then res:= [bsf,msf] else res:= `gmin/gmin`(fn, fnp, fnpp, a, b, fa, fb,bsf,msf) fi; Digits:= Digits-5; if nargs >= 3 then m:= evalf(res[2]) fi; evalf(res[1]); end; `gmin/infinite` := proc(f::algebraic, r::(name = realcons .. realcons), m::name) options `Maple Advisor Database 1.00 for Maple V Release 4 and 5`, `Copyright (c) 1998 by Robert B. Israel. All rights reserved`; local x,a,b,res,b1,m1,b2,m2,b3,m3,v,c; ## infinite interval x := op(1, r); a := evalf(op(1, op(2, r))); b := evalf(op(2, op(2, r))); if a >= b then ERROR(`invalid interval `, [a, b]) fi; Digits:= Digits-5; if a = -infinity then if b = infinity then # two-sided infinite # on (-infinity, -10] transform v = 1/x, v = -.1 .. 0 b2:= gmin(f, x = -10 .. 10,m2); b1:= gmin(subs(x=1/v, f), v = -.1 .. 0, m1,b2); m1:= map(proc(t) if t = 0 then -infinity else 1/t fi end, m1); b3:= gmin(subs(x=1/v, f), v = 0 .. 0.1, m3,min(b1,b2)); m3:= map(proc(t) if t = 0 then infinity else 1/t fi end, m3); res:= `gmin/newmin`(m3,b3,op(`gmin/newmin`(m1,b1,b2,m2))); else # infinite on left c:= -10-abs(b); b1:= gmin(subs(x=1/v, f), v = 1/c .. 0, m1); m1:= map(proc(t) if t = 0 then -infinity else 1/t fi end, m1); b2:= gmin(f, x = c .. b, m2,b1); res:= `gmin/newmin`(m1,b1,b2,m2); fi else # infinite on right c:= 10+abs(a); b1:= gmin(subs(x=1/v,f), v = 0 .. 1/c, m1); m1:= map(proc(t) if t = 0 then infinity else 1/t fi end, m1); b2:= gmin(f, x = a .. c, m2,b1); res:= `gmin/newmin`(m1,b1,b2,m2); fi; m:= res[2]; res[1]; end; `gmin/newmin`:= proc(a,fa,bsf,msf) ## add new candidate to best so far options `Maple Advisor Database 1.00 for Maple V Release 4 and 5`, `Copyright (c) 1998 by Robert B. Israel. All rights reserved`; local m,b, sa; if type(a,set) then sa:= a else sa:= {a} fi; if fa < bsf - `gmin/eps2` then # new best replaces previous if sa = {} then userinfo(5,gmin,`new best `,fa) else userinfo(5,gmin,`new best `,fa,` at `,a) fi; [fa, sa] elif fa < bsf + `gmin/eps2` then # new best tied with previous b:= min(bsf,fa); if sa <> {} then userinfo(5,gmin,`best `,fa,` at `,a); sa:= msf union sa; else sa:= msf fi; [ b, sa] else # retain old [bsf, msf] fi end; `gmin/gmin`:= proc(fn, fnp, fnpp, a, b, fa, fb, bsf, msf) ## return list of minimum value, set of points where achieved ## so far know best value bsf, attained at msf ## taking a and b into account options `Maple Advisor Database 1.00 for Maple V Release 4 and 5`, `Copyright (c) 1998 by Robert B. Israel. All rights reserved`; local bb, mb, fr, f1, f2, fp1, fp2, fa1, fb1, c, fc, fpp1, fpp2, a1, b1, fpa, fpb, q, w; userinfo(4, gmin, `Checking interval `, a .. b, `end values`, fa, fb); if b-a < `gmin/eps1` then RETURN(`gmin/newmin`(a,fa,bsf,msf)) fi; # interval reduced to a point q:= fn(INTERVAL(a .. b)); fr := traperror(evalf(readlib(evalr)(q))); if fr = lasterror then fr:= traperror(evalf(evalr(normal(q)))) fi; if fr = lasterror then assume(w>=a,w<=b); if signum(0,fn(w) - bsf,0) = 1 then ## cut off RETURN([bsf,msf]) else fr:= INTERVAL(-infinity..infinity) fi fi; if type(fr, {specfunc(range, INTERVAL), list(range)}) then f1 := op(1, op(1,fr)); f2:= op(2, op(-1,fr)); elif type(fr, realcons) then f1 := fr; f2:= fr; else ERROR(`evalr returned `, fr) fi; userinfo(5,gmin,`range on`,[a,b],` is `,[f1,f2]); ## now range of f on [a,b] is in [f1, f2] if f1 > bsf then ## cut off RETURN([bsf,msf]) elif f2 <= f1 then ## flat on interval RETURN(`gmin/newmin`(INTERVAL(a .. b), min(fa,fb),bsf,msf)) fi; ## now not flat or cut off q:= fnp(INTERVAL(a .. b)); fr:= traperror(evalf(evalr(q))); ## interval for derivative if fr = lasterror then fr:= traperror(evalf(evalr(normal(q)))) fi; if fr = lasterror then assume(w > a, w < b); q:= fnp(w); if signum(0,q,1) = 1 then fr:= INTERVAL(0 .. infinity) elif signum(0,q,-1) = -1 then fr:= INTERVAL(-infinity .. 0) else fr:= INTERVAL(-infinity .. infinity) fi fi; if type(fr, {specfunc(range, INTERVAL), list(range)}) then fp1 := op(1, op(1,fr)); fp2 := op(2, op(-1,fr)) elif type(fr, realcons) then fp1 := fr; fp2 := fr else ERROR(`evalr returned `, fr) fi; if -`gmin/eps3` <= fp1 then ## function is nondecreasing if fp2 <= `gmin/eps3` then ## function is flat RETURN(`gmin/newmin`(INTERVAL(a..b),min(fa,fb),bsf,msf)) else RETURN(`gmin/newmin`(a,fa,bsf,msf)) fi elif `gmin/eps3` >= fp2 then ## function is nonincreasing if fp1 >= -`gmin/eps3` then ## function is flat RETURN(`gmin/newmin`(INTERVAL(a..b),min(fa,fb),bsf,msf)) else RETURN(`gmin/newmin`(b,fb,bsf,msf)) fi fi; ## now not monotone, check for convex/concave q:= fnpp(INTERVAL(a .. b)); fr := traperror(evalf(evalr(q))); if fr = lasterror then fr:= traperror(evalf(evalr(normal(q)))) fi; if fr = lasterror then assume(w > a, w < b); q:= fnpp(w); if signum(0,q,1) = 1 then fr:= INTERVAL(0 .. infinity) elif signum(0,q,-1) = -1 then fr:= INTERVAL(-infinity .. 0) else fr:= INTERVAL(-infinity .. infinity) fi fi; if type(fr, {specfunc(range, INTERVAL),list(range)}) then fpp1 := op(1, op(1,fr)); fpp2 := op(2, op(-1,fr)) elif type(fr, realcons) then fpp1 := fr; fpp2 := fr fi; if fpp1 >= 0 then ## function is convex RETURN(`gmin/convex`(args, evalf(fnp(a)), evalf(fnp(b)))) elif fpp2 <= 0 then ## function is concave if fa <= fb then ## min at a if fb = fa then ## also b RETURN(`gmin/newmin`({a,b},min(fa,fb),bsf,msf)) else ## only a RETURN(`gmin/newmin`(a,fa,bsf,msf)) fi else ## only b RETURN(`gmin/newmin`(b,fb,bsf,msf)) fi; fi; ## not monotone, convex or concave fpa:= traperror(evalf(fnp(a))); if fpa = lasterror then fpa:= evalf(limit(fnp(x),x=a,right)) fi; fpb:= traperror(evalf(fnp(b))); if fpb = lasterror then fpb:= evalf(limit(fnp(x),x=b,left)) fi; if fp1 <> -infinity then a1:= a + (bsf - fa)/fp1; # f(x) >= f(a) + (x-a) fp1 >= bsf else a1:= a fi; # note fp1 < 0 and bsf <= fa if fpp1 <> -infinity then a1:= max(a1, a - (fpa + sqrt(fpa^2 - 2 * fpp1*(fa - bsf)))/fpp1); fi; # f(x) >= fa + (x-a)*fpa + (x-a)^2 fpp1/2 > bsf # note fpp1 < 0 if (fpa < 0) and (fpp2 <> infinity) then a1:= max(a1, a - fpa/fpp2*(1-Float(1,5-Digits))) fi; # f'(x) <= fpa + (x-a) fpp2 < 0 for x-a < -fpa/fpp2 # note fpp2 > 0 if fp2 <> infinity then b1:= b + (bsf - fb)/fp2 else b1:= b fi; # f(x) >= f(b) + (x-b) fp2 >= bsf # note fp2 > 0 if fpp1 <> -infinity then b1:= min(b1, b - (fpb - sqrt(fpb^2 - 2 * fpp1*(fb - bsf)))/fpp1) fi; if (fpb > 0) and (fpp2 <> infinity) then b1:= min(b1, b - fpb/fpp2*(1-Float(1,5-Digits))) fi; ## check if cut off if b1 < a1 - `gmin/eps1` then RETURN([bsf,msf]) elif b1 < a1 then a1:= (a1+b1)/2; if a1 = a then fa1:= fa else fa1:= evalf(fn(a1)) fi; RETURN(`gmin/newmin`(a1,fa1,bsf,msf)) fi; ## not cut off, calculate values at new points if a1 = a then fa1:= fa else fa1:= evalf(fn(a1)) fi; if b1 = b then fb1:= fb else fb1:= evalf(fn(b1)) fi; bb:= `gmin/newmin`({},fb1,op(`gmin/newmin`({},fa1,bsf,msf))); ## should we split? if 2*(b1 - a1) < b-a then ## interval cut down enough, don't split RETURN(`gmin/gmin`(fn, fnp, fnpp, a1, b1, fa1, fb1, op(bb))); else ## split in two c:= (a1+b1)/2; fc:= evalf(fn(c)); bb:= `gmin/newmin`({},fc,op(bb)); ## now evaluate most promising side first if fa1 < fb1 then bb:= `gmin/gmin`(fn, fnp, fnpp, a1, c, fa1, fc, op(bb)); RETURN(`gmin/gmin`(fn, fnp, fnpp, c, b1, fc, fb1, op(bb))) else bb:= `gmin/gmin`(fn, fnp, fnpp, c, b1, fc, fb1, op(bb)); RETURN(`gmin/gmin`(fn, fnp, fnpp, a1, c, fa1, fc, op(bb))) fi; fi; end; `gmin/convex`:= proc(fn, fnp, fnpp, a, b, fa, fb, bsf, msf, fpa, fpb) ## function known to be convex on [a,b] options `Maple Advisor Database 1.00 for Maple V Release 4 and 5`, `Copyright (c) 1998 by Robert B. Israel. All rights reserved`; local bb, mb, fr, f1, f2, fp1, fp2, fa1, fb1, c, fc, fpc, a1, b1, d, fd, fpd, fppd, fpa1, fpb1; userinfo(4, gmin, `Convex on interval `, a .. b, `end values`, fa, fb); if b-a < `gmin/eps1` then RETURN(`gmin/newmin`(a,fa,bsf,msf)) fi; # interval reduced to a point if fpa >= -`gmin/eps3` then # function is nondecreasing if fpb > `gmin/eps3` then # not flat RETURN(`gmin/newmin`(a,fa,bsf,msf)) elif fa > bsf then RETURN([bsf,msf]) ## flat and above min else RETURN([bsf,msf minus {a,b} union {INTERVAL(a..b)}]); # flat at min fi elif fpb <= `gmin/eps3` then # function is nonincreasing if fpa < -`gmin/eps3` then # not flat RETURN(`gmin/newmin`(b,fb,bsf,msf)) elif fb > bsf then RETURN([bsf,msf]) ## flat and above min else RETURN([bsf,msf minus {a,b} union {INTERVAL(a..b)}]); # flat at min fi fi; ## now f'(a) < 0, f'(b) > 0 ## Newton iteration for min if fpa + fpb >= 0 then d:= a; fpd:= fpa ## f'(a) closer to 0 else d:= b; fpd:= fpb ## f'(b) closer to 0 fi; fppd:= evalf(fnpp(d)); if fppd > 0 then c:= d - fpd/fppd; if (c < a) or (c > b) then ## Newton is bad, try interpolation c:= (fpb * a - fpa * b)/(fpb - fpa); fi else c:= (fpb * a - fpa * b)/(fpb - fpa); fi; fc:= evalf(fn(c)); fpc:= evalf(fnp(c)); if abs(fpc) <`gmin/eps3` then # hit it exactly RETURN(`gmin/newmin`(c,fc,bsf,msf)) elif fpc > 0 then # know it's on the left; try cutting off from a bb:= `gmin/newmin`({},fc,bsf,msf); a1:= a + (bsf - fa)/fpa; if a1 > c then RETURN(op(bb)) fi; # cut off d:= .5*(a1+c); # to make sure interval is halved fd:= evalf(fn(d)); bb:= `gmin/newmin`({},fd,op(bb)); fpd:= evalf(fnp(d)); if fpd >= 0 then # between a1 and d fa1:= evalf(fn(a1)); bb:= `gmin/newmin`({},fa1,op(bb)); fpa1:= evalf(fnp(a1)); RETURN(`gmin/convex`(fn,fnp,fnpp,a1,d,fa1,fd,op(bb), fpa1,fpd)) else # between d and c RETURN(`gmin/convex`(fn,fnp,fnpp,d,c,fd,fc,op(bb),fpd,fpc)) fi else # know it's on the right bb:= `gmin/newmin`({},fc,bsf,msf); b1:= b + (bsf - fb)/fpb; if b1 < c then RETURN(op(bb)) fi; # cut off d:= .5*(c+b1); fd:= evalf(fn(d)); bb:= `gmin/newmin`({},fd,op(bb)); fpd:= evalf(fnp(d)); if fpd <= 0 then # between d and b1 fb1:= evalf(fn(b1)); bb:= `gmin/newmin`({},fb1,op(bb)); fpb1:= evalf(fnp(b1)); RETURN(`gmin/convex`(fn,fnp,fnpp,d,b1,fd,fb1,op(bb), fpd,fpb1)) else # between c and d RETURN(`gmin/convex`(fn,fnp,fnpp,c,d,fc,fd,op(bb),fpc,fpd)) fi fi end;