Alternative Implicit Plotting

Alternative Implicit Plot Functions

This page will describe an alternative implicit plot algorithm written by Dr. Wilkinson.

Dr. Wilkinson's Mathematica Code:

implicitPlot[equation_, xRange : {x_, xmin_, xmax_, dx_}, 
  yRange : {y_, ymin_, ymax_, dy_}] :=
 
 Module[{values, f, processSquare, lineSegments = {}},
  f[s_, t_] = Apply[Subtract, equation] /. {x -> s, y -> t};
  values = Table[{x, y, f[x, y]}, xRange, yRange];
  
  (* Inner function processes the square with indices for given lower left corner. *)
  processSquare[{s_, t_}] := Module[{roots = {}},
    (* Check along the bottom. *)
    
    If[Last[values[[s, t]]]*Last[values[[s + 1, t]]] < 0,
     roots = 
      Join[roots, {{(values[[s, t]][[1]]*values[[s + 1, t]][[3]] - 
            values[[s + 1, t]][[1]]*
             values[[s, t]][[3]])/(values[[s + 1, t]][[3]] - 
            values[[s, t]][[3]]), values[[s, t]][[2]]}}]
     ];
    (* Check along the left side. *)
    
    If[Last[values[[s, t]]]*Last[values[[s, t + 1]]] < 0,
     roots = 
      Join[roots, {{values[[s, t]][[
          1]], (values[[s, t]][[2]]*values[[s, t + 1]][[3]] - 
            values[[s, t + 1]][[2]]*
             values[[s, t]][[3]])/(values[[s, t + 1]][[3]] - 
            values[[s, t]][[3]])}}]
     ];
    (* Check along the top. *)
    
    If[Last[values[[s, t + 1]]]*Last[values[[s + 1, t + 1]]] < 0,
     roots = 
      Join[roots, {{(values[[s, t + 1]][[1]]*
             values[[s + 1, t + 1]][[3]] - 
            values[[s + 1, t + 1]][[1]]*
             values[[s, t + 1]][[3]])/(values[[s + 1, t + 1]][[3]] - 
            values[[s, t + 1]][[3]]), values[[s, t + 1]][[2]]}}]
     ];
    (* Check along the right side. *)
    
    If[Last[values[[s + 1, t]]]*Last[values[[s + 1, t + 1]]] < 0,
     roots = 
      Join[roots, {{values[[s + 1, t]][[
          1]], (values[[s + 1, t]][[2]]*values[[s + 1, t + 1]][[3]] - 
            values[[s + 1, t + 1]][[2]]*
             values[[s + 1, t]][[3]])/(values[[s + 1, t + 1]][[3]] - 
            values[[s + 1, t]][[3]])}}]
     ];
    If[Length[roots] == 2, 
     lineSegments = Join[lineSegments, {Line[roots]}]
     ]
    ];
  Do[processSquare[{s, t}], {s, (xmax - xmin)/dx}, {t, (ymax - ymin)/
    dy}];
  Return[
   Graphics[lineSegments, Axes -> True, 
    PlotRange -> {{xmin, xmax}, {ymin, ymax}}]]
  ]

Dr. Wilkinson's Algorithm Rewritten in Sage

def table(a,b,f):
    assert len(a) == 3 and len(b) == 3
    ls = []
    for x in srange(a[0],a[1],a[2],include_endpoint=True):
        inner = []
        for y in srange(b[0],b[1],b[2],include_endpoint=True):
            inner.append( (float(x),float(y),float( f(x,y) ) ) )
        ls.append(inner)
    return ls

'''
  
  eq : symbolic equation of two variables
  xRange, yRange : dictionaries containg the following keys
    var: symbolic variable 
    min: minimum x/y value
    max: maximum x/y value
    step: dx/dy
    
'''
def myImplicitPlot( eq, xRange, yRange ):
    f( xRange['var'] , yRange['var'] ) = eq.subtract_from_both_sides( eq.lhs() ).right()
    values = table((xRange['min'],xRange['max'],xRange['step']),(yRange['min'],yRange['max'],yRange['step']),f)
    lineSegements = []
    g = Graphics()
    for s in range( 0, len(values) - 1):
      for t in range( 0, len(values[s]) - 1):
	roots = []
	#Check the bottom
	if values[s][t][2] * values[s+1][t][2] < 0:
	  x1 = values[s][t][0]   ; f1 = values[s][t][2]  ; y1 = values[s][t][1]
	  x2 = values[s+1][t][0] ; f2 = values[s+1][t][2] 
	  roots.append( [ (x1*f2 - x2*f1) / (f2 - f1) , y1 ] )
	#Check on the left side
	if values[s][t][2] * values[s][t+1][2] < 0:
	  y1 = values[s][t][1]   ; f1 = values[s][t][2]  ; x1 = values[s][t][0]
	  y2 = values[s][t+1][1] ; f2 = values[s][t+1][2]
	  roots.append( [ x1, (y1*f2 - y2*f1) / (f2 - f1) ] )
	#Check the top
	if values[s][t+1][2] * values[s+1][t+1][2] < 0:
	  x1 = values[s][t+1][0]   ; f1 = values[s][t+1][2]  ; y1 = values[s][t+1][1]
	  x2 = values[s+1][t+1][0] ; f2 = values[s+1][t+1][2]
	  roots.append( [ (x1*f2 - x2*f1) / (f2 - f1) , y1 ] )
	#Check on the right side
	if values[s+1][t][2] * values[s+1][t+1][2] < 0:
	  y1 = values[s+1][t][1]   ; f1 = values[s+1][t][2]  ; x1 = values[s+1][t][0]
	  y2 = values[s+1][t+1][1] ; f2 = values[s+1][t+1][2] 
	  roots.append( [ x1, (y1*f2 - y2*f1) / (f2 - f1) ] )
	if len(roots) == 2:
	  lineSegements.append( line( roots ) )
    g += plot( lineSegements, rgbcolor=(0,0,0) ) 
    return g

Notes on the Sage Version

In the table function you can see that x,y, and f(x,y) are all converted to float after they are calculated. The type of x,y and f(X,y) is initially sage.rings.real_mpfr.RealNumber. Though sage.rings.real_mpfr.RealNumber can use any specified precision, converting to float is necessary because the arithmetic performed in the ImplicitPlot function will take significantly longer with sage.rings.real_mpfr.RealNumber than with float. To give you an idea with a x and y step size of 0.1 with sage.rings.real_mpfr.RealNumber would take 10 to 15 minutes and with float it takes between 20 and 25 seconds.