Here's the code I'm working with:
# Computes the theta correspondence function for the boundary defined
by G w/ center a going to 0
def boundaryMapper(G,a,n):
B = 2*pi
tk = [x*B/n for x in range(n)] # the collocation points
Gprime = derivative(G) # generate various intermediate functions
Hconj(t) = 1/(2*pi*i)*(Gprime(t)/abs(Gprime(t))/(a-G(t))).conjugate
()
A(t,s) = 1/(2*pi*i)*(Gprime(s)/abs(Gprime(s))/(G(s)-G(t))-(Gprime
(t)/abs(Gprime(t))/(G(s)-G(t))).conjugate())
K(t,s) = sqrt(abs(Gprime(s))) * A(t,s) * sqrt(abs(Gprime(t)))
g(t) = sqrt(abs(Gprime(t)))*Hconj(t)
Ka = Matrix(CDF,n,n) # Generate the matrices for Nystrom's method
for k in range(n):
for j in range(n):
if k==j:
Ka[k,j] = 0
else:
Ka[k,j] = B/n*K(tk[k],tk[j])
G = matrix(CDF,n,1)
for k in range(n):
G[k,0] = g(tk[k])
F = identity_matrix(CDF,n) - Ka
theta_points = F.solve_right(G) # find the psi function values at
the collocation points
phi(t) = g(t) # find the psi and theta functions: This gets
ridiculous
for j in range(n):
phi(t) = phi(t) + K(t,tk[j])*theta_points[j,0]*2*pi/n
phi2(t) = -i*phi(t)^2*Gprime(t)
print "got this far"
# below: my workaround for the lack of a symbolic 'arg' function
theta(t) = acos(real(phi2(t))/abs(phi2(t))) # This (the correct
one) never completes.
#theta(t) = acos(real(phi(t))/abs(phi(t))) # This (not the correct
one) has no problem.
print "hey I made it"
GprimeA = map(Gprime,tk)
theta_array = range(n)
for j in range(n): # find the theta correspondence function at the
collocation points
a = N(arg(-i*theta_points[j,0]^2*GprimeA[j]))
if a < -.0000001:
theta_array[j] = (N(tk[j]),N(a + 2*pi))
else:
theta_array[j] = (N(tk[j]),a)
return theta_array, theta
ellipse(x) = e^(i*x)-.5*e^(-i*x)
#show(parametric_plot((real(ellipse),imag(ellipse)),0,2*pi))
theta_array , theta = boundaryMapper(ellipse,0,16) # args are figure,
center, collocation points.
# decrease the number of collocation points to run it faster.
print "done"
list_plot(theta_array)
This function is intended to be an implementation of the Szego kernel
method for computing a riemann map. (google it if you want, I couldn't
find a free site giving details). Essentially the 'theta_array'
contains points for the function I'm trying to compute. The theta
function is meant to interpolate those points using a special
interpolation formula (the interpolation happens in the declaration of
the psi function). My primary problem is that the program hangs when I
try to execute the final step of generating the theta function. Can
anyone tell me why/how to fix it?
As a side note, my code also runs much slower than I think it should.
If anyone spots innefficient things that I'm doing i.e. "a for loop?
Everyone knows map is faster" please point it out.
If this should be posted somewhere else, let me know.
Thanks,
Ethan
--~--~---------~--~----~------------~-------~--~----~
To post to this group, send email to [email protected]
To unsubscribe from this group, send email to
[email protected]
For more options, visit this group at
http://groups.google.com/group/sage-support
URLs: http://www.sagemath.org
-~----------~----~----~----~------~----~------~--~---
