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 20217 days, 23 hours, 40 minutes ago.