Notebook

1 Plots.plot による関数のプロット
- 1.1 一変数関数
  - 1.1.1 きれいに書ける例
  - 1.1.2 カクカクしてしまう例
- 1.2 パラメトリックな二次元曲線
  - 1.2.1 うまく行く例
  - 1.2.2 カクカクしてしまう例
2 Plots.plot のソースファイルでカクカクの原因を確認
- 2.1 一次元関数の場合
- 2.2 パラメトリックな２次元曲線について
3 より柔軟な関数プロット手順を考える
4 気になる点

Plots.plot による関数のプロット¶

Plots.plot に関数と範囲を渡すとグラフが書けるのだけれど、結果がカクカクしてしまうことがある。

In [1]:

versioninfo()

Julia Version 0.6.2
Commit d386e40 (2017-12-13 18:08 UTC)
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: Intel(R) Core(TM) i7-6600U CPU @ 2.60GHz
  WORD_SIZE: 64
  BLAS: libopenblas (USE64BITINT DYNAMIC_ARCH NO_AFFINITY Haswell)
  LAPACK: libopenblas64_
  LIBM: libopenlibm
  LLVM: libLLVM-3.9.1 (ORCJIT, skylake)

In [2]:

Pkg.status("Plots")

 - Plots                         0.16.0+            master

In [3]:

using Plots
pyplot()

Out[3]:

Plots.PyPlotBackend()

一変数関数¶

きれいに書ける例¶

In [4]:

plot(x->abs(sin(x)), 0, 2pi, marker=:circle, legend=:none)

Out[4]:

曲率の大きなところにサンプル点を多数配置して、カクカクしづらくなるよう工夫がされている

カクカクしてしまう例¶

In [5]:

plot(x->abs(sin(x)) + 10, 0, 4pi, marker=:circle, legend=:none)

Out[5]:

同じ形の関数でも、縦軸がゼロから離れているとカクカクしてしまう。

パラメトリックな二次元曲線¶

うまく行く例¶

In [6]:

plot(t->cos(t), t->sin(t), 0, 2pi, aspectratio=:equal, marker=:circle, legend=:none)

Out[6]:

カクカクしてしまう例¶

In [7]:

plot(t->cos(t), t->sin(t), 0, 2pi, marker=:circle, legend=:none, xlims=(-0.1,0.1), ylims=(0.995,1.001))

Out[7]:

一部を拡大するとカクカクしてしまう

In [8]:

plot(t->cos(14t)cos(t), t->cos(15t)sin(t), 0, 2pi, aspectratio=:equal, marker=:circle, legend=:none)

Out[8]:

複雑な関数でもサンプル点数が変わらないためカクカクになる

Plots.plot のソースファイルでカクカクの原因を確認¶

このあたりの @recipe で関数からサンプル点を求めていると twitter で黒木さんから教わった。https://twitter.com/genkuroki/status/984086197164060672

https://github.com/JuliaPlots/Plots.jl/blob/1ed7899296f668c644e41a5de02701db2fffc2ac/src/series.jl#L512-L521

一次元関数の場合¶

対応するソースコード：

@recipe function f(f::Function, xmin::Number, xmax::Number)
    xs = adapted_grid(f, (xmin, xmax))
    xs, f
end

adapted_grid でサンプル点の x 座標を求めて、それらの点の f の値をプロットしている。

adapted_grid のソースはこちら：

https://github.com/JuliaPlots/PlotUtils.jl/blob/fb8404578c21834a48d5984b0ea8e260b5052cd0/src/adapted_grid.jl#L15-L21

曲率の大きなところに多数のサンプル点を取るように、最適化されたサンプル点を求めている。

サンプル点数を増やしてカクカクしないようにするには max_recursions や max_curvature を Plots.plot から指定できたら良かったのだけれど、現時点ではそういうことはできない。

また、adapted_grid には以下の問題がある：

