CODEWITHSUNDEEP
pythonscipyode
Oscar Lingard
Oscar Lingard

Reputation: 49

In Python, excess work done problem for SciPY ode.int

I'm trying to model a pursuit problem of bugs chasing one another in a 2 dimensional plane and I'm using SciPY.odeint to help me with this. With the following code, the model works however as the bugs get closer together the model breaks down and emits the excess work done on this call (perhaps wrong Dfun type) error.

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

As I'm quite new with using the Scipy.odeint function, I was wondering if there is a solution to this excess work done? Thank-you for your time.

Upvotes: 1

Views: 182

Answers (1)

Lutz Lehmann
Lutz Lehmann

Reputation: 26040

Your problem is that in the exact solution the paths meet at time t=1.48 to t=1.5. In an exact solution, you would get a division by zero error, with floating point noise that is "degraded" to a stiff situation where the step size is regulated down until the output time step requires more than mxstep=500 internal steps.

You could add conditions so that the right side is set to zero if the positions get to close. One quick hack to achieve that is to modify the distance function abso to

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

so that you never divide by zero, and for visible distances the perturbation is negligible.

Upvotes: 2

Related Questions