1
$\begingroup$

I'm trying to combine plotting a Mandelbrot set with Monte Carlo randomization to plot an equation using random points for complex number z, for the function $z^3-2z+2=0$. Below is the code I have wrote for each separately, the Mandelbrot code is functional already.

f[z_] := z^3 - 2 z + 2; range = 3;
Estimate[z_] := 
  Module[{iteration = 20, error = 10^(-10), breakif = 10^(-10), i = 1,
     z0 = z}, 
   While[i <= iteration, (function = f[z0]; 
     If[Abs[f'[z0]] < breakif, Break[]]; z1 = z0 - function/f'[z0]; 
     If[Abs[z1 - z0] < error, Break[]]; z0 = z1); i++]; i];
DensityPlot[
 Estimate[xc + I yc], {xc, -range, range}, {yc, -range, range}, 
 PlotPoints -> 50, Mesh -> False, Frame -> False, 
 ColorFunction -> (If[# != 1, Hue[#], Hue[0, 0, 0]] &)]

I am unsure if the easiest way to do this would be to somehow make z0 different every time? or to define xc and xy as random numbers within the range? Or if I will need to completely rewrite the code to determine the convergence or divergence at each random point chosen.

Below is the Monte Carlo code I have previously written to choose random points and determine if they are inside or outside of an ellipse (in order to estimate the area). I was told that this code would be useful but I am not sure how to adapt it, especially considering there isn't an explicit function for "d". Any tips about this problem are much appreciated - combining the Modules especially - if I should nest them or be able to combine all of the variables into one Module.

MonteCarloPi[n0_] :=   
  Module[{d, i}, n = n0; Pin = Pout = {}; 
   For[i = 1, i <= n, i++, X = RandomReal[{-1, 1}]; 
    Y = RandomReal[{-1, 1}]; d = X^2 + Y^2/0.5^2; 
    If[d <= 1, Pin = Append[Pin, {X, Y}], 
     Pout = Append[Pout, {X, Y}];];]; m = Length[Pin]; 
   k = Length[Pout]; p = m/n; approx = p*4.0;
   Return[approx];];

MonteCarloResults[n0_] := 
 Module[{}, MonteCarloPi[n0]; Pin = Map[Point, Pin]; 
  DOTSin = Graphics[{Red, PointSize[0.01], Pin}]; 
  Pout = Map[Point, Pout];
  DOTSout = Graphics[{Green, PointSize[0.01], Pout}];
  ellipse = Ellipsoid[{0, 0}, {1, 0.5}];
  ellipseline = 
   Graphics[{EdgeForm[{Thick, Blue}], FaceForm[None], ellipse}];

  Print["Points Generated n = ", n];
  Print["  Inside Points m = ", m];
  Print["  Outside Points k = ", k];
  Print["Approximation to \[Pi] \[TildeTilde] ", approx];
  Print["Order of error \!\(\*FractionBox[\(\[Pi]\), \
SqrtBox[\(n\)]]\) \[TildeTilde] ", Pi/Sqrt[n]];
  Print["      |\[Pi]-approx| \[TildeTilde] ", Abs[Pi - approx]];

  Show[DOTSin, DOTSout, ellipseline, Axes -> True, 
       Ticks -> {Range[0, 1, 0.5], Range[0, 1, 0.5]}]]
MonteCarloResults[1000]
$\endgroup$
1
  • $\begingroup$ Dude, I'm having trouble with number 5 too. $\endgroup$
    – user27583
    Commented Apr 7, 2015 at 23:51

1 Answer 1

1
$\begingroup$

Here's my attempt to plot the Mandelbrot set with Monte Carlo randomization:

mand = Compile[{{z0, _Complex}, {imax, _Integer}}, 
   Module[{z = z0, i = 0}, 
    While[i < imax && Abs[z] <= 2, z = z^3 - 2 z + 2; i++]; i], 
   Parallelization -> True, RuntimeAttributes -> Listable(*,
   CompilationTarget->"C"*)];

n = 10^6; range = 2; resolution = 700; iteration = 500;
zlst = RandomComplex[range {-(1 + I), 1 + I}, n];
pos = 1 + Floor[resolution Rescale[(*ReIm*){Re@#, Im@#} &@zlst]]\[Transpose];
val = mand[zlst, iteration]; // AbsoluteTiming
System`SetSystemOptions[
  "SparseArrayOptions" -> {"TreatRepeatedEntries" -> 1}];
mat = Normal@SparseArray[pos -> val]; // AbsoluteTiming
System`SetSystemOptions[
  "SparseArrayOptions" -> {"TreatRepeatedEntries" -> 0}];
ArrayPlot[mat, ColorFunction -> "AvocadoColors"]

Remember to add CompilationTarget->"C" if you have a C compiler installed. And you can use ReIm instead of {Re@#, Im@#} & if you are in v10.1. As to the "TreatRepeatedEntries" -> 1 trick you may want to refer to this post. Finally here is the resulting picture, hmm, not that interesting, I think:

enter image description here

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.