Generar fractales de Newton

24

Todos ustedes conocen el método de Newton para aproximar las raíces de una función, ¿no? Mi objetivo en esta tarea es presentarle un aspecto interesante de este algoritmo.

El algoritmo de Newton converge solo para ciertos, pero sobre todo para valores de entrada complejos. Si imagina la convergencia del método para todos los valores de entrada en el plano complejo, generalmente obtiene un hermoso fractal como este:

Fractal de Newton para f (x) = x ^ 3-1 Imagen de Wikimedia commons

Presupuesto

El objetivo de esta tarea es generar tales fractales. Esto significa que obtienes un polinomio como entrada y tienes que imprimir el fractal correspondiente como una imagen en el formato que elijas como salida.

Entrada

La entrada es una lista de números complejos separados por espacios en blanco. Están escritas en el estilo <Real part><iImaginary part>, al igual que este número: 5.32i3.05. Puede suponer que el número de entrada no tiene más de 4 decimales y es menor que 1000. El primero de ellos no debe ser cero. Por ejemplo, esto podría ser una entrada a su programa:

1 -2i7.5 23.0004i-3.8 i12 0 5.1233i0.1

Los números se interpretan como los coeficientes de un polinomio, comenzando con la potencia más alta. Durante el resto de esta especificación, el polinomio de entrada se llama P . La entrada anterior es igual a este polinomio:

f (x) = x 5 + (-2 + 7.5 i ) x 4 + (23.0004 - 3.8 i ) x 3 + 12 i x 2 + 5.1233 + 0.1 i

La entrada puede venir de stdin, de un argumento pasado al programa o de un mensaje que se muestra a su programa. Puede suponer que la entrada no contiene ningún espacio en blanco inicial o final.

Representación

Tienes que renderizar el fractal de la siguiente manera:

  • Elija tantos colores como raíces de P más un color extra para divergencia
  • Para cada número en el plano visible, determine si el método converge y, en caso afirmativo, a qué raíz. Colorea el punto de acuerdo con el resultado.
  • No imprima reglas u otras cosas elegantes
  • Imprima un punto negro en los puntos, que son las raíces de los polinomios para la orientación. Puede imprimir hasta cuatro píxeles alrededor de cada raíz.
  • Encuentre una manera de elegir el plano visible de manera que todas las raíces sean distinguibles y se extiendan ampliamente, si es posible. Aunque no se requiere una colocación perfecta del marco de salida, me reservo el derecho de negarme a aceptar una respuesta que elija el marco de una manera inaceptable, por ejemplo. siempre en las mismas coordenadas, todas las raíces están en un punto, etc.
  • La imagen de salida debe tener un tamaño de 1024 * 1024 píxeles.
  • El tiempo de procesamiento es de 10 minutos como máximo.
  • Usar valores de punto flotante de precisión simple es suficiente

Salida

El resultado debe ser una imagen de gráficos de trama en un formato de archivo de su elección, legible por software estándar para un sistema operativo de marca X. Si desea utilizar un formato poco común, considere agregar un enlace a un sitio web donde se pueda descargar un visor.

Salida del archivo a stdout. Si su idioma no admite poner algo en stdout o si encuentra esta opción menos conveniente, busque otra forma. De cualquier manera, debe ser posible guardar la imagen generada.

Restricciones

  • No hay bibliotecas de procesamiento de imágenes
  • No hay bibliotecas generadoras de fractales
  • El código más corto gana

Extensiones

Si le gusta esta tarea, puede intentar colorear los puntos según la velocidad de convergencia o algún otro criterio. Me gustaría ver algunos resultados interesantes.

FUZxxl
fuente
66
No estoy completamente seguro de si esto es adecuado como código de golf. En mi opinión, la tarea es demasiado compleja. Sin embargo, puedo demostrar que estoy equivocado.
Joey
55
@Joey: De hecho. Me gustaría que esto sea un desafío de código para mí.
Joey Adams
2
... o PPM para el caso.
Joey
1
@Joey: Mi intención era crear una tarea bastante difícil, ya que a muchas personas no les gustan las tareas muy fáciles.
FUZxxl
1
Se divide fácilmente en tareas separadas, y si su idioma admite números complejos de coma flotante de forma nativa, puede guardar una gran parte. Tengo una versión no totalmente desarrollada con 1600 caracteres, de los cuales 340 son la clase de números complejos. Todavía no identifica las raíces y usa colores, pero estoy tratando de rastrear lo que presumo que es un error en el código NR. (¡Encontrar una raíz de x ^ 3-1 a partir de -0.5 + 0.866i seguramente no debería divergir!)
Peter Taylor

Respuestas:

13

Python, 827777 caracteres

