FORUM DE DISCUSSION SUR LE LANGAGE PANORAMIC

Développement d'applications avec le langage Panoramic
 
AccueilAccueil  FAQFAQ  RechercherRechercher  S'enregistrerS'enregistrer  MembresMembres  GroupesGroupes  Connexion  
Derniers sujets
» KGF_dll - nouvelles versions
par Klaus Aujourd'hui à 23:33

» KGF.dll - demandes ou suggestions de modifications ou ajouts
par Pedro Alvarez Aujourd'hui à 22:45

» StratégoV8
par Jean Claude Aujourd'hui à 19:55

» probleme de touche (verr num)
par pascal10000 Aujourd'hui à 16:20

» Mah-Jong européen new-look
par jjn4 Aujourd'hui à 15:58

» track_bar circulaire
par Klaus Hier à 21:41

» Gestionnaire de Projets Panoramic
par Jean Claude Hier à 20:41

» Gestion de l'Unicode
par jean_debord Hier à 11:21

» Button_picture
par pascal10000 Jeu 14 Déc 2017 - 11:41

» Pourquoi le compilateur stagne
par Minibug Jeu 14 Déc 2017 - 11:09

» 4 (en analyse): SYNEDIT_TARGET_IS_OBJECT devient inactif
par Jack Jeu 14 Déc 2017 - 10:09

» 3 (en analyse): Mauvaise interprétation du string "THEN"
par Jack Jeu 14 Déc 2017 - 10:03

» API Windows
par Klaus Mar 12 Déc 2017 - 3:21

» Cartes de voeux, menus, etc.
par JL35 Lun 11 Déc 2017 - 17:48

» a l'aide klaus
par Minibug Lun 11 Déc 2017 - 11:42

Navigation
 Portail
 Index
 Membres
 Profil
 FAQ
 Rechercher
Rechercher
 
 

Résultats par :
 
Rechercher Recherche avancée
Décembre 2017
LunMarMerJeuVenSamDim
    123
45678910
11121314151617
18192021222324
25262728293031
CalendrierCalendrier

Partagez | 
 

 Fractale de Newton

Voir le sujet précédent Voir le sujet suivant Aller en bas 
AuteurMessage
jean_debord

avatar

Nombre de messages : 766
Age : 63
Localisation : Limoges
Date d'inscription : 21/09/2008

MessageSujet: Fractale de Newton   Sam 24 Mai 2014 - 11:19

La méthode de Newton résout une équation par approximations successives. On colore les points selon le nombre d'itérations nécessaires pour atteindre l'une des racines de l'équation. Pour chaque racine, le "bassin d'atraction" est l'ensemble des points qui convergent vers cette racine. La frontière entre les différents bassins est une courbe fractale.

L'exemple montre l'équation classique z^3-1=0 qui a 3 racines dans l'ensemble des nombres complexes.

Le programme est prévu pour le compilateur. Il passera aussi avec l'interpréteur, mais beaucoup plus lentement.

Tout cela sera développé dans un prochain article ...

Code:

' **********************************************************************
' Fractale de Newton
' **********************************************************************

' Variables definies par l'utilisateur

dim PicWidth%  : PicWidth%  = 640   : ' Taille de l'image en pixels
dim PicHeight% : PicHeight% = 480

dim x0         : x0 = 0             : ' Coord. du centre de l'image
dim y0         : y0 = 0
dim Eps        : Eps = 0.000001     : ' Precision requise
dim MaxIter%   : MaxIter% = 100     : ' Nb maxi d'iterations
dim ZoomFact   : ZoomFact = 1       : ' Facteur de zoom
dim ColorFact  : ColorFact = -2     : ' Facteur de coloration
dim CstVal     : CstVal = 0.9       : ' Luminosite HSV (constante)

' Variables supplementaires

dim HalfPicWidth%  : HalfPicWidth%  = PicWidth% / 2
dim HalfPicHeight% : HalfPicHeight% = PicHeight% / 2

dim R%, G%, B%  : ' Parametres de coloration
dim AbsColor    : ' abs(ColorFact)
dim ScaleFact   : ' Facteur d'echelle = distance entre 2 pixels
dim Nx%, Ny%    : ' Coordonnees d'un point (pixels)
dim xt, yt      : ' Coordonnees d'un point (algebriques)
dim f_x, f_y    : ' Fonction complexe f(z)
dim fp_x, fp_y  : ' Derivee f'(z)

' ----------------------------------------------------------------------
' Description des objets
' ----------------------------------------------------------------------

' Fenetre principale

left    0, 50
top     0, 50
width   0, PicWidth% + 70
height  0, PicHeight% + 120
caption 0, "Fractale de Newton ... Veuillez patienter."

' Image

picture 1
left    1, 30
top     1, 40
width   1, PicWidth%
height  1, PicHeight%

2d_target_is 1

' ----------------------------------------------------------------------
' Programme principal
' ----------------------------------------------------------------------

ColorFact = 0.01 * ColorFact
AbsColor  = abs(ColorFact)
ScaleFact = 4 / (PicHeight% * ZoomFact)

for Ny% = 0 to PicHeight% - 1
 yt = y0 - ScaleFact * (Ny% - HalfPicHeight%)
 for Nx% = 0 to PicWidth% - 1
   xt = x0 + ScaleFact * (Nx% - HalfPicWidth%)
   Newton(xt, yt)
   2d_pen_color R%, G%, B%
   2d_point Nx%, Ny%
 next Nx%
next Ny%

file_save 1, "newton.bmp"

caption 0, "Fractale de Newton ... Terminé."

end

' ----------------------------------------------------------------------
' Sous-programmes
' ----------------------------------------------------------------------

