Oops! a piece of the code is missing. The complete code is below.
n=6 # polynomial order
listOmega=[]
for i in range(n+1):
listOmega.append(i/n)
listPhases=[0]
for i in range(1,n+1):
listPhases.append(listOmega[i]^2/2 + 9/2*listOmega[i])
Alpha0=listOmega[1]/tan(listPhases[1])
listAlpha=[Alpha0]
for i in range(1,n):
firstline=listOmega[i+1]^2-listOmega[i]^2
lastline=Alpha0-listOmega[i+1]/tan(listPhases[i+1])
ic=1
while ic<i:
lastline=listAlpha[ic]-((listOmega[i+1]^2 -
listOmega[ic]^2)/lastline)
ic+=1
listAlpha.append(firstline/lastline)
var('s,omega')
H=[1,s + listAlpha[0]]
for i in range(2,n+1):
H.append(listAlpha[i-1]*H[i-1] + (s^2 + listOmega[i-1]^2)*H[i-2])
f=sqrt(H[6](s=I*omega)*H[6](s=-I*omega));
g=[]
for i in range(n+1):
g.append(f*chebyshev_T(i,omega)/sqrt(1-omega^2))
c=[]
for i in range(n+1):
if i==0:
c.append((1/pi)*g[i].nintegral(omega,-1,1)[0])
else:
c.append((2/pi)*g[i].nintegral(omega,-1,1)[0])
up=0
for i in range(n+1):
up=up + c[i]*chebyshev_T(i,omega)
hmod(omega)=(H[6](s=I*omega)).polynomial(CC)
den=lambda omega: up/abs(hmod(omega))
norm_factor=find_local_maximum(den,-1,1,tol=1e-18); norm_factor
print("Point where local maximum occurs: " + str(norm_factor[1]))
print("Local maximum: " + str(norm_factor[0](omega=norm_factor[1]).n()))
plot(den(omega),-1,1)
On Monday, January 15, 2018 at 11:11:50 AM UTC+1, João Alberto de França
Ferreira wrote:
>
> Sorry for the delay, I took some time to put the code to work. However,
> there's something strange.
>
> Below is a "concise" version of my code. The strange thing that I noticed
> is that the *find_local_maximum* function returns a polynomial as the
> maximum value, and changing the tolerance (*tol*) argument of the
> function do not change the result. You can see in the graph that the
> maximum point is slightly away from the calculated point.
>
> Thank you for the help!
>
> n=6 # polynomial order
> listOmega=[]
> for i in range(n+1):
> listOmega.append(i/n)
>
> listPhases=[0]
> for i in range(1,n+1):
> listPhases.append(listOmega[i]^2/2 + 9/2*listOmega[i])
>
> Alpha0=listOmega[1]/tan(listPhases[1])
> listAlpha=[Alpha0]
> for i in range(1,n):
> firstline=listOmega[i+1]^2-listOmega[i]^2
> lastline=Alpha0-listOmega[i+1]/tan(listPhases[i+1])
> ic=1
> while ic<i:
> lastline=listAlpha[ic]-((listOmega[i+1]^2 -
> listOmega[ic]^2)/lastline)
> ic+=1
> listAlpha.append(firstline/lastline)
>
> var('s,omega')
> H=[1,s + listAlpha[0]]
> for i in range(2,n+1):
> H.append(listAlpha[i-1]*H[i-1] + (s^2 + listOmega[i-1]^2)*H[i-2])
>
> f=sqrt(H[6](s=I*omega)*H[6](s=-I*omega));
> g=[]
> for i in range(n+1):
> g.append(f*chebyshev_T(i,omega)/sqrt(1-omega^2))
>
> up=0
> for i in range(n+1):
> up=up + c[i]*chebyshev_T(i,omega)
>
> hmod(omega)=(H[6](s=I*omega)).polynomial(CC)
> den=lambda omega: up/abs(hmod(omega))
> norm_factor=find_local_maximum(den,-1,1,tol=1e-18); norm_factor
> print("Point where local maximum occurs: " + str(norm_factor[1]))
> print("Local maximum: " + str(norm_factor[0](omega=norm_factor[1]).n()))
> plot(den(omega),-1,1)
>
> On Friday, January 12, 2018 at 4:51:32 PM UTC+1, Vegard Lima wrote:
>>
>> Hi
>>
>> On Fri, Jan 12, 2018 at 4:10 PM, João Alberto de França Ferreira
>> <[email protected]> wrote:
>> > Is there a way to compute the local maximum of the absolute value of a
>> > polynomial with complex coefficients? the below snippet should
>> exemplify my
>> > problem.
>> ...
>>
>> Maybe a trick like this works:
>> f(x) = x^2 + 2*I*x + 2
>> g = lambda x: abs(f(x))
>> find_local_maximum(g,-1,1)
>>
>> this gives
>> (3.6055511727769516, 0.99999996297566163)
>>
>> while
>> find_local_maximum(abs(f(x)),-1,1)
>> gives
>> TypeError: unable to coerce to a real number
>>
>> The difference is that one is a python function while
>> the other one is a SymbolicExpression, not sure
>> why one works and the other doesn't.
>>
>>
>> Thanks,
>> --
>> Vegard
>>
>
--
You received this message because you are subscribed to the Google Groups
"sage-support" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to [email protected].
To post to this group, send email to [email protected].
Visit this group at https://groups.google.com/group/sage-support.
For more options, visit https://groups.google.com/d/optout.
