池田のアトラクタ

program-id. chaos01.

data division.
working-storage section.
01 i pic 9(6).
01 a pic s9v999999.
01 b pic s9v999999.
01 t comp-2.
01 x0 comp-2.
01 x1 comp-2.
01 y0 comp-2.
01 y1 comp-2.
01 ad pic x(20).

procedure division.
display "IKeda attractor" upon syserr
display 'a=' upon syserr no advancing accept ad compute a = function numval(ad)
display 'b=' upon syserr no advancing accept ad compute b = function numval(ad)
compute x0 = 0.1 compute y0 = 0
perform varying i from 1 by 1 until i > 50000
compute t = x0 * x0 + y0 * y0 - 2.0
compute x1 = a + b * (x0 * function cos(t)- y0 * function sin(t))
compute y1 = b * (x0 * function sin(t) + y0 * function cos(t))
display x1 ' ' y1
compute x0 = x1 compute y0 = y1
end-perform
exit program.

gnuplot -p -e "plot 'tmp' ps 0 lt rgb 'red'"

gnuplot -p -e "set title 'a=4.0 b=0.4'; plot 'tmp' ps 0 lt rgb 'blue'"

gnuplot -p -e "set title 'a=4.0 b=0.8'; plot 'tmp' ps 0 lt rgb 'green'"

gnuplot -p -e "set title 'a=-4.0 b=0.6'; plot 'tmp' ps 0 lt rgb 'gold'"