import re,random
N=1024
M=N*N
R=range
P=map(lambda x:eval(re.sub('i','+',x)+'j'if 'i'in x else x),raw_input().split())[::-1]
Q=[i*P[i]for i in R(len(P))][1:]
E=lambda p,x:sum(x**k*p[k]for k in R(len(p)))
def Z(x):
 for j in R(99):
  f=E(P,x);g=E(Q,x)
  if abs(f)<1e-9:return x,1
  if abs(x)>1e5or g==0:break
  x-=f/g
 return x,0
T=[]
a=9e9
b=-a
for i in R(999):
 x,f=Z((random.randrange(-9999,9999)+1j*random.randrange(-9999,9999))/99)
 if f:a=min(a,x.real,x.imag);b=max(b,x.real,x.imag);T+=[x]
s=b-a
a,b=a-s/2,b+s/2
s=b-a
C=[[255]*3]*M
H=lambda x,k:int(x.real*k)+87*int(x.imag*k)&255
for i in R(M):
 x,f=Z(a+i%N*s/N+(a+i/N*s/N)*1j)
 if f:C[i]=H(x,99),H(x,57),H(x,76)
for r in T:C[N*int(N*(r.imag-a)/s)+int(N*(r.real-a)/s)]=0,0,0
print'P3',N,N,255
for c in C:print'%d %d %d'%c

Encuentra límites de visualización (y raíces) al encontrar puntos de convergencia para un montón de muestras aleatorias. Luego dibuja el gráfico calculando los puntos de convergencia para cada punto de partida y utilizando una función hash para obtener colores aleatorios para cada punto de convergencia. Mire muy de cerca y podrá ver las raíces marcadas.

Aquí está el resultado para el polinomio de ejemplo.

resultado, por ejemplo, polinomio

Keith Randall
fuente
¡Bueno! Me gusta esto.
FUZxxl
14

Java, 1093 1058 1099 1077 caracteres

public class F{double r,i,a,b;F(double R,double I){r=R;i=I;}F a(F c){return
new F(r+c.r,i+c.i);}F m(F c){return new F(r*c.r-i*c.i,r*c.i+i*c.r);}F
r(){a=r*r+i*i;return new F(-r/a,i/a);}double l(F c){a=r-c.r;b=i-c.i;return
Math.sqrt(a*a+b*b);}public static void main(String[]a){int
n=a.length,i=0,j,x,K=1024,r[]=new int[n];String o="P3\n"+K+" "+K+"\n255 ",s[];F z=new
F(0,0),P[]=new F[n],R[]=new F[n],c,d,e,p,q;for(;i<n;)P[i]=new
F((s=a[i++].split("i"))[0].isEmpty()?0:Float.parseFloat(s[0]),s.length==1?0:Float.parseFloat(s[1]));double
B=Math.pow(P[n-1].m(P[0].r()).l(z)/2,1./n),b,S;for(i=1;i<n;){b=Math.pow(P[i].m(P[i-1].r()).l(z),1./i++);B=b>B?b:B;}S=6*B/K;for(x=0;x<K*K;){e=d=c=new
F(x%K*S-3*B,x++/K*S-3*B);for(j=51;j-->1;){p=P[0];q=p.m(new
F(n-1,0));for(i=1;i<n;){if(i<n-1)q=q.m(c).a(P[i].m(new
F(n-1-i,0)));p=p.m(c).a(P[i++]);}c=c.a(d=q.r().m(p));if(d.l(z)<S/2)break;}i=j>0?0:n;for(;i<n;i++){if(R[i]==null)R[i]=c;if(R[i].l(c)<S)break;}i=java.awt.Color.HSBtoRGB(i*1f/n,j<1||e.l(c)<S&&r[i]++<1?0:1,j*.02f);for(j=0;j++<3;){o+=(i&255)+" ";i>>=8;}System.out.println(o);o="";}}}

La entrada es argumentos de línea de comandos, por ejemplo, ejecutar java F 1 0 0 -1 . La salida es stdout en formato PPM (mapa de píxeles ASCII).

La escala se elige utilizando el límite de Fujiwara en el valor absoluto de las raíces complejas de un polinomio; Luego multiplico ese límite por 1.5. Realmente ajusto el brillo por la tasa de convergencia, por lo que las raíces estarán en los parches más brillantes. Por lo tanto, es lógico usar blanco en lugar de negro para marcar las ubicaciones aproximadas de las raíces (lo que me está costando 41 caracteres para algo que ni siquiera se puede hacer "correctamente". Si etiqueto todos los puntos que convergen dentro de 0.5 píxeles de sí mismos entonces algunas raíces salen sin etiquetar; si etiqueto todos los puntos que convergen dentro de 0.6 píxeles de sí mismos, algunas raíces salen etiquetadas en más de un píxel; entonces, para cada raíz, etiqueto el primer punto encontrado para converger dentro de 1 píxel de sí mismo )

