Générer des fractales de Newton

24

Vous connaissez tous la méthode de Newton pour approximer les racines d'une fonction, n'est-ce pas? Mon objectif dans cette tâche est de vous présenter un aspect intéressant de cet algorithme.

L'algorithme de Newton ne converge que pour certaines, mais surtout pour les valeurs d'entrée complexes. Si vous imaginez la convergence de la méthode pour toutes les valeurs d'entrée sur le plan complexe, vous obtenez généralement une belle fractale comme celle-ci:

Fractale de Newton pour f (x) = x ^ 3-1 Image de Wikimedia commons

Caractéristiques

Le but de cette tâche est de générer de telles fractales. Cela signifie que vous obtenez un polynôme en entrée et que vous devez imprimer la fractale correspondante en tant qu'image dans un format de votre choix en tant que sortie.

Contribution

L'entrée est une liste de nombres complexes séparés par des espaces. Ils sont écrits dans le style <Real part><iImaginary part>, comme ce numéro: 5.32i3.05. Vous pouvez supposer que le nombre d'entrée n'a pas plus de 4 décimales et est inférieur à 1000. Le premier d'entre eux ne doit pas être nul. Par exemple, cela pourrait être une entrée pour votre programme:

1 -2i7,5 23,0004i-3,8 i12 0 5,1233i0,1

Les nombres sont interprétés comme les coefficients d'un polynôme, en commençant par la puissance la plus élevée. Dans le reste de cette description, le polynôme d'entrée est appelé P . L'entrée ci-dessus est égale à ce polynôme:

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

L'entrée peut vous venir soit de stdin, soit d'un argument passé au programme, soit d'une invite affichée à votre programme. Vous pouvez supposer que l'entrée ne contient aucun espace blanc de début ou de fin.

Le rendu

Vous devez rendre la fractale de la manière suivante:

  • Choisissez autant de couleurs que de racines de P plus une couleur supplémentaire pour la divergence
  • Pour chaque nombre dans le plan visible, déterminez si la méthode converge et si oui vers quelle racine. Colorez le point en fonction du résultat.
  • N'imprimez pas de règles ou d'autres objets de fantaisie
  • Imprimez un point noir aux points, qui sont les racines des polynômes d'orientation. Vous pouvez imprimer jusqu'à quatre pixels autour de chaque racine.
  • Trouvez un moyen de choisir le plan visible de manière à ce que toutes les racines se distinguent et se propagent largement à travers celui-ci, si possible. Bien qu'un placement parfait du cadre de sortie ne soit pas requis, je me réserve le droit de refuser d'accepter une réponse qui choisit le cadre d'une manière inacceptable, par exemple. toujours aux mêmes coordonnées, toutes les racines sont en un seul point, etc.
  • L'image de sortie doit avoir une taille de 1024 * 1024 pixels.
  • Le temps de rendu est de 10 minutes maximum
  • L'utilisation de valeurs à virgule flottante simple précision suffit

Sortie

La sortie doit être une image graphique raster dans un format de fichier de votre choix, lisible par un logiciel standard pour un système d'exploitation de marque X. Si vous souhaitez utiliser un format rare, pensez à ajouter un lien vers un site Web où l'on peut télécharger une visionneuse pour celui-ci.

Sortez le fichier sur stdout. Si votre langue ne prend pas en charge la mise de quelque chose sur stdout ou si vous trouvez cette option moins pratique, trouvez un autre moyen. De toute façon, il doit être possible d'enregistrer l'image générée.

Restrictions

  • Aucune bibliothèque de traitement d'image
  • Aucune bibliothèque génératrice de fractales
  • Le code le plus court gagne

Extensions

Si vous aimez cette tâche, vous pouvez essayer de colorer les points en fonction de la vitesse de convergence ou d'autres critères. J'aimerais voir des résultats intéressants.

FUZxxl
la source
6
Je ne suis pas tout à fait sûr de savoir si cela convient comme un golf de code. À mes yeux, la tâche est beaucoup trop complexe. Mais je peux me tromper.
Joey
5
@Joey: En effet. J'aimerais que ce soit moi - même un défi de code .
Joey Adams
2
... ou PPM d'ailleurs.
Joey
1
@Joey: Mon intention était de créer une tâche plutôt difficile, car beaucoup de gens n'aiment pas les tâches très faciles.
FUZxxl
1
Il est facilement décomposé en tâches distinctes, et si votre langue prend en charge nativement les nombres à virgule flottante complexes, vous pouvez enregistrer un gros morceau. J'ai une version non entièrement jouée fonctionnant à 1600 caractères, dont 340 sont la classe des nombres complexes. Il n'identifie pas encore les racines et n'utilise pas les couleurs, mais j'essaie de retrouver ce que je suppose être un bogue dans le code NR. (Trouver une racine de x ^ 3-1 à partir de -0,5 + 0,866i ne devrait certainement pas diverger!)
Peter Taylor

