Fractal de Newton-Raphson

Este fractal se crea usando el Metodo de Newton-Raphson y su forma de generarlo es similar al de Julia, salvo su ecuación, el código esta en VisualC++ 5.0, su ecuación es la siguiente:

f(z) = z³ - 1 = 0

donde z es una varieble compleja z=x+yi. El Metodo de Newton-Raphson se define de la siguiente manera :

x_n = x_n-1 - f(x_n-1)/f '(x_n-1)

para todo n mayor o igual que 1, y f ' es la derivada de la ecuación, este metodo se puede generalizar para el campo de las variables complejas. Aqui les va la forma de obtener las ecuaciones para obtener el Fractal de Newton.

f(z) = z3 - 1
f'(z) = 3z2, derivada de f(z)
z = z - f(z)/f'(z)
z = z - (z3 - 1)/3z2
  = (3z3 - (z3 - 1))/3z2
  = (2z3 + 1)/3z2
  = (2/3)*z + 1/3z2

como z=x+yi, la ecuación anterior nos queda, operando y haciendo la division

z = (2/3)*(x+yi) + 1/3*(x+yi)2
  = (2/3)*(x+yi) + 1/3*((x2 - y2) + 2*x*yi)
  = (2/3)*(x+yi) + ((x2 - y2) - 2*x*yi)/3*((x2 - y2)2 + 4*x2*y2)

para simplificar un poco la ecuación llamo d=3*((x2 - y2)2 + 4*x2*y2)

z = (2/3)*(x+yi) + ((x2 - y2) - 2*x*yi)/d

separando la parte real e imaginaria de

z = [(2/3)*x + (x2 - y2)/d] + [(2/3)*y - 2*x*y/d]i

me quedan finalmente las ecuaciones, a tengo que cuidar que d no sea cero
en el código, le asigno a d 0.000001, cuando este vale cero, es un truco para
prevenir la divición por cero, a y tambien utilizo una variable temporal 
tmp para guardar el ultimo valor de x.

tmp = x
x = (2/3)*x + (x2 - y2)/d
y = (2/3)*y - 2*tmp*y/d

Esta imagen es la salida del programa, el rango es de [-2.0, -2.0], [2.0, 2.0] y el número de iteraciones es de 2048.

Código Fuente

Aqui tienen una parte del código, el resto se encuentar en newton.zip (3 Kb), esta en VisualC++ 5.0.

//=== Función del procedimiento de ventana WndProc() === LRESULT CALLBACK WndProc( HWND hwnd, UINT message, WPARAM wParam, LPARAM lParam ) { HDC hdc; PAINTSTRUCT ps; RECT rect; BYTE R, G, B; double deltax, deltay, x, y, tmp, xx, yy, d; float xmin = -2.0, ymin=-2.0, xmax=2.0, ymax=2.0; int maxiter = 2048; int maxcol, maxrow; int color, row, col, count; switch( message ) { // mensaje producido en la creación de la ventana case WM_CREATE: break; //dibuja el fractal de Newton-Raphson case WM_PAINT: CrearPaleta(); hdc = BeginPaint( hwnd, &ps ); GetClientRect(hwnd, &rect); maxcol = rect.right - rect.left ; maxrow = rect.bottom - rect.top; deltax = (xmax - xmin)/maxcol; deltay = (ymax - ymin)/maxrow; for (col=0; col<=maxcol; col++) { for (row=0; row<=maxrow; row++) { x = xmin + col * deltax; y = ymin + row * deltay; count=0; while (count<maxiter) { xx = x*x; yy = y*y; d = 3.0*((xx - yy)*(xx - yy) + 4.0*xx*yy); if (d == 0.0) d = 0.000001; tmp=x; x = (2.0/3.0)*x + (xx - yy)/d; y = (2.0/3.0)*y - 2.0*tmp*y/d; count+=1; } if (x>0.0) color = count%64; else { if ((x<-0.3) && (y>0.0)) color = (count%64) + 64; else color = (count%64) + 128; } R=ColorRGB[color][0]; G=ColorRGB[color][1]; B=ColorRGB[color][2]; SetPixel(hdc, col, row, RGB(R,G,B)); } } EndPaint( hwnd, &ps ); break; // mensaje producido al cerrar la ventana case WM_DESTROY: PostQuitMessage( 0 ); break; // resto de mensajes, dar una respuesta estándar. default: return( DefWindowProc( hwnd, message, wParam, lParam ) ); } return(0); }