2

Eu tentei adaptar o seguinte código em fortran para python:

PROGRAM Escape_tent
IMPLICIT NONE
! Declare local variables
INTEGER :: i,j,tfinal,P,k,rede
Real,dimension(0:2800) :: x,xn
real, external :: f
Real :: sigma,pi,soma,r
pi=4.*atan(1.0)
!
! Parametros e CIs
rede=100
do j=0,75
x(j)=0.3+0.1*sin((j*2*pi)/rede)
enddo
do j=50,rede
x(j)=1.1+0.1*sin((j*2*pi)/rede+pi)
enddo
tfinal=100
sigma=0.38
r=0.32
P=int(rede*r)

!
! Mapa logistico
!
do i=1,tfinal
!atualiza os vetores com o mapa
do j=0,rede
x(j)=f(x(j))
enddo
!
!
! Finite Range Coupling
!
do j=0,rede
soma=0.0
do k=-P,P
soma=soma+x(modulo(j+K,rede+1))-x(j)
enddo
xn(j)=x(j)+(sigma/(2.*P))*soma
enddo
do j = 0,rede
x(j)=xn(j)
enddo
! Write out values of the map
write(1,100) (x(k),k=0,rede)
enddo
END PROGRAM Escape_tent

real function f(x)
    implicit none
    real, intent(in) :: x
    real, parameter :: r=2.78   
f=r*x(1 - x)
end function f

Quando eu adaptei para o python ficou da seguinte maneira:

import numpy as np
from time import time

# Variáveis do modelo
x, xn = [], []
r = 0.32  # 'raio' de acoplamento r = P/N
pi = np.pi
soma = 0.0

# Parâmetros para condições iniciais
rede = 100  # tamanho da rede
P = int(rede * r)
sigma = 0.38
t = np.arange(0, 1000)  # tempo de iteração

#Condições iniciais da rede
for j in range(0, int(rede / 2)):
    # determina as condições iniciais de metade da rede
    x.insert(j, 0.3 + 0.1 * np.sin((j * 2 * pi) /     (rede)))

for j in range(int(rede / 2), rede):
    # determina as condições iniciais da outra metade da rede
    x.insert(j, 1.1 + 0.1 * np.sin((j * 2 * pi) / (rede + pi)))

# função usada na iteração
def f(x):
    return 3.8 * x * (1 - x)

t1 = time()
#iterando a rede
for i in t:
    for j in range(0, rede):
        x[j] = f(x[j])  # aplica a função em todos os     elementos da rede

    for j in range(0, rede):  # para cada elemento j da rede
        for k in range(0, rede):  # temos um elemento k que esta a uma distância do j
            soma = 0.0
            if k != j:  # k não interage com k
                # calcula a distância do k até o j
                soma = x[k] - x[j]
                xn.insert(j, x[j] + (sigma / 2 * P) * soma)  # add no vetor x

    for j in range(0, rede):
        x[j] = xn[j] % (1.0) #Aqui eu limito o valor de Xn[j] para não tender ao infinito.
    print("T {}\n X :{}".format(i,x))
t2 = time()
print("Tempo decorrido: {}".format(t2 - t1))

Depois da Nona iteração ocorre a exception "RuntimeWarning: overflow encountered in double_scalars" e "RuntimeWarning: invalid value encountered in double_scalars" no terminal do Pycharm e os dados adicionados ao vetor x são 'nan'. Gostaria de saber o que esta acontecendo, pois no fortran roda normalmente ja no python tenho esse problema. Obrigado.

1 Resposta 1

3

Not a number (nan)

Seu código está gerando nan na lista que vc nomeia de x (sugiro que assista esse video), e o que é um nan? De forma bem objetiva e suscinta, nan é o acrônimo para "Not a number", ou seja, qualquer operação matemática que resulte em um valor que não possa ser expressado em um número, a representação desse resultado será nan, geralmente uma exceção é gerada quando se tenta executar uma operação que resulte em um nan, por exempo se vc tentar fazer uma divisão por zero, a exceção ZeroDivisionError será levantada.

a = 5
b = 0
c = a/b
ZeroDivisionError: division by zero

Mas nem sempre a exceção é levantada, a menos que vc trate isso, é dificil criar, propositadamente, um nan, mas existe uma operação com inf que é clássica para um exemplo:

a = float('Inf')
b = float('Inf')
c = a - b
print(c)
nan

O que acontece é que seu codigo em python (não vou comparar com o codigo em fortran, pq não conheço a linguagem), está gerando valores nan na lista x, então quando voce envia o valor para a função f() o valor nan é retornado, pois qualqer operação que envolva um nan resultará um nan, a mensagem RuntimeWarning: overflow encountered in double_scalars na maioria das vezes está relacionada com operações que envolvem valores nan, descubra pq esses valores estão sendo gerados e vc resolverá o problema, vou fazer algumas sugestões para voce debugar, mas voce pode acrescentar mais procedimentos de acordo com seu contexto:

1) Altere a função f() para interceptar o valor nan em x, da seguinte forma:

# função usada na iteração
def f(x):
    result =  3.8 * x * (1 - x)
    if np.isnan(result):
        line = '='*22
        print(line,'x recebido na função f: ', x,line, sep='\n')
    return result

2) Comente a saida de x no loop de iteração, apenas para que a saida do debug fique mais limpa, além disso antes do procedimento que chama a função f() verifique se o elemento j em x que vc esta mandando para a função f não é nan, se for apresente as variaveis envolvidas, o loop de iteração deverá ficar assim:

t1 = time()
#iterando a rede
for i in t:
    for j in range(0, rede):

        if np.isnan(x[j]):
            print('valor de j na iteracao, quando x[j] é nan', j)
            print('valor de x na iteracao, quando x[j] é nan', x)

        x[j] = f(x[j])  # aplica a função em todos os elementos da rede

    for j in range(0, rede):  # para cada elemento j da rede
        for k in range(0, rede):  # temos um elemento k que esta a uma distância do j
            soma = 0.0
            if k != j:  # k não interage com k
                # calcula a distância do k até o j
                soma = x[k] - x[j]
                xn.insert(j, x[j] + (sigma / 2 * P) * soma)  # add no vetor x

    for j in range(0, rede):
        x[j] = xn[j]
    # print("T {}\n X :{}".format(i,x))
t2 = time()
print("Tempo decorrido: {}".format(t2 - t1))

Após essa alterações o resultado deverá ser algo próximo de:

valor de j na iteracao, quando x[j] é nan 0
valor de x na iteracao, quando x[j] é nan [nan, nan, nan,... 1.0350165585196508e+289,... ]
======================
x recebido na função f: 
nan
======================
...

A saida acima é apenas um fragmento, rode o codigo com as alterações para ver a saida real, voce pode (e deve) usar print nas variáveis envolvidas nos cálculos para gerar os elementos em x para descobrir o que está gerando os valores nan.

1
  • Excelente resposta, muito obrigado. Eu consegui perceber meu erro depois de um tempo, o fato de eu usar uma função logística e pra alguns casos ela tem um regime caótico logo ela converge ao infinito muito rápido e acho que é por isso que ela gera um NaN. Para corrigir esse erro eu limitei o 'xn[j]' entre 0 e 1. 12/02/2019 às 17:20

Você deve fazer log-in para responder a esta pergunta.

Esta não é a resposta que você está procurando? Pesquise outras perguntas com a tag .