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 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] */