curvatures を求める際に使っている max_f が「関数値の絶対値の最大値」なので、同じ形のグラフでも関数値がゼロから離れているだけでカクカクしやすくなってしまう
adapted_grid でサンプル点を決めるために関数値を求めているのに、グラフを書く際にはその値を使わずもう一度関数を呼ぶので効率が悪い
境界値である xmin, xmax にサンプル点が取られない

パラメトリックな２次元曲線について¶

対応するソースコード:

@recipe f(fx::FuncOrFuncs{F}, fy::FuncOrFuncs{G}, umin::Number, umax::Number, n = 200) where {F<:Function,G<:Function} = fx, fy, linspace(umin, umax, n)

umin から umax を n 等分し、それらに対応する点をプロットしている。

n のデフォルト値は 200 だが、これを増やせば、よりスムーズな曲線を描ける。

In [9]:

plot(t->cos(14t)cos(t), t->cos(15t)sin(t), 0, 2pi, 600, aspectratio=:equal, marker=:circle, legend=:none)

Out[9]:

したがって、単にスムーズな曲線を描きたければ n を増やすだけで良いのだが、現状のコードには以下の問題がある：

曲率の大きなところも小さなところも同じ間隔で計算するため効率が悪い
サンプル点数のデフォルト値が 200 固定なので、その都度 n を指定してやる必要がある

より柔軟な関数プロット手順を考える¶

プロットされる範囲を意識して曲率を評価する
パラメトリック曲線の描画でも曲率を元にしてサンプル点を決める
サンプル点を決める際に求めた関数値を、プロットする際に再利用する
境界点 xmin, xmax にサンプル点を配置する

を指針とした。

ソースコード¶

アルゴリズムは以下の通り：

始め initialdivision+1 個のサンプリング点を等間隔に配置する(initialdivision 個のセクションに分割される)
３つずつ点を見ていって、折れ曲がり量が maxbending を越えれば、それらの点に挟まれる２つのセクションを2等分するようにサンプル点を増やす
2.の操作を maxrecursion 回繰り返す

この方針により、

すべてのセクションで maxbending により指定されるより小さな折れ曲がりが実現されるか
セクションの幅が (pmax-pmin)/initialdivision/2^(maxrecursion) となるまで

サンプル点が増やされる

折れ曲がり量の評価は、隣り合う３点 $\boldsymbol{p}_1,\boldsymbol{p}_2,\boldsymbol{p}_3$ を繋ぐ２つのベクトル $\boldsymbol{v}_1 = \boldsymbol{p}_2-\boldsymbol{p}_1, \boldsymbol{v}_2 = \boldsymbol{p}_3 - \boldsymbol{p}_2$ をのなす角を $\delta$ とすると、

$\begin{align} \sqrt{2(\left|\boldsymbol{v}_1\right|\,\left|\boldsymbol{v}_2\right|-\boldsymbol{v}_1\cdot\boldsymbol{v}_2)} &=\sqrt{\left|\boldsymbol{v}_1\right|\,\left|\boldsymbol{v}_2\right|}\sqrt{2(1-\cos\delta)}\\ &=\sqrt{\left|\boldsymbol{v}_1\right|\,\left|\boldsymbol{v}_2\right|}\sqrt{2(\frac{1}{2}\delta^2-O(\delta^4))}\\ &\sim\underbrace{\sqrt{\left|\boldsymbol{v}_1\right|\,\left|\boldsymbol{v}_2\right|}}_\mathrm{geometric\ mean} \delta\\ \end{align}$

として評価される。折れ曲がり角に長さの相乗平均を掛けてあるのは、同じ曲がり角でも長いベクトル同士の折れ曲りはより目に付きやすいためだ。

また、 $\boldsymbol{p}_1,\boldsymbol{p}_2,\boldsymbol{p}_3$ は絶対座標ではなく、表示範囲に対する相対座標で取る。なぜなら、例えば縦軸の長さが1,000、横軸の長さが0.001のグラフでは、横方向の $10^{-4}$ の違いは重要だが、縦方向の $0.1$ の違いは無視できる。このような違いを考慮するため、絶対座標をグラフ表示範囲で割った相対座標で折れ曲りを評価する。

