В Python проблема с лишней работой для SciPY ode.int - PullRequest
Новая
       6

В Python проблема с лишней работой для SciPY ode.int

0 голосов
/ 22 ноября 2018

Я пытаюсь смоделировать проблему преследования ошибок, преследующих друг друга в двухмерной плоскости, и я использую SciPY.odeint, чтобы помочь мне в этом.С помощью следующего кода модель работает, однако, когда ошибки сближаются, модель ломается и выдает излишнюю работу, выполненную с этим вызовом (возможно, неправильный тип Dfun). 

import numpy as np
from scipy.integrate import odeint

def split_list(a_list):
    half = len(a_list)//2
    return a_list[:half], a_list[half:]

def diff(w, t):
    x_points, y_points = split_list(w)
    def abso(f, s):
        return np.sqrt((x_points[f] - x_points[s])**2 + (y_points[f] - y_points[s])**2)
    x_funct = [(x_points[i+1] - x_points[i]) / abso(i+1, i) for i in range(len(x_points) - 1)]
    x_funct.append((x_points[0] - x_points[-1]) / abso(0,-1))

    y_funct = [(y_points[i+1] - y_points[i]) / abso(i+1,i) for i in range(len(x_points) - 1)]
    y_funct.append((y_points[0] - y_points[-1]) / abso(0,-1))
    funct = x_funct + y_funct
    return funct

def ode(tstart, tend, init_cond):

    t = np.linspace(tstart, tend, step_size)

    wsol = odeint(diff, init_cond, t)
    sols = [wsol[:,i] for i in range(len(init_cond))]
    x_sols, y_sols = split_list(sols)
    return x_sols, y_sols, t, tend

bug_init_cond = [[-0.5, 1],
                 [0.5, 1],
                 [0.5, -1],
                 [-0.5, -1],]

amount_of_bugs = 4
step_size = 10000

x_sols, y_sols, t, tend = ode(0, 5, [bug_init_cond[i][j] for j in range(2) for i in range(amount_of_bugs)])

Как я вполнеМне было интересно, есть ли решение для этой лишней работы, если вы недавно начали использовать функцию Scipy.odeint?Спасибо за ваше время.

1 Ответ

0 голосов
/ 23 ноября 2018

Ваша проблема в том, что в точном решении пути встречаются во время от t=1.48 до t=1.5.В точном решении вы получите ошибку деления на ноль с шумом с плавающей запятой, который «ухудшается» до жесткой ситуации, когда размер шага регулируется вниз до тех пор, пока выходной временной шаг не потребует более mxstep=500 внутренних шагов.

Вы можете добавить условия, чтобы правая сторона была установлена ​​на ноль, если позиции закрываются.Один быстрый способ добиться этого - изменить функцию расстояния от abso до 

  def abso(f, s):
    return np.sqrt(1e-12+(x_points[f] - x_points[s])**2 + (y_points[f] - y_points[s])**2)

, чтобы вы никогда не делили на ноль, а для видимых расстояний возмущение пренебрежимо мало.

...