ASKSAGE: Sage Q&A Forum - Individual question feedhttps://ask.sagemath.org/questions/Q&A Forum for SageenCopyright Sage, 2010. Some rights reserved under creative commons license.Sun, 21 Feb 2016 22:17:23 -0600Double Integralhttps://ask.sagemath.org/question/7636/double-integral/So I have equations:<br />
B(x,d,h,G)=(a01+b01*x+c01*x^2+d01*x^3) + (a02+b02*x+c02*x^2+d02*x^3)*h + (a11+b11*x+c11*x^2+d11*x^3)*(G/d) + (a12+b12*x+c12*x^2+d12*x^3)*(G/d)*h + (a21+b21*x+c21*x^2+d21*x^3)*(G/d)^2 + (a22+b22*x+c22*x^2+d22*x^3)*((G/d)^2)*h
normal(x,av,sd)=(1/(sd*sqrt(2*pi)))*exp(-(x-av)^2/(2*sd^2))
brk(x,d,h,D,H,Dstd,Hstd,G)=B(x,d,h,G)*normal(d,D,Dstd)*normal(h,H,Hstd)
(a##-d## are all decimals)
And I want to double integrate over d and h. So just for example I can integrate over just d:
f(x) = break1(x,d,11.16,2.85,11.16,0.61,13.85,.5).integrate(d,0,infinity)
f(500).numerical_approx()
1.04299933381999
and it evaluates just fine. Likewise if I put in an value for d and integrate over h, it also produces a value.
So I want the numerical approximation of this double integral, but when I try for it using:
f(x) = break1(x,d,h,2.85,11.16,0.61,13.85,.5).integrate(d,0,infinity).integrate(h,0,infinity)
f(500).numerical_approx()
/home/wil/sage/local/lib/python2.6/site-packages/sage/symbolic/integration/integral.pyc in _evalf_(self, f, x, a, b, parent)
199 # The gsl routine returns a tuple, which also contains the error.
200 # We only return the result.
--> 201 return numerical_integral(f, a, b)[0]
202
203 def _tderivative_(self, f, x, a, b, diff_param=None):
/home/wil/sage/local/lib/python2.6/site-packages/sage/gsl/integration.so in sage.gsl.integration.numerical_integral (sage/gsl/integration.c:1551)()
ValueError: Integrand has wrong number of parameters
At any rate, I think what is happening is that it is trying to evaluate the inside integral numerically first perhaps, which it is not able to do as it has a variable?
I tried using the lambda in the first answer, and I was able to evaluate, but I had to set it to max_points=10 to get an answer that was even close to correct, plus I could find no way to plot that one.
Thanks for any help!
WilWed, 25 Aug 2010 09:49:41 -0500https://ask.sagemath.org/question/7636/double-integral/Answer by kkumer for <p>So I have equations:<br/></p>
<pre><code>B(x,d,h,G)=(a01+b01*x+c01*x^2+d01*x^3) + (a02+b02*x+c02*x^2+d02*x^3)*h + (a11+b11*x+c11*x^2+d11*x^3)*(G/d) + (a12+b12*x+c12*x^2+d12*x^3)*(G/d)*h + (a21+b21*x+c21*x^2+d21*x^3)*(G/d)^2 + (a22+b22*x+c22*x^2+d22*x^3)*((G/d)^2)*h
normal(x,av,sd)=(1/(sd*sqrt(2*pi)))*exp(-(x-av)^2/(2*sd^2))
brk(x,d,h,D,H,Dstd,Hstd,G)=B(x,d,h,G)*normal(d,D,Dstd)*normal(h,H,Hstd)
</code></pre>
<p>(a##-d## are all decimals)</p>
<p>And I want to double integrate over d and h. So just for example I can integrate over just d:</p>
<pre><code>f(x) = break1(x,d,11.16,2.85,11.16,0.61,13.85,.5).integrate(d,0,infinity)
f(500).numerical_approx()
1.04299933381999
</code></pre>
<p>and it evaluates just fine. Likewise if I put in an value for d and integrate over h, it also produces a value.</p>
<p>So I want the numerical approximation of this double integral, but when I try for it using:</p>
<pre><code>f(x) = break1(x,d,h,2.85,11.16,0.61,13.85,.5).integrate(d,0,infinity).integrate(h,0,infinity)
f(500).numerical_approx()
/home/wil/sage/local/lib/python2.6/site-packages/sage/symbolic/integration/integral.pyc in _evalf_(self, f, x, a, b, parent)
199 # The gsl routine returns a tuple, which also contains the error.
200 # We only return the result.
--> 201 return numerical_integral(f, a, b)[0]
202
203 def _tderivative_(self, f, x, a, b, diff_param=None):
/home/wil/sage/local/lib/python2.6/site-packages/sage/gsl/integration.so in sage.gsl.integration.numerical_integral (sage/gsl/integration.c:1551)()
ValueError: Integrand has wrong number of parameters
</code></pre>
<p>At any rate, I think what is happening is that it is trying to evaluate the inside integral numerically first perhaps, which it is not able to do as it has a variable?</p>
<p>I tried using the lambda in the first answer, and I was able to evaluate, but I had to set it to max_points=10 to get an answer that was even close to correct, plus I could find no way to plot that one.</p>
<p>Thanks for any help!
Wil</p>
https://ask.sagemath.org/question/7636/double-integral/?answer=11552#post-id-11552You could try with numerical_integral. You didn't provide definition of B(), so here is example with simpler function:
sage: var('d h')
(d, h)
sage: integral(integral(exp(-(h-1)^2-(d-2)^2), d, 0,Infinity), h, 0, Infinity)
1/4*pi + 1/4*(pi + pi*erf(1))*erf(2) + 1/4*pi*erf(1)
sage: n(_)
2.88773776713433
sage: numerical_integral(lambda h: numerical_integral(lambda d:exp(-(h-1)^2-(d-2)^2), 0,Infinity)[0], 0, Infinity)
(2.8877377671374997, 2.8275945874944256e-06)
BTW: I'd like to know how to do the above without 'lambda'.
BTW2: Are you sure your integral converges at d=0?
Wed, 25 Aug 2010 13:06:07 -0500https://ask.sagemath.org/question/7636/double-integral/?answer=11552#post-id-11552Comment by willmwade for <p>You could try with numerical_integral. You didn't provide definition of B(), so here is example with simpler function:</p>
<pre><code>sage: var('d h')
(d, h)
sage: integral(integral(exp(-(h-1)^2-(d-2)^2), d, 0,Infinity), h, 0, Infinity)
1/4*pi + 1/4*(pi + pi*erf(1))*erf(2) + 1/4*pi*erf(1)
sage: n(_)
2.88773776713433
sage: numerical_integral(lambda h: numerical_integral(lambda d:exp(-(h-1)^2-(d-2)^2), 0,Infinity)[0], 0, Infinity)
(2.8877377671374997, 2.8275945874944256e-06)
</code></pre>
<p>BTW: I'd like to know how to do the above without 'lambda'.</p>
<p>BTW2: Are you sure your integral converges at d=0?</p>
https://ask.sagemath.org/question/7636/double-integral/?comment=22840#post-id-22840I figured out what the [0] was! It was limiting the return to only the value and leaving out the error calculation.Thu, 26 Aug 2010 06:14:41 -0500https://ask.sagemath.org/question/7636/double-integral/?comment=22840#post-id-22840Comment by jmarcellopereira for <p>You could try with numerical_integral. You didn't provide definition of B(), so here is example with simpler function:</p>
<pre><code>sage: var('d h')
(d, h)
sage: integral(integral(exp(-(h-1)^2-(d-2)^2), d, 0,Infinity), h, 0, Infinity)
1/4*pi + 1/4*(pi + pi*erf(1))*erf(2) + 1/4*pi*erf(1)
sage: n(_)
2.88773776713433
sage: numerical_integral(lambda h: numerical_integral(lambda d:exp(-(h-1)^2-(d-2)^2), 0,Infinity)[0], 0, Infinity)
(2.8877377671374997, 2.8275945874944256e-06)
</code></pre>
<p>BTW: I'd like to know how to do the above without 'lambda'.</p>
<p>BTW2: Are you sure your integral converges at d=0?</p>
https://ask.sagemath.org/question/7636/double-integral/?comment=32622#post-id-32622hi, willmwade
how to discovery the function of item [0]?Sun, 21 Feb 2016 22:17:23 -0600https://ask.sagemath.org/question/7636/double-integral/?comment=32622#post-id-32622Comment by willmwade for <p>You could try with numerical_integral. You didn't provide definition of B(), so here is example with simpler function:</p>
<pre><code>sage: var('d h')
(d, h)
sage: integral(integral(exp(-(h-1)^2-(d-2)^2), d, 0,Infinity), h, 0, Infinity)
1/4*pi + 1/4*(pi + pi*erf(1))*erf(2) + 1/4*pi*erf(1)
sage: n(_)
2.88773776713433
sage: numerical_integral(lambda h: numerical_integral(lambda d:exp(-(h-1)^2-(d-2)^2), 0,Infinity)[0], 0, Infinity)
(2.8877377671374997, 2.8275945874944256e-06)
</code></pre>
<p>BTW: I'd like to know how to do the above without 'lambda'.</p>
<p>BTW2: Are you sure your integral converges at d=0?</p>
https://ask.sagemath.org/question/7636/double-integral/?comment=22845#post-id-22845Well that does work for me as well. What is the "[0]" after each integral?
I am trying this on my original equation, but it will take a little while. I will post the result when it finally completes.
Thanks for your help!Thu, 26 Aug 2010 03:07:07 -0500https://ask.sagemath.org/question/7636/double-integral/?comment=22845#post-id-22845Comment by willmwade for <p>You could try with numerical_integral. You didn't provide definition of B(), so here is example with simpler function:</p>
<pre><code>sage: var('d h')
(d, h)
sage: integral(integral(exp(-(h-1)^2-(d-2)^2), d, 0,Infinity), h, 0, Infinity)
1/4*pi + 1/4*(pi + pi*erf(1))*erf(2) + 1/4*pi*erf(1)
sage: n(_)
2.88773776713433
sage: numerical_integral(lambda h: numerical_integral(lambda d:exp(-(h-1)^2-(d-2)^2), 0,Infinity)[0], 0, Infinity)
(2.8877377671374997, 2.8275945874944256e-06)
</code></pre>
<p>BTW: I'd like to know how to do the above without 'lambda'.</p>
<p>BTW2: Are you sure your integral converges at d=0?</p>
https://ask.sagemath.org/question/7636/double-integral/?comment=22849#post-id-22849Yes, even running it through numerical_integral I get the response of "Integrand has wrong number of parameters"
I broke down the equations to the smallest I could get and still get the error:
integral(integral((1/d)*h*e^d*e^(h/2),d,4,5),h,4,5)
That will cause the same error. Ideas?Thu, 26 Aug 2010 01:47:11 -0500https://ask.sagemath.org/question/7636/double-integral/?comment=22849#post-id-22849Comment by kkumer for <p>You could try with numerical_integral. You didn't provide definition of B(), so here is example with simpler function:</p>
<pre><code>sage: var('d h')
(d, h)
sage: integral(integral(exp(-(h-1)^2-(d-2)^2), d, 0,Infinity), h, 0, Infinity)
1/4*pi + 1/4*(pi + pi*erf(1))*erf(2) + 1/4*pi*erf(1)
sage: n(_)
2.88773776713433
sage: numerical_integral(lambda h: numerical_integral(lambda d:exp(-(h-1)^2-(d-2)^2), 0,Infinity)[0], 0, Infinity)
(2.8877377671374997, 2.8275945874944256e-06)
</code></pre>
<p>BTW: I'd like to know how to do the above without 'lambda'.</p>
<p>BTW2: Are you sure your integral converges at d=0?</p>
https://ask.sagemath.org/question/7636/double-integral/?comment=22848#post-id-22848Hm. I am puzzled. This works for me:
sage: version()
'Sage Version 4.4.4, Release Date: 2010-06-23'
sage: numerical_integral(lambda h: numerical_integral(lambda d:(1/d)*h*e^d*e^(h/2),4,5)[0],4,5)[0]
894.91270304016609
Thu, 26 Aug 2010 02:41:15 -0500https://ask.sagemath.org/question/7636/double-integral/?comment=22848#post-id-22848Comment by willmwade for <p>You could try with numerical_integral. You didn't provide definition of B(), so here is example with simpler function:</p>
<pre><code>sage: var('d h')
(d, h)
sage: integral(integral(exp(-(h-1)^2-(d-2)^2), d, 0,Infinity), h, 0, Infinity)
1/4*pi + 1/4*(pi + pi*erf(1))*erf(2) + 1/4*pi*erf(1)
sage: n(_)
2.88773776713433
sage: numerical_integral(lambda h: numerical_integral(lambda d:exp(-(h-1)^2-(d-2)^2), 0,Infinity)[0], 0, Infinity)
(2.8877377671374997, 2.8275945874944256e-06)
</code></pre>
<p>BTW: I'd like to know how to do the above without 'lambda'.</p>
<p>BTW2: Are you sure your integral converges at d=0?</p>
https://ask.sagemath.org/question/7636/double-integral/?comment=22843#post-id-22843The full equation takes a really long time. I am guessing that it is just trying to be too precise. I don't need much. Anyway to decrease the precision?Thu, 26 Aug 2010 04:43:53 -0500https://ask.sagemath.org/question/7636/double-integral/?comment=22843#post-id-22843