sub Func(x, y)
' Calcul de la fonction complexe f(z) = (f_x, f_y)
' et de sa derivee f'(z) = (fp_x, fp_y)
' Ici f(z) = x^3 - 1 et f'(z) = 3 x^2

 dim_local x2, y2, x3, y3
 
 if x = 0 and y = 0
   
   f_x = 10
   f_y = 10
   
 else  
   
   x2 = x * x : x3 = x2 * x
   y2 = y * y : y3 = y2 * y
 
   f_x = x3 - 3 * x * y2 - 1
   f_y = 3 * x2 * y - y3
 
   fp_x = 3 * (x2 - y2)
   fp_y = 6 * x * y
 
 end_if    
end_sub

sub Newton(xt, yt)
' Iteration de la fonction complexe au point (xt, yt)

 dim_local Iter%, mf2, mq2, x, y, zn_x, zn_y, d, q_x, q_y

 zn_x = xt
 zn_y = yt
 
 Iter% = 0
 mf2 = 1
 mq2 = 1

 while Iter% < MaxIter% and mf2 > Eps and mq2 > Eps

   x = zn_x
   y = zn_y

   Func(x, y)

   d = fp_x * fp_x + fp_y * fp_y
 
   q_x = (f_x * fp_x + f_y * fp_y) / d  : ' q = f(z) / f'(z)  
   q_y = (f_y * fp_x - f_x * fp_y) / d
 
   zn_x = x - q_x
   zn_y = y - q_y

   mf2 = f_x * f_x + f_y * f_y
   mq2 = q_x * q_x + q_y * q_y
   
   Iter% = Iter% + 1

 end_while
 
 if Iter% > MaxIter%  : ' Pas de convergence
   R% = 0
   G% = 0
   B% = 0
 else
   MdbCol(Iter%)
 end_if
end_sub

sub MdbCol(Iter%)
' Determine la teinte (Hue, H) et la Saturation (S)
' d'apres le nombre d'iterations

 dim_local Q, Angle, Radius, H, S, V
 
 Q = log(Iter%) * AbsColor

 if Q < 0.5
   Q = 1 - 1.5 * Q
   Angle = 1 - Q
 else
   Q = 1.5 * Q - 0.5
   Angle = Q
 end_if

 Radius = sqr(Q)

' Si ColorFact > 0, assombrir une bande sur 2

 V = CstVal

 if (ColorFact > 0) and (odd(Iter%) > 0)
   V = 0.85 * V
   Radius = 0.667 * Radius
 end_if

 H = Angle * 10
 H = H - int(H)
 H = H * 360

 S = Radius - int(Radius)

' Convertir HSV en RGB

 HSVtoRGB(H, S, V)
end_sub

sub HSVtoRGB(H, S, V)
' Conversion HSV --> RGB. Resultats dans R%, G%, B%

 dim_local II%, ZZ, FF, PP, QQ, TT, RR, GG, BB
 
 if S = 0

   R% = int(V * 255)
   G% = R%
   B% = R%

 else

   ZZ  = H / 60
   II% = int(ZZ)
   FF  = ZZ - int(ZZ)
   PP  = V * (1 - S)
   QQ  = V * (1 - S * FF)
   TT  = V * (1 - S * (1 - FF))

   select II%
     case 0
       RR = V  : GG = TT : BB = PP
     case 1
       RR = QQ : GG = V  : BB = PP
     case 2
       RR = PP : GG = V  : BB = TT
     case 3
       RR = PP : GG = QQ : BB = V
     case 4
       RR = TT : GG = PP : BB = V
     case 5
       RR = V  : GG = PP : BB = QQ
   end_select

   R% = int(RR * 255)
   G% = int(GG * 255)
   B% = int(BB * 255)

 end_if
end_sub
Revenir en haut Aller en bas
Voir le profil de l'utilisateur http://www.unilim.fr/pages_perso/jean.debord/index.htm
papydall

avatar

Nombre de messages : 5613
Age : 67
Localisation : Moknine (Tunisie) Entre la chaise et le clavier
Date d'inscription : 03/03/2012

MessageSujet: Re: Fractale de Newton   Sam 24 Mai 2014 - 14:40

Une fois de plus, Bravo Jean_Debord pour cette belle fractale!
Revenir en haut Aller en bas
Voir le profil de l'utilisateur http://papydall-panoramic.forumarabia.com/
Jicehel

avatar

Nombre de messages : 5865
Age : 45
Localisation : 77500
Date d'inscription : 19/04/2011

MessageSujet: Re: Fractale de Newton   Sam 24 Mai 2014 - 18:38

Superbe et belle démo pour l'intérêt de compiler ses programmes.
Revenir en haut Aller en bas
Voir le profil de l'utilisateur
jean_debord

avatar

Nombre de messages : 766
Age : 63
Localisation : Limoges
Date d'inscription : 21/09/2008

MessageSujet: Re: Fractale de Newton   Mar 24 Juin 2014 - 10:08

Version améliorée. Fonctionne avec un exposant entier quelconque et utilise des SUBs pour les calculs sur les complexes.

En faisant varier l'exposant vous pouvez obtenir différentes images.

Le compilateur n'étant pas très à l'aise avec les longues séries de SUBs, je n'ai intégré que ceux dont j'avais besoin.

Code:

' **************************************************************
' Fractale de Newton : Fonction z^p - 1 (p entier)
' **************************************************************

' Variables definies par l'utilisateur

dim PicWidth%  : PicWidth%  = 640   : ' Taille de l'image
dim PicHeight% : PicHeight% = 480   : '   en pixels
dim p%         : p% = 5             : ' Exposant (z^p - 1)
dim x0         : x0 = 0             : ' Coordonnees du centre
dim y0         : y0 = 0             : '   de l'image
dim Eps        : Eps = 0.000001     : ' Precision requise
dim MaxIter%   : MaxIter% = 100     : ' Nb maxi d'iterations
dim ZoomFact   : ZoomFact = 1       : ' Facteur de zoom
dim ColorFact  : ColorFact = -2     : ' Facteur de coloration
dim CstVal     : CstVal = 0.9       : ' Luminosite HSV constante

' Variables supplementaires

dim HalfPicWidth%  : HalfPicWidth%  = PicWidth% / 2
dim HalfPicHeight% : HalfPicHeight% = PicHeight% / 2
dim Eps2           : Eps2 = Eps * Eps

dim R%, G%, B%  : ' Parametres de coloration
dim AbsColor    : ' abs(ColorFact)
dim ScaleFact   : ' Facteur d'echelle = distance entre 2 pixels
dim Nx%, Ny%    : ' Coordonnees d'un point (pixels)
dim xt, yt      : ' Coordonnees d'un point (algebriques)
dim f_x, f_y    : ' Fonction complexe f(z)
dim fp_x, fp_y  : ' Derivee f'(z)

' --------------------------------------------------------------
' Variables globales de la bibliotheque (nombres complexes)
' --------------------------------------------------------------

' Constantes mathematiques

dim MaxNum, MinNum, MaxLog, MinLog, Pi, PiDiv2

MaxLog =  709.78     : ' Argument max. pour EXP
MinLog = -708.39     : ' Argument min. pour EXP
MaxNum = exp(MaxLog) : ' Nb reel max. ~ 2^1024
MinNum = exp(MinLog) : ' Nb reel min. ~ 2^(-1022)
Pi     = 4 * atn(1)
PiDiv2 = Pi / 2

' Resultats des calculs
' Partie reelle, partie imaginaire, module, argument, signe

dim r_x, r_y, r_mod, r_arg, r_sgn

' Code d'erreur
'  0 = pas d'erreur
' -1 = argument hors bornes
' -2 = singularite
' -3 = overflow
' -4 = underflow

dim ErrCode%

' --------------------------------------------------------------
' Description des objets
' --------------------------------------------------------------

' Fenetre principale

left    0, 50
top     0, 50
width   0, PicWidth% + 70
height  0, PicHeight% + 120
caption 0, "Fractale de Newton ... Veuillez patienter."

' Image

picture 1
left    1, 30
top     1, 40
width   1, PicWidth%
height  1, PicHeight%

2d_target_is 1

' --------------------------------------------------------------
' Programme principal
' --------------------------------------------------------------

ColorFact = 0.01 * ColorFact
AbsColor  = abs(ColorFact)
ScaleFact = 4 / (PicHeight% * ZoomFact)

for Ny% = 0 to PicHeight% - 1
 yt = y0 - ScaleFact * (Ny% - HalfPicHeight%)
 for Nx% = 0 to PicWidth% - 1
   xt = x0 + ScaleFact * (Nx% - HalfPicWidth%)
   Newton(xt, yt)
   2d_pen_color R%, G%, B%
   2d_point Nx%, Ny%
 next Nx%
next Ny%

file_save 1, "newton.bmp"

caption 0, "Fractale de Newton ... Terminé."

end

' --------------------------------------------------------------
' Sous-programmes
' --------------------------------------------------------------

sub Func(x, y)
' Calcul de la fonction complexe f(z) = (f_x, f_y)
' et de sa derivee f'(z) = (fp_x, fp_y)
' Ici f(z) = z^p - 1 et f'(z) = p z^(p - 1)

 dim_local u_x, u_y     : ' u = z^(p - 1)

 CIntPower(x, y, p% - 1)
 
 u_x = r_x
 u_y = r_y
 
 CMul(u_x, u_y, x, y)

 f_x = r_x - 1
 f_y = r_y

 fp_x = p% * u_x
 fp_y = p% * u_y
end_sub

sub Newton(xt, yt)
' Iteration de la fonction complexe au point (xt, yt)

 dim_local Iter%, mf2, mq2, x, y, zn_x, zn_y, q_x, q_y

 zn_x = xt
 zn_y = yt

 Iter% = 0

 repeat
   x = zn_x
   y = zn_y

   Func(x, y)

   if fp_x = 0 and fp_y = 0
     Iter% = MaxIter%             : ' Arret du calcul si f'(z) = 0
   else
     CDiv(f_x, f_y, fp_x, fp_y)   : ' q = f(z) / f'(z)
     
     q_x = r_x
     q_y = r_y

     zn_x = x - q_x
     zn_y = y - q_y

     mf2 = f_x * f_x + f_y * f_y  : ' |f(z)|^2
     mq2 = q_x * q_x + q_y * q_y  : ' |q|^2

     Iter% = Iter% + 1
   end_if
 until (mf2 < Eps2 and mq2 < Eps2) or (Iter% > MaxIter%)

 if Iter% > MaxIter%  : ' Pas de convergence
   R% = 0
   G% = 0
   B% = 0
 else
   MdbCol(Iter%)
 end_if
end_sub

sub MdbCol(Iter%)
' Determine la teinte (Hue, H) et la Saturation (S)
' d'apres le nombre d'iterations

 dim_local Q, Angle, Radius, H, S, V

 Q = log(Iter%) * AbsColor

 if Q < 0.5
   Q = 1 - 1.5 * Q
   Angle = 1 - Q
 else
   Q = 1.5 * Q - 0.5
   Angle = Q
 end_if

 Radius = sqr(Q)

' Si ColorFact > 0, assombrir une bande sur 2

 V = CstVal

 if (ColorFact > 0) and (odd(Iter%) > 0)
   V = 0.85 * V
   Radius = 0.667 * Radius
 end_if

 H = Angle * 10
 H = H - int(H)
 H = H * 360

 S = Radius - int(Radius)

' Convertir HSV en RGB

 HSVtoRGB(H, S, V)
end_sub

sub HSVtoRGB(H, S, V)
' Conversion HSV --> RGB. Resultats dans R%, G%, B%

 dim_local II%, ZZ, FF, PP, QQ, TT, RR, GG, BB

 if S = 0
   R% = int(V * 255)
   G% = R%
   B% = R%
 else
   ZZ  = H / 60
   II% = int(ZZ)
   FF  = ZZ - int(ZZ)
   PP  = V * (1 - S)
   QQ  = V * (1 - S * FF)
   TT  = V * (1 - S * (1 - FF))

   select II%
     case 0
       RR = V  : GG = TT : BB = PP
     case 1
       RR = QQ : GG = V  : BB = PP
     case 2
       RR = PP : GG = V  : BB = TT
     case 3
       RR = PP : GG = QQ : BB = V
     case 4
       RR = TT : GG = PP : BB = V
     case 5
       RR = V  : GG = PP : BB = QQ
   end_select

   R% = int(RR * 255)
   G% = int(GG * 255)
   B% = int(BB * 255)
 end_if
end_sub

' --------------------------------------------------------------
' Procedures de la bibliotheque (nombres complexes)
' --------------------------------------------------------------

sub CMul(a_x, a_y, b_x, b_y)
' Multiplication : r_x + i r_y = (a_x + i a_y) * (b_x + i b_y)

 ErrCode% = 0

 r_x = a_x * b_x - a_y * b_y
 r_y = a_x * b_y + a_y * b_x
end_sub

sub CDiv(a_x, a_y, b_x, b_y)
' Division : r_x + i r_y = (a_x + i a_y) / (b_x + i b_y)
' Algorithme d'apres "Numerical Recipes"

 dim_local q, t

 if b_x = 0 and b_y = 0
   ErrCode% = -3
   r_x = MaxNum
   r_y = MaxNum
 else
   ErrCode% = 0
   if abs(b_x) >= abs(b_y)
     q = b_y / b_x
     t = b_x + b_y * q
     r_x = (a_x + a_y * q) / t
     r_y = (a_y - a_x * q) / t
   else
     q = b_x / b_y
     t = b_x * q + b_y
     r_x = (a_x * q + a_y) / t
     r_y = (a_y * q - a_x) / t
   end_if
 end_if
end_sub

sub CSquare(a_x, a_y)
' Carre : r_x + i r_y = (a_x + i a_y)^2

 ErrCode% = 0

 r_x = a_x * a_x - a_y * a_y
 r_y = 2 * a_x * a_y
end_sub

sub CIntPower(a_x, a_y, n%)
' Puissance entiere : r_x + i r_y = (a_x + i a_y)^n

 dim_local m%, b_x, b_y, res_x, res_y
 
 ErrCode% = 0
 
 if a_x = 0 and a_y = 0
   if n% = 0
     ' 0^0 = lim x^x quand x --> 0 = 1
     r_x = 1
     r_y = 0
   else
     if n% > 0
       ' 0^n = 0 si n > 0
       r_x = 0
       r_y = 0
     else
       ' 0^n indefini si n < 0
       ErrCode% = -2
       r_x = MaxNum
       r_y = MaxNum
     end_if
   end_if
 else
   if n% < 0
     m% = abs(n%)
     CInv(a_x, a_y)
     b_x = r_x
     b_y = r_y
   else
     m% = n%
     b_x = a_x
     b_y = a_y
   end_if

   res_x = 1 : res_y = 0

   while m% > 0
     if odd(m%) = 1
       CMul(b_x, b_y, res_x, res_y)
       res_x = r_x
       res_y = r_y
     end_if
     CSquare(b_x, b_y)
     b_x = r_x
     b_y = r_y
     m% = int(m% / 2)
   end_while
 
   r_x = res_x
   r_y = res_y
 end_if  
end_sub
Revenir en haut Aller en bas
Voir le profil de l'utilisateur http://www.unilim.fr/pages_perso/jean.debord/index.htm
jean_debord

avatar

Nombre de messages : 766
Age : 63
Localisation : Limoges
Date d'inscription : 21/09/2008

MessageSujet: Re: Fractale de Newton   Mar 24 Juin 2014 - 10:11

Version pour exposant réel :

Code:

' **************************************************************
' Fractale de Newton : Fonction z^p - 1 (p reel)
' **************************************************************

' Variables definies par l'utilisateur

dim PicWidth%  : PicWidth%  = 640   : ' Taille de l'image
dim PicHeight% : PicHeight% = 480   : '   en pixels
dim p          : p = 3.5            : ' Exposant (z^p - 1)
dim x0         : x0 = 0             : ' Coordonnees du centre
dim y0         : y0 = 0             : '   de l'image
dim Eps        : Eps = 0.000001     : ' Precision requise
dim MaxIter%   : MaxIter% = 100     : ' Nb maxi d'iterations
dim ZoomFact   : ZoomFact = 1       : ' Facteur de zoom
dim ColorFact  : ColorFact = -2     : ' Facteur de coloration
dim CstVal     : CstVal = 0.9       : ' Luminosite HSV constante

' Variables supplementaires

dim HalfPicWidth%  : HalfPicWidth%  = PicWidth% / 2
dim HalfPicHeight% : HalfPicHeight% = PicHeight% / 2
dim Eps2           : Eps2 = Eps * Eps

dim R%, G%, B%  : ' Parametres de coloration
dim AbsColor    : ' abs(ColorFact)
dim ScaleFact   : ' Facteur d'echelle = distance entre 2 pixels
dim Nx%, Ny%    : ' Coordonnees d'un point (pixels)
dim xt, yt      : ' Coordonnees d'un point (algebriques)
dim f_x, f_y    : ' Fonction complexe f(z)
dim fp_x, fp_y  : ' Derivee f'(z)

' --------------------------------------------------------------
' Variables globales de la bibliotheque (nombres complexes)
' --------------------------------------------------------------

' Constantes mathematiques

dim MaxNum, MinNum, MaxLog, MinLog, Pi, PiDiv2

MaxLog =  709.78     : ' Argument max. pour EXP
MinLog = -708.39     : ' Argument min. pour EXP
MaxNum = exp(MaxLog) : ' Nb reel max. ~ 2^1024
MinNum = exp(MinLog) : ' Nb reel min. ~ 2^(-1022)
Pi     = 4 * atn(1)
PiDiv2 = Pi / 2

' Resultats des calculs
' Partie reelle, partie imaginaire, module, argument, signe

dim r_x, r_y, r_mod, r_arg, r_sgn

' Code d'erreur
'  0 = pas d'erreur
' -1 = argument hors bornes
' -2 = singularite
' -3 = overflow
' -4 = underflow

dim ErrCode%

' --------------------------------------------------------------
' Description des objets
' --------------------------------------------------------------

' Fenetre principale

left    0, 50
top     0, 50
width   0, PicWidth% + 70
height  0, PicHeight% + 120
caption 0, "Fractale de Newton ... Veuillez patienter."

' Image

picture 1
left    1, 30
top     1, 40
width   1, PicWidth%
height  1, PicHeight%

2d_target_is 1

' --------------------------------------------------------------
' Programme principal
' --------------------------------------------------------------

ColorFact = 0.01 * ColorFact
AbsColor  = abs(ColorFact)
ScaleFact = 4 / (PicHeight% * ZoomFact)

for Ny% = 0 to PicHeight% - 1
 yt = y0 - ScaleFact * (Ny% - HalfPicHeight%)
 for Nx% = 0 to PicWidth% - 1
   xt = x0 + ScaleFact * (Nx% - HalfPicWidth%)
   Newton(xt, yt)
   2d_pen_color R%, G%, B%
   2d_point Nx%, Ny%
 next Nx%
next Ny%

file_save 1, "newton.bmp"

caption 0, "Fractale de Newton ... Terminé."

end

' --------------------------------------------------------------
' Sous-programmes
' --------------------------------------------------------------

sub Func(x, y)
' Calcul de la fonction complexe f(z) = (f_x, f_y)
' et de sa derivee f'(z) = (fp_x, fp_y)
' Ici f(z) = z^p - 1 et f'(z) = p z^(p - 1)

 dim_local u_x, u_y     : ' u = z^(p - 1)

 CRealPower(x, y, p - 1)
 
 u_x = r_x
 u_y = r_y
 
 CMul(u_x, u_y, x, y)

 f_x = r_x - 1
 f_y = r_y

 fp_x = p * u_x
 fp_y = p * u_y
end_sub

sub Newton(xt, yt)
' Iteration de la fonction complexe au point (xt, yt)

 dim_local Iter%, mf2, mq2, x, y, zn_x, zn_y, q_x, q_y

 zn_x = xt
 zn_y = yt

 Iter% = 0

 repeat
   x = zn_x
   y = zn_y

   Func(x, y)

   if fp_x = 0 and fp_y = 0
     Iter% = MaxIter%             : ' Arret du calcul si f'(z) = 0
   else
     CDiv(f_x, f_y, fp_x, fp_y)   : ' q = f(z) / f'(z)
     
     q_x = r_x
     q_y = r_y

     zn_x = x - q_x
     zn_y = y - q_y

     mf2 = f_x * f_x + f_y * f_y  : ' |f(z)|^2
     mq2 = q_x * q_x + q_y * q_y  : ' |q|^2

     Iter% = Iter% + 1
   end_if
 until (mf2 < Eps2 and mq2 < Eps2) or (Iter% > MaxIter%)

 if Iter% > MaxIter%  : ' Pas de convergence
   R% = 0
   G% = 0
   B% = 0
 else
   MdbCol(Iter%)
 end_if
end_sub

sub MdbCol(Iter%)
' Determine la teinte (Hue, H) et la Saturation (S)
' d'apres le nombre d'iterations

 dim_local Q, Angle, Radius, H, S, V

 Q = log(Iter%) * AbsColor

 if Q < 0.5
   Q = 1 - 1.5 * Q
   Angle = 1 - Q
 else
   Q = 1.5 * Q - 0.5
   Angle = Q
 end_if

 Radius = sqr(Q)

' Si ColorFact > 0, assombrir une bande sur 2

 V = CstVal

 if (ColorFact > 0) and (odd(Iter%) > 0)
   V = 0.85 * V
   Radius = 0.667 * Radius
 end_if

 H = Angle * 10
 H = H - int(H)
 H = H * 360

 S = Radius - int(Radius)

' Convertir HSV en RGB

 HSVtoRGB(H, S, V)
end_sub

