2366
Comment: Started section with "Winding number of a plane curve"
|
← Revision 7 as of 2020-06-01 18:42:11 ⇥
2538
|
Deletions are marked like this. | Additions are marked like this. |
Line 7: | Line 7: |
by Pablo Angulo. Computes winding number as an integral, and also as a intersection number with a half line through the origin. | by Pablo Angulo. Computes winding number (with respect to the origin!) as an integral, and also as a intersection number with a half line through the origin. |
Line 9: | Line 9: |
{{{ | {{{#!sagecell |
Line 13: | Line 13: |
r2 = x^2 + y^2 | r2 = x**2 + y**2 |
Line 16: | Line 16: |
integrando = (x*yp -y*xp)/r2 | integrando = (x*yp -y*xp) / r2 |
Line 18: | Line 18: |
return round(i/(2*pi)) | return round(i / (2 * pi)) |
Line 24: | Line 24: |
delta= (b-a)/N | delta = (b - a) / N |
Line 26: | Line 26: |
for t in srange(a, b ,delta): | for t in srange(a, b, delta): |
Line 28: | Line 28: |
zeros.append(find_root(f, t-epsilon, t+delta+epsilon)) | zeros.append(find_root(f, t - epsilon, t + delta + epsilon)) |
Line 32: | Line 32: |
if not zeros: return zeros if abs(zeros[0] + 2*pi - zeros[-1])<epsilon: |
if not zeros: return zeros if abs(zeros[0] + 2*pi - zeros[-1]) < epsilon: |
Line 37: | Line 38: |
if abs(c - zeros_cleaned[-1])>epsilon: | if abs(c - zeros_cleaned[-1]) > epsilon: |
Line 39: | Line 40: |
if abs(zeros[0] + 2*pi - zeros[-1])<epsilon: | if abs(zeros[0] + 2*pi - zeros[-1]) < epsilon: |
Line 46: | Line 47: |
x = x.function(t); y = y.function(t); | x = x.function(t) y = y.function(t) |
Line 48: | Line 50: |
raise ValueError, "Curve is not closed!" | raise ValueError("Curve is not closed!") |
Line 55: | Line 57: |
html(r'$\int \frac{1}{x^2 + y^2}(xdy-ydx)=%d$'%winding_number_integral(x,y,a,b)) | pretty_print(html(r'$\int \frac{1}{x^2 + y^2}(xdy-ydx)=%d$'%winding_number_integral(x,y,a,b))) |
Line 57: | Line 59: |
zeros = all_the_zeros(x,a, b) | zeros = all_the_zeros(x, a, b) |
Line 70: | Line 72: |
print 'Winding number = (# of red points) - (# of green points): ', wn | print('Winding number = (number of red points) - (number of green points): {}'.format(wn)) |
Line 73: | Line 75: |
arrow((x(0),y(0)), (x(0) + xp1(0), y(0) + yp1(0))) ) | arrow((x(0),y(0)), (x(0) + xp1(0), y(0) + yp1(0))) + point2d([(0,0)], color = 'black', pointsize = 70)) |
Sage Interactions - Topology
goto interact main page
Winding number of a plane curve
by Pablo Angulo. Computes winding number (with respect to the origin!) as an integral, and also as a intersection number with a half line through the origin.