# Alternative Implicit Plotting

(Difference between revisions)
 Revision as of 13:53, 17 October 2008 (edit)← Previous diff Revision as of 13:59, 17 October 2008 (edit) (undo)m Next diff → Line 134: Line 134:

var('x y')                                                                                                                     var('x y')
-                                                            xrange = { 'var' : x, 'min' : -1, 'max' : 5, 'step' : .05 }  +                                                                 xrange = { 'var' : x, 'min' : -3, 'max' : 3, 'step' : .05 }
yrange = { 'var' : y, 'min' : -3, 'max' : 3, 'step' : .05 }                                                                    yrange = { 'var' : y, 'min' : -3, 'max' : 3, 'step' : .05 }
myplot = myImplicitPlot( y^2 == x^3 - x, xrange, yrange )                                                                      myplot = myImplicitPlot( y^2 == x^3 - x, xrange, yrange )

## 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.

#### Example

Code:

```var('x y')
xrange = { 'var' : x, 'min' : -3, 'max' : 3, 'step' : .05 }
yrange = { 'var' : y, 'min' : -3, 'max' : 3, 'step' : .05 }
myplot = myImplicitPlot( y^2 == x^3 - x, xrange, yrange )
```

Output: