using Plots

the logistic map: first steps into chaos

read Sethna sec 4.3

$$ x_{n+1} = r x_n (1-x_n) $$

• stable when r in 1:3
• oscillation between 2 values when r in 3:4.5
• CRAZY? when r > 3.45

function iterate(x; r=0.5)
    # iterate once
    return r * x*(1-x)
end

iterate (generic function with 1 method)

let
    
    rr = 2.5
    N1 = 10
   
    # different starting points

    x1_arr = zeros(N1)
    x2_arr = zeros(N1)

    x1 = 0.01
    x2 = 0.9

    for i1 = 1:N1

        x1_arr[i1] = x1
        x2_arr[i1] = x2

        x1 = iterate(x1; r = rr)
        x2 = iterate(x2; r = rr)

    end

    scatter(0:N1-1, x1_arr, ylim = (0, 1.1), marker=:square)
    scatter!(0:N1-1, x2_arr, ylim = (0, 1.1), title="the logistic map @ r=$rr", xlabel="iterations")

    # analytic limit
    plot!(0:N1-1, x -> 1 - 1 / rr, label = "1-1/r", color="gray", line=(dash=:dash))

end

let
    
    rr = 3.5
    N1 = 20
   
    # different starting points

    x1_arr = zeros(N1)
    x2_arr = zeros(N1)

    x1 = 0.01
    x2 = 0.9

    for i1 = 1:N1

        x1_arr[i1] = x1
        x2_arr[i1] = x2

        x1 = iterate(x1; r = rr)
        x2 = iterate(x2; r = rr)

    end

    scatter(0:N1-1, x1_arr, ylim = (0, 1.1), marker=:square)
    scatter!(0:N1-1, x2_arr, ylim = (0, 1.1), title="the logistic map @ r=$rr", xlabel="iterations")

    # analytic limit
    plot!(0:N1-1, x -> 1 - 1 / rr, label = "1-1/r", color="gray", line=(dash=:dash))

end

let

    # visual aid to understand end points
    # r = 3.4 has 2 end points

    # solution is at x = f(x) = f(f(x)) = ...
    r = 3.4
    f(x) = r * x * (1 - x)

    rv = range(0.0, 1.0, length=250)

    Plots.plot(rv, x -> x, label = "x",legend=:bottomright)
    Plots.plot!(rv, x -> f(x), label = "f(x)")
    Plots.plot!(rv, x -> f(f(x)), label = "f(f(x))")
    Plots.plot!(rv, x -> f(f(f(x))), label = "f(f(f(x)))")
    Plots.plot!(rv, x -> f(f(f(f(x)))), label = "f(f(f(f(x))))")

end

invariant mass density

use a histogram to collect the set of final points at a given r

let

    # invariant density
    # just plot the number of occurances in a histogram
    rval = 3.64

    # random start
    x = rand()
    x_arr = [x]

    Nmax = 25000
    x_arr = zeros(Nmax)
    
    for i1 = 1:Nmax
        x = iterate(x; r = rval)
        x_arr[i1] = x
    end

    # analytic result for r -> 4
    f(x) = 1 / sqrt(x * (1 - x)) / pi

    histogram(x_arr, normalize = true, bins = 150, label="numerical", legend=:topleft)
    plot!(0:0.01:1, f,
        label="analytic (r=4)",
        title="invariant density @ r = $rval",
        color="gray", line=(dash=:dash)
    )

end

let

    # invariant density
    # just plot the number of occurances in a histogram
    rval = 4.0

    # random start
    x = rand()
    x_arr = [x]

    Nmax = 25000
    x_arr = zeros(Nmax)
    
    for i1 = 1:Nmax
        x = iterate(x; r = rval)
        x_arr[i1] = x
    end

    # analytic result for r -> 4
    f(x) = 1 / sqrt(x * (1 - x)) / pi

    histogram(x_arr, normalize = true, bins = 150, label="numerical", legend=:topleft)
    plot!(0:0.01:1, f,
        label="analytic (r=4)",
        title="invariant density @ r = $rval",
        color="gray", line=(dash=:dash)
    )

end

let

    # for each r, iterates and plot the last N1 points
    N1 = 50

    # given r, iterates and collects
    function track(r1; N=10*N1)
        xarr = zeros(N)

        x = rand()
        for i1 = 1:N
            x = iterate(x; r = r1)
            xarr[i1] = x
        end
        
        return xarr
    end

    r_arr = range(0.9, 4.01, length = 120)


    fig = plot(r_arr, r->1-1/r, label="1-1/r",legend=:topleft,
            color="red", line=(dash=:dash))

    for i1 = 1:N1
        
        # collect the end-i1 - th point from a r-track for each r
        fig = scatter!(r_arr, r -> track(r)[end-i1],
            markersize = 2.0,
            color="black",
            alpha=0.5,
            ylim = (0.0, 1),
            label = "",
            )
        
    end

    fig

end

other map

$$ x_{n+1} = r \sin(\pi x_n) $$

let

    # for each r, iterates and plot the last N1 points
    N1 = 50

    function sin_map(x; r=0.5)
        return r*sin(x*pi)
    end
    
    # given r, iterates and collects
    function track(r1; N=6*N1)
        xarr = zeros(N)

        x = rand()
        for i1 = 1:N
            x = sin_map(x; r = r1)
            xarr[i1] = x
        end
        
        return xarr
    end

    r_arr = range(0.0, 1.0, length = 120)

    
    fig = scatter(r_arr, r-> track(r)[end], 
        markersize=2.6, color=:red, alpha=1.0, label="")

    ana = r -> r > 1/pi ? sqrt((r*pi-1)*6/pi^3/r) : Inf

    fig = plot(r_arr, ana, label="small x",legend=:topleft,
        color="red", line=(dash=:dash))
    

    for i1 = 1:N1
        
        # collect the end-i1 - th point from a r-track for each r
        fig = scatter!(r_arr, r -> track(r)[end-i1],
            markersize = 2.0,
            color="black",
            alpha=0.5,
            ylim = (0.0, 1),
            label = "",
            )
        
    end

    fig

end

further reading

• runs this with more points offline
• read more in J.F. Boudreau and E.S. Swanson, Applied Computational Physics (ch.13)
• article: simple math models

iteration of matrix

let

    using LinearAlgebra

    # for each r, iterates and plot the last N1 points
    N1 = 70

    # iterate a matrix, collect its trace
    function sin_map(x; r = 0.5)
        # iterate once
        return r * sin(x * pi)
    end

    # given r, iterates and collects
    function track(r1; N=40*N1)
        xarr = zeros(N)

        x = rand(2,2)
        for i1 = 1:N
            x = sin_map(x; r = r1)
            # tr or det
            xarr[i1] = tr(x)
        end
        
        return xarr
    end

    r_arr = range(0.5, 4.25, length = 120)


    fig = scatter(r_arr, r-> track(r)[end], 
        markersize=2.6, color=:red, alpha=1.0,label="")
    
    for i1 = 1:N1
        
        # collect the end-i1 - th point from a r-track for each r
        fig = scatter!(r_arr, r -> track(r)[end-i1],
            markersize = 1.2,
            color="black",
            alpha=0.5,
            label = "",
            )
        
    end

    fig

end