Hello!
I need to find the local maximum of an expression in an interval. It seems that the find_local_maximum() function has a problem when the argument expression contains the abs() function. This problem was reported to the sage-support group and two tickets were opened (https://trac.sagemath.org/ticket/24536 and https://trac.sagemath.org/ticket/24537), but it seems that it will take a while to be solved.
As a workaround, I can plot the function in that interval and look for the point with the maximum value, but this is not precise and produces a slightly different value each time the code is run.
Another workaround is to interface with Mathematica, but this is something I would like to avoid.
So, does anyone knows an alternative way to find the local maximum in an interval for an expression containing abs() functions? A piece of code is shown below as an example.
s = polygen(CC, "s")
omega = polygen(CC, "omega")
pi = CC.pi()
I = CC.0
p_order = 6
ePoly = s^6 + 2.52349407705767*s^5 + 4.57149144245730*s^4 + 4.95095921397024*s^3 + 3.70209076479638*s^2 + 1.64521458323101*s + 0.361999139300727
##### Generate P(s) polynomial #####
pMod(omega) = sqrt(ePoly(s=I*omega)*ePoly(s=-I*omega))
integrand = []
for a in range(p_order+1):
integrand.append(pMod*chebyshev_T(a, omega)/sqrt(1-omega^2))
c = []
for a in range(p_order+1):
if a == 0:
c.append((1/pi)*Rational(integrand[a].nintegral(omega, -1, 1)[0]))
else:
c.append((2/pi)*Rational(integrand[a].nintegral(omega, -1, 1)[0]))
up = 0
for a in range(p_order+1):
up = up + c[a]*chebyshev_T(a, omega)
eMod = ePoly(s=I*omega)
### First workaround: Find the plot point with maximum value
den = lambda omega: up/abs(eMod)
pl = plot(den(omega), -1, 1)
plPoints = []
for pair in list(pl[0]):
plPoints.append(pair[1])
# normFactor is multiplied by a constant to ensure that P(s) lies below the abscissa axis
normFactor = 1.00001*Rational(max(plPoints))
### Second workaround: Interfacing with Mathematica to calculate normFactor
# den = up/abs(eMod)
# mden = mathematica(den)
# mNormFactor = mden.NMaximize(-1 < omega < 1, omega)[1]
# normFactor = mNormFactor.sage()
pPoly = up/normFactor
pPoly = pPoly(omega=s/I)
pPoly = pPoly.polynomial(CC)
pPoly = pPoly.subs(x=s)
print(pPoly)
If I declare my expression as below, it works
den = lambda omega: (0.141162317194057*omega^6 + 0.173051768044343*omega^4 + 0.140157830888798*omega^2 + 0.359575336011182)/abs(-1.00000000000000*omega^6 + 2.52349407705767*I*omega^5 + 4.57149144245730*omega^4 - 4.95095921397024*I*omega^3 - 3.70209076479638*omega^2 + 1.64521458323101*I*omega + 0.361999139300727)
print(find_local_maximum(den,-1,1))
the result is
(1.0053729181238868, -0.32744025779801017)
If I declare my expression as
den = lambda omega: up/abs(eMod)
print(find_local_maximum(den,-1,1))
it returns
((0.141162317194057*omega^6 + 0.173051768044343*omega^4 + 0.140157830888798*omega^2 + 0.359575336011182)/abs(-1.00000000000000*omega^6 + 2.52349407705767*I*omega^5 + 4.57149144245730*omega^4 - 4.95095921397024*I*omega^3 - 3.70209076479638*omega^2 + 1.64521458323101*I*omega + 0.361999139300727), 0.99999996297566163)
