Download or view maximize.frink in plain text format
use allTransforms.frink
use functionUtils.frink
/** This file contains routines to minimize or maximize a function over a
certain range using interval arithmetic techniques. And, maybe, even
better, using symbolic differentiation techniques! */
maximize[f, minX, maxX, minXDiff] := intervalOptimize[f, minX, maxX, minXDiff, 1]
minimize[f, minX, maxX, minXDiff] := intervalOptimize[f, minX, maxX, minXDiff, -1]
intervalOptimize[f, minX, maxX, minXDiff, factor=1] :=
{
widthStep = maxX - minX
fullX = new interval[minX, minX + widthStep/2, maxX]
// solutions is an array of arrays. Each array contained is a contiguous
// range of intervals of still-valid solutions.
solutions = [[minX, fullX, maxX]]
// solutions.push[[fullX]]
fullY = factor * f[fullX]
lowerBound = infimum[fullY]
hardBound = lowerBound
// println["lowerBound is $lowerBound"]
// Try the endpoints because most functions are monotonic
minY = factor f[minX]
ly = infimum[minY]
if ly < hardBound
hardBound = ly
maxY = factor f[maxX]
ly = infimum[maxY]
if ly < hardBound
hardBound = ly
while (widthStep > minXDiff)
{
widthStep = widthStep / 2
// println["input solutions are: $solutions"]
newSolutions = new array
for solutionBlock = solutions
{
newSolutionBlock = undef
SOLUTION:
for solution = solutionBlock
{
y = factor f[solution]
ly = infimum[y]
// println["ly is $ly"]
if ly > lowerBound
{
lowerBound = ly
// println["new lower bound is $lowerBound"]
}
// If y is possibly greater or equal to the lower bound, then this
// interval may contain a solution.
if y PGE lowerBound
{
// println["width is $width"]
// Bifurcate the x interval and add the 2 new intervals to the
// newSolutionBlock
lower = infimum[solution]
upper = supremum[solution]
width = (upper-lower)/2
halfWidth = width/2
mid = lower + width
midLower = lower + halfWidth
nx = new interval[lower, midLower, mid]
ny = factor f[nx]
// println["nx is $nx, ny is $ny"]
my = mainValue[ny]
if my == undef
println["No mid for $nx, $ny"]
if my != undef && my > hardBound
{
// println["hardBound was $hardBound"]
// println["hardBound getting set by $my"]
hardBound = my
lowerBound = hardBound
// println["new hard bound is $hardBound for $nx (1)"]
}
if ny PGE hardBound and ny PGE lowerBound
{
ly = infimum[ny]
if ly > lowerBound
{
lowerBound = ly
// println["new lower bound is $lowerBound (1)"]
}
// println["Pushing $nx (1)"]
// Construct a new solution block if it doesn't exist
if newSolutionBlock == undef
newSolutionBlock = new array
newSolutionBlock.push[nx]
} else
{
if newSolutionBlock != undef
{
newSolutions.push[newSolutionBlock]
newSolutionBlock = undef
}
}
if lower == upper
next SOLUTION
midUpper = mid + halfWidth
nx = new interval[mid, midUpper, upper]
ny = factor f[nx]
my = mainValue[ny]
if my != undef && my > hardBound
{
// println["hardBound was $hardBound"]
// println["hardBound getting set by $my"]
hardBound = my
lowerBound = hardBound
// println["new hard bound is $hardBound for $nx (2)"]
}
if ny PGE hardBound and ny PGE lowerBound
{
ly = infimum[ny]
if ly > lowerBound
{
lowerBound = ly
// println["new lower bound is $lowerBound (2)"]
}
// println["Pushing $nx (2)"]
// Construct a new solution block if it doesn't exist
if newSolutionBlock == undef
newSolutionBlock = new array
newSolutionBlock.push[nx]
} else
{
if newSolutionBlock != undef
{
newSolutions.push[newSolutionBlock]
newSolutionBlock = undef
}
}
}
}
if newSolutionBlock != undef
newSolutions.push[newSolutionBlock]
}
solutions = newSolutions
println["lowerBound is $lowerBound"]
println["There are " + length[flatten[solutions]] + " solutions"]
// println["Solutions: $solutions"]
println["width is $widthStep"]
}
// println[solutions]
s = new array
for sols = solutions
{
maxMain = undef
xMax = minX
s1 = undef
if (length[sols] > 0)
{
for i = 0 to length[sols] - 1
{
nx = sols@i
ny = factor f[nx]
// println["Prospective solution " + sols@i + " = " + ny]
my = mainValue[ny]
if my != undef && (maxMain == undef || my > maxMain)
{
maxMain = my
if mainValue[nx] != undef
{
xMax = mainValue[nx]
// println["New xMax = $xMax at $nx"]
}
}
if ny PGE lowerBound
if s1 == undef
s1 = sols@i
else
s1 = union[s1, sols@i]
}
}
if xMax != undef
s1 = new interval[infimum[s1], xMax, supremum[s1]]
if s1 != undef
s.push[[s1, f[s1]]]
}
return s
}
symbolicOptimize[f, minX, maxX, factor=1] :=
{
args = functionArgumentsAsSymbols[f]
if length[args] != 1
{
println["Wrong number of arguments"]
return undef
}
println["Arguments are $args"]
derivativeFunc2 = makeDerivative[f]
println["Before solve: $derivativeFunc2"]
derivative2 = transformExpression[derivativeFunc2]
println["Derivative is: $derivative2"]
solve = makeSolve[derivative2, 0, args@0]
println["Before solve: $solve"]
solved = transformExpression[solve]
println["After solve: $solved"]
type = type[solved]
println["Type is " + type]
if (type == "Solve") // Is this a solved expression with ===
{
right = getChild[solved, 1]
println["Right is " + right]
solution = eval[right]
println["Solution is " + solution]
return solution
} else
return undef
}
//setPrecision[15]
//println[maximize[{|x| -x^2+2x-3}, 0, 4 pi]]
//println[maximize[{|x| x}, 0, 10, 1e-20]]
//println[]
//println[join["\n",maximize[{|x| x}, 0, 10, 1e-20]]]
//println[]
//println[join["\n",maximize[{|x| x sin[x]}, -15, 15, 1e-8]]]
//println[]
//println[join["\n",maximize[{|x| x cos[x]}, -23, 23, 1e-7]]]
//println[]
//println[join["\n",minimize[{|x| x^x}, 0.01, 4 pi, 1e-8]]]
//println[]
//println[join["\n",maximize[{|x| x^2}, -4., 4., 1e-10]]]
//println[]
//println[join["\n",minimize[{|x| -x^2}, -4., 4., 1e-10]]]
//println[]
//println[join["\n",maximize[{|x| -x^2}, -4., 4., 1e-20]]]
//println[]
//println[join["\n",maximize[{|x| -abs[x^2]}, -4., 4., 1e-10]]]
//println[]
//println[join["\n",maximize[{|x| abs[x^2]}, -4., 4., 1e-10]]]
//println[]
//println[join["\n",maximize[{|x| x cos[x] + sin[x]}, -10, 10, 1e-8]]]
/*
use sun.frink">sun.frink
println[]
println[maximize[{|date| sunMoonAngularSeparation[date, 40 degrees North, 105 degrees West]}, #2019-01-20 7:30 PM #, #2019-01-20 11:59 PM#, .1 s]]
*/
/*
showApproximations[false]
symbolicMode[true]
a = symbolicOptimize[{|x| sin[2 x]}, -pi, pi]
*/
Download or view maximize.frink in plain text format
This is a program written in the programming language Frink.
For more information, view the Frink
Documentation or see More Sample Frink Programs.
Alan Eliasen was born 19972 days, 7 hours, 21 minutes ago.