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
» un nouveau editeur panobasic
par Jean Claude Aujourd'hui à 10:18

» Compilateur FBPano
par jean_debord Aujourd'hui à 10:07

» Le compilateur.
par Pedro Alvarez Aujourd'hui à 8:36

» COMPILATEUR V 0.9 beta 7 du 10 aout 2017
par Jack Hier à 20:23

» Pb 16 (en analyse): ON_CLOSE plante à l'exécution
par Jack Hier à 20:00

» Pb 15 (en analyse): TIMER_ON plante à l'exécution
par Jack Hier à 19:58

» KGF_dll - nouvelles versions
par Yannick Dim 13 Aoû 2017 - 17:35

» probleme d'outil
par Yannick Dim 13 Aoû 2017 - 17:32

» Carte de France des régions
par Yannick Sam 12 Aoû 2017 - 21:33

» Pb 14 (en analyse): PRINT_LOCATE plante à l'exécution
par Jack Ven 11 Aoû 2017 - 22:37

» Petit avertissement [Mots réservés]
par papydall Ven 11 Aoû 2017 - 13:45

» Distances sur plan
par JL35 Jeu 10 Aoû 2017 - 21:29

» Tracé : Triangle, Carrée, Dents de scie, Sinusoïde redressée
par papydall Jeu 10 Aoû 2017 - 14:52

» Troncature dans une image
par JL35 Mer 9 Aoû 2017 - 13:45

» A chacun son point de vue
par papydall Mar 8 Aoû 2017 - 17:20

Navigation
 Portail
 Index
 Membres
 Profil
 FAQ
 Rechercher
Rechercher
 
 

Résultats par :
 
Rechercher Recherche avancée
Août 2017
LunMarMerJeuVenSamDim
 123456
78910111213
14151617181920
21222324252627
28293031   
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 : 743
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 : 5501
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 : 5849
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 : 743
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 : 743
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 : 743
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
» Super Fractale : L'arbre de Pythagore
» Le binome de Newton (développer (a+b)^n)
» (m) Harvey Newton-Haydon + I hate everything about you.
» Apophysis de septembre

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: