混沌与分形

Published on Jun 27, 2020

2010年的某个夏天，第一次看到混沌的痕迹，人生也忽然发生了奇怪的变化。

宇宙万物的开始到结束，也许都是注定的。

随机只是一种幻觉。

真如如一，具有不变异性。

Lorenz System   ATTACH

微分方程

\begin{eqnarray} \frac{\Delta {X}}{\Delta {t}} &=& \sigma(Y-X) \label{eq:X}\\ \frac{\Delta {Y}}{\Delta {t}} &=& X (r - Z) -Y \label{eq:Y}\\ \frac{\Delta {Z}}{\Delta {t}} &=& XY - bZ \label{eq:Z} \end{eqnarray}

计算数据

## From wikipedia lorenz system

fileID = fopen("/tmp/fractal/lorenz_system.dat", "w");

function dx = lorenzatt(X)
    sigma = 10;beta = 8/3;rho = 28;
    dx = zeros(3,1);
    dx(1) = sigma*(X(2) - X(1));
    dx(2) = X(1)*(rho - X(3)) - X(2);
    dx(3) = X(1)*X(2) - beta*X(3);
    return
end

lsode_options("absolute tolerance", 1e-3)
lsode_options("relative tolerance",1e-4)
t = linspace(0,25,1e3); X0 = [1,1,1];
[X,T,MSG]=lsode(@lorenzatt,X0,t);

for i=1:length(X(:,1));
    fprintf(fileID, "%d %d %d\n", X(i,1), X(i,2), X(i,3))
end

绘图

set view 60, 30, 1, 1
unset pm3d
splot "/tmp/fractal/lorenz_system.dat" using 1:2:3 w l notitle
lorenz_attractor.png

Chua Circuit   ATTACH

计算数据

## https://www.chuacircuits.com/matlabsim.php

fileID = fopen("/tmp/fractal/chua_circuit.dat", "w");

function dx = chua(X)
    alpha = 15.6;beta = 28;m0 = -1.143;m1 = -0.714;
    dx = zeros(3,1);
    h = m1*X(1)+0.5*(m0-m1)*(abs(X(1)+1)-abs(X(1)-1));
    dx(1) = alpha*(X(2)-X(1)-h);
    dx(2) = X(1) - X(2) + X(3);
    dx(3) = - beta*X(2);
    return
end

lsode_options("absolute tolerance", 1e-3)
lsode_options("relative tolerance",1e-4)
t = linspace(0,100,1e3); X0 = [0.7,0,0];
[X,T,MSG]=lsode(@chua,X0,t);

for i=1:length(X(:,1));
    fprintf(fileID, "%d %d %d\n", X(i,1), X(i,2), X(i,3))
end

绘图

reset
set view 60, 30, 1, 1
unset pm3d
splot "/tmp/fractal/chua_circuit.dat" using 1:2:3 w l notitle
chua_circuit.pngt.png

Logistic Map   ATTACH

迭代公式

简单的 \(x(n)=rx(n-1)(1-x(n-1))\)

f = fopen("/tmp/fractal/logistic_map.dat", "w")
nn=1000;
r=linspace(2.0,3.9,1000);
for k=1:length(r);
    x(1)=0.4;
    for n=1:nn-1;
        x(n+1)=r(k)*x(n)*(1-x(n));
        if(n>0.9*nn)
            fprintf(f, "%d %d\n", r(k), x(n))
        end
    end
end

plot [2.8:] "/tmp/fractal/logistic_map.dat" using 1:2 w p pt 7 ps 0.1 notitle
logistic_map.png

Mandelbrot Set   ATTACH

rmax = 2
nmax = 100
complex (x, y) = x * {1, 0} + y * {0, 1}
mandelbrot (z, z0, n) = n == nmax || abs (z) > rmax ? n : mandelbrot (z ** 2 + z0, z0, n + 1)
set samples 200
set isosamples 200
set pm3d map
set size square
splot [-2 : .8] [-1.4 : 1.4] mandelbrot (complex (0, 0), complex (x, y), 0) notitle
mandelbrot.png

Julia Set   ATTACH

a = -0.79
b = 0.150
N = 100 # recusion limit
complex (x, y) = x * {1, 0} + y * {0, 1}
julia(x, y, z, n) = (abs(z)>2 || n<=0) ? 0 : julia(x, y, z*z+complex(x,y), n-1) + 1
set samples 200
set isosamples 200
set pm3d map
set size square
splot [-1.5 : 1.5] [-1.5 : 1.5] julia (a, b, complex (x, y), N) notitle
julia.png

