La théorie du chaos

' ********************************************************************
' Itérateur quadratique : x(n + 1) = a * x(n) * [1 - x(n)]
' ********************************************************************
' Ce programme montre le comportement d'un système dynamique :
' l'itérateur quadratique (ou suite logistique).
'
' Ce système est défini par l'équation :
'
'                      x(n+1) = a.x(n).[1 - x(n)]
'
' Pour a <= 3 la suite tend vers une valeur unique
'
' Pour a > 3 la suite devient périodique, avec une cascade
' de doublements de période :
'
' a        3.1    3.5      3.55      3.566
' Période  2      4        8         16
'
' Le régime chaotique débute à partir de a ~ 3.57, avec toutefois des
' zones de comportement périodiques à l'intérieur de ce régime, p. ex.
'
' a        3.63   3.74     3.83
' Période  6      5        3
'
' On peut aussi observer une alternance de phases périodiques et
' chaotiques, phénomène connu sous le nom d'intermittence (p. ex.
' pour a = 3.82812)
'
' Pour plus d'informations, voir la page de Wikipédia :
' https://fr.wikipedia.org/wiki/Suite_logistique
' **********************************************************************


' --------------------------------------------------------------------
' Constantes et variables globales
' --------------------------------------------------------------------

const Nmax = 100, Nstep = 10          ' Echelle sur Ox (itérations)
const Xmin = 0, Xmax = 1, Xstep = 0.1 ' Echelle sur Oy (valeurs x(n))

const Npixels = 5                      ' Nb de pixels par itération

const WCAN = Nmax * Npixels            ' Largeur du CANVAS
const HCAN = 400                       ' Hauteur du CANVAS

const Xfact = HCAN / (Xmax - Xmin)     ' Facteur d'échelle

const Font = FL_HELVETICA, Size = 15   ' Police de caract. et taille

dim X(Nmax)                            ' Suite x(n)

' --------------------------------------------------------------------
' Construction de l'interface graphique
' --------------------------------------------------------------------

dim_widget frm, alf, edit1, edit2, btn, can, canx, cany

widget frm,   T_Form,       100, 100, WCAN + 110, HCAN + 210, "Itérateur quadratique"
widget alf,   T_Alpha,       10,  10, WCAN + 90, 30, "x(n + 1) = a * x(n) * [1 - x(n)]"
widget edit1, T_Edit_Float,  80,  50, 100, 30, "a [1..4] = "
widget edit2, T_Edit_Float, 330,  50, 100, 30, "x(0) ]0..1[ = "
widget btn,   T_Button,     500,  50, 100, 30, "Tracer", CLICK_btn
widget can,   T_Canvas,      90, 120, WCAN, HCAN, "", DRAW_can
widget canx,  T_Canvas,      80, 130 + HCAN, WCAN + 30, 70, "", DRAW_canx, TIMER_canx
widget cany,  T_Canvas,      10, 110, 70, HCAN + 20, "", DRAW_cany, TIMER_cany

set_box_type alf, FL_FLAT_BOX
color_fl alf, FL_DARK_BLUE
color_caption_fl alf, FL_CYAN
font_caption alf, FL_HELVETICA_BOLD, 18

color_text_fl edit1, FL_RED
color_text_fl edit2, FL_RED

text edit1, "3.5"
text edit2, "0.5"

' --------------------------------------------------------------------
' Lancement de l'application
' --------------------------------------------------------------------

start frm

' --------------------------------------------------------------------
' Sous-programme auxiliaire
' --------------------------------------------------------------------

sub plot_point (i%, x)
' Marque le point pour l'itération i

  dim ip%, xp%

  ip = i * Npixels         ' Abscisse du point en pixels
  xp = Xfact * (x - Xmin)  ' Ordonnée du point en pixels

  if i = 0 then
    fl_move ip, xp
  else
    fl_draw ip, xp
  end_if

  fl_pie ip, xp, 3
end_sub

' --------------------------------------------------------------------
' Sous-programmes associes aux evenements
' --------------------------------------------------------------------

sub DRAW_canx ()
'Tracé de l'axe Ox

   dim x1%, xt%, yt%, wt%, ht%
   dim i%, n%, stp, txt$

   fl_pen_color FL_GRAY
   fl_rectangle_fill 0, FL_CANVAS_H, FL_CANVAS_W, FL_CANVAS_H

   fl_pen_color FL_BLACK
   fl_move 10, FL_CANVAS_H
   fl_draw FL_CANVAS_W - 20, FL_CANVAS_H

   fl_font font, size

   n = Nmax / Nstep
   stp = WCAN / n

   x1 = 10
   for i = 0 to n
     txt = str(i * Nstep)
     fl_move x1, FL_CANVAS_H
     fl_draw x1, FL_CANVAS_H - 10
     fl_text_extension txt, xt, yt, wt, ht
     fl_text_print txt, x1 - wt \ 2, 30 + ht
     x1 = x1 + stp
   next i

   txt = "Itérations"
   fl_text_extension txt, xt, yt, wt, ht
   fl_text_print txt, (WCAN - wt) \ 2, 20
end_sub

sub TIMER_canx
  widget_redraw canx
end_sub

sub DRAW_cany
' Tracé de l'axe Oy

   dim y1%, xt%, yt%, wt%, ht%
   dim i%, n%, stp, txt$

   fl_pen_color FL_GRAY
   fl_rectangle_fill 0, FL_CANVAS_H, FL_CANVAS_W, FL_CANVAS_H

   fl_pen_color FL_BLACK
   fl_move FL_CANVAS_W - 1, 10
   fl_draw FL_CANVAS_W - 1, FL_CANVAS_H - 10

   fl_font font, size

   n = (Xmax - Xmin) / Xstep
   stp = HCAN / n

   y1 = 10
   for i = 0 to n
     txt = dec(Xmin + i * Xstep, "0.0")
     fl_move FL_CANVAS_W - 1, y1
     fl_draw FL_CANVAS_W - 10, y1
     fl_text_extension txt, xt, yt, wt, ht
     fl_text_print txt, FL_CANVAS_W - 35, y1 - ht \ 2
     y1 = y1 + stp
   next i

   txt = "x(n)"
   fl_text_extension txt, xt, yt, wt, ht
   fl_text_print txt, 20, (HCAN - wt) \ 2, 90
end_sub

sub TIMER_cany
  widget_redraw cany
end_sub

sub DRAW_can()
' Tracé du diagramme

  dim i%

  fl_pen_color FL_BLACK
  fl_rectangle_fill 0, FL_CANVAS_H, FL_CANVAS_W, FL_CANVAS_H

  fl_pen_color FL_GREEN

  for i = 0 to Nmax
    plot_point i, X(i)
  next i
end_sub


sub CLICK_btn ()
' Calcul des itérations

  dim a, i%

  a = val(get_text(edit1))

  if a < 1 or a > 4 then
    message "a doit être entre 1 et 4"
    exit_sub
  end_if

  X(0) = val(get_text(edit2))

  if X(0) <= 0 or X(0) >= 1 then
    message "x(0) doit être entre 0 et 1 (non compris)"
    exit_sub
  end_if

  for i = 1 to Nmax
    X(i) = a * X(i - 1) * (1 - X(i - 1))
  next i

  widget_redraw can
end_sub

' ********************************************************************
' Itérateur quadratique : Calcul du spectre de puissance
' ********************************************************************

' --------------------------------------------------------------------
' Constantes et variables globales
' --------------------------------------------------------------------

const Tmin = 2, Tmax = 20, Tstep = 1       ' Echelle sur Ox (périodes)
const dBmin = 0, dBmax = 200, dBstep = 20  ' Echelle sur Oy (décibels)

const Tfact = 30                           ' Nb de pixels par période

const WCAN = (Tmax - Tmin) * Tfact         ' Largeur du CANVAS
const HCAN = 400                           ' Hauteur du CANVAS

const dBfact = HCAN / (dBmax - dBmin)      ' Facteur d'échelle sur Oy

const font = FL_HELVETICA, size = 15       ' Police de caract. et taille

const Ncoef = 50                           ' Nb de coefficients pour le calcul du spectre
dim P(WCAN)                                ' Puissances

' --------------------------------------------------------------------
' Construction de l'interface graphique
' --------------------------------------------------------------------

dim_widget frm, alf, edit1, edit2, spin1, spin2, btn, can, canx, cany

widget frm,   T_Form,       100, 100, WCAN + 100, HCAN + 250, "Itérateur quadratique : Spectre de puissance"
widget alf,   T_Alpha,       10,  10, WCAN + 80, 30, "x(n + 1) = a * x(n) * [1 - x(n)]"
widget edit1, T_Edit_Float,  90,  50, 100, 30, "a ]3..4[ = "
widget edit2, T_Edit_Float,  90,  90, 100, 30, "x(0) ]0..1[ = "
widget spin1, T_Spin,       305,  50,  75, 30, "Calculer "
widget spin2, T_Spin,       555,  50,  75, 30, " points à partir du point n° "
widget btn,   T_Button,     245,  90, 385, 30, "Calculer", CLICK_btn
widget can,   T_Canvas,      90, 160, WCAN, HCAN, "", DRAW_can
widget canx,  T_Canvas,      80, 170 + HCAN, WCAN + 30, 70, "", DRAW_canx, TIMER_canx
widget cany,  T_Canvas,      10, 150, 70, HCAN + 20, "", DRAW_cany, TIMER_cany

set_box_type alf, FL_FLAT_BOX
color_fl alf, FL_DARK_BLUE
color_caption_fl alf, FL_CYAN
font_caption alf, FL_HELVETICA_BOLD, 18

color_text_fl edit1, FL_RED
color_text_fl edit2, FL_RED
color_text_fl spin1, FL_RED
color_text_fl spin2, FL_RED

text edit1, "3.566"
text edit2, "0.5"

position_range spin1, 1000, 10000, 1000
position_range spin2, 100, 1000, 100

position spin1, 1000
position spin2, 100

' --------------------------------------------------------------------
' Lancement de l'application
' --------------------------------------------------------------------

start frm

' --------------------------------------------------------------------
' Sous-programme auxiliaire
' --------------------------------------------------------------------

function y_pixel% (P)
' Calcule l'ordonnée pour une valeur de puissance

  dim dB = 4.342944819032518 * log(P) ' 10 ln(P) / ln(10)
  return dBfact * (dB - dBmin)
end_function

' --------------------------------------------------------------------
' Sous-programmes associes aux evenements
' --------------------------------------------------------------------

sub DRAW_canx ()
' Tracé de l'axe Ox

  dim x1%, xt%, yt%, wt%, ht%
  dim i%, n%, stp, txt$

  fl_pen_color FL_GRAY
  fl_rectangle_fill 0, FL_CANVAS_H, FL_CANVAS_W, FL_CANVAS_H

  fl_pen_color FL_BLACK
  fl_move 10, FL_CANVAS_H
  fl_draw FL_CANVAS_W - 20, FL_CANVAS_H

  fl_font font, size

  n = (Tmax - Tmin) / Tstep
  stp = WCAN / n

  x1 = 10
  for i = 0 to n
   txt = str(Tmin + i * Tstep)
   fl_move x1, FL_CANVAS_H
   fl_draw x1, FL_CANVAS_H - 10
   fl_text_extension txt, xt, yt, wt, ht
   fl_text_print txt, x1 - wt \ 2, 30 + ht
   x1 = x1 + stp  
  next i

  txt = "Période"
  fl_text_extension txt, xt, yt, wt, ht
  fl_text_print txt, (WCAN - wt) \ 2, 15
end_sub

sub TIMER_canx
  widget_redraw canx
end_sub

sub DRAW_cany ()
' Tracé de l'axe Oy

  dim y1%, xt%, yt%, wt%, ht%
  dim i%, n%, stp, txt$

  fl_pen_color FL_GRAY
  fl_rectangle_fill 0, FL_CANVAS_H, FL_CANVAS_W, FL_CANVAS_H

  fl_pen_color FL_BLACK
  fl_move FL_CANVAS_W - 1, 10
  fl_draw FL_CANVAS_W - 1, FL_CANVAS_H - 10

  fl_font font, size

  n = (dBmax - dBmin) / dBstep
  stp = HCAN / n

  y1 = 10
  for i = 0 to n
    txt = dec(i * dBstep, "###")
    fl_move FL_CANVAS_W - 1, y1
    fl_draw FL_CANVAS_W - 10, y1
    fl_text_extension txt, xt, yt, wt, ht
    fl_text_print txt, FL_CANVAS_W - 35, y1 - ht \ 2
    y1 = y1 + stp
  next i

  txt = "Puissance relative (dB)"
  fl_text_extension txt, xt, yt, wt, ht
  fl_text_print txt, 20, (HCAN - wt) \ 2, 90
end_sub

sub TIMER_cany
  widget_redraw cany
end_sub

sub draw_can()
' Tracé du spectre

  dim i%

  fl_pen_color FL_BLACK
  fl_rectangle_fill 0, FL_CANVAS_H, FL_CANVAS_W, FL_CANVAS_H

  fl_pen_color FL_GREEN
  fl_move 1, y_pixel(P(1))
  
  for i = 2 to WCAN
    fl_draw i, y_pixel(P(i))
  next i
end_sub

sub click_btn ()
' Calcul du spectre de puissance

  dim a, x0, n1%, n2%, i%, T

  a = val(get_text(edit1))

  if a <= 3 or a >= 4 then
    message "a doit être entre 3 et 4 (non compris)"
    exit_sub
  end_if

  x0 = val (get_text(edit2))

  if x0 <= 0 or x0 >= 1 then
    message "x(0) doit être entre 0 et 1 (non compris)"
    exit_sub
  end_if

  n1 = get_position(spin1)
  n2 = get_position(spin2)

  for i = 1 to n2
    x0 = a * x0 * (1 - x0)
  next i

  dim x(n1), coef(Ncoef)

  x(0) = x0

  for i = 1 to n1
    x(i) = a * x(i - 1) * (1 - x(i - 1))
  next i

  spectrum_coef x(), coef()

  for i = 1 to WCAN
    T = Tmin + i / Tfact
    P(i) = spectrum_power(1 / T, 1, coef())
  next i

  P(0) = P(1)
  for i = 2 to WCAN
    if P(i) < P(0) then P(0) = P(i)
  next i

  for i = 1 to WCAN
    P(i) = P(i) / P(0)  ' Puissance relative
  next i

  widget_redraw can
end_sub

' ------------------------------------------------------------
' Spectre de puissance
' D'après "Numerical Recipes"
' http://www.phys.uri.edu/nigh/NumRec/bookfpdf/f13-7.pdf
' ------------------------------------------------------------

sub spectrum_coef (x(), coef())
' Calcul des coefficients du spectre
' Entrée : x(1..np) = tableau des valeurs du signal (points équidistants)
' Sortie : coef(0..nc) = tableau des coefficients

  dim np% = ubound(x)
  dim nc% = ubound(coef)

  dim wk1(np), wk2(np), wkm(nc)

  dim i%, j%, k%, num, denom, s

  s = 0
  for i = 1 to np
    s = s + x(i) * x(i)
  next i

  coef(0) = s / np

  wk1(1) = x(1)
  wk2(np - 1) = x(np)

  for i = 2 to np - 1
    wk1(i) = x(i)
    wk2(i - 1) = x(i)
  next i

  for k = 1 to nc
    num = 0
    denom = 0
    for i = 1 to np - k
      num = num + wk1(i) * wk2(i)
      denom = denom + wk1(i) * wk1(i) + wk2(i) * wk2(i)
    next i

    coef(k) = 2.0 * num / denom
    coef(0) = coef(0) * (1.0 - coef(k) * coef(k))

    if k > 1 then
      for j = 1 to k - 1
        coef(j) = wkm(j) - coef(k) * wkm(k - j)
      next j
    end_if

    if k = ncoef then exit_for

    for j = 1 to k
      wkm(j) = coef(j)
    next j

    for i = 1 to np - k - 1
      wk1(i) = wk1(i) - wkm(k) * wk2(i)
      wk2(i) = wk2(i + 1) - wkm(k) * wk1(i + 1)
    next i
  next k
end_sub

function spectrum_power (f, dt, coef())
' Retourne la valeur de la puissance pour une fréquence f
' La fréquence est comprise entre 0 et N/2 où N est la fréquence d'échantillonnage
' dt est l'unité de temps = 1/N

  const TwoPi = 6.28318530717959

  dim k%, theta, wpr, wpi, wr, wi, sumi, sumr, wtemp

  theta = TwoPi * f * dt
  wpr = cos(theta)
  wpi = sin(theta)
  wr = 1.0
  wi = 0.0
  sumr = 1.0
  sumi = 0.0

  for k = 1 to ubound(coef)
    wtemp = wr
    wr = wr * wpr - wi * wpi
    wi = wi * wpr + wtemp * wpi
    sumr = sumr - coef(k) * wr
    sumi = sumi - coef(k) * wi
  next k

  spectrum_power = coef(0) / (sumr * sumr + sumi * sumi)
end_function