scipy dblquad providing the wrong result in simple double integral
I am trying to calculate a straightforward doble definite integral in Python: function Max(0, (4-12x) + (6-12y)) in the square [0,1] x [0,1].
We can do it with Mathematica and get the exact result:
Integrate[Max[0, 4-12*u1 + 6-12*u2], {u1, 0, 1}, {u2, 0,1}] = 125/108.
With a simple Monte Carlo simulation I can confirm this result. However, using scipy.integrate.dblquad I am getting a value of 0.0005772072907971, with error 0.0000000000031299
from scipy.integrate import dblquad
def integ(u1, u2):
return max(0, (4 - 12*u1) + (6 - 12*u2))
sol_int, err = dblquad(integ, 0, 1, lambda _:0, lambda _:1, epsabs=1E-12, epsrel=1E-12)
print("dblquad: %0.16f. Error: %0.16f" % (sol_int, err) )
Agreed that the function is not derivable, but it is continuous, I see no reason for this particular integral to be problematic.
I thought maybe dblquad has an 'options' argument where I can try different numerical methods, but I found nothing like that.
So, what am I doing wrong?
scipy integration
asked Nov 12 at 14:59
zeycus
186113
|
show 1 more comment
I am trying to calculate a straightforward doble definite integral in Python: function Max(0, (4-12x) + (6-12y)) in the square [0,1] x [0,1].
We can do it with Mathematica and get the exact result:
Integrate[Max[0, 4-12*u1 + 6-12*u2], {u1, 0, 1}, {u2, 0,1}] = 125/108.
With a simple Monte Carlo simulation I can confirm this result. However, using scipy.integrate.dblquad I am getting a value of 0.0005772072907971, with error 0.0000000000031299
from scipy.integrate import dblquad
def integ(u1, u2):
return max(0, (4 - 12*u1) + (6 - 12*u2))
sol_int, err = dblquad(integ, 0, 1, lambda _:0, lambda _:1, epsabs=1E-12, epsrel=1E-12)
print("dblquad: %0.16f. Error: %0.16f" % (sol_int, err) )
Agreed that the function is not derivable, but it is continuous, I see no reason for this particular integral to be problematic.
I thought maybe dblquad has an 'options' argument where I can try different numerical methods, but I found nothing like that.
So, what am I doing wrong?
scipy integration
asked Nov 12 at 14:59
zeycus
186113
1
When I run this using scipy 1.1.0 and numpy 1.15.4, the output isdblquad: 1.1574073869306833. Error: 0.0000000000031299
– Warren Weckesser
Nov 12 at 15:16
@Warren thank you, very interesting, that would be correct. I see I am using numpy 1.14.3, scipy 1.1.0. Python 3.6.5
– zeycus
Nov 12 at 16:00
Are you sure the code that you show in the question is exactly the same as the code that gave you0.0005772072907971?
– Warren Weckesser
Nov 12 at 16:01
Yes, I just recopied from this screen to a new file, and executed it. Same 0.000577 result.
– zeycus
Nov 12 at 16:36
2
Are you on Windows? How did you install scipy? You might be hitting this bug: github.com/scipy/scipy/issues/6882; see also stackoverflow.com/questions/27270044/…
– Warren Weckesser
Nov 12 at 17:30
|
show 1 more comment
I am trying to calculate a straightforward doble definite integral in Python: function Max(0, (4-12x) + (6-12y)) in the square [0,1] x [0,1].
We can do it with Mathematica and get the exact result:
Integrate[Max[0, 4-12*u1 + 6-12*u2], {u1, 0, 1}, {u2, 0,1}] = 125/108.
With a simple Monte Carlo simulation I can confirm this result. However, using scipy.integrate.dblquad I am getting a value of 0.0005772072907971, with error 0.0000000000031299
from scipy.integrate import dblquad
def integ(u1, u2):
return max(0, (4 - 12*u1) + (6 - 12*u2))
sol_int, err = dblquad(integ, 0, 1, lambda _:0, lambda _:1, epsabs=1E-12, epsrel=1E-12)
print("dblquad: %0.16f. Error: %0.16f" % (sol_int, err) )
Agreed that the function is not derivable, but it is continuous, I see no reason for this particular integral to be problematic.
I thought maybe dblquad has an 'options' argument where I can try different numerical methods, but I found nothing like that.
So, what am I doing wrong?
scipy integration
asked Nov 12 at 14:59
zeycus
186113
I am trying to calculate a straightforward doble definite integral in Python: function Max(0, (4-12x) + (6-12y)) in the square [0,1] x [0,1].
We can do it with Mathematica and get the exact result:
Integrate[Max[0, 4-12*u1 + 6-12*u2], {u1, 0, 1}, {u2, 0,1}] = 125/108.
With a simple Monte Carlo simulation I can confirm this result. However, using scipy.integrate.dblquad I am getting a value of 0.0005772072907971, with error 0.0000000000031299
from scipy.integrate import dblquad
def integ(u1, u2):
return max(0, (4 - 12*u1) + (6 - 12*u2))
sol_int, err = dblquad(integ, 0, 1, lambda _:0, lambda _:1, epsabs=1E-12, epsrel=1E-12)
print("dblquad: %0.16f. Error: %0.16f" % (sol_int, err) )
Agreed that the function is not derivable, but it is continuous, I see no reason for this particular integral to be problematic.
I thought maybe dblquad has an 'options' argument where I can try different numerical methods, but I found nothing like that.
So, what am I doing wrong?
scipy integration
scipy integration
asked Nov 12 at 14:59
zeycus
186113
asked Nov 12 at 14:59
zeycus
186113
asked Nov 12 at 14:59
zeycus
186113
asked Nov 12 at 14:59
zeycus
186113
asked Nov 12 at 14:59
zeycus
186113
186113
1
When I run this using scipy 1.1.0 and numpy 1.15.4, the output isdblquad: 1.1574073869306833. Error: 0.0000000000031299
– Warren Weckesser
Nov 12 at 15:16
@Warren thank you, very interesting, that would be correct. I see I am using numpy 1.14.3, scipy 1.1.0. Python 3.6.5
– zeycus
Nov 12 at 16:00
Are you sure the code that you show in the question is exactly the same as the code that gave you0.0005772072907971?
– Warren Weckesser
Nov 12 at 16:01
Yes, I just recopied from this screen to a new file, and executed it. Same 0.000577 result.
– zeycus
Nov 12 at 16:36
2
Are you on Windows? How did you install scipy? You might be hitting this bug: github.com/scipy/scipy/issues/6882; see also stackoverflow.com/questions/27270044/…
– Warren Weckesser
Nov 12 at 17:30
|
show 1 more comment
1
When I run this using scipy 1.1.0 and numpy 1.15.4, the output isdblquad: 1.1574073869306833. Error: 0.0000000000031299
– Warren Weckesser
Nov 12 at 15:16
@Warren thank you, very interesting, that would be correct. I see I am using numpy 1.14.3, scipy 1.1.0. Python 3.6.5
– zeycus
Nov 12 at 16:00
Are you sure the code that you show in the question is exactly the same as the code that gave you0.0005772072907971?
– Warren Weckesser
Nov 12 at 16:01
Yes, I just recopied from this screen to a new file, and executed it. Same 0.000577 result.
– zeycus
Nov 12 at 16:36
2
Are you on Windows? How did you install scipy? You might be hitting this bug: github.com/scipy/scipy/issues/6882; see also stackoverflow.com/questions/27270044/…
– Warren Weckesser
Nov 12 at 17:30
1
1
When I run this using scipy 1.1.0 and numpy 1.15.4, the output is
dblquad: 1.1574073869306833. Error: 0.0000000000031299– Warren Weckesser
Nov 12 at 15:16
When I run this using scipy 1.1.0 and numpy 1.15.4, the output is
dblquad: 1.1574073869306833. Error: 0.0000000000031299– Warren Weckesser
Nov 12 at 15:16
@Warren thank you, very interesting, that would be correct. I see I am using numpy 1.14.3, scipy 1.1.0. Python 3.6.5
– zeycus
Nov 12 at 16:00
@Warren thank you, very interesting, that would be correct. I see I am using numpy 1.14.3, scipy 1.1.0. Python 3.6.5
– zeycus
Nov 12 at 16:00
Are you sure the code that you show in the question is exactly the same as the code that gave you
0.0005772072907971?– Warren Weckesser
Nov 12 at 16:01
Are you sure the code that you show in the question is exactly the same as the code that gave you
0.0005772072907971?– Warren Weckesser
Nov 12 at 16:01
Yes, I just recopied from this screen to a new file, and executed it. Same 0.000577 result.
– zeycus
Nov 12 at 16:36
Yes, I just recopied from this screen to a new file, and executed it. Same 0.000577 result.
– zeycus
Nov 12 at 16:36
2
2
Are you on Windows? How did you install scipy? You might be hitting this bug: github.com/scipy/scipy/issues/6882; see also stackoverflow.com/questions/27270044/…
– Warren Weckesser
Nov 12 at 17:30
Are you on Windows? How did you install scipy? You might be hitting this bug: github.com/scipy/scipy/issues/6882; see also stackoverflow.com/questions/27270044/…
– Warren Weckesser
Nov 12 at 17:30
|
show 1 more comment
1 Answer
1
active
oldest
votes
try different numerical methods
That's what I would suggest, given the trouble that iterated quad has on Windows. After changing it to an explicit two-step process, you can replace one of quad with another method, romberg seems the best alternative to me.
from scipy.integrate import quad, romberg
def integ(u1, u2):
return max(0, (4 - 12*u1) + (6 - 12*u2))
sol_int = romberg(lambda u1: quad(lambda u2: integ(u1, u2), 0, 1)[0], 0, 1)
print("romberg-quad: %0.16f " % sol_int)
This prints 1.1574073959987758 on my computer, and hopefully you will get the same.
This does solve it, thanks! I did not know about therombergfunction.
– zeycus
Nov 13 at 8:14
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%2f53264803%2fscipy-dblquad-providing-the-wrong-result-in-simple-double-integral%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
try different numerical methods
That's what I would suggest, given the trouble that iterated quad has on Windows. After changing it to an explicit two-step process, you can replace one of quad with another method, romberg seems the best alternative to me.
from scipy.integrate import quad, romberg
def integ(u1, u2):
return max(0, (4 - 12*u1) + (6 - 12*u2))
sol_int = romberg(lambda u1: quad(lambda u2: integ(u1, u2), 0, 1)[0], 0, 1)
print("romberg-quad: %0.16f " % sol_int)
This prints 1.1574073959987758 on my computer, and hopefully you will get the same.
This does solve it, thanks! I did not know about therombergfunction.
– zeycus
Nov 13 at 8:14
add a comment |
try different numerical methods
That's what I would suggest, given the trouble that iterated quad has on Windows. After changing it to an explicit two-step process, you can replace one of quad with another method, romberg seems the best alternative to me.
from scipy.integrate import quad, romberg
def integ(u1, u2):
return max(0, (4 - 12*u1) + (6 - 12*u2))
sol_int = romberg(lambda u1: quad(lambda u2: integ(u1, u2), 0, 1)[0], 0, 1)
print("romberg-quad: %0.16f " % sol_int)
This prints 1.1574073959987758 on my computer, and hopefully you will get the same.
This does solve it, thanks! I did not know about therombergfunction.
– zeycus
Nov 13 at 8:14
add a comment |
try different numerical methods
That's what I would suggest, given the trouble that iterated quad has on Windows. After changing it to an explicit two-step process, you can replace one of quad with another method, romberg seems the best alternative to me.
from scipy.integrate import quad, romberg
def integ(u1, u2):
return max(0, (4 - 12*u1) + (6 - 12*u2))
sol_int = romberg(lambda u1: quad(lambda u2: integ(u1, u2), 0, 1)[0], 0, 1)
print("romberg-quad: %0.16f " % sol_int)
This prints 1.1574073959987758 on my computer, and hopefully you will get the same.
try different numerical methods
That's what I would suggest, given the trouble that iterated quad has on Windows. After changing it to an explicit two-step process, you can replace one of quad with another method, romberg seems the best alternative to me.
from scipy.integrate import quad, romberg
def integ(u1, u2):
return max(0, (4 - 12*u1) + (6 - 12*u2))
sol_int = romberg(lambda u1: quad(lambda u2: integ(u1, u2), 0, 1)[0], 0, 1)
print("romberg-quad: %0.16f " % sol_int)
This prints 1.1574073959987758 on my computer, and hopefully you will get the same.
answered Nov 13 at 4:43
user6655984
This does solve it, thanks! I did not know about therombergfunction.
– zeycus
Nov 13 at 8:14
add a comment |
This does solve it, thanks! I did not know about therombergfunction.
– zeycus
Nov 13 at 8:14
This does solve it, thanks! I did not know about the
romberg function.– zeycus
Nov 13 at 8:14
This does solve it, thanks! I did not know about the
romberg function.– zeycus
Nov 13 at 8:14
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.
Some of your past answers have not been well-received, and you're in danger of being blocked from answering.
Please pay close attention to the following guidance:
- 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%2f53264803%2fscipy-dblquad-providing-the-wrong-result-in-simple-double-integral%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
1
When I run this using scipy 1.1.0 and numpy 1.15.4, the output is
dblquad: 1.1574073869306833. Error: 0.0000000000031299– Warren Weckesser
Nov 12 at 15:16
@Warren thank you, very interesting, that would be correct. I see I am using numpy 1.14.3, scipy 1.1.0. Python 3.6.5
– zeycus
Nov 12 at 16:00
Are you sure the code that you show in the question is exactly the same as the code that gave you
0.0005772072907971?– Warren Weckesser
Nov 12 at 16:01
Yes, I just recopied from this screen to a new file, and executed it. Same 0.000577 result.
– zeycus
Nov 12 at 16:36
2
Are you on Windows? How did you install scipy? You might be hitting this bug: github.com/scipy/scipy/issues/6882; see also stackoverflow.com/questions/27270044/…
– Warren Weckesser
Nov 12 at 17:30