In [10]:

# Copyright (c) 2018 Osamu Takeuchi <osamu@big.jp>
# Released under the MIT license

module FuncPlot

    using Plots

    ######## structure to represent a sample point

    mutable struct Point
        p::Number                                   # parameter
        v                                           # value corresponds to the parameter p
        next::Int                                   # index of the next point in the points array
        shoulddivide::Bool                          # sould divide this segment (from this point to next point)
        function Point(p::Number, value, next::Int)
            return new(p, value, next, false)
        end
    end
    hasnext(point::Point) = point.next > 0
    next(points::Array{Point, 1}, point::Point) = points[point.next]


    ######## calculate sample points

    # calcvalue(p)
    #   should return a single point from the parameter p
    #
    # pmin, pmax
    #   defines parameter range
    #
    # initialdivision
    #   the region from pmin to pmax is initially divided into (initialdivision) segments
    #
    # maxrecursion
    #   the curve is scanned by the function `shoulddivide` (maxrecursion) times
    #
    # shoulddivide(p1::Point, p2::Point, p3::Point)
    #   should return if the segment from p1 to p3 should be divided
    #   if it returns true, the two segements [p1, p2] and [p2, p3] are divided into halves
    #
    # returns calculated points in the form `Array{Point, 1}`

    function adaptedpoints(calcvalue::Function, pmin::Number, pmax::Number, 
                            initialdivision::Int, maxrecursion::Int, shoulddivide::Function)

        # prepare for equally separated (initialdivision+1) points

        points = Point[]
        pdelta = (pmax-pmin)/initialdivision
        for i in 0:initialdivision
            p = (i == initialdivision) ? pmax : pmin + pdelta * i # 末尾を pmax に合わせる
            push!(points, Point(p, calcvalue(p), length(points)+2))
        end
        last(points).next = -1

        for recursion in 1:maxrecursion
            shouldbreak = true

            p1 = first(points)
            p2 = next(points, p1)
            while hasnext(p2) > 0
                p3 = next(points, p2)
                if shoulddivide(p1, p2, p3)
                    p1.shoulddivide = true
                    p2.shoulddivide = true
                    shouldbreak = false
                end
                p1 = p2
                p2 = p3
            end

            if shouldbreak
                break
            end

            p1 = first(points)
            while hasnext(p1) > 0
                p2 = next(points, p1)
                if p1.shoulddivide
                    p = ( p1.p + p2.p ) / 2
                    push!(points, Point(p, calcvalue(p), p1.next))
                    p1.next = length(points)
                    p1.shoulddivide = false
                end
                p1 = p2
            end
        end

        points
    end


    ######## plot a function with a single argument

    function funcplot!(f::Function, xmin::Number, xmax::Number; 
                        initialdivision = 21, maxrecursion = 10, maxbending = 0.01, 
                        xlims=false, ylims=false, other...)

        # we need plot region to evaluate amount of bending
        xwidth = abs( xlims==false ? xmax-xmin : xlims[2]-xlims[1] )
        ymin = ylims==false ?  Inf : ylims[1]
        ymax = ylims==false ? -Inf : ylims[2]

        function calcvalue(x)
            try
                y = f(x)        # may cause exception
                if ylims==false # we need plot region to evaluate amount of bending
                    ymin = min(ymin, y)
                    ymax = max(ymax, y)
                end
            catch
                y = NaN
            end
            return y
        end

        function shoulddivide(p1, p2, p3)
            if !isfinite(p1.v) || !isfinite(p2.v) || !isfinite(p3.v)
                return true
            end
        
            # amount of bending is evaluated relatively to the plot region
            ywidth = abs(ymax-ymin)
            xw = xwidth != 0 ? xwidth :
                 ywidth != 0 ? ywidth/10 : 1
            yw = ywidth != 0 ? ywidth :
                 xwidth != 0 ? xwidth/10 : 1
            dx1 = (p2.p-p1.p)/xw
            dy1 = (p2.v-p1.v)/yw
            d1 = sqrt( dx1^2 + dy1^2 )
            dx2 = (p3.p-p2.p)/xw
            dy2 = (p3.v-p2.v)/yw
            d2 = sqrt( dx2^2 + dy2^2 )

            # sqrt(2(|v1| |v2| - v1⋅v2)) = sqrt(|v1| |v2|) sqrt( 2(1 - cosδ) ) ∼ gmean(|v1|, |v2|)δ
            # here `gmena` represents geometric mean
            # `abs` is needed to avoid domain error caused by calculation errors
            sqrt(2abs(d1*d2-(dx1*dx2+dy1*dy2))) > maxbending
        end

        points = adaptedpoints(calcvalue, xmin, xmax, 
                                initialdivision, maxrecursion, shoulddivide)

        # extract coordinates from points
        x = Float64[]
        y = Float64[]
        point = first(points)
        while true
            push!(x, point.p)
            push!(y, point.v)
            if !hasnext(point)
                break
            end
            point = next(points, point)
        end
    
        if xlims!=false
            push!(other, (:xlims, xlims))
        end
        if ylims!=false
            push!(other, (:ylims, ylims))
        end
        plot!(x, y; other...)
    end

    function funcplot(f::Function, xmin::Number, xmax::Number; options...)
        plot()
        funcplot!(f, xmin, xmax; options...)
    end


    ######## plot a 2D parametric curve

    function funcplot!(fx::Function, fy::Function, pmin::Number, pmax::Number; 
                        initialdivision = 21, maxrecursion = 10, maxbending = 0.01, 
                        xlims=false, ylims=false, other...)

        # we need plot region to evaluate amount of bending
        xmin = xlims==false ?  Inf : xlims[1]
        xmax = xlims==false ? -Inf : xlims[2]
        ymin = ylims==false ?  Inf : ylims[1]
        ymax = ylims==false ? -Inf : ylims[2]

        function calcvalue(p)
            try
                x = fx(p)       # may cause exception
                y = fy(p)
                if xlims==false # we need plot region to evaluate amount of bending
                    xmin = min(xmin, x)
                    xmax = max(xmax, x)
                end
                if ylims==false
                    ymin = min(ymin, y)
                    ymax = max(ymax, y)
                end
            catch
                x = NaN
                y = NaN
            end
            return (x, y)
        end

        function shoulddivide(p1, p2, p3)
            if !isfinite(p1.v[1]) || !isfinite(p1.v[2]) ||
                    !isfinite(p2.v[1]) || !isfinite(p2.v[2]) ||
                    !isfinite(p3.v[1]) || !isfinite(p3.v[2]) 
                return true
            end
        
            # amount of bending is evaluated relatively to the plot region
            xwidth = abs(xmax-xmin)
            ywidth = abs(ymax-ymin)
            xw = xwidth != 0 ? xwidth :
                 ywidth != 0 ? ywidth/10 : 1
            yw = ywidth != 0 ? ywidth :
                 xwidth != 0 ? xwidth/10 : 1
            dx1 = (p2.v[1]-p1.v[1])/xw
            dy1 = (p2.v[2]-p1.v[2])/yw
            d1 = sqrt( dx1^2 + dy1^2 )
            dx2 = (p3.v[1]-p2.v[1])/xw
            dy2 = (p3.v[2]-p2.v[2])/yw
            d2 = sqrt( dx2^2 + dy2^2 )

            # sqrt(2(|v1| |v2| - v1⋅v2)) = sqrt(|v1| |v2|) sqrt( 2(1 - cosδ) ) ∼ gmean(|v1|, |v2|)δ
            # here `gmena` represents geometric mean
            # `abs` is needed to avoid domain error caused by calculation errors
            sqrt(2abs(d1*d2-(dx1*dx2+dy1*dy2))) > maxbending
        end

        points = adaptedpoints(calcvalue, pmin, pmax, 
                                initialdivision, maxrecursion, shoulddivide)

        # extract coordinates from points
        x = Float64[]
        y = Float64[]
        point = first(points)
        while true
            push!(x, point.v[1])
            push!(y, point.v[2])
            if !hasnext(point)
                break
            end
            point = next(points, point)
        end
    
        if xlims!=false
            push!(other, (:xlims, xlims))
        end
        if ylims!=false
            push!(other, (:ylims, ylims))
        end
        plot!(x, y; other...)
    end

    function funcplot(fx::Function, fy::Function, pmin::Number, pmax::Number; options...)
        plot()
        funcplot!(fx, fy, pmin, pmax; options...)
    end


    ######## plot a 3D parametric curve

    function funcplot!(fx::Function, fy::Function, fz::Function, pmin::Number, pmax::Number; 
                        initialdivision = 21, maxrecursion = 10, maxbending = 0.01, 
                        xlims=false, ylims=false, zlims=false, other...)

        # we need plot region to evaluate amount of bending
        xmin = xlims==false ?  Inf : xlims[1]
        xmax = xlims==false ? -Inf : xlims[2]
        ymin = ylims==false ?  Inf : ylims[1]
        ymax = ylims==false ? -Inf : ylims[2]
        zmin = zlims==false ?  Inf : zlims[1]
        zmax = zlims==false ? -Inf : zlims[2]

        function calcvalue(p)
            try
                x = fx(p)       # may cause exception
                y = fy(p)
                z = fz(p)
                if xlims==false # we need plot region to evaluate amount of bending
                    xmin = min(xmin, x)
                    xmax = max(xmax, x)
                end
                if ylims==false
                    ymin = min(ymin, y)
                    ymax = max(ymax, y)
                end
                if zlims==false
                    zmin = min(zmin, z)
                    zmax = max(zmax, z)
                end
            catch
                x = NaN
                y = NaN
                z = NaN
            end
            return (x, y, z)
        end

        function shoulddivide(p1, p2, p3)
            if !isfinite(p1.v[1]) || !isfinite(p1.v[2]) || !isfinite(p1.v[3]) ||
                    !isfinite(p2.v[1]) || !isfinite(p2.v[2]) || !isfinite(p2.v[3]) ||
                    !isfinite(p3.v[1]) || !isfinite(p3.v[2]) || !isfinite(p3.v[3]) 
                return true
            end
        
            # amount of bending is evaluated relatively to the plot region
            xwidth = abs(xmax-xmin)
            ywidth = abs(ymax-ymin)
            zwidth = abs(zmax-zmin)
            xw = xwidth != 0 ? xwidth :
                 ywidth != 0 ? ywidth/10 :
                 zwidth != 0 ? zwidth/10 : 1
            yw = ywidth != 0 ? ywidth :
                 zwidth != 0 ? zwidth/10 : 
                 xwidth != 0 ? xwidth/10 : 1
            zw = zwidth != 0 ? zwidth :
                 ywidth != 0 ? ywidth/10 : 
                 xwidth != 0 ? xwidth/10 : 1
            dx1 = (p2.v[1]-p1.v[1])/xw
            dy1 = (p2.v[2]-p1.v[2])/yw
            dz1 = (p2.v[3]-p1.v[3])/zw
            d1 = sqrt( dx1^2 + dy1^2 + dz1^2 )
            dx2 = (p3.v[1]-p2.v[1])/xw
            dy2 = (p3.v[2]-p2.v[2])/yw
            dz2 = (p3.v[3]-p2.v[3])/zw
            d2 = sqrt( dx2^2 + dy2^2 + dz2^2 )

            # sqrt(2(|v1| |v2| - v1⋅v2)) = sqrt(|v1| |v2|) sqrt( 2(1 - cosδ) ) ∼ gmean(|v1|, |v2|)δ
            # here `gmena` represents geometric mean
            # `abs` is needed to avoid domain error caused by calculation errors
            sqrt(2abs(d1*d2-(dx1*dx2+dy1*dy2+dz1*dz2))) > maxbending
        end

        points = adaptedpoints(calcvalue, pmin, pmax, 
                                initialdivision, maxrecursion, shoulddivide)

        # extract coordinates from points
        x = Float64[]
        y = Float64[]
        z = Float64[]
        point = first(points)
        while true
            push!(x, point.v[1])
            push!(y, point.v[2])
            push!(z, point.v[3])
            if !hasnext(point)
                break
            end
            point = next(points, point)
        end
    
        if xlims!=false
            push!(other, (:xlims, xlims))
        end
        if ylims!=false
            push!(other, (:ylims, ylims))
        end
        if zlims!=false
            push!(other, (:zlims, zlims))
        end
        plot!(x, y, z; other...)
    end

    function funcplot(fx::Function, fy::Function, fz::Function, pmin::Number, pmax::Number; options...)
        plot()
        funcplot!(fx, fy, fz, pmin, pmax; options...)
    end

end

Out[10]:

FuncPlot

一変数関数¶

In [11]:

plot(
FuncPlot.funcplot(x->abs(sin(x)), 0, 2pi, title="funcplot: |sin(x)|"),
plot(x->abs(sin(x)), 0, 2pi, title="Plots.plot: |sin(x)|"),
size=(800,300), marker=:circle, legend=:none
)

Out[11]:

Plots.plot に比べて場所によるサンプル点間隔の変化が大きいみたいだけれど、仕上がりは悪くない感じ？

In [12]:

plot(
FuncPlot.funcplot(x->abs(sin(x)) + 10, 0, 2pi, title="funcplot: |sin(x)|+10"),
plot(x->abs(sin(x)) + 10, 0, 2pi, title="Plots.plot: |sin(x)|+10"),
size=(800,300), marker=:circle, legend=:none
)

Out[12]:

縦軸がゼロから離れていてもサンプル点間隔が変わらないことを確認できる

In [13]:

plot(
    FuncPlot.funcplot(x->1/(1+exp(-40x))-1, -1, 1, title="funcplot"),
    plot(x->1/(1+exp(-40x))-1, -1, 1, title="Plots.plot"),
    size=(800,300), marker=:circle, legend=:none
)

Out[13]:

曲率の小さな部分ではかなり手を抜く

In [14]:

plot(
    FuncPlot.funcplot(x->abs(4-abs(x)), -10, 10, title="funcplot"),
    plot(x->abs(4-abs(x)), -10, 10, title="Plots.plot"),
    size=(800,300), marker=:circle, legend=:none
)

Out[14]:

尖点を正しく表示できる

In [15]:

plot(
    FuncPlot.funcplot(x->log(x), 0, 1, title="funcplot: log(x)"),
    plot(x->log(x), 0, 1, title="Plots.plot: log(x)"),
    size=(800,300), marker=:circle, legend=:none
)

Out[15]:

領域の端で発散してしまう場合でも、それなりの結果が得られている

In [16]:

plot(
    FuncPlot.funcplot(x->floor(sqrt(x)), 0, 100, title="funcplot: floor(sqrt(x))"),
    plot(x->floor(sqrt(x)), 0, 100, title="Plots.plot: floor(sqrt(x))"),
    size=(800,300), marker=:circle, legend=:none
)

Out[16]:

不連続点の周りもまともに表示される。最終点が正しく表示される。

In [17]:

plot(
    FuncPlot.funcplot(x->gamma(x), -4, 1, ylims=(-100, 100), title="funcplot: gamma(x)"),
    plot(x->gamma(x), -4, 1, ylims=(-100, 100), title="Plots.plot: gamma(x)"),
    size=(800,300), marker=:circle, legend=:none
)

Out[17]:

無限大に飛んでいても大丈夫

In [18]:

plot(
    FuncPlot.funcplot(x->sin(cos(x)*exp(-x/2)), -10, 10, title="funcplot"),
    plot(x->sin(cos(x)*exp(-x/2)), -10, 10, title="Plots.plot"),
    size=(800,300), marker=:circle, legend=:none
)

Out[18]:

In [19]:

plot(
    FuncPlot.funcplot(x->sin(cos(x)*exp(-x/2)), -10, 10, title="funcplot"),
    plot(x->sin(cos(x)*exp(-x/2)), -10, 10, title="Plots.plot"),
    size=(800,300), legend=:none
)

Out[19]:

ちゃんと上下端が -1 や 1 に揃っている

2Dパラメトリック曲線¶

In [20]:

plot(
FuncPlot.funcplot(t->cos(t)t, t->sin(t)t, 0, 2pi, title="funcplot"),
plot(t->cos(t)t, t->sin(t)t, 0, 2pi, title="Plots.plot"),
aspect_ratio=:equal,legend=:none, marker=:circle,
size=(1000,300)
)

Out[20]:

曲率に合わせて点数を決めるため、単純な形状のグラフでは点数を少なくできる

In [21]:

plot(
FuncPlot.funcplot(t->cos(t), t->sin(4t), 0, 2pi, title="funcplot"),
plot(t->cos(t), t->sin(4t), 0, 2pi, title="Plots.plot"),
aspect_ratio=:equal,legend=:none, marker=:circle,
size=(1000,500)
)

