Sie alle kennen die Newton-Methode, um die Wurzeln einer Funktion zu approximieren, nicht wahr? Mein Ziel bei dieser Aufgabe ist es, Ihnen einen interessanten Aspekt dieses Algorithmus vorzustellen.

Newtons Algorithmus konvergiert nur für bestimmte, vor allem aber für komplexe Eingabewerte. Wenn Sie sich die Konvergenz der Methode für alle Eingabewerte über die komplexe Ebene vorstellen, erhalten Sie normalerweise ein schönes Fraktal wie das folgende:

Bild aus Wikimedia Commons

Spezifikationen

Ziel dieser Aufgabe ist es, solche Fraktale zu generieren. Dies bedeutet, dass Sie ein Polynom als Eingabe erhalten und das entsprechende Fraktal als Bild in einem Format Ihrer Wahl als Ausgabe ausdrucken müssen.

Eingang

Die Eingabe ist eine durch Leerzeichen getrennte Liste komplexer Zahlen. Sie sind im Stil niedergeschrieben <Real part><iImaginary part>, wie diese Zahl ein : 5.32i3.05. Sie können davon ausgehen, dass die eingegebene Zahl nicht mehr als 4 Dezimalstellen hat und kleiner als 1000 ist. Die erste von ihnen darf nicht Null sein. Dies könnte beispielsweise eine Eingabe für Ihr Programm sein:

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

Die Zahlen werden als Koeffizienten eines Polynoms interpretiert, beginnend mit der höchsten Potenz. Im Rest dieser Beschreibung wird das Eingabepolynom P genannt . Die obige Eingabe entspricht diesem Polynom:

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

Die Eingabe kann entweder von stdin, von einem an das Programm übergebenen Argument oder von einer an Ihr Programm angezeigten Eingabeaufforderung kommen. Sie können davon ausgehen, dass die Eingabe keine führenden oder nachfolgenden Leerzeichen enthält.

Rendern

Sie müssen das Fraktal folgendermaßen rendern:

• Wählen Sie so viele Farben wie die Wurzeln von P und eine zusätzliche Farbe für die Divergenz
• Bestimmen Sie für jede Zahl in der sichtbaren Ebene, ob die Methode konvergiert und wenn ja, zu welcher Wurzel. Färben Sie den Punkt entsprechend dem Ergebnis.
• Drucken Sie keine Lineale oder andere ausgefallene Dinge
• Drucken Sie einen schwarzen Punkt an den Punkten, die die Polynomwurzeln zur Orientierung darstellen. Sie können bis zu vier Pixel um jede Wurzel drucken.
• Suchen Sie nach einer Möglichkeit, die sichtbare Ebene so zu wählen, dass alle Wurzeln unterscheidbar und möglichst weit verteilt sind. Obwohl eine perfekte Platzierung des Ausgaberahmens nicht erforderlich ist, behalte ich mir das Recht vor, die Annahme einer Antwort zu verweigern, die den Rahmen in nicht akzeptabler Weise auswählt, z. immer bei den gleichen Koordinaten, alle Wurzeln sind in einem Punkt usw.
• Das Ausgabebild sollte eine Größe von 1024 * 1024 Pixel haben.
• Die Renderzeit beträgt maximal 10 Minuten
• Die Verwendung von Gleitkommawerten mit einfacher Genauigkeit ist ausreichend

Ausgabe

Die Ausgabe sollte ein Raster-Grafikbild in einem Dateiformat Ihrer Wahl sein, das von Standardsoftware für ein Brand X-Betriebssystem gelesen werden kann. Wenn Sie ein seltenes Format verwenden möchten, sollten Sie einen Link zu einer Website hinzufügen, auf der Sie einen Viewer herunterladen können.

Gib die Datei auf stdout aus. Wenn Ihre Sprache es nicht unterstützt, etwas auf Standard zu setzen, oder wenn Sie diese Option als weniger praktisch empfinden, suchen Sie einen anderen Weg. In jedem Fall muss es möglich sein, das generierte Bild zu speichern.

Beschränkungen

• Keine Bildverarbeitungsbibliotheken
• Keine fraktalerzeugenden Bibliotheken
• Der kürzeste Code gewinnt

Erweiterungen

Wenn Ihnen diese Aufgabe gefällt, können Sie versuchen, die Punkte nach der Konvergenzgeschwindigkeit oder nach anderen Kriterien einzufärben. Ich würde gerne interessante Ergebnisse sehen.

Python, 827 777 Zeichen

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

Findet Anzeigegrenzen (und -wurzeln), indem Konvergenzpunkte für eine Reihe von Zufallsstichproben ermittelt werden. Anschließend wird das Diagramm gezeichnet, indem die Konvergenzpunkte für jeden Startpunkt berechnet und mit einer Hash-Funktion zufällige Farben für jeden Konvergenzpunkt ermittelt werden. Schauen Sie genau hin und Sie können die markierten Wurzeln sehen.

Hier ist das Ergebnis für das Beispielpolynom.

Java, 1093 1058 1099 1077 Zeichen

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="";}}}

Eingabe sind Kommandozeilenargumente - zB run java F 1 0 0 -1. Die Ausgabe erfolgt im PPM-Format (ASCII-Pixmap).

Die Skala wird unter Verwendung des Fujiwara gewählt, der an den Absolutwert der komplexen Wurzeln eines Polynoms gebunden ist; Ich multipliziere dann die Grenze mit 1,5. Ich passe die Helligkeit durch die Konvergenzrate an, so dass die Wurzeln in den hellsten Flecken sind. Daher ist es logisch, Weiß anstelle von Schwarz zu verwenden, um die ungefähren Positionen der Wurzeln zu markieren (das kostet mich 41 Zeichen für etwas, was nicht einmal "richtig" gemacht werden kann. Wenn ich alle Punkte beschrifte, die innerhalb von 0,5 Pixel von sich selbst konvergieren dann kommen einige Wurzeln unbeschriftet heraus; wenn ich alle Punkte beschrifte, die innerhalb von 0,6 Pixeln von sich selbst konvergieren, dann kommen einige Wurzeln über mehr als ein Pixel hinaus beschriftet, so beschrifte ich für jede Wurzel den ersten Punkt, der angetroffen wird, um innerhalb von 1 Pixel von sich selbst zu konvergieren ).

Bild für das Beispielpolynom (mit GIMP in PNG konvertiert):

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()

Nach Beschleunigung und Verschönerung (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()

Das folgende Diagramm ist für Newton Fractal der log (z) -Funktion.