ASKSAGE: Sage Q&A Forum - RSS feedhttps://ask.sagemath.org/questions/Q&A Forum for SageenCopyright Sage, 2010. Some rights reserved under creative commons license.Fri, 06 Jul 2018 13:35:18 +0200Problem with find_local_maximum when expression has an abs() functionhttps://ask.sagemath.org/question/42844/problem-with-find_local_maximum-when-expression-has-an-abs-function/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)Fri, 06 Jul 2018 12:27:02 +0200https://ask.sagemath.org/question/42844/problem-with-find_local_maximum-when-expression-has-an-abs-function/Answer by dan_fulea for <p>Hello!</p>
<p>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 (<a href="https://trac.sagemath.org/ticket/24536">https://trac.sagemath.org/ticket/24536</a> and <a href="https://trac.sagemath.org/ticket/24537">https://trac.sagemath.org/ticket/24537</a>), but it seems that it will take a while to be solved.</p>
<p>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.</p>
<p>Another workaround is to interface with Mathematica, but this is something I would like to avoid.</p>
<p>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.</p>
<pre><code>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)
</code></pre>
<p>If I declare my expression as below, it works</p>
<pre><code>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))
</code></pre>
<p>the result is</p>
<pre><code>(1.0053729181238868, -0.32744025779801017)
</code></pre>
<p>If I declare my expression as</p>
<pre><code>den = lambda omega: up/abs(eMod)
print(find_local_maximum(den,-1,1))
</code></pre>
<p>it returns</p>
<pre><code>((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)
</code></pre>
https://ask.sagemath.org/question/42844/problem-with-find_local_maximum-when-expression-has-an-abs-function/?answer=42845#post-id-42845Insist to have a number in the lambda function, the two "omega variables" involved are not really two omega variables...
Compare with:
sage: find_local_maximum( lambda u: up(omega=u) / abs(eMod(omega=u)), -1, 1)
(1.00331600582447, 0.23606797749978958)
which works.
Above, `u` takes some values in the given interval. Fix one of them. Then `eMod(omega=u)` is the value of `eMod` in this fixed `u`-value. One can apply the `abs` on it, getting a number. In contrast, in `lambda u: up/abs(eMod)` we have for each `u` as result a symbolic expression. And the `abs` is applied also symbolically.
Fri, 06 Jul 2018 12:44:10 +0200https://ask.sagemath.org/question/42844/problem-with-find_local_maximum-when-expression-has-an-abs-function/?answer=42845#post-id-42845Comment by joaoff for <p>Insist to have a number in the lambda function, the two "omega variables" involved are not really two omega variables...
Compare with:</p>
<pre><code>sage: find_local_maximum( lambda u: up(omega=u) / abs(eMod(omega=u)), -1, 1)
(1.00331600582447, 0.23606797749978958)
</code></pre>
<p>which works.</p>
<p>Above, <code>u</code> takes some values in the given interval. Fix one of them. Then <code>eMod(omega=u)</code> is the value of <code>eMod</code> in this fixed <code>u</code>-value. One can apply the <code>abs</code> on it, getting a number. In contrast, in <code>lambda u: up/abs(eMod)</code> we have for each <code>u</code> as result a symbolic expression. And the <code>abs</code> is applied also symbolically. </p>
https://ask.sagemath.org/question/42844/problem-with-find_local_maximum-when-expression-has-an-abs-function/?comment=42846#post-id-42846Ok, thanks for explaining!
However, this doesn't solve the problem. I cannot discretize the interval. This is what I do in the first workaround.
As you can see, the correct result is
(1.0053729181238868, -0.32744025779801017)
and you have obtained
(1.00331600582447, 0.23606797749978958)
I can obtain a more precise result with the first workaround mentioned above.
With your answer, I discovered that I need to upgrade my Sage system. I am using version 7.5.1 and in this version, your solution doesn't work, it returns
RuntimeError: ECL says: THROW: The catch MACSYMA-QUIT is undefined.
but it works in Sage Cell.Fri, 06 Jul 2018 13:35:18 +0200https://ask.sagemath.org/question/42844/problem-with-find_local_maximum-when-expression-has-an-abs-function/?comment=42846#post-id-42846