sub HSVtoRGB(H, S, V)
' Conversion HSV --> RGB. Resultats dans R%, G%, B%

 dim_local II%, ZZ, FF, PP, QQ, TT, RR, GG, BB

 if S = 0
   R% = int(V * 255)
   G% = R%
   B% = R%
 else
   ZZ  = H / 60
   II% = int(ZZ)
   FF  = ZZ - int(ZZ)
   PP  = V * (1 - S)
   QQ  = V * (1 - S * FF)
   TT  = V * (1 - S * (1 - FF))

   select II%
     case 0
       RR = V  : GG = TT : BB = PP
     case 1
       RR = QQ : GG = V  : BB = PP
     case 2
       RR = PP : GG = V  : BB = TT
     case 3
       RR = PP : GG = QQ : BB = V
     case 4
       RR = TT : GG = PP : BB = V
     case 5
       RR = V  : GG = PP : BB = QQ
   end_select

   R% = int(RR * 255)
   G% = int(GG * 255)
   B% = int(BB * 255)
 end_if
end_sub

' --------------------------------------------------------------
' Procedures de la bibliotheque (nombres complexes)
' --------------------------------------------------------------

sub CMul(a_x, a_y, b_x, b_y)
' Multiplication : r_x + i r_y = (a_x + i a_y) * (b_x + i b_y)

 ErrCode% = 0

 r_x = a_x * b_x - a_y * b_y
 r_y = a_x * b_y + a_y * b_x
end_sub

sub CDiv(a_x, a_y, b_x, b_y)
' Division : r_x + i r_y = (a_x + i a_y) / (b_x + i b_y)
' Algorithme d'apres "Numerical Recipes"

 dim_local q, t

 if b_x = 0 and b_y = 0
   ErrCode% = -3
   r_x = MaxNum
   r_y = MaxNum
 else
   ErrCode% = 0
   if abs(b_x) >= abs(b_y)
     q = b_y / b_x
     t = b_x + b_y * q
     r_x = (a_x + a_y * q) / t
     r_y = (a_y - a_x * q) / t
   else
     q = b_x / b_y
     t = b_x * q + b_y
     r_x = (a_x * q + a_y) / t
     r_y = (a_y * q - a_x) / t  
   end_if    
 end_if
end_sub

sub CAbs(a_x, a_y)
' Module : r_mod = |a_x + i a_y|
' Algorithme d'apres "Numerical Recipes"

 ErrCode% = 0

 dim_local AbsX, AbsY, R, C

 AbsX = abs(a_x)
 AbsY = abs(a_y)

 if a_x = 0
   r_mod = abs(a_y)
 else
   if a_y = 0
     r_mod = abs(a_x)
   else
     if AbsX > AbsY
       R = AbsY / AbsX
       C = AbsX
     else
       R = AbsX / AbsY
       C = AbsY
     end_if
     r_mod = C * sqr(1 + R * R)
   end_if
 end_if
end_sub

sub CArg(a_x, a_y)
' Argument : r_arg = arg(a_x + i a_y)
' Resultat dans [-Pi, Pi)
' Equivaut a atan2(a_y, a_x)

 ErrCode% = 0

 if a_x = 0
   r_arg = sgn(a_y) * PiDiv2
 else
   ' 4e / 1er quadrant : -Pi/2..Pi/2
   r_arg = atn(a_y / a_x)
   if a_x < 0
     if a_y > 0
       ' 2e quadrant : Pi/2..Pi
       r_arg = r_arg + Pi
     else
       ' 3e quadrant : -Pi..-Pi/2
       r_arg = r_arg - Pi
     end_if
   end_if
 end_if
end_sub

sub CRealPower(a_x, a_y, p)
' Puissance (exposant reel) : (a_x + i a_y)^p
' Resultat dans r_x, r_y
' Resultat aussi dans r_mod, r_arg si a <> 0

 ErrCode% = 0

 if a_x = 0 and a_y = 0
   if p = 0
     ' 0^0 = lim x^x quand x --> 0 = 1
     r_x = 1
     r_y = 0
   else
     if p > 0
       ' 0^p = 0 si p > 0
       r_x = 0
       r_y = 0
     else
       ' 0^p indefini si p < 0
       ErrCode% = -2
       r_x = MaxNum
       r_y = MaxNum
     end_if
   end_if
 else
   CAbs(a_x, a_y)
   CArg(a_x, a_y)
   r_mod = power(r_mod, p)
   r_arg = r_arg * p
   r_x = r_mod * cos(r_arg)
   r_y = r_mod * sin(r_arg)
 end_if
end_sub
Revenir en haut Aller en bas
Voir le profil de l'utilisateur http://www.unilim.fr/pages_perso/jean.debord/index.htm
jean_debord

avatar

Nombre de messages : 766
Age : 63
Localisation : Limoges
Date d'inscription : 21/09/2008

MessageSujet: Re: Fractale de Newton   Mar 24 Juin 2014 - 10:15

Version pour exposant complexe :

Code:

' **************************************************************
' Fractale de Newton : Fonction z^p - 1 (p complexe)
' **************************************************************

' Variables definies par l'utilisateur

dim PicWidth%  : PicWidth%  = 640   : ' Taille de l'image
dim PicHeight% : PicHeight% = 480   : '   en pixels
dim p_x, p_y   : p_x = 3 : p_y = 2  : ' Exposant (z^p - 1)
dim x0         : x0 = 0             : ' Coordonnees du centre
dim y0         : y0 = 0             : '   de l'image
dim Eps        : Eps = 0.000001     : ' Precision requise
dim MaxIter%   : MaxIter% = 100     : ' Nb maxi d'iterations
dim ZoomFact   : ZoomFact = 1       : ' Facteur de zoom
dim ColorFact  : ColorFact = -2     : ' Facteur de coloration
dim CstVal     : CstVal = 0.9       : ' Luminosite HSV constante

' Variables supplementaires

