FAKULTET ELEKTROTEHNIKE I RACUNARSTVA
"Matematicko modeliranje pomocu 
                                                                             Wolframove matematike"                                             ©  (Mario Matijevic, ENE-ET )

Fraktalni oblici u numerickim aproksimacijama 

(Newton-Raphsonova metoda)



(Sir Isaac Newton, 1643-1727)

● Struktura programa (!)

Ulazni paramteri ⟶ IZVRSNI KOD - SEKCIJE ⟶ Izlazni parametri

Napomena: za ispravnu komunikaciju izmedju pojedinih sekcija (potprograma) potrebno je inicijalno izvrsiti prvu sekciju!
		(Prva sekcija-"Generiranje osnovnog fraktalnog skupa"- sadrzi pozive dodatnih matematickih paketa bez kojih ostale sekcije ne rade ispravno.)

● Generiranje osnovnog fraktalnog skupa			 			             	           Osnovni skup

<< LinearAlgebra`MatrixManipulation`
<< Graphics`Animation`
<< Graphics`ComplexMap`
<< Graphics`Arrow`
Off[General::spell1]
Off[General::spell]

f[z_] = 
        Input["Ucitaj kompleksni polinom u varijabli z stupnja \
⩾3 (npr. z^3-1 ):" ];
derivacija[z_] = f'[z];
nule = NSolve[f[z] == 0, z];
polovi = NSolve[derivacija[z] == 0, z];

\(\(\( (*reset*) \)\(
\)\(\(xgrid = 
        300;\)\(\ \ \ \ \ \ \ \ \)\( (*gustoca\ na\ uzorku\ ravnine*) \)\(
\)
    \(ygrid = 300;\)\(
\)
    \(xMin = \(-2.0\);\)\(
\)
    \(xMax = 2.0;\)\(
\)
    \(yMin = \(-2.0\);\)\(
\)
    \(yMax = 2.0;\)\(
\)
    \(ε = 
      10\^\(-4\); \  (*zadana\ tocnost\ numerickog\ algoritma*) 
dx = 
      N[\((xMax - xMin)\)/xgrid];\)\(\ \ \ \ \ \ \)\(
\)
    \(dy = 
        N[\((yMax - yMin)\)/
            ygrid];\)\(\ \ \ \ \ \ \ \ \)\(
\)
    \(data = \(zoom = ZeroMatrix[xgrid, ygrid]\);\)\(\ \)\)\)\)

dim1 = 
        Part[Dimensions[nule], 1];    
dim2 = Part[Dimensions[polovi], 1];

(* k  -  vektor nultocaka*) 
 (* 
      p  -  vektor polova*) 
lista1 = niz1 = k = tocke1 = 
              Range[dim1] ;   (*alokacija memorije*) 
lista2 = niz2 = p = tocke2 = 
              Range[dim2]  ;

\(For[i = 1, 
      i ≤ 
        dim1, \(i++\), 
	\(k[\([i]\)] = 
          nule[\([i, 1, 2]\)];\)]\)
\(For[i = 1, 
      i ≤ 
        dim2, \(i++\), 
	\(p[\([i]\)] = 
          polovi[\([i, 1, 2]\)];\)]\)

For[i = 1, 
      i ≤ 
        dim1, i++, 
lista1[[i]] = {Re[k[[i]]], 
          Im[k[[i]]]}]
For[i = 1, 
      i ≤ 
        dim2, i++, 
lista2[[i]] = {Re[p[[i]]], 
          Im[p[[i]]]}]

For[i = 1, 
        i ≤ dim1, i++, 
tocke1[[i]] = 
          Point[lista1[[i]]]];      

    For[i = 1, 
        i ≤ dim2, i++, 
tocke2[[i]] = 
          Point[lista2[[i]]]];