Réponses:

13

Python, 827 777 caractères

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

Les trouvailles affichent les limites (et les racines) en trouvant des points de convergence pour un groupe d'échantillons aléatoires. Il trace ensuite le graphique en calculant les points de convergence pour chaque point de départ et en utilisant une fonction de hachage pour obtenir des couleurs aléatoires pour chaque point de convergence. Regardez de très près et vous pouvez voir les racines marquées.

Voici le résultat de l'exemple de polynôme.

résultat par exemple polynôme

Keith Randall
la source
Bien! J'aime ça.
FUZxxl
14

Java, 1093 1058 1099 1077 caractères

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

L'entrée est des arguments de ligne de commande - par exemple, run java F 1 0 0 -1. La sortie est à stdout au format PPM (pixmap ASCII).

L'échelle est choisie en utilisant la borne de Fujiwara sur la valeur absolue des racines complexes d'un polynôme; Je multiplie ensuite celui lié par 1,5. J'ajuste la luminosité par le taux de convergence, donc les racines seront dans les patchs les plus brillants. Par conséquent, il est logique d'utiliser du blanc plutôt que du noir pour marquer les emplacements approximatifs des racines (ce qui me coûte 41 caractères pour quelque chose qui ne peut même pas être fait "correctement". Si j'étiquette tous les points qui convergent à moins de 0,5 pixel d'eux-mêmes) alors certaines racines sortent sans étiquette; si j'étiquette tous les points qui convergent à moins de 0,6 pixels d'eux-mêmes, alors certaines racines sortent étiquetées sur plus d'un pixel; donc pour chaque racine j'étiquette le premier point rencontré pour converger à moins de 1 pixel de lui-même ).

Image pour l'exemple de polynôme (converti en png avec le GIMP): Racines de x ^ 5 + (- 2 + 7.5i) x ^ 4 + (23.0004-3.8i) x ^ 3 + 12i x ^ 2 + (5.1233 + 0.1i)

Peter Taylor
la source
@FUZxxl, l'image provient de l'ancienne version. J'en téléchargerai un avec le taux de convergence plus tard. Mais le problème avec le marquage des racines est de déterminer quel pixel marquer. C'est le problème classique qu'avec la virgule flottante, vous ne pouvez pas utiliser des tests d'égalité exacts, vous devez donc comparer avec un epsilon. Par conséquent, je n'ai pas de valeurs "canoniques" pour les racines. Je pourrais marquer des pixels qui convergent en une seule étape, mais cela ne garantit rien, et je pourrais marquer un bloc de 4 pixels pour une seule racine.
Peter Taylor
@Peter Taylor: Comme vous le voyez, Keith Randall a également trouvé une solution à ce problème. J'ai ajouté cette exigence comme une difficulté supplémentaire. Une approche pour le faire serait de calculer le pixel le plus proche pour chaque racine, puis de vérifier que chaque pixel est égal à lui.
FUZxxl
@FUZxxl, vous n'avez pas compris mon point. "Le pixel le plus proche" d'une racine n'est pas bien défini. Cependant, je peux pirater quelque chose qui pourrait fonctionner pour tous les cas de test que vous lui lancez, et j'ai l'impression que cela vous rendra heureux. Je vais le colorer en blanc, pas en noir, car c'est plus logique.
Peter Taylor
@Peter Taylor: D'accord.
FUZxxl
6
Ma photo de profil devrait bientôt changer x^6-9x^3+8, soigneusement conçue en choisissant les racines puis en utilisant Wolfram Alpha pour simplifier le polynôme. Ok, j'ai triché en échangeant les teintes après dans le GIMP.
Peter Taylor
3

Python, 633 octets

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

Après les accélérations et l'embellissement (756 octets)

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

Le graphique ci-dessous concerne la fractale de Newton de la fonction log (z).

Fractale de Newton pour log (z)

Pain au fromage
la source
Vous pouvez utiliser des noms plus courts (1 caractère) et supprimer les espaces en combinant plusieurs lignes à l'aide de ;. Supprimez également tous les espaces possibles.
mbomb007
Certains golfs réguliers réduisent cela à seulement 353 octets ! Je ne l'ai pas testé (non matplotlibici), donc aucune garantie qu'il fonctionne toujours.
Khuldraeseth na'Barya