Simplify conditions that integrate returns

89 views
Skip to first unread message

Ondřej Čertík

unread,
Sep 5, 2014, 11:43:51 AM9/5/14
to sympy
Hi,

I needed these two integrals for my work:

In [23]: integrate(erf(alpha*r)*exp(-alpha**2 * r**2), (r, 0, oo))
Out[23]:
⎧ ___
⎪ ╲╱ π │ ⎛ 2 ⎞│ π
⎪ ───── for │periodic_argument⎝polar_lift (α), ∞⎠│ ≤ ─
⎪ 4⋅α 2

⎪∞
⎨⌠
⎪⎮ 2 2
⎪⎮ -α ⋅r
⎪⎮ ℯ ⋅erf(α⋅r) dr otherwise
⎪⌡
⎪0


In [24]: integrate(erf(alpha*r)*exp(-alpha**2 * r**2)*r, (r, 0, oo))
Out[24]:
⎧ ___
⎪ ╲╱ 2 ⎛│ ⎛ 2 ⎞│ π │
⎪ ───── for ⎜│periodic_argument⎝polar_lift (α), ∞⎠│ ≤ ─ ∧ │p
⎪ 2 ⎝ 2
⎪ 4⋅α

⎨∞
⎪⌠
⎪⎮ 2 2
⎪⎮ -α ⋅r
⎪⎮ r⋅ℯ ⋅erf(α⋅r) dr
⎪⌡
⎩0


⎛ 2 ⎞│ π⎞ │ ⎛ 2
eriodic_argument⎝polar_lift (α), ∞⎠│ < ─⎟ ∨ │periodic_argument⎝polar_lift (α),
2⎠






otherwise




⎞│ π
∞⎠│ < ─
2




What does periodic_argument and polar_lift mean? It seems really
complicated. For comparison, Mathematica returns:

In [1]: Integrate[Erf[alpha*r]*Exp[-alpha^2*r^2], {r, 0, Infinity}]

Out [1]: ConditionalExpression[Sqrt[\[Pi]]/(4 alpha),
Re[alpha^2] > 0 && Re[alpha] > 0]

In [2]: Integrate[Erf[alpha*r]*Exp[-alpha^2*r^2]*r, {r, 0, Infinity}]

Out [2]: ConditionalExpression[1/(2 Sqrt[2] alpha^2),
Re[alpha^2] > 0 && Re[alpha] > 0]

Which is much simpler. So I think we should return something similar.

Ondrej

Aaron Meurer

unread,
Sep 5, 2014, 1:16:33 PM9/5/14
to sy...@googlegroups.com
If you set alpha as positive the arguments reduce.

These come from Tom Bachmann's GSoC project. The docstrings of those
functions explain what they are. Basically, in order to get correct
results, we have to keep track of complex numbers on the domain of the
complex logarithm, rather than just the complex plane.

Unfortunately, there has been no work done on simplifying these
conditions that are returned by the algorithm, even though in many
(most?) cases they can be simplified a lot.

There is als a conds flag that you can pass to integrate as
'piecewise', 'separate', or 'none'.

Aaron Meurer
> --
> You received this message because you are subscribed to the Google Groups "sympy" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to sympy+un...@googlegroups.com.
> To post to this group, send email to sy...@googlegroups.com.
> Visit this group at http://groups.google.com/group/sympy.
> To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/CADDwiVD%3DVsR5rcqDx%3D4w%3D0gYAXpx0ZT8NbiGc0Kta5xOk_%2Btew%40mail.gmail.com.
> For more options, visit https://groups.google.com/d/optout.

Aaron Meurer

unread,
Sep 5, 2014, 2:28:11 PM9/5/14
to sy...@googlegroups.com, Tom Bachmann
So looking closer, something that I don't fully understand is what
periodic_argument(x, oo) means. According to the docstring (which
should be made much clearer), periodic_argument(x, P) is a value in
(-P/2, P/2] via exp(P*I) = 1 (and it is the argument of x on some
branch of the logarithm).

Some things I am not clear on:

- Does the second argument of periodic_argument need to be a multiple
of pi? I tried a non-multiple and got an answer, but is it just
nonsense?

- What does periodic_argument(x, oo) mean? Is it the same as
unbranched_argument (which has no docstring BTW)?

Moving to polar_lift, this is easier to understand. It lifts the value
to the Riemann surface of the logarithm. So polar_lift(I) gives
exp_polar(I*pi/2).

So, if I understand the condition correctly, where alpha is a complex
number, does it mean that alpha**2 (should polar_lift(alpha)**2 be the
same as polar_lift(alpha**2)?) should have a nonnegative real part?

There are some more docs on this at
http://docs.sympy.org/latest/modules/integrals/g-functions.html if you
want to dig in further (by the way, is it just me or is the math not
rendering on that page?).

Aaron Meurer

Aaron Meurer

unread,
Sep 5, 2014, 3:05:18 PM9/5/14
to sy...@googlegroups.com, Tom Bachmann
I have fixed the math in the docs. See https://github.com/sympy/sympy/pull/7984.

Aaron Meurer

Ondřej Čertík

unread,
Sep 14, 2014, 5:17:01 PM9/14/14
to sympy, Tom Bachmann
So I think that refine() should be able to simplify things, but it
doesn't, e.g.:

In [5]: s = integrate(exp(-alpha**2*r**2)*r**2, (r, 0, oo))

In [6]: refine(s, Q.positive(alpha))
Out[6]:
⎧ ___
⎪ ╲╱ π │ ⎛ 2 ⎞│ π
⎪ ───── for │periodic_argument⎝polar_lift (α), ∞⎠│ < ─
⎪ 3 2
⎪ 4⋅α

⎨∞
⎪⌠
⎪⎮ 2 2
⎪⎮ 2 -α ⋅r
⎪⎮ r ⋅ℯ dr otherwise
⎪⌡
⎩0
> To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/CAKgW%3D6J28JQRJOwnzC7KYoq0D0BUG1sAsnQRArAK1fH5g165dg%40mail.gmail.com.

Aaron Meurer

unread,
Sep 14, 2014, 10:11:35 PM9/14/14
to sy...@googlegroups.com, Tom Bachmann
I don't think the new assumptions know about these functions.

Aaron Meurer
> To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/CADDwiVBzKB7anYM6LhMnmkYKmctZPmGsMvcAqRt72Ufcqwn4mg%40mail.gmail.com.
Reply all
Reply to author
Forward
0 new messages