Out[21]:

曲率の大きなところにより多くの点を配置する

In [22]:

a, b=7,5
plot(
FuncPlot.funcplot(t->cos(a*t)cos(t), t->cos(b*t)sin(t), 0, 2pi, title="funcplot: ("*string(Int(a))*", "*string(Int(b))*")"),
plot(t->cos(a*t)cos(t), t->cos(b*t)sin(t), 0, 2pi, title="Plots.plot: ("*string(Int(a))*", "*string(Int(b))*")"),
aspect_ratio=:equal,legend=:none, marker=:circle,
size=(1000,500)
)

Out[22]:

曲率の小さなところは積極的に手抜きする

In [23]:

plot(
    FuncPlot.funcplot(t->cos(t), t->sin(t), 0, 2pi, xlims=(-0.1,0.1), ylims=(0.995,1.001)),
    plot(t->cos(t), t->sin(t), 0, 2pi, xlims=(-0.1,0.1), ylims=(0.995,1.001)),
    marker=:circle, legend=:none, size=(800,300)
)

Out[23]:

一部を切り出してもキレイ

In [24]:

plot(
    FuncPlot.funcplot(t->cos(t), t->sin(t), 0, 2pi),
    plot(t->cos(t), t->sin(t), 0, 2pi),
    marker=:circle, legend=:none, size=(800,300), xlims=(-0.1,0.1), ylims=(0.995,1.001)
)

Out[24]:

とはいえ、後から拡大されると対応できないので、xlims, ylims は始めから指定してやって下さい。

In [25]:

plot(
    FuncPlot.funcplot(t->cos(t), t->sin(t), 0, 2pi, maxbending=0.0001),
    plot(t->cos(t), t->sin(t), 0, 2pi),
    marker=:circle, legend=:none, size=(800,300), xlims=(-0.1,0.1), ylims=(0.995,1.001)
)

Out[25]:

あるいは、initialdivision, maxrecursion, maxbending などのキーワード引数を指定して、うまく行くよう調整します。

In [26]:

a, b=14,8
plot(
FuncPlot.funcplot(t->cos(a*t)cos(t), t->cos(b*t)sin(t), 0, 2pi, title="funcplot: ("*string(Int(a))*", "*string(Int(b))*")"),
plot(t->cos(a*t)cos(t), t->cos(b*t)sin(t), 0, 2pi, title="Plots.plot: ("*string(Int(a))*", "*string(Int(b))*")"),
aspect_ratio=:equal,legend=:none, marker=:circle,
size=(1000,500)
)

Out[26]:

In [27]:

a, b=14,8
plot(
FuncPlot.funcplot(t->cos(a*t)cos(t), t->cos(b*t)sin(t), 0, 2pi, title="funcplot: ("*string(Int(a))*", "*string(Int(b))*")"),
plot(t->cos(a*t)cos(t), t->cos(b*t)sin(t), 0, 2pi, title="Plots.plot: ("*string(Int(a))*", "*string(Int(b))*")"),
aspect_ratio=:equal, legend=:none,
size=(1000,500)
)

Out[27]:

複雑な図形になる場合もカクカクしないよう自動的に点数が増やされる

In [28]:

plot([
        FuncPlot.funcplot(t->cos(a*t)cos(t), t->cos(b*t)sin(t), 0, 2pi, 
            aspect_ratio=:equal, size=(800,800), legend=:none, 
            xshowaxis=false, yshowaxis=false, grid=false
        ) for a in 1:10, b in 1:10
]...)

Out[28]:

楽しい！

3Dパラメトリック曲線¶

In [29]:

plot([
    FuncPlot.funcplot(t->cos(t)t,t->sin(t)t,t->t,0,10pi,title="funcplot"),
    plot(t->cos(t)t,t->sin(t)t,t->t,0,10pi,title="Plots.plot")
]..., size=(800,400), marker=:circle)

Out[29]:

In [30]:

a=8; b=7; c=5
plot([
    FuncPlot.funcplot(t->cos(a*t)cos(t),t->cos(b*t)sin(t),t->sin(c*t),0,2pi,title="funcplot"),
    plot(t->cos(a*t)cos(t),t->cos(b*t)sin(t),t->sin(c*t),0,2pi,title="Plots.plot")
]..., size=(800,400), marker=:circle, aspect_ratio=:equal, legend=:none)

Out[30]:

In [31]:

a=8; b=7; c=5
plot([
    FuncPlot.funcplot(t->cos(a*t)cos(t),t->cos(b*t)sin(t),t->sin(c*t),0,2pi,title="funcplot", zlims=(0.75,1)),
    plot(t->cos(a*t)cos(t),t->cos(b*t)sin(t),t->sin(c*t),0,2pi,title="Plots.plot", zlims=(0.75,1))
]..., size=(800,400), marker=:circle, aspect_ratio=:equal, legend=:none)

Out[31]:

拡大すると自動的にサンプリング点を増やす。

気になる点¶

現状では、定数値関数に対して表示範囲がうまく選ばれず、プロットが画面に表示されない。

ylims を手動で指定すれば現れるけれど、これは不親切。

In [32]:

plot([
        FuncPlot.funcplot(x->1, 0, 1, title="funcplot"),
        plot(x->1, 0, 1,title="Plots.plot")
]..., size=(600,300), legend=:none)

Out[32]:

In [33]:

ylims!(0,2)

Out[33]:

これは関数値のプロットに限った問題ではないので、ここれは対策しない。

In [34]:

plot([1,1], legend=:none)

Out[34]:

In [35]:

ylims!(0,2)

Out[35]:

In [ ]:

Table of Contents

Plots.plot による関数のプロット¶

一変数関数¶

きれいに書ける例¶

カクカクしてしまう例¶

パラメトリックな二次元曲線¶

うまく行く例¶

カクカクしてしまう例¶

Plots.plot のソースファイルでカクカクの原因を確認¶

一次元関数の場合¶

パラメトリックな２次元曲線について¶

より柔軟な関数プロット手順を考える¶

ソースコード¶

一変数関数¶

2Dパラメトリック曲線¶

3Dパラメトリック曲線¶

気になる点¶