Fortran code adaptation for python [RuntimeWarning: overflow encountered in double_scalars]

0

I tried to adapt the following code in fortran for 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

When I adapted to python it was as follows:

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]
    print("T {}\n X :{}".format(i,x))
t2 = time()
print("Tempo decorrido: {}".format(t2 - t1))

After the ninth iteration occurs the "RuntimeWarning: overflow encountered in double_scalars" and "RuntimeWarning: invalid value encountered in double_scalars" exception in the Pycharm terminal, and the data added to the vector x is 'nan'. I would like to know what is happening, because in the fortran wheel normally already in python I have this problem. Thank you.

    
asked by anonymous 30.12.2018 / 20:08

1 answer

0

Not a number (nan)

Your code is generating nan in the list that you name from x (I suggest you watch this video ), and what is a nan ? In a very objective and succinct way, nan is the acronym for "Not a number", that is, any mathematical operation that results in a value that can not be expressed in a number, the representation of that result will be nan , usually an exception is generated when trying to execute an operation that results in a nan , for example if you try to divide by zero, the ZeroDivisionError exception will be raised.

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

But the exception is not always raised, unless you deal with this, it is difficult to purposely create a nan , but there is an operation with inf that is classic for an example:

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

What happens is that your code in python (I will not compare with the code in fortran, because I do not know the language), is generating nan values in the x list, so when you send the value to the function f() the nan value is returned, since any operation that involves a nan will result in a nan , the RuntimeWarning: overflow encountered in double_scalars message is most often related to operations involving nan values, find out these values are being generated and you will solve the problem, I will make some suggestions for you to debug, but you can add more procedures according to your context:

1) Change the f() function to intercept the nan value in x , as follows:

# 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) Comment the output of x in the iteration loop, just so that the output of the debug is cleaner, in addition before the procedure that calls the function f() , check that the element j in x that you are sending to the function f is not nan , if you present the variables involved, the iteration loop should look like this:

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

After these changes the result should be something close to:

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
======================
...

The output above is just a fragment, run the code with the changes to see the actual output, you can (and should) print in the variables involved in the calculations to generate the elements in x to find out what is generating nan values.

    
31.12.2018 / 15:00