Logistic Map
May 29, 2019
There’s nothing tricky about this:
; logistic map
(define r #f)
(define x #f)
(define (reset z)
(set! r z)
(set! x 0.5))
(define (logistic)
(let ((new-x (* r x (- 1 x))))
(set! x new-x)
x))
(define (drop n)
(do ((n n (- n 1))) ((zero? n))
(logistic)))
(define (take n)
(let loop ((n n) (zs (list)))
(if (zero? n) (reverse zs)
(loop (- n 1) (cons (logistic) zs)))))
The tricky part is looking at the results. With r = 0.75, the sequence tends to zero:
> (reset 0.75) > (drop 1000) > (take 10) (1.1801375113691321e-126 8.851031335268491e-127 6.638273501451368e-127 4.978705126088526e-127 3.734028844566395e-127 2.8005216334247964e-127 2.1003912250685973e-127 1.575293418801448e-127 1.181470064101086e-127 8.861025480758144e-128)
With r = 1.75, the sequence tends to a single non-zero value:
> (reset 1.75) > (drop 1000) > (take 10) (0.42857142857142855 0.42857142857142855 0.42857142857142855 0.42857142857142855 0.42857142857142855 0.42857142857142855 0.42857142857142855 0.42857142857142855 0.42857142857142855 0.42857142857142855)
Likewise with r = 2.75:
> (reset 2.75) > (drop 1000) > (take 10) (0.6363636363636365 0.6363636363636362 0.6363636363636365 0.6363636363636362 0.6363636363636365 0.6363636363636362 0.6363636363636365 0.6363636363636362 0.6363636363636365 0.6363636363636362)
But with r = 3.75, the results are chaotic:
(reset 3.75) > (drop 1000) > (0.8591179385506535 0.45387864829173385 0.9295230784372591 0.24566221908667554 0.6949210995003217 0.7950216186359463 0.6111084170153485 0.8912059487562879 0.36359214621634744 0.8677233653480164)
I’ll be reading some more about chaos. You can run the program at https://ideone.com/N5zQ93.
Here’s a Haskell program that prints the bifurcation diagram for the logistic map. I’ve included output for a couple of ranges of r values.
import Data.List (nub, sort) import System.Environment (getArgs, getProgName) import Text.Read (readMaybe) logistic :: Double -> Double -> [Double] logistic r = iterate (\x -> r*x*(1-x)) -- A string with dots at the 1-based positions in ns. dotsAt :: [Int] -> String dotsAt ns = concat $ zipWith dotAt (0:ns) ns where dotAt i n = replicate (n-i-1) ' ' ++ "." -- One row of dots corresponding to the x values that were visited for one value -- of r. oneRow :: Int -> [Double] -> String oneRow cols = let w = 1/(fromIntegral cols - 1) in dotsAt . nub . sort . map (\x -> floor (x/w) + 1) -- One row of the bifurcation diagram, corresponding to one r value. Drop some -- initial values of the logistic map to let it converge to a value (if it in -- fact will converge). showRow :: Int -> Double -> Double -> String showRow cols r = oneRow cols . take 100 . drop 1000 . logistic r -- Plot the bifurcation diagram of the logistic map, for rmin ≤ r ≤ rmax, on a -- graph having the given number of columns and rows. plot :: Int -> Int -> Double -> Double -> Double -> IO () plot cols rows rmin rmax x = mapM_ (\r -> putStrLn $ showRow cols r x) $ rs rows rmin rmax -- N evenly spaced values of r, where rmin ≤ r ≤ rmax. rs :: Int -> Double -> Double -> [Double] rs n rmin rmax = let rb = (rmax - rmin)/(fromIntegral n - 1) in [i*rb + rmin | j <- [0..n-1], let i = fromIntegral j] data Args = Args Int Int Double Double deriving (Read) parseArgs :: [String] -> Maybe Args parseArgs args = readMaybe $ "Args " ++ unwords args main :: IO () main = do prog <- getProgName args <- getArgs let x = 0.4 -- starting x value case parseArgs args of Just (Args width height rmin rmax) -> plot width height rmin rmax x _ -> error $ "Usage: " ++ prog ++ " width height rmin rmax"$ ./logistic 120 50 2.8 3.8 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . .. ...... ... .. ....... ...... ... ................ . .. .... .. .. ............. ....... . . ..... .. .. . .... .... ........ ....... ....... ... ... ... . ... . . . ..... ............. ..... ........... . ... . .. . ..... ... . . .. .... ..... ..... .... . .... . .. ... .. .. .. ... .. ..... ... .. ............. .......... ........ . . .. . ... .. .... .. .. . .... . ...... ... . ......... ..... ..... . .. .. . .. . ... . ... . ................... ... . .... . . . . . ...... .. .. . .. ... .. ... . . . ........ . .. . ..... ...... ..... ... . .. .. ... .. . .. .. .. . . ... .. .. ... . .. ..... .. ... ... ... . .... .. . . . . . .. . . ..... .. . .. .......... .... .. .......... .Zooming in on the region 3.4 <= r <= 3.7.
$ ./logistic 120 50 3.4 3.7 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . .. . . . . . . .. . . . . . . .. . . . . .. .. . ... .. ... ...... .... .. ... ... .... . . . . . . ... . . ... . .. .. . ... . .. .. . .... ......... . .... .......... ..... ....... .......... .. ..... .. . .... .... ........ .... ..... . ..... .. ....... . ............... ...... ............. . ... . ..... ............... ........ .......... ... . ..... . .. ................ .... .......... ... . . . . .... .. . .... ......... . . . . . . . .. . .. . .. .. ... . . . .... .... . ..... ... . ...... ... .... . ... .... .... .. . . .. ..... . . ........ .......... .... ... . . . .. ... . ... . ..... ..... .... . .. .... ....... . . .. . .. .. ...... . ... . . . .......... ......... ........ ...... . ..... . . . .. ... .... . ... ..... ............... .. .. .. . . ..... . .. . . ..... ... .... .. ..... ...... . ........ .... ... .... . .... . . . . . . .. . . ...... .. . ... ................. .... .. . .. ... .... . . . ... ... .. . .... .................. ... ...... .. .. . . . . ... .. . . . ... ..... ...... . ..... . .. .. . ... .. .. .... . ... ... . . ....... .. ......... ..... . ........... ... . . . . . .. . .. ...... ... .. .... ............... ....... ......Here’s a solution in Python that shows a bifurcation diagram, inspired by @Globules’s solution above.
import math def logistic_map(x, r): while True: yield x x = r * x * (1 - x) CHAR = 'x' INIT = 0.5 R_ITERS = 50 T_SKIP = 1000 T_KEEP = 1000 WIDTH = 80 R_MIN = 2.4 # inclusive R_MAX = 4.0 # exclusive r_range = R_MAX - R_MIN for iter in range(R_ITERS): r = R_MIN + iter * (R_MAX - R_MIN) / R_ITERS gen = logistic_map(INIT, r) chars = [' '] * WIDTH for _ in range(T_SKIP): next(gen) for _ in range(T_KEEP): x = next(gen) idx = math.floor(x * WIDTH) chars[idx] = CHAR print(''.join(chars))Output:
x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x xx x x x x x x x x xx x x x x xx xxxxxx xxxxxxxxxxxx xxxxxxxxx xxxxxxxxxxxxxxxxxxxxxxxxx xxxxxxxxxxx xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx xxxxxxxxxxxxx xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx xx xxxx xx xxx x xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx x x x xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx xxxxxxxxxxxxxxxxxx xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx