using Plots, LinearAlgebra

the logistic map

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 = 0.01
    x2 = 0.9

    x1_arr = [x1]
    x2_arr = [x2]

    for i1 = 1:N1-1

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

        push!(x1_arr, x1)
        push!(x2_arr, x2)
        
    end

    scatter(0:N1-1, x1_arr, ylim = (0, 1.1))
    scatter!(0:N1-1, x2_arr, ylim = (0, 1.1), title="the logistic map @ r=$rr", xlabel="no. of iter.")

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



end

let
    
    rr = 3.5
    N1 = 20
   
    # different starting points
    x1 = 0.01
    x2 = 0.9

    x1_arr = [x1]
    x2_arr = [x2]

    for i1 = 1:N1-1

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

        push!(x1_arr, x1)
        push!(x2_arr, x2)
        
    end

    scatter(0:N1-1, x1_arr, ylim = (0, 1.1))
    scatter!(0:N1-1, x2_arr, ylim = (0, 1.1), title="the logistic map @ r=$rr", xlabel="no. of iter.")

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



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)

    Plots.plot(range(0.0, 1.0, length = 250), x -> x, label = "x",legend=:bottomright)
    Plots.plot!(range(0.0, 1.0, length = 250), x -> f(x), label = "f(x)")
    Plots.plot!(range(0.0, 1.0, length = 250), x -> f(f(x)), label = "f(f(x))")
    Plots.plot!(range(0.0, 1.0, length = 250), x -> f(f(f(x))), label = "f(f(f(x)))")
    Plots.plot!(range(0.0, 1.0, length = 250), x -> f(f(f(f(x)))), label = "f(f(f(f(x))))")

end

invariant mass density

A histogram of values in the trajectory for 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]
    for i1 = 1:25000

        x = iterate(x; r = rval)

        push!(x_arr, 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"
    )

end

let
    
    r_arr = range(0.5, 4.1, length = 140)

    function track(r1)
        xarr = []
        x = rand()
        for i1 = 1:250
            x = iterate(x; r = r1)
            push!(xarr, x)
        end
        return xarr
    end

    fig = plot(r_arr, r->1-1/r, label="analytic (small r)",legend=:topleft)

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

    fig

end

other map

$$ \begin{aligned} x_{n+1} = r \, sin(\pi \, x_n) \end{aligned} $$

or

$$ \begin{aligned} x_{n+1} = x_n \, e^{\, r \, (1-x_n)} \end{aligned} $$

have similar feature.

let
    
    function sin_map(x; r = 0.5)
        # iterate once
        return r * sin(x * pi)
    end
    
    r_arr = range(0.0, 1.0, length = 140)

    function track(r1)
        xarr = []
        x = rand()
        for i1 = 1:250
            x = sin_map(x; r = r1)
            push!(xarr, x)
        end
        return xarr
    end

    fig = scatter(r_arr, r->track(r)[end], label="", markersize=0.8)

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

    fig

end

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

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

    function track(r1)
        xarr = []
        x = rand()
        for i1 = 1:250
            x = exp_map(x; r = r1)
            push!(xarr, x)
        end
        return xarr
    end

    fig = scatter(r_arr, r->track(r)[end], label="", markersize=0.8)

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

    fig

end

further reading

• runs this with more points offline
• https://www.google.com/search?q=logistic+map

• J.F. Boudreau and E.S. Swanson, Applied Computational Physics (ch.13)

• simple math models

• what about matrix?

let
    # more non-sense:
    
    # iterate a matrix, collect its trace
    
    function sin_map(x; r = 0.5)
        # iterate once
        return r * sin(x * pi)
    end

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

    function track(r1)
        xarr = []
        x = rand(2, 2)
        for i1 = 1:250
            x = sin_map(x; r = r1)
            push!(xarr, tr(x))
        end
        return xarr
    end

    fig = scatter(r_arr, r->track(r)[end], label="", markersize=0.8)

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

    fig

end