Kalles Fraktaler
|
|
« on: November 26, 2013, 06:55:53 PM » |
|
I was going to implement the pertubation method to render 3rd degree Mandelbrot zooms in previously unbelievable speeds. However I found out that any reference outside <0,0i> distort the image to beyond recognition - even if it is still inside the main Mandelbrot "body". Also for second degree the main Mandelbrot get distorted - but if the reference is outside the Mandelbrot, but never inside it. Maybe the power of three just makes values unreachable with double float precision.
Here is what I came up with, inserting the 3rd degree function in master K.I.Martin's equations:
The function for a pixel is X[n+1] = X[n]^3 + X[0]
The function for a nearby pixel is then Y[n+1] = Y[n]^3 + Y[0]
The difference for any iteration is D[n] = Y[n]-X[n] D[n+1] = Y[n+1] - X[n+1]
Y for any iteration can be calculated from Y[n] = X[n] + D[n]
Replace Y[n+1] with how it was calculated from iteration n D[n+1] = (Y[n]^3 + Y[0]) - (X[n]^3 + X[0])
Replace Y with X and D D[n+1] = ((X[n]+D[n])^3 + X[0]+D[0]) - (X[n]^3 + X[0])
Expand the parentheses D[n+1] = X[n]^3 + 3*X[n]^2*D[n] + 3*X[n]*D[n]^2 + D[n]^3 + X[0] + D[0] - X[n]^3 - X[0]
Reduce D[n+1] = 3*X[n]^2*D[n] + 3*X[n]*D[n]^2 + D[n]^3 + D[0]
So the pertubation function for the real and imaginary parts would be Dr[n+1] = 3*Xr[n]*Xr[n]*Dr[n] - 6*Xr[n]*Xi[n]*Di[n] - 3*Xi[n]*Xi[n]*Dr[n] + 3*Xr[n]*Dr[n]*Dr[n] - 3*Xr[n]*Di[n]*Di[n] - 3*Xi[n]*2*Dr[n]*Di[n] + Dr[n]*Dr[n]*Dr[n] - 3*Dr[n]*Di[n]*Di[n] + Dr[0] Di[n+1] = 3*Xr[n]*Xr[n]*Di[n] + 6*Xr[n]*Xi[n]*Dr[n] - 3*Xi[n]*Xi[n]*Di[n] + 3*Xr[n]*2*Dr[n]*Di[n] + 3*Xi[n]*Dr[n]*Dr[n] - 3*Xi[n]*Di[n]*Di[n] + 3*Dr[n]*Dr[n]*Di[n] - Di[n]*Di[n]*Di[n] + Di[0]
And it works fine with the reference in <0,0i> see first image. The second image has it's reference very near but not at <0,0i> and get distorted. The third image has it's reference farther away but still inside the main Mandelbrot's body...
So... if anyone have any suggestions on how to do this in a workable way I would appreciate that a lot. Is there anything in the final functions that can be reduced to avoid double float bailout?
|
|
|
Logged
|
|
|
|
knighty
Fractal Iambus
Posts: 819
|
|
« Reply #1 on: November 26, 2013, 08:38:03 PM » |
|
I've checked your furmula and couldn't find any error. Maybe there is a mistyping somewhere in your code. With formulas as long as those it may happen. You could also do some factorisations: shorter, thus faster and less error prone. Judging from the equations it is kind of obvious that a reference point of 0 gives the right result because its orbit points are all zeroes so the perturbation formula reduces to the original one: D[n+1]=D[n]^3+D[0].
|
|
|
Logged
|
|
|
|
Kalles Fraktaler
|
|
« Reply #2 on: November 26, 2013, 10:48:37 PM » |
|
Thanks for checking. I hoped for a typo somewhere, but unfortunately it doesn't seem there is. Also Martin's equation (1) reduce to the standard Mandelbrot formula for reference point 0, so I think it should do so also here http://www.superfractalthing.co.nf/sft_maths.pdf
|
|
|
Logged
|
|
|
|
knighty
Fractal Iambus
Posts: 819
|
|
« Reply #3 on: November 27, 2013, 01:01:05 AM » |
|
That is interesting! You are right. I've done a quick experiment (without multiprecision) and got the same results. It seems that the perturbation formula is itself chaotic so very sensitive to rounding errors. //evladraw script static frm=-1; static sx=0., sy=1.01; (x,y) { if(frm<numframes){ frm=numframes; sy=mousx*0.001; calcseed(100); } /*cx=x;cy=y; for(i=0;i<50 && x*x+y*y<4;i++){ step(x,y,x,y,cx,cy); } i/50*/ perturbation(x,y,100) } step(x,y,&x1,&y1,cx,cy){ //(x+i*y)^3=x^3-3*x*y^2+i*(3*x^2*y-y^3); x1=x^3-3*x*y^2+cx; y1=3*x^2*y-y^3+cy; } static Seed[1024][2]; calcseed(n){ x=sx;y=sy; for(i=0;i<n;i++){ step(x,y,x,y,sx,sy); Seed[i][0]=x; Seed[i][1]=y; } } perturbation(x,y,n){ d0r=x-sx; d0i=y-sy; dr=d0r; di=d0i; for(i=0;i<n && x*x+y*y<4;i++){ x=Seed[i][0]; y=Seed[i][1]; //Dr1 = 3*x*x*Dr - 6*x*y*Di - 3*y*y*Dr + 3*x*Dr*Dr - 3*x*Di*Di - 6*y*Dr*Di + Dr*Dr*Dr - 3*Dr*Di*Di + D0r; //Di1 = 3*x*x*Di + 6*x*y*Dr - 3*y*y*Di + 6*x*Dr*Di + 3*y*Dr*Dr - 3*y*Di*Di + 3*Dr*Dr*Di - Di*Di*Di + D0i; dr1=3*((x^2-y^2)*dr+x*(dr^2-di^2)-di*(2*y*(x+dr)+dr*di))+dr^3+d0r; di1=3*((x^2-y^2)*di+y*(dr^2-di^2)+dr*(2*x*(y+di)+dr*di))-di^3+d0i; dr=dr1; di=di1; x+=dr; y+=di; } i/n }
|
|
|
Logged
|
|
|
|
knighty
Fractal Iambus
Posts: 819
|
|
« Reply #4 on: May 07, 2014, 07:54:44 PM » |
|
Ahem! That's really embarrassing It seems we made the same mistake. The error wasn't in the formula but here: static Seed[1024][2]; calcseed(n){ x=sx;y=sy; for(i=0;i<n;i++){ step(x,y,x,y,sx,sy); Seed[i][0]=x; Seed[i][1]=y; } } ... Seed[0] will contain the second iterate instead of the first. It should be like this: static Seed[1024][2]; calcseed(n){ x=0;y=0; for(i=0;i<n;i++){ step(x,y,x,y,sx,sy); Seed[i][0]=x; Seed[i][1]=y; } } ...or this: static Seed[1024][2]; calcseed(n){ x=sx;y=sy; for(i=0;i<n;i++){ Seed[i][0]=x; Seed[i][1]=y; step(x,y,x,y,sx,sy); } } Yeah, silly mistakes are the most difficult to find. It WORKS!... hope there won't be other issues. //evladraw script static frm=-1; static sx=0., sy=1.01; static bailout=2^(2/(3-1));//for f(z)=z^p+c minimum bailout radius is 2^(1/(p-1)); (x,y) { if(frm<numframes){ frm=numframes; sx=(mousx-0.5*xres)*0.001; sy=(0.5*yres-mousy)*0.001; calcseed(50); } /*cx=x;cy=y; for(i=0;i<50 && x*x+y*y<4;i++){ step(x,y,x,y,cx,cy); } i/50*/ perturbation(x,y,50) } step(x,y,&x1,&y1,cx,cy){ //(x+i*y)^3=x^3-3*x*y^2+i*(3*x^2*y-y^3); x1=x^3-3*x*y^2+cx; y1=3*x^2*y-y^3+cy; } static Seed[1024][2]; calcseed(n){ x=0;y=0; for(i=0;i<n;i++){ step(x,y,x,y,sx,sy); Seed[i][0]=x; Seed[i][1]=y; } } perturbation(x,y,n){ d0r=x-sx; d0i=y-sy; dr=d0r; di=d0i; for(i=0;i<n && x*x+y*y<bailout;i++){ x=Seed[i][0]; y=Seed[i][1]; //Dr1 = 3*x*x*Dr - 6*x*y*Di - 3*y*y*Dr + 3*x*Dr*Dr - 3*x*Di*Di - 6*y*Dr*Di + Dr*Dr*Dr - 3*Dr*Di*Di + D0r; //Di1 = 3*x*x*Di + 6*x*y*Dr - 3*y*y*Di + 6*x*Dr*Di + 3*y*Dr*Dr - 3*y*Di*Di + 3*Dr*Dr*Di - Di*Di*Di + D0i; dr1=3*((x^2-y^2)*dr+x*(dr^2-di^2)-di*(2*y*(x+dr)+dr*di))+dr^3+d0r; di1=3*((x^2-y^2)*di+y*(dr^2-di^2)+dr*(2*x*(y+di)+dr*di))-di^3+d0i; x+=dr; y+=di; dr=dr1; di=di1; //x+=dr; y+=di;//oops another mistake :D } i/n }
|
|
« Last Edit: May 08, 2014, 02:09:27 PM by knighty, Reason: anoter silly mistake »
|
Logged
|
|
|
|
cKleinhuis
|
|
« Reply #5 on: May 07, 2014, 08:33:16 PM » |
|
nice!
|
|
|
Logged
|
---
divide and conquer - iterate and rule - chaos is No random!
|
|
|
|
|
|
Kalles Fraktaler
|
|
« Reply #9 on: May 08, 2014, 09:19:45 PM » |
|
Next challenge - construct the Series Approximation functions!
Anyone?
|
|
|
Logged
|
|
|
|
|
Kalles Fraktaler
|
|
« Reply #11 on: May 09, 2014, 08:16:25 PM » |
|
But it is very difficult if you don't notice that Thanks a lot I will try it
|
|
|
Logged
|
|
|
|
Kalles Fraktaler
|
|
« Reply #12 on: May 12, 2014, 08:50:14 PM » |
|
Knighty, as you know by now, the formulas you assembled works super. On my tests the first term yielded for example 2000 iterations. The second - none! The third - 2. That is ok I suspect the impact of each term decrease much faster than on the 2nd degree Mandelbrot. However I have trouble understanding how you came up with the 'running derivates'. The derivata of x 3 = 3x 2 I follow but then I am lost. Can you explain please? Because if I could understand this I could apply it on other formulas as well
|
|
|
Logged
|
|
|
|
knighty
Fractal Iambus
Posts: 819
|
|
« Reply #13 on: May 12, 2014, 09:39:54 PM » |
|
That's just derivative rule for products of functions: (f(x)*g(x))' = f'(x)*g(x) + f(x)*g'(x) And the chain rule: have to be seen as a function of so we have in realtiy: and the derivative is w.r.t. . ... ; because 1 is a constant so it's derivative is 0 ; Product rule ; Chain rule applied to the term. On my tests the first term yielded for example 2000 iterations. The second - none! The third - 2.
I can't say what's exactly happening there... But it works! Just keep in mind that it will become more and more difficult to do higher exponents both for the perturbation formula and the series approximation: The formulas will become quite long when developping them to real an imaginary parts. (Maybe using a "Complex" class (the C++ complex template class?) would be useful at some point...)
|
|
|
Logged
|
|
|
|
Kalles Fraktaler
|
|
« Reply #14 on: May 12, 2014, 11:59:59 PM » |
|
Yeah, I have already start playing with higher exponents. And I gave up expanding the perturbation functions already at 4th degree. I'm only using the first term on SA, the one I understand... Because since the 2nd and 3rd terms' impact decreased so fast I don't think they will help at all on higher exponents. It needs more testing, but it works really well, beyond expectation, much slower of course but hopefully much faster than the traditional way Here is a 7th degree closeup at e19
|
|
|
Logged
|
|
|
|
|