\(\(\( (*definicija\ kompleksnog\ dinamickog\ sustava*) \)\(
\)\(Clear[iterator]

    \(iterator[z_] = N[Simplify[z - \((f[z]/\(f'\)[z])\)]];\)\)\)\)

matrica = 
          Table[FixedPoint[iterator, (x + y I), 150], {y, yMin, 
              yMax - dy, dy}, {x, xMin, xMax - dx, dx}]; // 
      Timing

(*transformacija konvergencijske matrice u boju*) 
For[x = 1, 
        x ≤ xgrid, x++, 
For[y = 1, 
          y ≤ 
            ygrid, y++, {
boja = 
              0.123; 
For[j = 1, 
              j ≤ dim1, j++, 
If[
                Abs[matrica[[y, x]] - 
                      k[[j]]] < ε, (boja = 
                    j/dim1)]]; 
data[[y, x]] = 
              boja}
]] // Timing

eksponent = Exponent[f[z], z];
faktor = 
        Coefficient[f[z], (z^eksponent)];
provjera = f[z] - faktor*(z^eksponent);

    
If[
        NumberQ[provjera], {

linije = 
            Range[dim1]; 
For[i = 1, 
            i ≤ dim1, i++, 
            linije[[i]] = 
              Line[{{0, 0}, lista1[[i]]}]]; 
For[i = 1, 
            i ≤ dim1, i++, 
niz1[[i]] = 
              ListPlot[{{Re[k[[i]]], Im[k[[i]]]}}, 
                DisplayFunction -> Identity, 
                PlotStyle -> {PointSize[0.025], Hue[0]}, 
                Epilog -> {Hue[0.1], linije, 
                    Circle[{0, 0}, 
                      Abs[k[[1]]]]}]]}, 

{For[i = 1, 
            i ≤ dim1, i++, 
niz1[[i]] = 
              ListPlot[{{Re[k[[i]]], Im[k[[i]]]}}, 
                DisplayFunction -> Identity, 
                PlotStyle -> {PointSize[0.025], 
                    Hue[0]}]]}];

For[i = 1, 
          i ≤ dim2, i++, 
niz2[[i]] = 
            ListPlot[{{Re[p[[i]]], Im[p[[i]]]}}, 
              DisplayFunction -> Identity, 
              PlotStyle -> {PointSize[0.025], Hue[0.7]}]];

pregled=Show[{niz1, niz2},DisplayFunction -> Identity,Background -> GrayLevel[0.9],AspectRatio -> Automatic,ImageSize -> 300,AxesStyle -> Thickness[0.006],AxesLabel->{StyleForm["Re(z)  ⇒",FontSize -> 12,FontWeight -> "Plain",FontFamily -> "Times"],StyleForm["Im(z)", FontSize -> 12, 
                  FontWeight -> "Plain", 
                  FontFamily -> "Times"]},PlotRange -> {{-2, 2}, {-2, 2}}]

\(fraktal = 
      ListDensityPlot[data, Mesh -> False, ImageSize -> 300, 
        ColorFunction -> Hue, DisplayFunction -> Identity, 
        PlotLabel -> "polinom: f(z)=" f[z], 
        FrameStyle -> Thickness[0.008], 
        FrameTicks -> {{{0, "-2"}, {xgrid/4, "-1"}, {xgrid/
                  2, "0"}, {3 
                   xgrid/4, "1"}, {xgrid, "2"}}, {{0, "-2i"}, {ygrid/4, "-i"}, {ygrid/2, "0"}, {3 
                   ygrid/4, "i"}, {ygrid, "2i"}}}]\)

Print[""];
Print[
      StyleForm["Parametri fraktalnog skupa:", FontColor -> Hue[0.6], 
        FontSize -> 20, FontFamily -> "Times", 
        FontWeight -> "Bold"]]
Print[StyleForm["●  kompleksni polinom: f(z) = ",FontSize -> 14,FontFamily -> "Times",FontWeight -> "Plain"],TraditionalForm[f[z]]]
Print[StyleForm["●  derivacija polinoma: ((df (z))/dz) = ",FontSize -> 14,FontFamily -> "Times",FontWeight -> "Plain"],TraditionalForm[derivacija[z]]]
Print[StyleForm["●  kompleksni iterator: (Subscript[(z), n + 1]) = ",FontSize -> 14,FontFamily -> "Times",FontWeight -> "Plain"],TraditionalForm[iterator[z]],StyleForm[",  z = (Subscript[(z), n])",FontSize -> 12,FontFamily -> "Times",FontWeight -> "Plain"]]
Print[StyleForm["●  nultocke polinoma: ",FontSize -> 14,FontFamily -> "Times",FontWeight -> "Plain"],TraditionalForm[MatrixForm[nule]]]
Print[StyleForm["●  polovi iteratora: ",FontSize -> 14,FontFamily -> "Times",FontWeight -> "Plain"],TraditionalForm[MatrixForm[polovi]] ]
Print[""]
Show[GraphicsArray[{pregled, fraktal},ImageSize -> 850,GraphicsSpacing -> 0.1,Frame -> False,
Epilog->{Text[StyleForm["● ⟶ nultočke polinoma",FontSize -> 12,FontWeight -> "Plain",FontFamily -> "Times"],{0.122, 0.95}],
Text[StyleForm["● ⟶ polovi iteratora",FontSize -> 12,FontWeight -> "Plain",FontFamily -> "Times"],{0.1, 0.91}]},
PlotLabel -> 
              StyleForm["Newton basins for: z" → 
                  TraditionalForm[f[z]], FontSize -> 15, 
                FontWeight -> "Bold", 
                FontFamily -> "Arial"],Background -> GrayLevel[0.9]]];

● Orbita kompleksne tocke								    	     KruznaAnimacija

Clear[tocka];
tocka = 
        Input["Ucitaj kompleksni broj (paste): "];
pixel = (xMin + 
              tocka[[1]]*dx) + (yMin + tocka[[2]]*dy)*
            I;
iteracija = 
        0;
maxiter = 
        150;
niz = 
      Range[maxiter];        
z0 = 
      pixel;
rez = 
        True;
broj = 0;
While[
        rez && (iteracija < maxiter), (g = iterator[z0]; 
          z0 = g; iteracija++; 
          For[i = 1, i ≤ dim1, i++, 
            rez = rez && (Abs[z0 - k[[i]]] > ε)]; 
          niz[[iteracija]] = 
            Graphics[{Hue[0.1], 
                Circle[{Re[g], Im[g]}, 1/((iteracija^2))]}, 
              Epilog -> {{PointSize[0.02], Hue[1], 
                    tocke1}, {PointSize[0.02], Hue[0.7], tocke2}}])]

For[i = 1, 
        i ≤ maxiter, i++, 
If[
          NumberQ[niz[[i]]], 
          Null, broj++]];
pom = 
        Range[broj];
For[i = 1, 
      i ≤ broj, i++, 
pom[[i]] = 
        niz[[i]]]
Print[StyleForm["→odabrana tocka: z = ",FontSize -> 14,FontWeight -> "Plain",FontFamily -> "Times"],TraditionalForm[pixel]]
Show[pom, 
        AspectRatio -> Automatic, AxesOrigin -> {0, 0}, 
        Axes -> True, PlotRange -> {{-2, 2}, {-2, 2}}, 
        Background -> GrayLevel[0.9], Frame -> False, 
        AxesLabel -> {StyleForm["Re(z)", FontSize -> 14, 
              FontWeight -> "Bold", FontFamily -> "Times"], 
            StyleForm["Im(z)", FontSize -> 14, 
              FontWeight -> "Bold", 
              FontFamily -> "Times"]}, 
        ImageSize -> 450];

Clear[tocka];
tocka = 
        Input["Ucitaj kompleksni broj (paste): "];
pixel = (xMin + 
              tocka[[1]]*dx) + (yMin + tocka[[2]]*dy)*
            I;
iteracija = 0;
maxiter = 150;
vektor = 
      Range[maxiter];        
z0 = pixel;
rez = True;
broj = 0;
\(\(While[
        rez && \((iteracija < maxiter)\), \((g = iterator[z0]; iza = z0; 
          z0 = g; \(iteracija++\); 
          For[i = 1, i ≤ dim1, \(i++\), 
            rez = rez && \((Abs[
                      z0 - k[\([i]\)]] > ε)\)]; 
vektor[\([iteracija]\)] = 
            ListPlot[{{Re[g], Im[g]}}, DisplayFunction -> Identity, 
              ImageSize -> 450, AspectRatio -> Automatic, 
              AxesOrigin -> {0, 0}, 
              PlotStyle -> {Hue[0.1], PointSize[0.015]}, 
              Background -> GrayLevel[0.9], Frame -> False, 
              PlotLabel -> 
                StyleForm["iteracija"⟶iteracija, 
                  FontSize -> 14, FontWeight -> "Bold", 
                  FontFamily -> "Times"], 
              PlotRange -> {{\(-2\), 2}, {\(-2\), 2}}, 
              Frame -> False, 
              Epilog -> {{PointSize[0.015], Hue[1], 
                    tocke1}, {PointSize[0.015], Hue[0.7], 
                    tocke2}, {Arrow[{Re[iza], Im[iza]}, {Re[g], 
                        Im[g]}]}}])\)];\)\)

For[i = 1, 
        i ≤ maxiter, i++, 
If[
          NumberQ[vektor[[i]]], 
          Null, broj++]];
pomocni = Range[broj];

For[i = 1, 
          i ≤ broj, i++, 
pomocni[[i]] = 
            vektor[[i]]];
  
Print[""];
Print[
      StyleForm["Parametri orbite tocke:", FontColor -> Hue[0.6], 
        FontSize -> 20, FontFamily -> "Times", 
        FontWeight -> "Bold"]]
Print[StyleForm["●  kompleksni polinom: f(z) = ",FontSize -> 14,FontFamily -> "Times",FontWeight -> "Plain"],TraditionalForm[f[z]]]
Print[StyleForm["●  odabrana tocka: z = ",FontSize -> 14,FontWeight -> "Plain",FontFamily -> "Times"],TraditionalForm[pixel]]
Print[""];
ListPlot[{{Re[pixel], Im[pixel]}}, 
      AspectRatio -> Automatic, AxesOrigin -> {0, 0}, 
      PlotStyle -> {Hue[0.1], PointSize[0.04]}, PlotLabel -> pixel, 
      Background -> GrayLevel[0.9], Frame -> False, 
      PlotRange -> {{-2, 2}, {-2, 2}}, 
      Epilog -> {Arrow[{0, 0}, {Re[pixel], Im[pixel]}]}, 
      ImageSize -> 200]
ShowAnimation[
        pomocni];

● Zoom explorer									           	     MagnifyReset

Clear[tocka, x, y];
tocka = 
        Input["Unesi centar uvecane slike (paste): "];
faktor = 
        Input["Unesi faktor uvecanja (npr. 70):"];
centar = (xMin + 
              tocka[[1]]*dx) + (yMin + tocka[[2]]*dy)*
            I;
pomak = faktor*dx;
T1 = (Re[centar] - pomak) + (Im[centar] + pomak)*
            I;
T2 = (Re[centar] + pomak) + (Im[centar] - pomak)*
            I;
  
xMax = Re[T2];
xMin = Re[T1];
yMax = Im[T1];
yMin = Im[T2];
dx = N[(xMax - xMin)/xgrid];
dy = N[(yMax - yMin)/ygrid];

zoom = 
          Table[FixedPoint[iterator, (x + y I), 150], {y, yMin, 
              yMax - dy, dy}, {x, xMin, xMax - dx, dx}]; // 
      Timing

      
For[x = 1, 
        x ≤ xgrid, x++, 
For[y = 1, 
          y ≤ 
            ygrid, y++, {
boja = 
              0.123; 
For[j = 1, 
              j ≤ dim1, j++, 
If[
                Abs[zoom[[y, x]] - 
                      k[[j]]] < ε, (boja = 
                    j/dim1)]]; 
data[[y, x]] = 
              boja}
]] // Timing

Print[""]
Print[
      StyleForm["Zoom explorer:", FontColor -> Hue[0.6], 
        FontSize -> 20, FontFamily -> "Times", 
        FontWeight -> "Bold"]]
Print[
      StyleForm["Faktor uvecanja = ", FontSize -> 14, 
        FontFamily -> "Times", FontWeight -> "Plain"], 
      TraditionalForm[faktor]]
Print[
      StyleForm["Centar uvecanja = ", FontSize -> 14, 
        FontFamily -> "Times", FontWeight -> "Plain"], 
      TraditionalForm[centar]]
Print[
      StyleForm["Velicina uzorka ℂ-ravnine:", 
        FontSize -> 14, FontFamily -> "Times", 
        FontWeight -> "Bold"]]
Print[
      StyleForm["X-max = ", FontSize -> 14, 
        FontFamily -> "Times", FontWeight -> "Plain"], 
      TraditionalForm[xMax]]
Print[
      StyleForm["X-min = ", FontSize -> 14, 
        FontFamily -> "Times", FontWeight -> "Plain"], 
      TraditionalForm[xMin]]
Print[
      StyleForm["Y-max = ", FontSize -> 14, 
        FontFamily -> "Times", FontWeight -> "Plain"], 
      TraditionalForm[yMax], 
      StyleForm["\[ImaginaryI]", FontFamily -> "Times", 
        FontWeight -> "Plain"]]
Print[
      StyleForm["Y-min = ", FontSize -> 14, 
        FontFamily -> "Times", FontWeight -> "Plain"], 
      TraditionalForm[yMin], 
      StyleForm["\[ImaginaryI]", FontFamily -> "Times", 
        FontWeight -> "Plain"]]
Print[""]
uvecanje = 
      ListDensityPlot[data, Mesh -> False, ColorFunction -> Hue, 
        ImageSize -> 400, Frame -> True, 
        Background -> GrayLevel[0.9], 
        PlotLabel -> 
          StyleForm["Newton zoom", FontSize -> 14, 
            FontWeight -> "Bold", FontFamily -> "Times"], 
        FrameStyle -> Thickness[0.01], 
        FrameTicks -> {{{0, xMin}, {xgrid, xMax}}, {{0, yMin*I}, {ygrid, 
                yMax*I}}}]
Null

● Kompleksno preslikavanje						            	 	           	           Preslikavanje

SetOptions[Graphics, 
        Frame -> True];
cartesian=CartesianMap[f,{-2, 2},{-2, 2},DisplayFunction -> Identity,PlotStyle -> {Hue[0.7]},Lines -> 30,AxesStyle -> Thickness[0.008],FrameStyle -> Thickness[0.008],PlotLabel->"CARTESIAN  (x, y)"];
polar=PolarMap[f,{0, (Sqrt[8])},{0, 2  Pi},DisplayFunction -> Identity,PlotStyle -> {Hue[0.7]},Lines -> 30,AxesStyle -> Thickness[0.008],FrameStyle -> Thickness[0.008],PlotLabel->"POLAR  (r, φ)"]
zajedno = 
        GraphicsArray[{cartesian, polar}];
Show[
      zajedno, ImageSize -> 800, GraphicsSpacing -> 0.2, 
      PlotLabel -> 
        StyleForm["Kompleksno preskivanje: z" → 
            TraditionalForm[f[z]], FontSize -> 15, 
          FontWeight -> "Bold", FontFamily -> "Arial"], 
      Background -> GrayLevel[0.9]]

● Korisni linkovi i literatura ⊳⊳⊳

MakeButtonForURL[u_String?(StringMatchQ[#, "http://*"] &)] := 
      DisplayForm[
        ButtonBox[u, ButtonData :> {URL[u], None}, 
          ButtonStyle -> "Hyperlink"]]
MakeButtonForURL["http://www.wolfram.com/webresources/"]
MakeButtonForURL["http://mathworld.wolfram.com/"]
MakeButtonForURL["http://scienceworld.wolfram.com/"]
MakeButtonForURL["http://staff.aes.rmit.edu.au/peter/FE/"]
MakeButtonForURL["http://www.fer.hr"]
\(MakeButtonForURL["http://spanky.triumf.ca/wwww/fractint/fractint/"]\)
MakeButtonForURL["http://www.astronomy.swin.edu.au/~pbourke/fractals/\
"]