dim HalfPicWidth%  : HalfPicWidth%  = PicWidth% / 2
dim HalfPicHeight% : HalfPicHeight% = PicHeight% / 2
dim Eps2           : Eps2 = Eps * Eps

dim R%, G%, B%  : ' Parametres de coloration
dim AbsColor    : ' abs(ColorFact)
dim ScaleFact   : ' Facteur d'echelle = distance entre 2 pixels
dim Nx%, Ny%    : ' Coordonnees d'un point (pixels)
dim xt, yt      : ' Coordonnees d'un point (algebriques)
dim f_x, f_y    : ' Fonction complexe f(z)
dim fp_x, fp_y  : ' Derivee f'(z)

' --------------------------------------------------------------
' Variables globales de la bibliotheque (nombres complexes)
' --------------------------------------------------------------

' Constantes mathematiques

dim MaxNum, MinNum, MaxLog, MinLog, Pi, PiDiv2

MaxLog =  709.78     : ' Argument max. pour EXP
MinLog = -708.39     : ' Argument min. pour EXP
MaxNum = exp(MaxLog) : ' Nb reel max. ~ 2^1024
MinNum = exp(MinLog) : ' Nb reel min. ~ 2^(-1022)
Pi     = 4 * atn(1)
PiDiv2 = Pi / 2

' Resultats des calculs
' Partie reelle, partie imaginaire, module, argument, signe

dim r_x, r_y, r_mod, r_arg, r_sgn

' Code d'erreur
'  0 = pas d'erreur
' -1 = argument hors bornes
' -2 = singularite
' -3 = overflow
' -4 = underflow

dim ErrCode%

' --------------------------------------------------------------
' Description des objets
' --------------------------------------------------------------

' Fenetre principale

left    0, 50
top     0, 50
width   0, PicWidth% + 70
height  0, PicHeight% + 120
caption 0, "Fractale de Newton ... Veuillez patienter."

' Image

picture 1
left    1, 30
top     1, 40
width   1, PicWidth%
height  1, PicHeight%

2d_target_is 1

' --------------------------------------------------------------
' Programme principal
' --------------------------------------------------------------

ColorFact = 0.01 * ColorFact
AbsColor  = abs(ColorFact)
ScaleFact = 4 / (PicHeight% * ZoomFact)

for Ny% = 0 to PicHeight% - 1
 yt = y0 - ScaleFact * (Ny% - HalfPicHeight%)
 for Nx% = 0 to PicWidth% - 1
   xt = x0 + ScaleFact * (Nx% - HalfPicWidth%)
   Newton(xt, yt)
   2d_pen_color R%, G%, B%
   2d_point Nx%, Ny%
 next Nx%
next Ny%

file_save 1, "newton.bmp"

caption 0, "Fractale de Newton ... Terminé."

end

' --------------------------------------------------------------
' Sous-programmes
' --------------------------------------------------------------

sub Func(x, y)
' Calcul de la fonction complexe f(z) = (f_x, f_y)
' et de sa derivee f'(z) = (fp_x, fp_y)
' Ici f(z) = z^p - 1 et f'(z) = p z^(p - 1)

 dim_local u_x, u_y        : ' u = z^(p - 1)

 CPower(x, y, p_x - 1, p_y)
 
 u_x = r_x
 u_y = r_y
 
 CMul(u_x, u_y, x, y)      : ' z^p

 f_x = r_x - 1
 f_y = r_y
 
 CMul(p_x, p_y, u_x, u_y)  : ' p z^(p - 1)

 fp_x = r_x
 fp_y = r_y
end_sub

sub Newton(xt, yt)
' Iteration de la fonction complexe au point (xt, yt)

 dim_local Iter%, mf2, mq2, x, y, zn_x, zn_y, q_x, q_y

 zn_x = xt
 zn_y = yt

 Iter% = 0

 repeat
   x = zn_x
   y = zn_y

   Func(x, y)

   if fp_x = 0 and fp_y = 0
     Iter% = MaxIter%             : ' Arret du calcul si f'(z) = 0
   else
     CDiv(f_x, f_y, fp_x, fp_y)   : ' q = f(z) / f'(z)
     
     q_x = r_x
     q_y = r_y

     zn_x = x - q_x
     zn_y = y - q_y

     mf2 = f_x * f_x + f_y * f_y  : ' |f(z)|^2
     mq2 = q_x * q_x + q_y * q_y  : ' |q|^2

     Iter% = Iter% + 1
   end_if
 until (mf2 < Eps2 and mq2 < Eps2) or (Iter% > MaxIter%)

 if Iter% > MaxIter%  : ' Pas de convergence
   R% = 0
   G% = 0
   B% = 0
 else
   MdbCol(Iter%)
 end_if
end_sub

sub MdbCol(Iter%)
' Determine la teinte (Hue, H) et la Saturation (S)
' d'apres le nombre d'iterations

 dim_local Q, Angle, Radius, H, S, V

 Q = log(Iter%) * AbsColor

 if Q < 0.5
   Q = 1 - 1.5 * Q
   Angle = 1 - Q
 else
   Q = 1.5 * Q - 0.5
   Angle = Q
 end_if

 Radius = sqr(Q)

' Si ColorFact > 0, assombrir une bande sur 2

 V = CstVal

 if (ColorFact > 0) and (odd(Iter%) > 0)
   V = 0.85 * V
   Radius = 0.667 * Radius
 end_if

 H = Angle * 10
 H = H - int(H)
 H = H * 360

 S = Radius - int(Radius)

' Convertir HSV en RGB

 HSVtoRGB(H, S, V)
end_sub

