Primeira vez aqui? Seja bem vindo e cheque o FAQ!
x

Como Implementar os chamados julia-sets, a partir da classe de números complexos ?

+2 votos
90 visitas
perguntada Jun 16, 2018 em Ciência da Computação por vanderson Delapedra (41 pontos)  
editado Jun 16, 2018 por vanderson Delapedra

b) Entenda a classe de números complexos desenvolvida nessas notas para lidar (explicitamente) com produtos de números complexos e o que mais for necessário para implementar os julia-sets.

Compartilhe

3 Respostas

+2 votos
respondida Jun 23, 2018 por MarcioGama (106 pontos)  

Segue adaptação da resposta usando orientação a objetos.

import numpy as np
from numba import autojit
import matplotlib.cm as cm
import matplotlib.pyplot as plt


im_width, im_height = 500, 500

zabs_max = 10000
nit_max = 1000

xmin, xmax = -1.5, 1.5
xwidth = xmax - xmin
ymin, ymax = -1.5, 1.5
yheight = ymax - ymin

fig, ax = plt.subplots()



class JuliaSet(object):

    def __init__(self, zmax, nmax, re, im):
        self.c = complex(re, im);
        self.zmax = zmax;
        self.nmax = nmax;
        self.julia = np.zeros((im_width, im_height));


    def calc(self):

        for ix in range(im_width):
            for iy in range(im_height):
                nit = 0

                z = complex(ix / im_width * xwidth + xmin,iy / im_height * yheight + ymin);

                while abs(z) <= self.zmax and nit < self.nmax:
                    z = z**2 + self.c
                    nit += 1
                ratio = nit / self.nmax
                self.julia[ix,iy] = ratio


julia = JuliaSet(zabs_max, nit_max, -0.70176, -0.365)
julia.calc();
ax.imshow(julia.julia, interpolation='nearest', cmap=cm.hot)

xtick_labels = np.linspace(xmin, xmax, xwidth / 0.5)
ax.set_xticks([(x-xmin) / xwidth * im_width for x in xtick_labels])
ax.set_xticklabels(['{:.1f}'.format(xtick) for xtick in xtick_labels])
ytick_labels = np.linspace(ymin, ymax, yheight / 0.5)
ax.set_yticks([(y-ymin) / yheight * im_height for y in ytick_labels])
ax.set_yticklabels(['{:.1f}'.format(ytick) for ytick in ytick_labels])


plt.show()
+1 voto
respondida Jun 16, 2018 por vanderson Delapedra (41 pontos)  

Os conjuntos denominados Julia-sets são uma extensão do Mandelbrot-set, demonstrado aqui . Na realidade eles não podem ser representados como uma figura única, tal como o Mandelbrot set; ao contrário disso, o conjunto Julia possui uma dimensão não fixa na parte dos números complexos da equação.

Ou seja, se no Mandelbrot-set, a cada iteração, a parte quadrática mudava enquanto a parte de números complexos se mantinha fixa em um ponto no plano complexo, para saber se ele convergia ou se dispersava para o infinito. No Julia set, ao invés de se substituir cada ponto do plano complexo pelo C, essa substituição ocorre apenas na parte inicial da equação, ou seja, no Z0.

Sendo assim, se inicialmente Z0= C, a próxima iteração será dada pela seguinte sequência: Z1= C^2 + C, a próxima, Z_2= [(C^(2 )+ C)^2 + C], e assim sucessivamente.

Uma vez que diferentes valores para C vão sendo introduzidos, diferentes formas vão surgindo, isso explica o motivo de uma aparente aleatoriedade inicial nas figuras e ao mesmo tempo a existência de uma lógica estruturada entre elas, preservando certa similaridade de replicação na medida em que se aumenta o zoom.

Para gerar um julia-set, será utilizada uma lógica aproximada da geração do mandelbrot-set.

A primeira coisa a fazer é importação do pacote numpy, pacote básico da linguagem Python que nos permite trabalhar com matrizes de N dimensões.

A função @autojit do pacote Numba é utilizado para dar agilidade no cálculo, uma vez que inúmeras simulações são realizadas com infinitas soluções. Além disso, a biblioteca matplotlib.pyplot possibilita a plotagem do resultado.

import numpy as np
from numba import autojit
import matplotlib.cm as cm
import matplotlib.pyplot as plt

O próximo passo necessário é o estabelecimento dos parâmetros da imagem onde será impressa a figura gerada pela equação.

im_width, im_height = 500, 500

Em seguida, cria-se valores, aleatoriamente, que irão compor o campo real e o imaginário dentro do Plano complexo. Note que cada ponto do plano complexo gerará um formato de figura.

c = complex(-0.70176, -0.365)

Em seguda, o número de iterações é informado

zabs_max = 10000
nit_max = 1000

Abaixo, ocorre a delimitação dos pixels dentro do plano complexo a serem utilizados

xmin, xmax = -1.5, 1.5
xwidth = xmax - xmin
ymin, ymax = -1.5, 1.5
yheight = ymax - ymin

A função que indicará o looping em torno da equação, é criada a seguir:

julia = np.zeros((im_width, im_height))
for ix in range(im_width):
    for iy in range(im_height):
        nit = 0

O ponto onde será impresso no plano complexo é identificado a partir da determinação do z

z = complex(ix / im_width * xwidth + xmin,
            iy / im_height * yheight + ymin)

Um comando while irá delimitar as iterações

while abs(z) <= zabs_max and nit < nit_max:
    z = z**2 + c
    nit += 1
ratio = nit / nit_max
julia[ix,iy] = ratio

fig, ax = plt.subplots()

ax.imshow(julia, interpolation='nearest', cmap=cm.hot)

Em seguida, são definidos os rótulos para as coordenadas de Zo no Plano complexo

xtick_labels = np.linspace(xmin, xmax, xwidth / 0.5)
ax.set_xticks([(x-xmin) / xwidth * im_width for x in xtick_labels])
ax.set_xticklabels(['{:.1f}'.format(xtick) for xtick in xtick_labels])
ytick_labels = np.linspace(ymin, ymax, yheight / 0.5)
ax.set_yticks([(y-ymin) / yheight * im_height for y in ytick_labels])
ax.set_yticklabels(['{:.1f}'.format(ytick) for ytick in ytick_labels])

Por fim, plota-se a imagem:

plt.show()

O código completo está expresso a seguir:

import numpy as np
from numba import autojit
import matplotlib.cm as cm
import matplotlib.pyplot as plt


im_width, im_height = 500, 500

c = complex(-0.70176, -0.365)

zabs_max = 10000
nit_max = 1000

xmin, xmax = -1.5, 1.5
xwidth = xmax - xmin
ymin, ymax = -1.5, 1.5
yheight = ymax - ymin


julia = np.zeros((im_width, im_height))
for ix in range(im_width):
    for iy in range(im_height):
        nit = 0

        z = complex(ix / im_width * xwidth + xmin,
                    iy / im_height * yheight + ymin)

        while abs(z) <= zabs_max and nit < nit_max:
            z = z**2 + c
            nit += 1
        ratio = nit / nit_max
        julia[ix,iy] = ratio

fig, ax = plt.subplots()
ax.imshow(julia, interpolation='nearest', cmap=cm.hot)


xtick_labels = np.linspace(xmin, xmax, xwidth / 0.5)
ax.set_xticks([(x-xmin) / xwidth * im_width for x in xtick_labels])
ax.set_xticklabels(['{:.1f}'.format(xtick) for xtick in xtick_labels])
ytick_labels = np.linspace(ymin, ymax, yheight / 0.5)
ax.set_yticks([(y-ymin) / yheight * im_height for y in ytick_labels])
ax.set_yticklabels(['{:.1f}'.format(ytick) for ytick in ytick_labels])


plt.show()

Para esse conjunto de pontos específicos escolhidos, o resultado será o seguinte:

A imagem será apresentada aqui.

comentou Jun 23, 2018 por MarcioGama (106 pontos)  
Vanderson, ficou bem legal a implementação. Demorei um pouco para entender o código todo. Mas, a proposta do exercício era fazer a implementação usando orientação a objetos. Vou postar uma adaptação do seu exercício para oop.
0 votos
respondida Out 24 por Stuart Mill (1,159 pontos)  

Apenas uma adição, ao definir a classe JuliaSets, coloquei como ficaria o cálculo dos Julia Sets conforme o algoritmo descrito no capítulo 8 do livro do Flake (The Computational Beauty of Nature) em que, em vez de se usar essa variávei ratio associada a cada elemento da matriz original, APENAS colorimos os valores que sobrevivem completamente à recursão. Isto é, apenas se um \( z \) passa pelo número máximo definido de iterações, associamos 1 a esse ponto no plano complexo e, caso contrário, associamos 0. A modificação no código seria (a modificação está justamente na definição do método calc_simple() ):

# Definindo a classe com o atributo de números complexos
class JuliaSet(object):

    def __init__(self, zmax, nmax, real, imaginary):
        self.c = complex(real, imaginary); # Método para acessar o número
        self.zmax = zmax; # Método para acessar o valor absoluto máximo de z (complexo)
        self.nmax = nmax; # Método para acessar o número máximo de iterações
        self.julia = np.zeros((im_width, im_height)); # Começa com uma matriz nula de 
                                    #im_width linhas e im_height colunas para a imagem.


    def calc(self): # Método para calcular a dinâmica da recursão de Mandelbrot.
                    # Quanto mais o ponto "sobrevive" (não explode), mais clara sua cor.

        for ix in range(im_width):
            for iy in range(im_height):
                nit = 0
                # O domínio são os valores inicias da recursão. Queremos inserir
                # uma série de z iniciais na recursão, mantendo a constante c (complexa) fixa.
                z = complex(ix / im_width * xwidth + xmin,iy / im_height * yheight + ymin);

                while abs(z) <= self.zmax and nit < self.nmax: #While com os parâmetros de parada definidos anteriormente
                    z = z**2 + self.c #Recursão de Mandelbrot. Veja que c está parado, estamos mapeando vários z's.
                    nit += 1 # Aumentando a profundidade da recursão.
                ratio = nit / self.nmax # Quanto mais próximo do nmax, maior o ratio e mais clara será a cor mapeada (isto é, mais o ponto "sobrevive" à recursão)
                self.julia[ix,iy] = ratio

    def calc_simple(self): #Versão mais simples do cálculo (presente no livro)
         for ix in range(im_width):
            for iy in range(im_height):
                nit = 0
                # O domínio são os valores inicias da recursão. Queremos inserir
                # uma série de z iniciais na recursão, mantendo a constante c (complexa) fixa.
                z = complex(ix / im_width * xwidth + xmin,iy / im_height * yheight + ymin);

                while abs(z) <= self.zmax and nit < self.nmax: #While com os parâmetros de parada definidos anteriormente
                    z = z**2 + self.c #Recursão de Mandelbrot. Veja que c está parado, estamos mapeando vários z's.
                    nit += 1 # Aumentando a profundidade da recursão.            
                    if nit == self.nmax:
                        self.julia[ix,iy] = 1 # Colorimos SOMENTE se a recursão continua até nosso limite mínimo. Usar 50, por exemplo.
                    else:
                       self.julia[ix,iy] = 0
...