RuntimeWarning: invalid value encountered in double_scalars

Question

При реализации связанной системы с одинаковыми параметрами х, у, z - все прекрасно работает. Но если использовать различные параметры - возникает следующая ошибка для переменных х и z. С чем это связано? Понятно что происходит переполнение после 10000 итераций, но не совсем понятно как это исправить. Возможно, есть како-то способ ограничить, чтобы дальше программа не считала и построила по первым 10000 значениям? Или решить как-то иначе? Искал подобные ошибки в интернете, но аналогичной не нашел в принципе. Сама ошибка:

z1_n1=z1_n+h*(b+z1_n*(x1_n-c))
RuntimeWarning: overflow encountered in double_scalars
z2_n1=z2_n+h*(b+z2_n*(x2_n-c))
RuntimeWarning: overflow encountered in double_scalars

Код программы:

from numpy import *
import numpy as np
from matplotlib import *
from scipy import *
from pylab import figure, show, setp
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
import scipy.io as sio
from scipy.signal import find_peaks
 
#We define a function which is going to be the recursive function.
def num_rossler(x1_n, y1_n, z1_n, x2_n, y2_n, z2_n, h, a, b, cz, k):
    x1_n1=x1_n+h*(-y1_n-z1_n+k*(x2_n-x1_n)) # dx/dt = x1 + h*( - y1 - z1 + k (x2 - x1 ))
    y1_n1=y1_n+h*(x1_n+a*y1_n)
    z1_n1=z1_n+h*(b+z1_n*(x1_n-c))   
 
    x2_n1=x2_n+h*(-y2_n-z2_n+k*(x1_n-x2_n))
    y2_n1=y2_n+h*(x2_n+a*y2_n)
    z2_n1=z2_n+h*(b+z2_n*(x2_n-c))
 
 
    return x1_n1, y1_n1, z1_n1, x2_n1, y2_n1, z2_n1
 
#Now we prepare some variables
#First the parameters
a=0.165
b=0.2
c=10
k=0.1
 
#Them the time interval and the step size
t_ini=0
t_fin=32*pi
h=0.0001
numsteps=int((t_fin-t_ini)/h)
 
#using this parameters we build the time.
t=numpy.linspace(t_ini,t_fin,numsteps)
#And the vectors for the solutions
x1=numpy.zeros(numsteps)
y1=numpy.zeros(numsteps)
z1=numpy.zeros(numsteps)
 
x2=numpy.zeros(numsteps)
y2=numpy.zeros(numsteps)
z2=numpy.zeros(numsteps)
 
 
#We set the initial conditions
x1[0]=0.001
y1[0]=0.001
z1[0]=0.001
 
x2[0]=0.002
y2[0]=0.002
z2[0]=0.002
 
#This is the main loop where we use the recursive system to obtain the solution
for k in range(x1.size-1):
    #We use the previous point to generate the new point using the recursion
    [x1[k+1],y1[k+1],z1[k+1],x2[k+1],y2[k+1],z2[k+1]]=num_rossler(x1[k],y1[k],z1[k],x2[k],y2[k],z2[k],t[k+1]-t[k],a,b,c,k)
 
#Now that we have the solution in vectors t,x,y,z is time to plot them.
 
#We create a figure and 4 axes on it. 3 of the axes are going to be 2D and the fourth one is a 3D plot.
fig = figure()
ax1 = fig.add_axes([0.30, 0.7, 0.4, 0.2])
ax3 = fig.add_axes([0.05, 0.1, 0.4, 0.5],projection='3d')
ax4 = fig.add_axes([0.55, 0.1, 0.4, 0.5],projection='3d')
 
 
#And we add vectors to each plot
ax1.plot(t, x1,color='red',lw=1,label='x1(t)')
ax1.set_xlabel('t')
ax1.set_ylabel('x(t)')
ax1.legend()
ax1.axis((t_ini,t_fin,min(x1),max(x1)))
 
 
ax1.plot(t, x2,color='green',lw=1,label='x2(t)')
ax1.set_xlabel('t')
ax1.set_ylabel('x(t)')
ax1.legend()
ax1.axis((t_ini,t_fin,min(x2),max(x2)))
 
ax4.plot(x2, y2,z2,color='green',lw=1)
ax4.set_xlabel('x2(t)')
ax4.set_ylabel('y2(t)')
ax4.set_zlabel('z2(t)')
 
ax3.plot(x1, y1,z1,color='red',lw=1)
ax3.set_xlabel('x1(t)')
ax3.set_ylabel('y1(t)')
ax3.set_zlabel('z1(t)')
 
show()

Answer 1

У вас ошибка вот в этой строчке:

[x1[k+1],y1[k+1],z1[k+1],x2[k+1],y2[k+1],z2[k+1]]=num_rossler(x1[k],y1[k],z1[k],x2[k],y2[k],z2[k],t[k+1]-t[k],a,b,c,k)