sub HSVtoRGB(H, S, V)
' Conversion HSV --> RGB. Resultats dans R%, G%, B%

 dim_local II%, ZZ, FF, PP, QQ, TT, RR, GG, BB

 if S = 0
   R% = int(V * 255)
   G% = R%
   B% = R%
 else
   ZZ  = H / 60
   II% = int(ZZ)
   FF  = ZZ - int(ZZ)
   PP  = V * (1 - S)
   QQ  = V * (1 - S * FF)
   TT  = V * (1 - S * (1 - FF))

   select II%
     case 0
       RR = V  : GG = TT : BB = PP
     case 1
       RR = QQ : GG = V  : BB = PP
     case 2
       RR = PP : GG = V  : BB = TT
     case 3
       RR = PP : GG = QQ : BB = V
     case 4
       RR = TT : GG = PP : BB = V
     case 5
       RR = V  : GG = PP : BB = QQ
   end_select

   R% = int(RR * 255)
   G% = int(GG * 255)
   B% = int(BB * 255)
 end_if
end_sub

' --------------------------------------------------------------
' Procedures de la bibliotheque (nombres complexes)
' --------------------------------------------------------------

sub CMul(a_x, a_y, b_x, b_y)
' Multiplication : r_x + i r_y = (a_x + i a_y) * (b_x + i b_y)

 ErrCode% = 0

 r_x = a_x * b_x - a_y * b_y
 r_y = a_x * b_y + a_y * b_x
end_sub

sub CDiv(a_x, a_y, b_x, b_y)
' Division : r_x + i r_y = (a_x + i a_y) / (b_x + i b_y)
' Algorithme d'apres "Numerical Recipes"

 dim_local q, t

 if b_x = 0 and b_y = 0
   ErrCode% = -3
   r_x = MaxNum
   r_y = MaxNum
 else
   ErrCode% = 0
   if abs(b_x) >= abs(b_y)
     q = b_y / b_x
     t = b_x + b_y * q
     r_x = (a_x + a_y * q) / t
     r_y = (a_y - a_x * q) / t
   else
     q = b_x / b_y
     t = b_x * q + b_y
     r_x = (a_x * q + a_y) / t
     r_y = (a_y * q - a_x) / t  
   end_if    
 end_if
end_sub

sub CAbs(a_x, a_y)
' Module : r_mod = |a_x + i a_y|
' Algorithme d'apres "Numerical Recipes"

 ErrCode% = 0

 dim_local AbsX, AbsY, R, C

 AbsX = abs(a_x)
 AbsY = abs(a_y)

 if a_x = 0
   r_mod = abs(a_y)
 else
   if a_y = 0
     r_mod = abs(a_x)
   else
     if AbsX > AbsY
       R = AbsY / AbsX
       C = AbsX
     else
       R = AbsX / AbsY
       C = AbsY
     end_if
     r_mod = C * sqr(1 + R * R)
   end_if
 end_if
end_sub

sub CArg(a_x, a_y)
' Argument : r_arg = arg(a_x + i a_y)
' Resultat dans [-Pi, Pi)
' Equivaut a atan2(a_y, a_x)

 ErrCode% = 0

 if a_x = 0
   r_arg = sgn(a_y) * PiDiv2
 else
   ' 4e / 1er quadrant : -Pi/2..Pi/2
   r_arg = atn(a_y / a_x)
   if a_x < 0
     if a_y > 0
       ' 2e quadrant : Pi/2..Pi
       r_arg = r_arg + Pi
     else
       ' 3e quadrant : -Pi..-Pi/2
       r_arg = r_arg - Pi
     end_if
   end_if
 end_if
end_sub

sub CExp(a_x, a_y)
' Exponentielle complexe : r_x + i r_y = exp(a_x + i a_y)

 dim_local ExpX

 if a_x < MinLog
   ErrCode% = -4
   r_x = 0
   r_y = 0
 else
   if a_x > MaxLog
     ErrCode = -3
     ExpX = MaxNum
   else
     ErrCode% = 0
     ExpX = exp(a_x)
   end_if
   r_x = ExpX * cos(a_y)
   r_y = ExpX * sin(a_y)
 end_if
end_sub

sub CPower(a_x, a_y, b_x, b_y)
' Puissance (exposant complexe) : (a_x + i a_y)^(b_x + i b_y)
' Resultat dans r_x, r_y

 ErrCode% = 0

 if a_x = 0 and a_y = 0
   if b_x = 0 and b_y = 0
     ' 0^0 = lim x^x quand x --> 0 = 1
     r_x = 1
     r_y = 0
   else
     ' 0^p = 0 si p > 0
     r_x = 0
     r_y = 0
   end_if
 else
   ' exp(b ln(a))
   CAbs(a_x, a_y)
   CArg(a_x, a_y)
   CMul(b_x, b_y, log(r_mod), r_arg)
   CExp(r_x, r_y)
 end_if
end_sub
Revenir en haut Aller en bas
Voir le profil de l'utilisateur http://www.unilim.fr/pages_perso/jean.debord/index.htm
Contenu sponsorisé




MessageSujet: Re: Fractale de Newton   

Revenir en haut Aller en bas
 
Fractale de Newton
Voir le sujet précédent Voir le sujet suivant Revenir en haut 
Page 1 sur 1
 Sujets similaires
-
» Fractale de Newton
» (m) Harvey Newton-Haydon + I hate everything about you.
» Apophysis de septembre
» fractales de Mai
» fractales d'Avril

Permission de ce forum:Vous ne pouvez pas répondre aux sujets dans ce forum
FORUM DE DISCUSSION SUR LE LANGAGE PANORAMIC :: PANORAMIC Le compilateur :: Le Compilateur-
Sauter vers: