Ask Your Question
0

Problem with find_local_maximum when expression has an abs() function

asked 2018-07-06 05:27:02 -0500

joaoff gravatar image

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)
edit retag flag offensive close merge delete

1 answer

Sort by ยป oldest newest most voted
1

answered 2018-07-06 05:44:10 -0500

dan_fulea gravatar image

Insist 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.

edit flag offensive delete link more

Comments

Ok, 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.

joaoff gravatar imagejoaoff ( 2018-07-06 06:35:18 -0500 )edit

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools

1 follower

Stats

Asked: 2018-07-06 05:27:02 -0500

Seen: 39 times

Last updated: Jul 06