Skip to content

Commit 66233e9

Browse files
committed
Added plot recipe
1 parent 0f85564 commit 66233e9

File tree

2 files changed

+47
-1
lines changed

2 files changed

+47
-1
lines changed

src/PlotRecipes.jl

Lines changed: 10 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,10 @@
1+
using RecipesBase
2+
3+
@recipe function poly_recipe(p::Poly, xs = polyinterval(p))
4+
buff = IOBuffer()
5+
printpoly(buff, p)
6+
label --> String(take!(buff))
7+
8+
ys = map(x -> polyval(p, x), xs)
9+
xs, ys
10+
end

src/Polynomials.jl

Lines changed: 37 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -5,7 +5,7 @@ module Polynomials
55

66
export Poly, poly
77
export degree, coeffs, variable, printpoly
8-
export polyval, polyint, polyder, roots, polyfit
8+
export polyval, polyint, polyder, roots, polyfit, polyinterval
99
export Pade, padeval
1010

1111
import Base: length, size, eltype, collect, eachindex, lastindex
@@ -716,7 +716,43 @@ function polyfit(x, y, n::Int=length(x)-1, sym::Symbol=:x)
716716
end
717717
polyfit(x,y,sym::Symbol) = polyfit(x,y,length(x)-1, sym)
718718

719+
720+
"""
721+
polyinterval(p)
722+
723+
Choose a sensible interval containing points of interest for the polynomial `p`.
724+
725+
The interval contains some margin either side to allow for plotting.
726+
727+
# Examples
728+
729+
```julia-repl
730+
julia> polyinterval(Poly([0, 0, 1])) # contains single critical point at x=0
731+
-1.0:0.01:1.0
732+
733+
julia> polyinterval(Poly([1, 0, 3, -1])) # points of interest between 0 and 3.1
734+
-3.103803402735535:0.031038034027355346:6.20760680547107
735+
```
736+
"""
737+
function polyinterval(p :: Poly)
738+
# Find points of interest
739+
zero_pts = roots(p)
740+
crit_pts = roots(polyder(p, 1))
741+
infl_pts = roots(polyder(p, 2))
742+
pts = sort([ real(pt) for pt in [zero_pts; crit_pts; infl_pts] if isreal(pt) ])
743+
744+
# Choose a range that shows all interesting points with some margin
745+
min_x, max_x = length(pts) > 0 ? (pts[1], pts[end]) : (-1, 1)
746+
diff = max(max_x - min_x, 1)
747+
a = min_x - diff
748+
b = max_x + diff
749+
750+
return a : diff/100 : b
751+
end
752+
753+
719754
### Pull in others
720755
include("pade.jl")
756+
include("PlotRecipes.jl")
721757

722758
end # module Poly

0 commit comments

Comments
 (0)