Imagen para el polinomio de ejemplo (convertido a png con el GIMP): Raíces de x ^ 5 + (- 2 + 7.5i) x ^ 4 + (23.0004-3.8i) x ^ 3 + 12i x ^ 2 + (5.1233 + 0.1i)

Peter Taylor
fuente
@FUZxxl, la imagen es de la versión anterior. Subiré uno con la tasa de convergencia más tarde. Pero el problema con marcar las raíces es determinar qué píxel marcar. Es el problema clásico que con el punto flotante no se pueden usar pruebas de igualdad exactas, por lo que se debe comparar con un épsilon. Como resultado, no tengo valores "canónicos" para las raíces. Podría marcar píxeles que convergen en un solo paso, pero eso no garantiza marcar nada, y podría marcar un bloque de 4 píxeles para una sola raíz.
Peter Taylor
@ Peter Taylor: Como puede ver, Keith Randall también encontró una solución a ese problema. Agregué este requisito como una dificultad adicional. Un enfoque para hacerlo sería calcular el píxel más cercano para cada raíz y luego verificar que cada píxel sea igual a él.
FUZxxl
@FUZxxl, no has entendido mi punto. "El píxel más cercano" de una raíz no está bien definido. Sin embargo, puedo hackear algo que podría funcionar para todos los casos de prueba que le arrojes, y tengo la impresión de que eso te hará feliz. Voy a colorearlo de blanco, no de negro, porque eso es más lógico.
Peter Taylor
@ Peter Taylor: De acuerdo.
FUZxxl
66
Mi foto de perfil pronto debería cambiar a x^6-9x^3+8, cuidadosamente diseñada eligiendo las raíces y luego usando Wolfram Alpha para simplificar el polinomio. Ok, hice trampa intercambiando los tonos después en el GIMP.
Peter Taylor
3

Python, 633 bytes

import numpy as np
import matplotlib.pyplot as plt
from colorsys import hls_to_rgb
def f(z):
    return (z**4 - 1)
def df(z):
    return (4*z**3) 
def cz(z):
    r = np.abs(z)
    arg = np.angle(z)   
    h = (arg + np.pi)  / (3 * np.pi)
    l = 1.0 - 1.0/(1.0 + r**0.1)
    s = 0.8 
    c = np.vectorize(hls_to_rgb) (h,l,s)
    c = np.array(c)
    c = c.swapaxes(0,2) 
    return c    
x, y = np.ogrid[-1.5:1.5:2001j, -1.5:1.5:2001j]
z = x + 1j*y    
for i in range(10):
    z -= (f(z) / df(z))
zz = z
zz[np.isnan(zz)]=0
zz=cz(zz)
plt.figure()
plt.imshow(zz, interpolation='nearest')
plt.axis('off')
plt.savefig('plots/nf.svg')
plt.close()

Después de Speed ​​Ups y Embellecimiento (756 bytes)

import numpy as np
from numba import jit
import matplotlib.pyplot as plt
from colorsys import hls_to_rgb 

@jit(nopython=True, parallel=True, nogil=True)
def f(z):
    return (z**4 - 1)   

@jit(nopython=True, parallel=True, nogil=True)
def df(z):
    return (4*z**3) 

def cz(z):
    r = np.abs(z)
    arg = np.angle(z)   

    h = (arg + np.pi)  / (3 * np.pi)
    l = 1.0 - 1.0/(1.0 + r**0.1)
    s = 0.8 

    c = np.vectorize(hls_to_rgb) (h,l,s)
    c = np.array(c)
    c = c.swapaxes(0,2) 
    return c    

x, y = np.ogrid[-1.5:1.5:2001j, -1.5:1.5:2001j]
z = x + 1j*y    

for i in range(10):
    z -= (f(z) / df(z))

zz = z
zz[np.isnan(zz)]=0
zz=cz(zz)
plt.figure()
plt.imshow(zz, interpolation='nearest')
plt.axis('off')
plt.savefig('plots/nf.svg')
plt.close()

El siguiente diagrama es para Newton Fractal de la función log (z).

Newton Fractal para log (z)

Pan de queso
fuente
Puede usar nombres más cortos (1 carácter) y eliminar espacios en blanco combinando varias líneas usando ;. Además, elimine todos los espacios posibles.
mbomb007
¡Algunos campos de golf regulares reducen esto a solo 353 bytes ! No lo he probado (no matplotlibaquí), así que no hay garantía de que todavía funcione.
Khuldraeseth na'Barya