Alternative Implicit Plotting

From www.norsemathology.org

Revision as of 14:57, 16 October 2008

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( (x,y,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 ImplicitPlot( 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()
    xIters = floor( (xRange['max'] - xRange['min']) / xRange['step'] )
    yIters = floor( (yRange['max'] - yRange['min']) / yRange['step'] )
    for s in range( 0, xIters ):
      for t in range( 0, yIters ):
	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