Вы обозначили константу связи между двумя аттракторами буквой k, и переменную итерации тоже буквой k. Поэтому в вашей системе x1 и x2 начинают бешено нарастать, ведь k растёт линейно по времени, и довольно быстро - на интервале от 0 до 30pi она вырастает до миллиона.

Обозначьте связь другой буквой, хоть K заглавной, и всё получится.

UPDATE

Если позволите, пара рекомендаций как ускорить счёт на Python.

Поставьте numba. Это JIT компилятор для Python, который умеет транслировать небольшое подмножество Python в машинный код. Ваша счётная функция оперирует только вещественной арифметикой, поэтому она после обработки numba будет считаться быстрее. Я сделал замеры - ускорение более, чем в два раза.

@numba.njit
def num_rossler(x1_n, y1_n, z1_n, x2_n, y2_n, z2_n, h, a, b, c, k):
    # Тело функции без изменений

Замените вызов функции на генератор, который строит траекторию. Фактически, этот итератор ведёт интеграцию методом Эйлера, выдавая каждый раз следующую точку.

@numba.njit
def gen_rossler(x1, y1, z1, x2,y2,z2, dt, a,b,c,k, numsteps):
    # Выдать начальную точку
    yield (x1,y1,z1, x2,y2,z2)
    for _ in range(numsteps):
        dx1 = dt*(-y1-z1 + k*(x2-x1))
        dy1 = dt*(x1 + a*y1)
        dz1 = dt*(b + z1*(x1 - c))
    
        dx2 = dt*(-y2-z2 + k*(x1-x2))
        dy2 = dt*(x2 + a*y2)
        dz2 = dt*(b + z2*(x2 - c))
    
        x1,y1,z1 = x1+dx1, y1+dy1, z1+dz1
        x2,y2,z2 = x2+dx2, y2+dy2, z2+dz2
        
        yield (x1,y1,z1, x2,y2,z2)

Перебор итератора вот таким циклом:

for i, (_x1,_y1,_z1,_x2,_y2,_z2) in zip(range(numsteps), gen_rossler(x1[0],y1[0],z1[0],x2[0],y2[0],z2[0],t[1]-t[0],a,b,c,K, numsteps)):
    #We use the previous point to generate the new point using the recursion
    x1[i],y1[i],z1[i],x2[i],y2[i],z2[i]=_x1,_y1,_z1,_x2,_y2,_z2

Ускорение по сравнению с первоначальным вариантом в девять раз.

UPDATE 2

По предложению @CrazyElf засунул основной цикл счёта в отдельную функцию и поместил её под декоратор numba.njit. Итоговое ускорение 70 раз.

@numba.njit(inline='always')
def step_rossler(x1,y1,z1,x2,y2,z2, dt, a, b, c, k):
    dx1 = dt*(-y1-z1 + k*(x2-x1))
    dy1 = dt*(x1 + a*y1)
    dz1 = dt*(b + z1*(x1 - c))

    dx2 = dt*(-y2-z2 + k*(x1-x2))
    dy2 = dt*(x2 + a*y2)
    dz2 = dt*(b + z2*(x2 - c))
    
    x1,y1,z1 = x1+dx1, y1+dy1, z1+dz1
    x2,y2,z2 = x2+dx2, y2+dy2, z2+dz2
    
    return (x1,y1,z1, x2,y2,z2)

@numba.njit
def orbit_rossler(x1,y1,z1,x2,y2,z2, dt, a, b, c, k, numsteps):
    XYZ = np.zeros((numsteps+1)*6)
    XYZ = np.reshape(XYZ, (numsteps+1, 6))
    XYZ[0,:] = x1,y1,z1,x2,y2,z2
    for i in range(1, numsteps+1):
        x1,y1,z1,x2,y2,z2 = step_rossler(x1,y1,z1,x2,y2,z2, dt, a,b,c,k)

        XYZ[i,:] = x1,y1,z1,x2,y2,z2
    return XYZ

Функция orbit_rossler возвращает массив из numsteps+1 шестёрок чисел. Этот массив разбирается на массивы для x1, y1 и остальных переменных следующим образом:

XYZ = orbit_rossler(x1_0, y1_0, z1_0, x2_0, y2_0, z2_0,t[1]-t[0],a,b,c,K, numsteps)

x1 = XYZ[:,0]
y1 = XYZ[:,1]
z1 = XYZ[:,2]
# Или все сразу в цикле
x1,y1,z1,x2,y2,z2 = [ XYZ[:,i] for i in range(6) ]

Неоптимизированная версия считала 4.5 секунд. Оптимизированная - 66 миллисекунд. Ускорение 70 раз.

Jupyter Notebook с реализацией можно скачать из Гитхаба

БЛОГ НА HUSL

RuntimeWarning: invalid value encountered in double_scalars

Ответы (1 шт):