Brownian Tree   ATTACH

DLA 过程

pkg load statistics
fileID = fopen("/tmp/browniantree.dat", "w")
function browniantree(fileID, xsize, ysize = xsize, numparticle = 1500)
  r = false(xsize, ysize);

  startx = unidrnd(xsize)
  starty = unidrnd(ysize)
  fprintf(fileID, "%d %d\n", startx, starty);
  r(startx, starty) = true;

  for i = 1:numparticle
    px = unidrnd(xsize-1)+1;
    py = unidrnd(ysize-1)+1;
    while(1)
      dx = unidrnd(3) - 2;
      dy = unidrnd(3) - 2;
      if ( (dx+px < 1) || (dx+px > xsize) || (dy+py < 1) || (dy+py > ysize) )
        px = unidrnd(xsize-1)+1;
        py = unidrnd(ysize-1)+1;
      elseif ( r(px+dx, py+dy) == true )
        px
        py
        fprintf(fileID, "%d %d\n", px, py);
        r(px, py) = 1;
        break;
      else
        px += dx;
        py += dy;
      endif
    endwhile
  endfor
endfunction
browniantree(fileID, 200, 200, 2000);

绘图

reset
set size square
plot "/tmp/browniantree.dat" pt 5 ps 1.0 notitle
browniantree.png

Sirepinski Triangle   ATTACH

一种随机的sierpinski 生成方法

每次都随🐔选择中点

fileID = fopen("/tmp/fractal/sierpinski.dat", "w")
x1 = [0,0];
x2 = [1,0];
x3 = [1/2, sqrt(3)/2];
start = rand(1,2)
for i=1:100000
  c = rand();
  if c >2/3
    start = (start + x3) / 2;
  elseif c > 1/3
    start = (start + x2) / 2;
  else
    start = (start + x1) / 2;
    fprintf(fileID, "%d %d\n", start(:,1),start(:,2))
  endif
endfor

unset xlabel
unset ylabel
plot "/tmp/fractal/sierpinski.dat" w p pt 7 ps 0.1 notitle
sierpinski.png

H-Tree   ATTACH

一个遍历平面所有点的曲线，具有奇异的特性

reset
set style arrow 1 nohead linewidth 1
set xrange [-1:1]
set yrange [0:sqrt(2)]
set notitle
unset xlabel
unset ylabel
unset xtick
unset ytick

number = 10
angle = pi / 2
s = sqrt(2)

rotx(x,y,t) = cos(t) * x - sin(t) * y
roty(x,y,t) = sin(t) * x + cos(t) * y

tree(n, x, y, dx, dy) = n >= number ? \
        sprintf("set arrow from %f,%f to %f,%f as 1;", x, y, x+dx/s, y+dy/s) : \
        sprintf("set arrow from %f,%f to %f,%f as 1;", x, y, x+dx/s, y+dy/s). \
        tree(n+1, x + dx/s, y + dy/s, rotx(dx/s, dy/s, angle), roty(dx/s, dy/s, angle )). \
        tree(n+1, x + dx/s, y + dy/s, rotx(dx/s, dy/s, -angle), roty(dx/s, dy/s, -angle ));

startx = 0.0
starty = 0.0
len = 1.0
eval(tree(0, startx, starty, 0, len))
plot -100 notitle
h-tree.png

稍微变换角度就成为了生命之树

reset
set style arrow 1 nohead linewidth 1
set xrange [-2:2]
set yrange [0:2.5]
set notitle
unset xlabel
unset ylabel
unset xtick
unset ytick

number = 10
angle = pi / 6
s = sqrt(2)

rotx(x,y,t) = cos(t) * x - sin(t) * y
roty(x,y,t) = sin(t) * x + cos(t) * y

tree(n, x, y, dx, dy) = n >= number ? \
        sprintf("set arrow from %f,%f to %f,%f as 1;", x, y, x+dx/s, y+dy/s) : \
        sprintf("set arrow from %f,%f to %f,%f as 1;", x, y, x+dx/s, y+dy/s). \
        tree(n+1, x + dx/s, y + dy/s, rotx(dx/s, dy/s, angle), roty(dx/s, dy/s, angle )). \
        tree(n+1, x + dx/s, y + dy/s, rotx(dx/s, dy/s, -angle), roty(dx/s, dy/s, -angle ));

startx = 0.0
starty = 0.0
len = 1.0
eval(tree(0, startx, starty, 0, len))
plot -100 notitle
tree.png

Built with in copyleft@2012-2022