Go to Qiita Advent Calendar 2024 Top
LoginSignup
5
3

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 5 years have passed since last update.

@yamasushi

Juliaで「グモウスキー・ミラの写像」を描いた

Last updated at Posted at 2019-04-22

グモウスキー・ミラの写像をJulia(version 1.1.x)で描きます。

初版

イテレータについては[お気楽 Julia プログラミング超入門 Julia の基礎知識]
(http://www.nct9.ne.jp/m_hiroi/light/juliaa02.html)を参考にしました。

gm.jl
struct DynamicalSystem
    initial::Array{Float64,1}
    map::Function
end
Base.iterate(ds::DynamicalSystem, pt=ds.initial ) = (pt, ds.map(pt) )

function push_points(initial::Array{Float64,1} ,
                    map::Function ,
                    drop::Int64  ,
                    num::Int64 ,
                    plot_x::Array{Float16,1} , plot_y::Array{Float16,1})
    for pt in Iterators.take(
                Iterators.drop(DynamicalSystem(initial,map ),drop ) , num)
        push!(plot_x,pt[1])
        push!(plot_y,pt[2])
    end
end

@enum GMType f1g1 f1g2 f2g1 f2g2

struct GumowskiMira
    gmtype::GMType
    α::Float64
    σ::Float64
    μ::Float64

    # g1,g2 : 有界非線形性 (bounded nonlinerity)
    g1::Function
    g2::Function

    f1::Function
    f2::Function

    function GumowskiMira(gmtype::GMType ;α::Float64=0. ,σ::Float64=0. ,μ::Float64)
        g1=x->μ*x + (2.0(1.0-μ)x^2 )/(1.0 + x^2)
        g2=x->μ*x + (1.0-μ)x^2*exp((1.0 - x^2)/4.0)
        f1=g->function(pt::Array{Float64,1})
            xx =  pt[2] + g(pt[1])
            yy = -pt[1] + g(xx)
            [xx,yy]
        end
        f2=g->function(pt::Array{Float64,1})
            xx =  pt[2] + α*pt[2]*( 1.0 - σ*pt[2]^2 ) + g(pt[1])
            yy = -pt[1] + g(xx)
            [xx,yy]
        end
        new(gmtype,α,σ,μ,g1,g2,f1,f2)
    end

end

gmmap(gm::GumowskiMira)=
                    if      (gm.gmtype==f1g1)
                        gm.f1(gm.g1)
                    elseif (gm.gmtype==f1g2)
                        gm.f1(gm.g2)
                    elseif (gm.gmtype==f2g1)
                        gm.f2(gm.g1)
                    elseif (gm.gmtype==f2g2)
                        gm.f2(gm.g2)
                    else
                        error("gmtypeが不正")
                    end


#gm = GumowskiMira(f2g1,α=0.008 , σ=0.05, μ=-0.496) # 「神話の鳥 (mythic bird)」3枚羽の翼
gm = GumowskiMira(f2g1,α=0.009 , σ=0.05, μ=-0.801) # 「神話の鳥 (mythic bird)」5枚羽の翼
#gm = GumowskiMira(f2g2,α=0.0083, σ=0.1 , μ=-0.38 ) #
#gm = GumowskiMira(f2g2,α=0.01  , σ=0.1 , μ=0.8 ) #
#gm = GumowskiMira(f1g1,μ=0.39 ) #
#gm = GumowskiMira(f1g1,μ=0.365 ) #
#gm = GumowskiMira(f1g1,μ=0.34 ) #

initial = [1.1,1.1]
plot_x=Float16[]
plot_y=Float16[]
push_points(initial,gmmap(gm),20000,200000,plot_x,plot_y)

using Plots
pyplot()
scatter(plot_x,plot_y,markersize=1,markercolor=:black)

Figure_1.png

第二版(色をつける、など)

  • moduleに分けました
  • 型パラメータを使用しました
  • パラメータの指定に名前付きタプルを使用しました
  • 色つきとモノクロの二種類に対応しました(with_rgbパラメータ)
DynamicalSystem.jl
module DynamicalSystem

using ColorTypes

export make_seq , make_initial_seq , push_points!
struct Seq{P}
    initial::P
    map::Function
end
Base.iterate(ds::Seq, pt=ds.initial ) = (pt, ds.map(pt) )

function make_seq(initial::P ,
                    map::Function ;
                    ndrop::Int64 = 0,
                    ntake::Int64 = 0) where{P} 
    Iterators.take(
        Iterators.drop(Seq{P}(initial,map ),ndrop ) , ntake)
end

function make_initial_seq(;
            density::Float64,radius::Float64=1.0) ::Base.Generator
    length =2*Int(round(density*radius))
    range_x=range(-radius,radius,length=length)
    range_y=range(-radius,radius,length=length)
    ((x,y) for x in range_x,y in range_y)
end

# push 2d points
function push_points!(plot_x::Array{T,1} ,
                       plot_y::Array{T,1} ;
                       initial_seq::Base.Generator ,
                       map::Function ,
                       ndrop::Int64  ,
                       ntake::Int64  ) where {T,U}
    for initial in initial_seq
        for pt in make_seq(initial,map,ndrop=ndrop,ntake=ntake)
            push!(plot_x,pt[1])
            push!(plot_y,pt[2])
        end
    end
end

# push 2d points with color
function push_points!(plot_x::Array{T,1} ,
                       plot_y::Array{T,1} ,
                       rgbs  ::Array{RGB{U},1} ;
                       make_color::Function ,
                       initial_seq::Base.Generator ,
                       map::Function ,
                       ndrop::Int64  ,
                       ntake::Int64  ) where {T,U}
    for initial in initial_seq
        rgb::RGB{U}=make_color(initial)
        for pt in make_seq(initial,map,ndrop=ndrop,ntake=ntake)
            push!(plot_x,pt[1])
            push!(plot_y,pt[2])
            push!(rgbs,rgb)
        end
    end
end


end # module DynamicalSystem
GumowskiMira.jl
module GumowskiMira
export gmmap,GMType

@enum GMType f1g1 f1g2 f2g1 f2g2

function g1(x::Float64 ; μ::Float64) ::Float64
    μ*x + (2.0(1.0-μ)x^2 )/(1.0 + x^2)
end

function g2(x::Float64 ; μ::Float64) ::Float64
    μ*x + (1.0-μ)x^2*exp((1.0 - x^2)/4.0)
end

function f1(g::Function,pt::Tuple{Float64,Float64};
            μ::Float64) ::Tuple{Float64,Float64}
    xx =  pt[2] + g(pt[1],μ=μ)
    yy = -pt[1] + g(xx   ,μ=μ)
    (xx,yy)
end

function f2(g::Function,pt::Tuple{Float64,Float64};
            α::Float64,σ::Float64,μ::Float64)::Tuple{Float64,Float64}
    xx =  pt[2] + α*pt[2]*( 1.0 - σ*pt[2]^2 ) + g(pt[1],μ=μ)
    yy = -pt[1] + g(xx,μ=μ)
    (xx,yy)
end

function gmmap(gmtype::GMType;
                α::Float64=0.,σ::Float64=0.,μ::Float64) ::Function
    if      (gmtype==f1g1) pt->f1(g1,pt,μ=μ)
    elseif (gmtype==f1g2)  pt->f1(g2,pt,μ=μ)
    elseif (gmtype==f2g1)  pt->f2(g1,pt,α=α,σ=σ,μ=μ)
    elseif (gmtype==f2g2)  pt->f2(g2,pt,α=α,σ=σ,μ=μ)
    else
        error("gmtypeが不正")
    end
end

end # module GumowskiMira
gm.jl
using Colors
using FixedPointNumbers
using Plots

# my modules
import DynamicalSystem
import GumowskiMira
# 

function calc_points(param::NamedTuple)
    
    function make_hsv(initial::Tuple{T,T}) where{T}
        hue= ( hypot(initial...) * 360 ) % 360
        HSV(hue,1,1)
    end
    
    initial_seq=DynamicalSystem.make_initial_seq(
                    density=param.initial_density,
                    radius =param.initial_radius)

    plot_x=Float16[]
    plot_y=Float16[]
    rgbs=RGB{N0f8}[]
    if param.with_rgb
        DynamicalSystem.push_points!(plot_x,plot_y,rgbs,
                        make_color  =make_hsv     ,
                        initial_seq =initial_seq  ,
                        map         =param.map    ,
                        ndrop       =param.ndrop  ,
                        ntake       =param.ntake        )
        param,plot_x,plot_y,rgbs
    else
        DynamicalSystem.push_points!(plot_x,plot_y,
                        initial_seq =initial_seq  ,
                        map         =param.map    ,
                        ndrop       =param.ndrop  ,
                        ntake       =param.ntake        )
        param,plot_x,plot_y,nothing
    end
end

function plot_points(param::NamedTuple ,
                plot_x::Array{T,1}  ,
                plot_y::Array{T,1}  ,
                rgbs::Union{Array{RGB{U},1},Nothing} 
                ) where{ T<:AbstractFloat , U }
    param.backend()
    if !isnothing(rgbs)
        scatter(plot_x,plot_y
                ,background_color=:black
                ,aspect_ratio=:equal
                ,size=param.size
                ,markerstrokewidth=0
                ,markeralpha=param.markeralpha
                ,markersize=param.markersize
                ,markercolor=rgbs)

    else
        scatter(plot_x,plot_y
                ,background_color=:black
                ,aspect_ratio=:equal
                ,size=param.size
                ,markerstrokewidth=0
                ,markeralpha=param.markeralpha
                ,markersize =param.markersize
                ,markercolor=:white)
    end
end

gm_map_1=GumowskiMira.gmmap(GumowskiMira.f2g1,α=0.008 , σ=0.05, μ=-0.496) # 「神話の鳥 (mythic bird)」3枚羽の翼
gm_map_2= GumowskiMira.gmmap(GumowskiMira.f2g1,α=0.009 , σ=0.05, μ=-0.80) #「神話の鳥 (mythic bird)」5枚羽の翼
gm_map_3= GumowskiMira.gmmap(GumowskiMira.f2g2,α=0.0083, σ=0.1 , μ=-0.38 ) #
gm_map_4= GumowskiMira.gmmap(GumowskiMira.f2g2,α=0.01  , σ=0.1 , μ=0.8 ) #
gm_map_5= GumowskiMira.gmmap(GumowskiMira.f1g1,μ=0.39 ) #
gm_map_6= GumowskiMira.gmmap(GumowskiMira.f1g1,μ=0.365) #
gm_map_7= GumowskiMira.gmmap(GumowskiMira.f1g1,μ=0.34 ) #

gm_map=gm_map_2

param_0_0=(
    backend        = pyplot ,
    map            = gm_map ,
    size           = (700,700),
    with_rgb       = true  ,markeralpha    =1 , markersize =1 ,
    initial_density= 5.0   ,initial_radius =3.0 , 
    ndrop=0 , ntake=1 )

param_0_1=(
    backend        = unicodeplots ,
    map            = gm_map ,
    size           = (700,700),
    with_rgb       = false ,markeralpha    =1 , markersize =1 ,
    initial_density= 5.0   ,initial_radius =3.0 , 
    ndrop=1000 , ntake=1 )

param_1=(
    backend        = pyplot ,
    map            = gm_map ,
    size           = (700,700),
    with_rgb       = false  ,markeralpha    =0.5 , markersize =1 ,
    initial_density= 5.0   ,initial_radius =3.0 , 
    ndrop=1000 , ntake=100 )

param_2=(
    backend        = pyplot ,
    map            = gm_map ,
    size           = (700,700),
    with_rgb       = true  ,markeralpha    =0.1 , markersize =1 ,
    initial_density= 5.0   ,initial_radius =3.0 , 
    ndrop=5000 , ntake=1000 )

plot_points(calc_points( param_2 )...)

Figure_1.png

with_rgbをfalseにするとモノクロで描画します
Figure_2.png

反復回数の最低値ndropを0に、取得数ntakeを1にすると初期値を描画します。
Figure_2.png

#第三版（複数描画

gm.jl
using Colors
using FixedPointNumbers
using Plots

# my modules
import DynamicalSystem
import GumowskiMira
#

function calc_points(param::NamedTuple)

    function make_hsv(initial::Tuple{T,T}) where{T}
        hue= ( hypot(initial...) * 360 ) % 360
        HSV(hue,1,1)
    end

    initial_seq=DynamicalSystem.make_initial_seq(
                    density=param.initial_density,
                    radius =param.initial_radius)

    plot_x=Float16[]
    plot_y=Float16[]
    rgbs=RGB{N0f8}[]
    if param.with_rgb
        DynamicalSystem.push_points!(plot_x,plot_y,rgbs,
                        make_color  =make_hsv     ,
                        initial_seq =initial_seq  ,
                        map         =param.map    ,
                        ndrop       =param.ndrop  ,
                        ntake       =param.ntake        )
        param,plot_x,plot_y,rgbs
    else
        DynamicalSystem.push_points!(plot_x,plot_y,
                        initial_seq =initial_seq  ,
                        map         =param.map    ,
                        ndrop       =param.ndrop  ,
                        ntake       =param.ntake        )
        param,plot_x,plot_y,nothing
    end
end

function plot_points(param::NamedTuple ,
                plot_x::Array{T,1}  ,
                plot_y::Array{T,1}  ,
                rgbs::Union{Array{RGB{U},1},Nothing}
                )::Plots.Plot where{ T<:AbstractFloat , U }
    if !isnothing(rgbs)
        scatter(plot_x,plot_y , markercolor=rgbs ,
                markeralpha=param.markeralpha    ,
                title =param.title)

    else
        scatter(plot_x,plot_y ,
                markeralpha=param.markeralpha ,
                title =param.title)
    end
end

param_f2g1_1=(
    title="f2g1:α=0.008,σ=0.05,μ=-0.496" ,
    map=GumowskiMira.gmmap(GumowskiMira.f2g1,
                        α=0.008 , σ=0.05, μ=-0.496) ,
                        # 「神話の鳥 (mythic bird)」3枚羽の翼
    with_rgb       = false  ,markeralpha   =0.3 ,
    initial_density= 10.   ,initial_radius =1. ,
    ndrop=1000 , ntake=50 )
param_f2g1_2=(
    title="f2g1:α=0.009,σ=0.05,μ=-0.80" ,
    map=GumowskiMira.gmmap(GumowskiMira.f2g1 ,
                        α=0.009 , σ=0.05, μ=-0.80) ,
                        #「神話の鳥 (mythic bird)」5枚羽の翼
    with_rgb       = false ,markeralpha    =0.3 ,
    initial_density= 10.   ,initial_radius =1. ,
    ndrop=1000 , ntake=50 )
param_f2g2_1=(
    title="f2g2:α=0.0083,σ=0.1,μ=-0.38" ,
    map=GumowskiMira.gmmap(GumowskiMira.f2g2,
                        α=0.0083, σ=0.1 , μ=-0.38 )  ,
    with_rgb       = false ,markeralpha   =0.3 ,
    initial_density= 10.   ,initial_radius =1. ,
    ndrop=1000 , ntake=50 )
param_f2g2_2=(
    title="f2g2:α=0.01,σ=0.1,μ=0.8" ,
    map=GumowskiMira.gmmap(GumowskiMira.f2g2,
                        α=0.01  , σ=0.1 , μ=0.8 )  ,
    with_rgb       = false ,markeralpha    =0.3 ,
    initial_density= 10.   ,initial_radius =1. ,
    ndrop=1000 , ntake=50 )
param_f1g1_1=(
    title="f1g1:μ=0.39" ,
    map=GumowskiMira.gmmap(GumowskiMira.f1g1,μ=0.39 )  ,
    with_rgb       = false   ,markeralpha   =0.1 ,
    initial_density= 10.0   ,initial_radius =3.0 ,
    ndrop=1000 , ntake=10 )
param_f1g1_2=(
    title="f1g1:μ=0.365" ,
    map=GumowskiMira.gmmap(GumowskiMira.f1g1,μ=0.365)  ,
    with_rgb       = false  ,markeralpha    =0.1 ,
    initial_density= 10.0   ,initial_radius =3.0 ,
    ndrop=1000 , ntake=10 )
param_f1g1_3=(
    title="f1g1:μ=0.34" ,
    map=GumowskiMira.gmmap(GumowskiMira.f1g1,μ=0.34 )  ,
    with_rgb       = false  ,markeralpha    =0.1 ,
    initial_density= 10.0   ,initial_radius =3.0 ,
    ndrop=1000 , ntake=10 )


pyplot()
plot(
    plot(
        plot_points(calc_points( param_f2g1_1 )...) ,
        plot_points(calc_points( param_f2g1_2 )...) ,
        layout=(1,2) 
    ),
    plot(
        plot_points(calc_points( param_f2g2_1 )...) ,
        plot_points(calc_points( param_f2g2_2 )...) ,
        layout=(1,2) 
    ),
    plot(
        plot_points(calc_points( param_f1g1_1 )...) ,
        plot_points(calc_points( param_f1g1_2 )...) ,
        plot_points(calc_points( param_f1g1_3 )...) ,
        layout=(1,3) 
    ),
    layout=(3,1),
    aspect_ratio=:equal,
    markercolor=:black ,
    markerstrokewidth=0 ,
    markersize =1 ,
    size=(700,800)
)

Figure_1.png

参考

5
3
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
5
3

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?