# Alternative Implicit Plotting

(Difference between revisions)
 Revision as of 14:40, 16 October 2008 (edit) (New page: == Alternative Implicit Plot Functions == This page will describe an alternative implicit plot algorithm written by Dr. Wilkinson. === Dr. Wilkinson's Mathematica Code: ===
impl...)← Previous diff                                                                                                               Current revision (15:17, 23 October 2008) (edit) (undo)m
(6 intermediate revisions not shown.)
Line 76:                                                                                                                                                                                                                                                                                                                                                                              Line 76:
inner = []                                                                                                                                                             inner = []
for y in srange(b[0],b[1],b[2],include_endpoint=True):                                                                                                                 for y in srange(b[0],b[1],b[2],include_endpoint=True):
-                                                                                                                                                                                                                                                                        inner.append( (x,y,f(x,y)) )                                                                                 +                                                         inner.append( (float(x),float(y),float( f(x,y) ) ) )
ls.append(inner)                                                                                                                                                       ls.append(inner)
return ls                                                                                                                                                              return ls
Line 90:                                                                                                                                                                                                                                                                                                                                                                              Line 90:

'''                                                                                                                                                                    '''
-                                                                                                                                                                                                                                                                        def ImplicitPlot( eq, xRange, yRange ):                                                                      +                                                         def myImplicitPlot( eq, xRange, yRange ):
f( xRange['var'] , yRange['var'] ) = eq.subtract_from_both_sides( eq.lhs() ).right()                                                                                   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)                                                            values = table((xRange['min'],xRange['max'],xRange['step']),(yRange['min'],yRange['max'],yRange['step']),f)
lineSegements = []                                                                                                                                                     lineSegements = []
g = Graphics()                                                                                                                                                         g = Graphics()
-                                                                                                                                                                                                                                                                        xIters = floor( (xRange['max'] - xRange['min']) / xRange['step'] )                                           +                                                         for s in range( 0, len(values) - 1):
-                                                                                                                                                                                                                                                                        yIters = floor( (yRange['max'] - yRange['min']) / yRange['step'] )                                           +                                                         for t in range( 0, len(values[s]) - 1):
-                                                                                                                                                                                                                                                                        for s in range( 0, xIters ):                                                                                 +
-                                                                                                                                                                                                                                                                        for t in range( 0, gIters ):                                                                                 +
roots = []                                                                                                                                                             roots = []
#Check the bottom                                                                                                                                                      #Check the bottom
Line 124:                                                                                                                                                                                                                                                                                                                                                                             Line 122:
g += plot( lineSegements, rgbcolor=(0,0,0) )                                                                                                                           g += plot( lineSegements, rgbcolor=(0,0,0) )
return g                                                                                                                                                               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: + + [[Image: Mod_implicit_plot.png]] + + === Benchmarks === + + Function: [itex]-y^2 + x^3 - x[/itex] + + XRange: (-3,3) + + YRange: (-3,3) + + Step Size: 0.1 + + In Sage the average execution time was '''10.563545155525208''' seconds + + In Mathematica the average execution time was '''0.359723''' seconds

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

### Benchmarks

Function: - y2 + x3 - x

XRange: (-3,3)

YRange: (-3,3)

Step Size: 0.1

In Sage the average execution time was 10.563545155525208 seconds

In Mathematica the average execution time was 0.359723 seconds