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.
python scipy ode
asked Nov 22 '18 at 15:48
Oscar LingardOscar Lingard
226
add a comment |
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.
python scipy ode
asked Nov 22 '18 at 15:48
Oscar LingardOscar Lingard
226
add a comment |
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.
python scipy ode
asked Nov 22 '18 at 15:48
Oscar LingardOscar Lingard
226
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.
python scipy ode
python scipy ode
asked Nov 22 '18 at 15:48
Oscar LingardOscar Lingard
226
asked Nov 22 '18 at 15:48
Oscar LingardOscar Lingard
226
asked Nov 22 '18 at 15:48
Oscar LingardOscar Lingard
226
asked Nov 22 '18 at 15:48
Oscar LingardOscar Lingard
226
asked Nov 22 '18 at 15:48
Oscar LingardOscar Lingard
226
226
add a comment |
add a comment |
1 Answer
1
active
oldest
votes
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.
answered Nov 23 '18 at 11:33
LutzLLutzL
13.9k21426
add a comment |
Your Answer
StackExchange.ifUsing("editor", function () {
StackExchange.using("externalEditor", function () {
StackExchange.using("snippets", function () {
StackExchange.snippets.init();
});
});
}, "code-snippets");
StackExchange.ready(function() {
var channelOptions = {
tags: "".split(" "),
id: "1"
};
initTagRenderer("".split(" "), "".split(" "), channelOptions);
StackExchange.using("externalEditor", function() {
// Have to fire editor after snippets, if snippets enabled
if (StackExchange.settings.snippets.snippetsEnabled) {
StackExchange.using("snippets", function() {
createEditor();
});
}
else {
createEditor();
}
});
function createEditor() {
StackExchange.prepareEditor({
heartbeatType: 'answer',
autoActivateHeartbeat: false,
convertImagesToLinks: true,
noModals: true,
showLowRepImageUploadWarning: true,
reputationToPostImages: 10,
bindNavPrevention: true,
postfix: "",
imageUploader: {
brandingHtml: "Powered by u003ca class="icon-imgur-white" href="https://imgur.com/"u003eu003c/au003e",
contentPolicyHtml: "User contributions licensed under u003ca href="https://creativecommons.org/licenses/by-sa/3.0/"u003ecc by-sa 3.0 with attribution requiredu003c/au003e u003ca href="https://stackoverflow.com/legal/content-policy"u003e(content policy)u003c/au003e",
allowUrls: true
},
onDemand: true,
discardSelector: ".discard-answer"
,immediatelyShowMarkdownHelp:true
});
}
});
Sign up or log in
StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Post as a guest
Required, but never shown
StackExchange.ready(
function () {
StackExchange.openid.initPostLogin('.new-post-login', 'https%3a%2f%2fstackoverflow.com%2fquestions%2f53434453%2fin-python-excess-work-done-problem-for-scipy-ode-int%23new-answer', 'question_page');
}
);
Post as a guest
Required, but never shown
1 Answer
1
active
oldest
votes
1 Answer
1
active
oldest
votes
active
oldest
votes
active
oldest
votes
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.
answered Nov 23 '18 at 11:33
LutzLLutzL
13.9k21426
add a comment |
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.
answered Nov 23 '18 at 11:33
LutzLLutzL
13.9k21426
add a comment |
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.
answered Nov 23 '18 at 11:33
LutzLLutzL
13.9k21426
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.
answered Nov 23 '18 at 11:33
LutzLLutzL
13.9k21426
answered Nov 23 '18 at 11:33
LutzLLutzL
13.9k21426
answered Nov 23 '18 at 11:33
LutzLLutzL
13.9k21426
answered Nov 23 '18 at 11:33
LutzLLutzL
13.9k21426
13.9k21426
add a comment |
add a comment |
Thanks for contributing an answer to Stack Overflow!
- Please be sure to answer the question. Provide details and share your research!
But avoid …
- Asking for help, clarification, or responding to other answers.
- Making statements based on opinion; back them up with references or personal experience.
To learn more, see our tips on writing great answers.
Sign up or log in
StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Post as a guest
Required, but never shown
StackExchange.ready(
function () {
StackExchange.openid.initPostLogin('.new-post-login', 'https%3a%2f%2fstackoverflow.com%2fquestions%2f53434453%2fin-python-excess-work-done-problem-for-scipy-ode-int%23new-answer', 'question_page');
}
);
Post as a guest
Required, but never shown
Sign up or log in
StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Post as a guest
Required, but never shown
Sign up or log in
StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Post as a guest
Required, but never shown
Sign up or log in
StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Post as a guest
Required, but never shown
Required, but never shown
Required, but never shown
Required, but never shown
Required, but never shown
Required, but never shown
Required, but never shown
Required, but never shown
Required, but never shown
