Revision history [back]

This is more a comment and code review than an answer, but since it's a bit long I'm using the answer field.

Below I'll use r for sqrt(x^2+y^2) and C_R for the circle of radius R centered at the origin.

1. In the definition of C you don't need parametric plots. Use Sage's built-in circle command.

2. It looks like c1, c3, c4 are meant to be discretized circles of radius 1, 3, and 4.

   In their definition, it seems that there should be no *u inside the cos and sin.

3. Are you using discretized circles on purpose, or as a workaround for continuous circles?

   If it's a workaround for continuous circles, note that dist((x,y),C_R) is just abs(R-r).

   To see that, draw the circles C_R and the circle of radius r centered at (x,y).

4. The quantity d = setdist(x,S) is always non-negative, so min(1,1-d) is always 1-d.

   So your functions h0 and h1 only differ if 1 <= r <= 3, and in this case h0 = 0 and h1 = 1.

   Did you mean something else in the definition of h1?

Assuming discretized circles were a workaround for continuous circles, I would rewrite your code as follows.

# Use "real double field" floating-point values for efficiency
zero = RDF.zero()
one = RDF.one()
four = 4 * one
twopi = 2 * RDF.pi()

# Circles as graphic objects
C1 = circle((0,0), 1)
C3 = circle((0,0), 3)
C4 = circle((0,0), 4)
C = C1 + C3 + C4

# useful functions

def norm2((x,y)):
    """
    Norm squared of `(x,y)`
    """
    return x^2 + y^2

def norm((x,y)):
    """
    Norm of `(x,y)`
    """
    return norm2((x,y)).sqrt()

# h0 and h1

def h0(x,y):
    """
    Return piecewise affine function of `r = norm(x,y)`

    Let C_n be the circle of radius n, for n = 1, 3, 4.

    With r = norm((x,y)) = sqrt(x^2+y^2), h0(x,y) is:
    * r inside C_1
    * 1 between C_1 and C_3
    * 4 - r between C_3 and C_4
    * 0 outside C_4
    """
    r = norm((x,y))
    if r < 1:
        return r
    if 1 <= r <= 3:
        return one
    if 3 <= r <= 4:
        return four - r
    return zero

def h1(x,y):
    """
    Return piecewise affine function of `r = norm(x,y)`

    Let C_n be the circle of radius n, for n = 1, 3, 4.

    With r = norm((x,y)) = sqrt(x^2+y^2), h0(x,y) is:
    * r inside C_1
    * 0 between C_1 and C_3
    * 4 - r between C_3 and C_4
    * 0 outside C_4
    """
    r = norm((x,y))
    if r < 1:
        return r
    if 1 <= r <= 3:
        return zero
    if 3 <= r <= 4:
        return four - r
    return zero

def dx(x,y):
    return y * h0(x,y)

def dy(x,y):
    return -x * h0(x,y) + y * h1(x,y)

If discretized circles are relevant, varying the number of points may help illustrate your point, so use a parameter.

# discretized circles
n = 99
DC1 = [(cos(i*twopi/n), sin(i*twopi/n)) for i in xrange(n)]
DC3 = [(3*x, 3*y) for (x,y) in DC1]
DC4 = [(4*x, 4*y) for (x,y) in DC1]

We would also need to rewrite the functions h0 and h1 using these discretized circles.

For that, I would rewrite your auxiliary functions as follows:

# more useful functions

def dist2((x,y),(xx,yy)):
    """
    Distance squared from `(x,y)` to `(xx,yy)`
    """
    return norm2((x - xx, y - yy))

def dist((x,y),(xx,yy)):
    """
    Distance squared from `(x,y)` to `(xx,yy)`
    """
    return norm((x - xx, y - yy))

def setdist2(u,S):
    """
    Distance squared from point to set
    """
    return min(dist2(u,v) for v in S)

def setdist(u,S):
    """
    Distance from point to set
    """
    return setdist2(u,S).sqrt()

Once you answer some of the questions above, we can investigate the vector-field-plotting part.

x, y = var('x y')
VF = plot_vector_field((dx(x,y),dy(x,y)), (x,-4,4), (y,-4,4))
(C + VF).show(aspect_ratio=1)