Factorial Precision

[tobia...@gmail.com]

Hello,

I need to integrate sinc(x) twice. To do this I have used the series expansion for sinc, which is a summation that includes a factorial term. I have attached my current code.

If I use the factorial function in sympy, this sequence does not appear to converge. If I use mpmath.factorial, it appears to converge after 8 terms (the number of terms can be changed on line 14). However this "convergence" is very abrupt - there is a large relative change between 7 & 8 terms, and then absolutely no change for 9 terms. Hence I am concerned I am hitting some precision limit and not actual convergence.

I would be very grateful if someone can sanity check my code. I am unfamiliar with the mpmath module, I only stumbled across the mpmath.factorial function through Google.

Thanks in advance,

Toby Wood

[Aaron Meurer]

For 9 terms, mpmath.fac gives me 2.14373016357422 and sympy.factorial gives me 1.88592361267614 (also mpmath.fac is much slower; SymPy in general handles exact integers better than floats). Given the two, I'm inclined to believe the SymPy version, since it deals with exact numbers, so there is no chance for numerical instability until the end (and the intermediate T_rfe has a lot of very small floats, so I think there is indeed a lot of instability going on).

Some other notes on your code:

- Si is already implemented in SymPy. You can use Si(x).series(x, 0, n).removeO() to get the series of Si up to x**n. 

- Just use sympy.pi. It will know how to evaluate itself in the end.

- Use T_rfe.evalf(subs={N: 8}) at the end to evaluate the expression, or if you want to use numpy, use lambdify. If you are worried about numerical instability you can increase the precision (at the end) by passing an argument to evalf(), like T_rfe.evalf(100, subs={N: 8}).  

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/05839053-a5ed-4c34-b02b-fe49773fa5c6%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[tobia...@gmail.com]

Hello Aaron,

Thank you so much for the help. I only dig out SymPy about once a year so am not as familiar with it as I would like.

I have attached updated code with alterations that you suggested. The formulation is slightly different (and simpler). The values I get for odd values of N look right now - they converge in a sensible fashion and match what I expect for this problem (You can see as N increases they tend towards ~0.1). However the even values of N do not look right - they oscillate about 0 as I increase the number of terms in the series. They should all be positive and follow the even values. Adding extra precision with evalf() did not seem to help.

I guess at this point it's actually to do with the properties of the sinc integral, but thought I should check that I'm not doing something stupid (again)!

Thanks again,

Toby

[Aaron Meurer]

My only suggestion here is to check the expressions to see if they are correct, either by working them out by hand, or checking against some other CAS or WolframAlpha. That will tell you if there is a problem in your formulation or a bug in SymPy.

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/2e5b13a0-bbf9-4334-80b5-6d6e4ca8b56c%40googlegroups.com.

[tobia...@gmail.com]

Thanks, I will.

(I think there might a bug in Wolfram Alpha! It claims the series expansion of integrate(integrate(sinc(x))) is 1 + x^2/2 - x^4/72. SymPy and doing it by hand do not have the leading 1)

[Aaron Meurer]

On Thu, Jul 16, 2015 at 10:33 AM, <tobia...@gmail.com> wrote:

Thanks, I will.

(I think there might a bug in Wolfram Alpha! It claims the series expansion of integrate(integrate(sinc(x))) is 1 + x^2/2 - x^4/72. SymPy and doing it by hand do not have the leading 1)

Remember that integrate() is allowed to add an arbitrary constant to an expression. The constant and linear term of that expansion are meaningless without integration limits or some kind of boundary or initial conditions, and can be anything.

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/40258d7e-4153-4f74-b8fa-253e7b10d053%40googlegroups.com.

[tobia...@gmail.com]

For completeness sake, there does not appear to be a bug in SymPy. A colleague ran some numerical simulations of this problem and we observed the same oscillatory behaviour for odd/even values of N.

Thanks again for the help,

Toby