Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Already on GitHub? Sign in to your account

Fresnel Integral Functions #938

Merged
merged 37 commits into from May 12, 2012

Conversation

Projects
None yet
6 participants
Contributor

raoulb commented Jan 6, 2012

This is a first implementation of the Fresnel integral functions S(x) and C(x).
It's the first part of Issue 2959: "Implement fresnel integrals".

The tests are still missing but I would like to get some feedback
before I proceed. (For example, are the function names ok?)

Owner

asmeurer commented Jan 7, 2012

Does Sage implement these? If so, what does it call them?

Contributor

raoulb commented Jan 7, 2012

In other CAS the names are:

Maxima: fresnel_s, fresnel_c
Maple: FresnelS, FresnelC
MMA: FresnelS, FresnelC
mpmath: fresnels, fresnelc

Sage probably calls maxima but I don't know.

Maybe we should use lower case only names too.

Contributor

raoulb commented Jan 8, 2012

Just added the testcases. Please review the branch thoroughly.

Contributor

raoulb commented Jan 8, 2012

Hmm, still having trouble with the limit function

In [2]: z = Symbol("z")

In [3]: limit(erf(z), z, oo)
Out[3]: 1

In [4]: limit(fresnels(z), z, oo)
Out[4]: 0

Why does it work for erf but not for fresnels?

Contributor

ness01 commented Jan 9, 2012

(reviewing this ASAP)

@ness01 ness01 commented on the diff Jan 9, 2012

sympy/functions/special/error_functions.py
+ def _eval_expand_complex(self, deep=True, **hints):
+ re_part, im_part = self.as_real_imag(deep=deep, **hints)
+ return re_part + im_part*S.ImaginaryUnit
+
+
+class fresnels(FresnelIntegral):
+ r"""
+ Fresnel integral S.
+
+ This function is defined by
+
+ .. math:: \operatorname{S}(z) = \int_0^z \sin{\frac{\pi}{2} t^2} \mathrm{d}t.
+
+ It is an entire function.
+
+ Examples
@asmeurer

asmeurer Jan 9, 2012

Owner

Actually, I like the order used here. Maybe I'm just used to Wikipedia, but to me, See Also and References should go last.

@ness01

ness01 Jan 9, 2012

Contributor

Actually, that's my personal preference and how I used to have it in all
my docstrings before it was suggested I change it g.

(Obviously you're the one to decide, not me ^^.)

On 09.01.2012 15:29, Aaron Meurer wrote:

  • def _eval_expand_complex(self, deep=True, **hints):
  •    re_part, im_part = self.as_real_imag(deep=deep, **hints)
    
  •    return re_part + im_part*S.ImaginaryUnit
    
    +class fresnels(FresnelIntegral):
  • r"""
  • Fresnel integral S.
  • This function is defined by
  • .. math:: \operatorname{S}(z) = \int_0^z \sin{\frac{\pi}{2} t^2} \mathrm{d}t.
  • It is an entire function.
  • Examples

Actually, I like the order used here. Maybe I'm just used to Wikipedia, but to me, See Also and References should go last.


Reply to this email directly or view it on GitHub:
https://github.com/sympy/sympy/pull/938/files#r336634

@raoulb

raoulb Jan 9, 2012

Contributor

I would also expect that the references are the last point. This is
the case for any research paper too. The referred material is probably
used less that the examples. But if you like I can exchange the order.

@ness01

ness01 Jan 9, 2012

Contributor

It's really not my call. But I agree with both you and aaron. Perhaps
the wiki should be changed?

On 09.01.2012 17:12, raoulb wrote:

  • def _eval_expand_complex(self, deep=True, **hints):
  •    re_part, im_part = self.as_real_imag(deep=deep, **hints)
    
  •    return re_part + im_part*S.ImaginaryUnit
    
    +class fresnels(FresnelIntegral):
  • r"""
  • Fresnel integral S.
  • This function is defined by
  • .. math:: \operatorname{S}(z) = \int_0^z \sin{\frac{\pi}{2} t^2} \mathrm{d}t.
  • It is an entire function.
  • Examples

I would also expect that the references are the last point. This is
the case for any research paper too. The referred material is probably
used less that the examples. But if you like I can exchange the order.


Reply to this email directly or view it on GitHub:
https://github.com/sympy/sympy/pull/938/files#r337005

@vperic

vperic Jan 9, 2012

Contributor

+1 on Examples/See Also/References as the preferred ordering. Should we bring this up on the mailing list? @miham, you wrote most of the file originally, why is it so? I guess just to follow NumPy standards?

@miham

miham Jan 9, 2012

Contributor

+1 on Examples/See Also/References as the preferred ordering. Should we bring this up on the mailing list? @miham, you wrote most of the file
originally, why is it so? I guess just to follow NumPy standards?

Yes, I just copied the order that is used in NumPy/SciPy. I already posted this question to the mailing list [1], but only Aaron and Ondrej voted. What do others think?

[1] http://groups.google.com/group/sympy/browse_thread/thread/7e958ed5edf61e86#

@ness01 ness01 and 1 other commented on an outdated diff Jan 9, 2012

sympy/functions/special/error_functions.py
+ ==========
+
+ .. [1] http://en.wikipedia.org/wiki/Fresnel_integral
+ .. [2] http://dlmf.nist.gov/7
+ .. [3] http://mathworld.wolfram.com/FresnelIntegrals.html
+ .. [4] http://functions.wolfram.com/GammaBetaErf/FresnelC
+ """
+
+ _trigfunc = C.cos
+ _sign = S.One
+
+ def _eval_rewrite_as_erf(self, z):
+ return (S.One-I)/4 * (erf((S.One+I)/2*sqrt(pi)*z) + I*erf((S.One-I)/2*sqrt(pi)*z))
+
+ def _eval_nseries(self, x, n, logx):
+ return x*(-x**4)**n*(2**(-2*n)*pi**(2*n))/((4*n+1)*C.factorial(2*n))
@ness01

ness01 Jan 9, 2012

Contributor

This shouldn't be necessary (and is not right as it stands anyway - this really should return a series of several terms, and an order. It is usually not necessary if differentiation is implemented, unless there are poles.).

@raoulb

raoulb Jan 9, 2012

Contributor

You are right, it is not necessary. However using diff is probably much slower, isn't it?
BTW: Maybe we should add symbolic order differentiation to sympy?

Which is the "best" example for this function where I can study how to do this right?

@raoulb

raoulb Jan 9, 2012

Contributor

BTW: The order term is not tight, I get:

In [5]: fresnels(z).series(n=12)
Out[5]: pi*z**3/6 - pi**3*z**7/336 + pi**5*z**11/42240 + O(z**12)

but it is really O(z**15) which is much better than O(z**12).

@ness01

ness01 Jan 9, 2012

Contributor

You can look at e.g. the trig functions. The way to go is apparently to
implement taylor_term. Personally I don't know if this is worth the
hassle (but it probably cannot hurt either :-)).

On 09.01.2012 17:19, raoulb wrote:

  • ==========
  • _trigfunc = C.cos
  • _sign = S.One
  • def _eval_rewrite_as_erf(self, z):
  •    return (S.One-I)/4 \* (erf((S.One+I)/2_sqrt(pi)_z) + I_erf((S.One-I)/2_sqrt(pi)*z))
    
  • def _eval_nseries(self, x, n, logx):
  •    return x_(-x__4)__n_(2**(-2_n)_pi**(2_n))/((4_n+1)_C.factorial(2_n))
    

You are right, it is not necessary. However using diff is probably much slower, isn't it?
BTW: Maybe we should add symbolic order differentiation to sympy?

Which is the "best" example for this function where I can study how to do this right?


Reply to this email directly or view it on GitHub:
https://github.com/sympy/sympy/pull/938/files#r337032

@ness01

ness01 Jan 9, 2012

Contributor

Yeah series code is a mess.

On 09.01.2012 17:23, raoulb wrote:

  • ==========
  • _trigfunc = C.cos
  • _sign = S.One
  • def _eval_rewrite_as_erf(self, z):
  •    return (S.One-I)/4 \* (erf((S.One+I)/2_sqrt(pi)_z) + I_erf((S.One-I)/2_sqrt(pi)*z))
    
  • def _eval_nseries(self, x, n, logx):
  •    return x_(-x__4)__n_(2**(-2_n)_pi**(2_n))/((4_n+1)_C.factorial(2_n))
    

BTW: The order term is not tight, I get:

In [5]: fresnels(z).series(n=12)
Out[5]: pi*z**3/6 - pi**3*z**7/336 + pi**5*z**11/42240 + O(z**12)

but it is really O(z**15) which is much better than O(z**12).


Reply to this email directly or view it on GitHub:
https://github.com/sympy/sympy/pull/938/files#r337040

@raoulb

raoulb Jan 9, 2012

Contributor

Maybe this is a gsoc task to clean up this mess.

@ness01 ness01 and 1 other commented on an outdated diff Jan 9, 2012

sympy/functions/special/error_functions.py
+
+ def _eval_rewrite_as_erf(self, z):
+ return (S.One-I)/4 * (erf((S.One+I)/2*sqrt(pi)*z) + I*erf((S.One-I)/2*sqrt(pi)*z))
+
+ def _eval_nseries(self, x, n, logx):
+ return x*(-x**4)**n*(2**(-2*n)*pi**(2*n))/((4*n+1)*C.factorial(2*n))
+
+ def _eval_aseries(self, n, args0, x, logx):
+ z = self.args[0]
+ #e = S.Half*I*pi*z**2
+ #h1 = C.hyper([S.One,S.Half],[],2*I/(pi*z**2))
+ #h2 = C.hyper([S.One,S.Half],[],-2*I/(pi*z**2))
+ #return (z**4)**C.Rational(3,4)/(2*z**3) + I/(2*pi*z)*(C.exp(-e)*h1 - C.exp(e)*h2)
+ return S.Half + 1/(pi*z)*C.sin(S.Half*pi*z**2)
+
+ def _eval_as_leading_term(self, x):
@ness01

ness01 Jan 9, 2012

Contributor

I believe this function is not needed either. (In fact I think it should be removed all over sympy, but that might need some more investigation ^^.)

@raoulb

raoulb Jan 9, 2012

Contributor

What is the sympy way to compute leading order behaviour of a function?
Just use the first time of series?

So you are right and we can remove this function too.

If we continue like that there will be not much left in the derived classes ^^

@ness01

ness01 Jan 9, 2012

Contributor

There is actually expr.compute_leading_term to do exactly this.

You are right the code size might shrink, but that's a good thing :-).

On 09.01.2012 17:25, raoulb wrote:

  • def _eval_rewrite_as_erf(self, z):
  •    return (S.One-I)/4 \* (erf((S.One+I)/2_sqrt(pi)_z) + I_erf((S.One-I)/2_sqrt(pi)*z))
    
  • def _eval_nseries(self, x, n, logx):
  •    return x_(-x__4)__n_(2**(-2_n)_pi**(2_n))/((4_n+1)_C.factorial(2_n))
    
  • def _eval_aseries(self, n, args0, x, logx):
  •    z = self.args[0]
    
  •    #e = S.Half_I_pi_z_*2
    
  •    #h1 = C.hyper([S.One,S.Half],[],2_I/(pi_z**2))
    
  •    #h2 = C.hyper([S.One,S.Half],[],-2_I/(pi_z**2))
    
  •    #return (z**4)**C.Rational(3,4)/(2_z__3) + I/(2_pi_z)_(C.exp(-e)_h1 - C.exp(e)_h2)
    
  •    return S.Half + 1/(pi_z)_C.sin(S.Half_pi_z**2)
    
  • def _eval_as_leading_term(self, x):

What is the sympy way to compute leading order behaviour of a function?
Just use the first time of series?

So you are right and we can remove this function too.

If we continue like that there will be not much left in the derived classes ^^


Reply to this email directly or view it on GitHub:
https://github.com/sympy/sympy/pull/938/files#r337048

@ness01 ness01 and 1 other commented on an outdated diff Jan 9, 2012

sympy/functions/special/error_functions.py
+ .. [1] http://en.wikipedia.org/wiki/Fresnel_integral
+ .. [2] http://dlmf.nist.gov/7
+ .. [3] http://mathworld.wolfram.com/FresnelIntegrals.html
+ .. [4] http://functions.wolfram.com/GammaBetaErf/FresnelC
+ """
+
+ _trigfunc = C.cos
+ _sign = S.One
+
+ def _eval_rewrite_as_erf(self, z):
+ return (S.One-I)/4 * (erf((S.One+I)/2*sqrt(pi)*z) + I*erf((S.One-I)/2*sqrt(pi)*z))
+
+ def _eval_nseries(self, x, n, logx):
+ return x*(-x**4)**n*(2**(-2*n)*pi**(2*n))/((4*n+1)*C.factorial(2*n))
+
+ def _eval_aseries(self, n, args0, x, logx):
@ness01

ness01 Jan 9, 2012

Contributor

Look at the gamma function implementation for what this function should do. (You need to return something oscillatory as you do, which probably means that this is not going to be much good for limits since gruntz cannot handle it.)

@raoulb

raoulb Jan 9, 2012

Contributor

You mean I should look at polygamma or loggamma? Gamma itself does not have this function.
Let me see what I can figure out.

@ness01

ness01 Jan 9, 2012

Contributor

Oh you are right. Either is fine (I think...).

On 09.01.2012 17:30, raoulb wrote:

  • _trigfunc = C.cos
  • _sign = S.One
  • def _eval_rewrite_as_erf(self, z):
  •    return (S.One-I)/4 \* (erf((S.One+I)/2_sqrt(pi)_z) + I_erf((S.One-I)/2_sqrt(pi)*z))
    
  • def _eval_nseries(self, x, n, logx):
  •    return x_(-x__4)__n_(2**(-2_n)_pi**(2_n))/((4_n+1)_C.factorial(2_n))
    
  • def _eval_aseries(self, n, args0, x, logx):

You mean I should look at polygamma or loggamma? Gamma itself does not have this function.
Let me see what I can figure out.


Reply to this email directly or view it on GitHub:
https://github.com/sympy/sympy/pull/938/files#r337073

@raoulb

raoulb Jan 9, 2012

Contributor

Hmm, for fresnel we can express the asymptotic expansion in 'closed' form, see:

http://functions.wolfram.com/GammaBetaErf/FresnelS/06/02/

@raoulb

raoulb Jan 9, 2012

Contributor

Btw, trying something like:

        z = self.args[0]
        result = S.Half + 1/(pi*z)*C.sin(S.Half*pi*z**2)
        o = C.Order(1/z**3, x)
        return result._eval_nseries(x, n, logx) + o 

to mimic the code of loggamma does only return an O term.
I must admit that I do not fully understand what's going on on the
sympy side here :-/

@ness01

ness01 Jan 9, 2012

Contributor

This is slightly tricky. The correct way to use aseries is a follows:

loggamma(1/x).series(x)

SymPy does not have a concept of "expansion at infinity". So
_eval_aseries should find an asymptotic series expansion as x -> 0, with
the understanding that self.args[0] will tend to infinity.

On 09.01.2012 18:02, raoulb wrote:

  • _trigfunc = C.cos
  • _sign = S.One
  • def _eval_rewrite_as_erf(self, z):
  •    return (S.One-I)/4 \* (erf((S.One+I)/2_sqrt(pi)_z) + I_erf((S.One-I)/2_sqrt(pi)*z))
    
  • def _eval_nseries(self, x, n, logx):
  •    return x_(-x__4)__n_(2**(-2_n)_pi**(2_n))/((4_n+1)_C.factorial(2_n))
    
  • def _eval_aseries(self, n, args0, x, logx):

Btw, trying something like:

         z = self.args[0]
         result = S.Half + 1/(pi*z)*C.sin(S.Half*pi*z**2)
         o = C.Order(1/z**3, x)
         return result._eval_nseries(x, n, logx) + o

to mimic the code of loggamma does only return an O term.
I must admit that I do not fully understand what's going on on the
sympy side here :-/


Reply to this email directly or view it on GitHub:
https://github.com/sympy/sympy/pull/938/files#r337207

@ness01 ness01 and 1 other commented on an outdated diff Jan 9, 2012

sympy/utilities/lambdify.py
@@ -44,7 +44,7 @@
#"uppergamma":"upper_gamma",
"LambertW":"lambertw",
"Matrix":"matrix",
- "conjugate":"conj",
@ness01

ness01 Jan 9, 2012

Contributor

I think you should avoid this sort of unrelated change.

@raoulb

raoulb Jan 9, 2012

Contributor

Oops, this is a relict from an earlier commit.

Contributor

ness01 commented Jan 9, 2012

Please add entries for fresnels and fresnelc to hyperexpand() table, so that e.g. integrate(sin(x**2), x, meijerg=True) returns a nice answer. There is an explanation on how to do this in the hyperexpand documentation (sphinx), if unclear ask here.

Regarding limits: there are two functions, "limit" and "gruntz". Limit is basically guesswork/heuristics whereas gruntz implements a "real" algorithm (guess which function I worked on a lot g). gruntz uses series expansions and cannot handle erf [simply because the necessary support code is not written. one needs to implement rewrite_as_tractable by defining a new helper function object, and then implement asymptotic series on the new object. ask if interested]. Making fresnel integrals work with gruntz would probably be non-trivial.

I'm not quite sure how limit works (except that it sometimes calls gruntz). Look there (simpy/series/limit.py) to figure out how to make the "standard" limits work without having to work with gruntz.

@ness01 ness01 commented on the diff Jan 9, 2012

sympy/functions/special/tests/test_error_functions.py
+
+ assert fresnelc(z).diff(z) == cos(pi*z**2/2)
+
+ assert fresnelc(z).as_leading_term(z) == z
+
+ assert fresnelc(z).rewrite(erf) == (S.One-I)/4 * (erf((S.One+I)/2*sqrt(pi)*z) + I*erf((S.One-I)/2*sqrt(pi)*z))
+
+ assert fresnelc(z)._eval_nseries(z, n, None) == z*(-z**4)**n*(2**(-2*n)*pi**(2*n))/((4*n+1)*factorial(2*n))
+
+ assert fresnelc(z)._eval_aseries(z, oo, 0, 0) == S.Half + sin(pi*z**2/2)/(pi*z)
+
+ assert fresnelc(w).is_real is True
+
+ assert fresnelc(z).as_real_imag() == ((fresnelc(re(z) - I*re(z)*Abs(im(z))/Abs(re(z)))/2 + fresnelc(re(z) + I*re(z)*Abs(im(z))/Abs(re(z)))/2,
+ I*(fresnelc(re(z) - I*re(z)*Abs(im(z))/Abs(re(z))) - fresnelc(re(z) + I*re(z)*Abs(im(z))/Abs(re(z))))*
+ re(z)*Abs(im(z))/(2*im(z)*Abs(re(z)))))
@ness01

ness01 Jan 9, 2012

Contributor

Please add tests for numerical evaluation (e.g. test against the integral definition by numerically evaluating an Integral).

@raoulb

raoulb Jan 9, 2012

Contributor

Shouldn't this go into the mpmath testsuite? Because the only thing I do is calling fresnel{s,c} from mpmath.

@ness01

ness01 Jan 9, 2012

Contributor

Yes. You want to test that numerically evaluating a sympy object
correctly calls the mpmath functions. (Or at least that's why I add
these sorts of tests. I'm not "authoritive" ;) )

On 09.01.2012 18:09, raoulb wrote:

  • assert fresnelc(z).diff(z) == cos(pi_z_*2/2)
  • assert fresnelc(z).as_leading_term(z) == z
  • assert fresnelc(z).rewrite(erf) == (S.One-I)/4 * (erf((S.One+I)/2_sqrt(pi)_z) + I_erf((S.One-I)/2_sqrt(pi)*z))
  • assert fresnelc(z).eval_nseries(z, n, None) == z(-z__4)_n(2**(-2_n)_pi**(2_n))/((4_n+1)_factorial(2_n))
  • assert fresnelc(z)._eval_aseries(z, oo, 0, 0) == S.Half + sin(pi_z__2/2)/(pi_z)
  • assert fresnelc(w).is_real is True
  • assert fresnelc(z).as_real_imag() == ((fresnelc(re(z) - I_re(z)_Abs(im(z))/Abs(re(z)))/2 + fresnelc(re(z) + I_re(z)_Abs(im(z))/Abs(re(z)))/2,
  •                                       I_(fresnelc(re(z) - I_re(z)_Abs(im(z))/Abs(re(z))) - fresnelc(re(z) + I_re(z)_Abs(im(z))/Abs(re(z))))_
    
  •                                       re(z)_Abs(im(z))/(2_im(z)*Abs(re(z)))))
    

Shouldn't this go into the mpmath testsuite? Because the only thing I do is calling fresnel{s,c} from mpmath.


Reply to this email directly or view it on GitHub:
https://github.com/sympy/sympy/pull/938/files#r337223

@raoulb

raoulb Jan 9, 2012

Contributor

Ah ok so testing if we use mpmath correctly. But in our case here we rely on the methods in Function
and do not intoduce anything new. For this reason I think there is nothing to test.

@ness01

ness01 Apr 7, 2012

Contributor

I still think this should be tested numerically.

@ness01

ness01 May 11, 2012

Contributor

Please add

from sympy.utilities.randtest import test_numerically

test_numerically(re(fresnelc(z)), fresnelc(z).as_real_imag()[0], z)
test_numerically(im(fresnelc(z)), fresnelc(z).as_real_imag()[1], z)
test_numerically(fresnelc(z), fresnelc(z).rewrite(meijerg), z)

etc

@raoulb

raoulb May 11, 2012

Contributor

Done. Do we want even more? (If yes, which ones)

@ness01

ness01 May 11, 2012

Contributor

On 11.05.2012 13:59, raoulb wrote:

  • assert fresnelc(z).diff(z) == cos(pi_z_*2/2)
  • assert fresnelc(z).as_leading_term(z) == z
  • assert fresnelc(z).rewrite(erf) == (S.One-I)/4 * (erf((S.One+I)/2_sqrt(pi)_z) + I_erf((S.One-I)/2_sqrt(pi)*z))
  • assert fresnelc(z).eval_nseries(z, n, None) == z(-z__4)_n(2**(-2_n)_pi**(2_n))/((4_n+1)_factorial(2_n))
  • assert fresnelc(z)._eval_aseries(z, oo, 0, 0) == S.Half + sin(pi_z__2/2)/(pi_z)
  • assert fresnelc(w).is_real is True
  • assert fresnelc(z).as_real_imag() == ((fresnelc(re(z) - I_re(z)_Abs(im(z))/Abs(re(z)))/2 + fresnelc(re(z) + I_re(z)_Abs(im(z))/Abs(re(z)))/2,
  •                                       I_(fresnelc(re(z) - I_re(z)_Abs(im(z))/Abs(re(z))) - fresnelc(re(z) + I_re(z)_Abs(im(z))/Abs(re(z))))_
    
  •                                       re(z)_Abs(im(z))/(2_im(z)*Abs(re(z)))))
    

Done. Do we want even more?

No I think that's enough.

Contributor

raoulb commented Jan 9, 2012

I think a good this to do is some documentation how to add new special functions and what to take care of.
Maybe I will write down some points once I understand this better.

Same applies for the extension of gruntz etc.

Contributor

ness01 commented Jan 9, 2012

I agree. If you can do this it would be awesome. I'll answer all
questions I can (but again, I'm not the one who wrote most of this code.
Though I modified quite a bit of it).

On 09.01.2012 18:05, raoulb wrote:

I think a good this to do is some documentation how to add new special functions and what to take care of.
Maybe I will write down some points once I understand this better.

Same applies for the extension of gruntz etc.


Reply to this email directly or view it on GitHub:
#938 (comment)

Contributor

raoulb commented Jan 9, 2012

About the hyperexpand stuff, I tried but my problem is that I only know fresnels(z) in terms of a 1F2 with argument -pi**2*z**4/16 and not z only:
http://functions.wolfram.com/GammaBetaErf/FresnelS/26/01/ShowAll.html

Hence I tried to transform this to match the form required by add(...) and came up with the rule:

add((S(3)/4,), (S(3)/2,S(7)/4), fresnels(root(-16*z/pi**2, 4)) * 6/(pi*(root(-16*z/pi**2, 4))**3)

The result of the integration is something A*fresnels(B*z) which looks promising
but A and B are big messy constants.

Contributor

raoulb commented Jan 9, 2012

This is a first try. The results seem to be ok IFF these ploar factors can simplify to 1.

If we take principal roots then in the case of fresnels we get a prefactor of I and
the same for the argument. Hence they combine into I*fresnels(I*z) which is fresnels(z).
For the case of fresnelc we get factors of I and -I as -I*fresnelc(I*z) which is fresnelc(z).

So I think the rules are correct. But it would be nice if we can tell sympy to fully simplify.

Contributor

raoulb commented Jan 9, 2012

I started a wiki page (I hope this is the right place and title.):
https://github.com/sympy/sympy/wiki/About-implementing-special-functions

I plan to fill in the gap while I'm lerning how to do all the parts.

Contributor

raoulb commented Jan 10, 2012

gruntz uses series expansions and cannot handle erf [simply because the necessary support code is not written. one needs to
implement rewrite_as_tractable by defining a new helper function object, and then implement asymptotic series on the new
object. ask if interested].

I tried to make erf traktable by Gruntz, please see my "erf_tractable" branch.
It seems to work ... quite astonishing given my current knowledge about sympy.

But I'm not 100% sure if it really takes the correct codepath, even though I removed
everything from the erf class including direct evaluation at oo and -oo.

I'm not quite sure how limit works (except that it sometimes calls gruntz). Look there (simpy/series/limit.py) to figure out how to
make the "standard" limits work without having to work with gruntz.

Hmm, I don't like these heuristics. (Hardcoded tan/cot ... will cause work to extend when updating my trig branch.)

Contributor

ness01 commented Jan 10, 2012

Your ideas are right, your implementation is poor ;). Comments:

  1. Your _eval_aseries does not work. Try erfs(1/x).series(x)...

This is for formal reasons (typos, missing imports) etc and I fixed this
for you in my branch erf.

  1. Your _eval_aseries is wrong.

Try (1-erf(1/x)).rewrite('tractable').replace(erfs, lambda t:erfs(t).series(x)). This will show you the asymptotic expansion of
erfc which you should compare to e.g. wikipedia.

A "realistic" way to see this is to run gruntz((1-erf(x))_exp(x__2)_x,
x, oo) - this should return 1/sqrt(pi) but it does not, because your
series starts only at 1/x**3. Wolfram alpha is a neat place to test limits.

  1. Why the limit still works with your broken code.

I was somewhat befuddled by this, since your aseries cannot even be
called. But all is well: rewrite_tractable must return a function with
only algebraic singularities (except for explicit exp/log terms), and
this enough knowledge to evaluate the limit.

  1. Cosmetic change.

Implement erfs._eval_rewrite_as_intractable (converting back to erf).
This is needed so that e.g. gruntz(erf(x), x, 0) does not return the
"obscure" erfs. [This is not the most elegant solution to the
tractable/intractable rewriting, but the way we currently do it.]

  1. Keep up the great work!

If you have the time, please brush up your code and submit it for
review. And if you are bored, similar extensions are possible for zeta,
Ei, and some bessel functions. You may wish to consult Gruntz' thesis
(see top of gruntz.py file) for extensive lists of examples.

On 10.01.2012 03:14, raoulb wrote:

gruntz uses series expansions and cannot handle erf [simply because the necessary support code is not written. one needs to
implement rewrite_as_tractable by defining a new helper function object, and then implement asymptotic series on the new
object. ask if interested].

I tried to make erf traktable by Gruntz, please see my "erf_tractable" branch.
It seems to work ... quite astonishing given my current knowledge about sympy.

But I'm not 100% sure if it really takes the correct codepath, even though I removed
everything from the erf class including direct evaluation at oo and -oo.

I'm not quite sure how limit works (except that it sometimes calls gruntz). Look there (simpy/series/limit.py) to figure out how to
make the "standard" limits work without having to work with gruntz.

Hmm, I don't like these heuristics. (Hardcoded tan/cot ... will cause work to extend when updating my trig branch.)


Reply to this email directly or view it on GitHub:
#938 (comment)

Contributor

ness01 commented Jan 10, 2012

I changed your fresnels entry slightly in my fresnel branch. I also
added unbranched = True to FresnelIntegral (this latter variable is used
by the function unpolarify). With these changes the integral comes out
nicely except for some gamma factor. Combsimp should simplify this but
it doesn't. One might argue whether this is a failure in combsimp or
gamma._eval_expand_func. (But it does not matter for this pull request.)

You might be worried about branch factors in the table entry I added
(since I nonchalantly did (-z)**(1/4) = exp(I*pi/4)*z**(1/4) etc). But
notice that the hypergeometric function is unbranched, hence so is the
right hand side ~ S(z**(1/4))/z**(3/4). It is thus fine to evaluate
z**(1/4) on the principal branch. Indeed you may verify that the
expanded version of hyper(.., x) agrees for all complex numbers x.

I suggest you change the fresnelc entry accordingly. Moreover, you
should compute better bases and use addb instead of add. A good way to
see why this is necessary is to try hyperexpand(hyper([S(3)/4-1], [S(3)/2, S(7)/4], x)).

"Computing a basis" basically means that you choose three functions B = [F1, F2, F3] such that this messy fresnels expression is a C(z)-linear
combination of them. You then find C = [C1, C2, C3] such that C.B = the
fresnel mess. Next you compute the matrix M such that z_B.diff(z) = M_B.
All of this should be chosen in such a way that the answers are "nice".

My suggestion for your basis would be F1 = (this fresnel mess), F2 = sqrt(z)*sinh((some constant)*sqrt(z) and F3 = cosh((same constant)*sqrt(z)).

On 09.01.2012 21:34, raoulb wrote:

This is a first try. The results seem to be ok IFF these ploar factors can simplify to 1.


Reply to this email directly or view it on GitHub:
#938 (comment)

@ness01 ness01 commented on an outdated diff Jan 10, 2012

sympy/functions/special/error_functions.py
@@ -82,3 +82,281 @@ def _eval_as_leading_term(self, x):
def _eval_is_real(self):
return self.args[0].is_real
+
+###############################################################################
+#################### FRESNEL INTEGRALS ########################################
+###############################################################################
+
+class FresnelIntegral(Function):
+ """ Base class for the Fresnel integrals."""
+
+ nargs = 1
+
+ _trigfunc = None
+ _sign = None
+
@ness01

ness01 Jan 10, 2012

Contributor

Remove these two lines. We decided in the gsoc-3 review that it is more helpful to leave this undefined.

Contributor

raoulb commented Jan 10, 2012

Making fresnel integrals work with gruntz would probably be non-trivial.

Maybe, but maybe not. We can rewrite the 'fresnel{s,c}' into 'erf'.
And if erf is tractable meybe we can do the fresnels this way.
It depends if we can do things like erf((1{+,-}I)/2*sqrt(pi)*z).
Currently we can not, Gruntz fails with an

NotImplementedError: Result depends on the sign of -sign((1/2 + I/2)**2)
Contributor

raoulb commented Jan 10, 2012

Your ideas are right, your implementation is poor ;).

I fully agree on the second part.

Comments: 1. Your _eval_aseries does not work. Try erfs(1/x).series(x)...
This is for formal reasons (typos, missing imports) etc
and I fixed this for you in my branch erf.

Some very stupid typos ... maybe it was a little bit too late last night.
This should be fixed now. The method was never called during limit
computation and I tested only erf(z), z->oo.

  1. Your _eval_aseries is wrong. Try (1-erf(1/x)).rewrite('tractable').replace(erfs, lambda t:erfs(t).series(x)). This
    will show you the asymptotic expansion of erfc which you should compare to e.g. wikipedia. A "realistic" way to
    see this is to run gruntz((1-erf(x))_exp(x__2)_x, x, oo) - this should return 1/sqrt(pi) but it does not, because your series
    starts only at 1/x**3. Wolfram alpha is a neat place to test limits.

I wanted to implement 5.17 of Gruntz, which is what I'm supposed to do
if I understand the text correctly. But obviously failed with the formula here too :-/

Fixed this (at least think so) and the realistic example now returns 1/sqrt(pi).
I'm still not sure about the big-O term.

  1. Why the limit still works with your broken code. I was somewhat befuddled by this, since your aseries cannot even be called.
    But all is well: rewrite_tractable must return a function with only algebraic singularities (except for explicit exp/log terms),
    and this enough knowledge to evaluate the limit.

Yes, the only function ever called was the rewrite tractable of erf.

  1. Cosmetic change. Implement erfs._eval_rewrite_as_intractable (converting back to erf). This is needed so that e.g.
    gruntz(erf(x), x, 0) does not return the "obscure" erfs.
    [This is not the most elegant solution to the tractable/intractable rewriting, but the way we currently do it.]

Done.

  1. Keep up the great work! If you have the time, please brush up your code and submit it for review.

Thanks! Yes I plan to do further work here, always as time permits.

Maybe I should merge this work into the erf improvement branch now.

And if you are bored, similar extensions are possible for zeta, Ei, and some bessel functions.

Yep :-)
Next to come after erf and fresnel. I plan to do the Airy functions too.

You may wish to consult Gruntz' thesis (see top of gruntz.py file) for extensive lists of examples.

I printed his thesis more than half a year ago but never finished reading it by now.
Maybe I should mention that I once had a lecture given by Dr. D. Gruntz (about software design).
It's a pitty that I never used this chance to talk about his algorithm.

Contributor

ness01 commented Jan 11, 2012

This is great. Please always notify me (@ness01) in the relevant pull
requests and CC me in bug reports (here "relevant" means related to
special functions / gruntz / hyperexpand / meijerint) since I'm usually
too busy to look through them.

On 10.01.2012 20:23, raoulb wrote:

Your ideas are right, your implementation is poor ;).

I fully agree on the second part.

Comments: 1. Your _eval_aseries does not work. Try erfs(1/x).series(x)...
This is for formal reasons (typos, missing imports) etc
and I fixed this for you in my branch erf.

Some very stupid typos ... maybe it was a little bit too late last night.
This should be fixed now. The method was never called during limit
computation and I tested only erf(z), z->oo.

  1. Your _eval_aseries is wrong. Try (1-erf(1/x)).rewrite('tractable').replace(erfs, lambda t:erfs(t).series(x)). This
    will show you the asymptotic expansion of erfc which you should compare to e.g. wikipedia. A "realistic" way to
    see this is to run gruntz((1-erf(x))_exp(x__2)_x, x, oo) - this should return 1/sqrt(pi) but it does not, because your series
    starts only at 1/x**3. Wolfram alpha is a neat place to test limits.

I wanted to implement 5.17 of Gruntz, which is what I'm supposed to do
if I understand the text correctly. But obviously failed with the formula here too :-/

Fixed this (at least think so) and the realistic example now return 1/sqrt(pi).
I'm still not sure about the big-O term.

  1. Why the limit still works with your broken code. I was somewhat befuddled by this, since your aseries cannot even be called.
    But all is well: rewrite_tractable must return a function with only algebraic singularities (except for explicit exp/log terms),
    and this enough knowledge to evaluate the limit.

Yes, the only function ever calles was the rewrite tractable of erf.

  1. Cosmetic change. Implement erfs._eval_rewrite_as_intractable (converting back to erf). This is needed so that e.g.
    gruntz(erf(x), x, 0) does not return the "obscure" erfs.
    [This is not the most elegant solution to the tractable/intractable rewriting, but the way we currently do it.]

Done.

  1. Keep up the great work! If you have the time, please brush up your code and submit it for review.

Thanks! Yes I plan to do further work here, always as time permits.

Maybe I should merge this work into the erf improvement branch now.

And if you are bored, similar extensions are possible for zeta, Ei, and some bessel functions.

Yep :-)
Next to come after erf and fresnel. I plan to do the Airy functions too.

You may wish to consult Gruntz' thesis (see top of gruntz.py file) for extensive lists of examples.

I printed his thesis more than half a year ago but never finished reading it by now.
Maybe I should mention that I once had a lecture given by Dr. D. Gruntz (about software design).
It's a pitty that I never used this chance to talk about his algorithm.


Reply to this email directly or view it on GitHub:
#938 (comment)

Owner

asmeurer commented Jan 11, 2012

The wiki is fine. When it gets fleshed out, we can migrate it to Sphinx.

Contributor

raoulb commented Jan 14, 2012

Rebased over master. Now we should review the series code and fix the remaining issues.

BTW: There ecists nice results for the Laplace transform. But sympy can not compute then yet.
What's missing?

Contributor

ness01 commented Jan 15, 2012

Oh yes, I almost forgot this. You should also extend the function _create_lookup_table at the top of sympy/integrals/meijerint.py. I'm afraid there is no documentation on how to do this, but I think it should be self-explanatory.

And then of course add tests for the transforms, and perhaps also some integrals involving fresenel functions (assuming they can be done, of course ^^).

@ness01 ness01 and 1 other commented on an outdated diff Jan 15, 2012

sympy/functions/special/error_functions.py
+ changed = True
+
+ nz = newarg.extract_multiplicatively(I)
+ if nz is not None:
+ prefact = cls._sign*I*prefact
+ newarg = nz
+ changed = True
+
+ if changed:
+ return prefact*cls(newarg)
+
+ # Values at infinities
+ if z is S.Infinity:
+ return S.Half
+ #elif z is S.NegativeInfinity:
+ # return -S.Half
@ness01

ness01 Jan 15, 2012

Contributor

Why is this commented out?

@ness01

ness01 Jan 15, 2012

Contributor

Ah I guess you extract signs automatically ... if so just delete this.

@raoulb

raoulb Jan 15, 2012

Contributor

Ok. I left a comment explaining why these cases are missing.

@ness01 ness01 and 1 other commented on an outdated diff Jan 15, 2012

sympy/functions/special/error_functions.py
+ ==========
+
+ .. [1] http://en.wikipedia.org/wiki/Fresnel_integral
+ .. [2] http://dlmf.nist.gov/7
+ .. [3] http://mathworld.wolfram.com/FresnelIntegrals.html
+ .. [4] http://functions.wolfram.com/GammaBetaErf/FresnelS
+ """
+
+ _trigfunc = C.sin
+ _sign = -S.One
+
+ def _eval_rewrite_as_erf(self, z):
+ return (S.One+I)/4 * (erf((S.One+I)/2*sqrt(pi)*z) - I*erf((S.One-I)/2*sqrt(pi)*z))
+
+ def _eval_nseries(self, x, n, logx):
+ return x**3*(-x**4)**n*(2**(-2*n-1)*pi**(2*n+1))/((4*n+3)*C.factorial(2*n+1))
@ness01

ness01 Jan 15, 2012

Contributor

If you really want to keep this code, implement it using taylor_term. See e.g. the sin function.

@raoulb

raoulb Jan 15, 2012

Contributor

Done.

I think is helpful for expressing the whole Taylor sum at once:

z = Symbol("z")
n = Symbol("n", integer=True)
Sum(fresnelc(z).taylor_term(n, z), (n,0,oo))

@ness01 ness01 and 1 other commented on an outdated diff Jan 15, 2012

sympy/functions/special/error_functions.py
+
+ .. [1] http://en.wikipedia.org/wiki/Fresnel_integral
+ .. [2] http://dlmf.nist.gov/7
+ .. [3] http://mathworld.wolfram.com/FresnelIntegrals.html
+ .. [4] http://functions.wolfram.com/GammaBetaErf/FresnelS
+ """
+
+ _trigfunc = C.sin
+ _sign = -S.One
+
+ def _eval_rewrite_as_erf(self, z):
+ return (S.One+I)/4 * (erf((S.One+I)/2*sqrt(pi)*z) - I*erf((S.One-I)/2*sqrt(pi)*z))
+
+ def _eval_nseries(self, x, n, logx):
+ return x**3*(-x**4)**n*(2**(-2*n-1)*pi**(2*n+1))/((4*n+3)*C.factorial(2*n+1))
+
@ness01

ness01 Jan 15, 2012

Contributor

I think the remaining two functions should be deleted.

(fresnels is intractable at infinity, but not even introducing a helper function is going to work here since gruntz (afaik) just cannot work with oscillatory limits. It cannot even do gruntz(sin(x)/x,x,oo))

@raoulb

raoulb Jan 15, 2012

Contributor

Done

@ness01 ness01 commented on the diff Jan 15, 2012

sympy/functions/special/tests/test_error_functions.py
+
+ assert fresnels(z).diff(z) == sin(pi*z**2/2)
+
+ assert fresnels(z).as_leading_term(z) == pi*z**3/6
+
+ assert fresnels(z).rewrite(erf) == (S.One+I)/4 * (erf((S.One+I)/2*sqrt(pi)*z) - I*erf((S.One-I)/2*sqrt(pi)*z))
+
+ assert fresnels(z)._eval_nseries(z, n, None) == z**3*(-z**4)**n*(2**(-2*n-1)*pi**(2*n+1))/((4*n+3)*factorial(2*n+1))
+
+ assert fresnels(z)._eval_aseries(z, oo, 0, 0) == S.Half - cos(pi*z**2/2)/(pi*z)
+
+ assert fresnels(w).is_real is True
+
+ assert fresnels(z).as_real_imag() == ((fresnels(re(z) - I*re(z)*Abs(im(z))/Abs(re(z)))/2 + fresnels(re(z) + I*re(z)*Abs(im(z))/Abs(re(z)))/2,
+ I*(fresnels(re(z) - I*re(z)*Abs(im(z))/Abs(re(z))) - fresnels(re(z) + I*re(z)*Abs(im(z))/Abs(re(z))))*
+ re(z)*Abs(im(z))/(2*im(z)*Abs(re(z)))))
@ness01

ness01 Jan 15, 2012

Contributor

consider testing this numerically

@raoulb

raoulb Jan 15, 2012

Contributor

Added a test with concrete numbers. Should I remove the symbolc expression test?

@ness01

ness01 Jan 15, 2012

Contributor

No, definitely keep the symbolic test.

What I meant was something like test_numerically(re(fresnels(z)), fresnels(z).as_real_imag()[0]). [Or, if you don't like random tests, test both sides on some specific value you choose.]

@ness01 ness01 and 1 other commented on an outdated diff Jan 15, 2012

sympy/functions/special/tests/test_error_functions.py
+ assert fresnels(-oo) == -S.Half
+
+ assert fresnels(z) == fresnels(z)
+ assert fresnels(-z) == -fresnels(z)
+ assert fresnels(I*z) == -I*fresnels(z)
+ assert fresnels(-I*z) == I*fresnels(z)
+
+ assert conjugate(fresnels(z)) == fresnels(conjugate(z))
+
+ assert fresnels(z).diff(z) == sin(pi*z**2/2)
+
+ assert fresnels(z).as_leading_term(z) == pi*z**3/6
+
+ assert fresnels(z).rewrite(erf) == (S.One+I)/4 * (erf((S.One+I)/2*sqrt(pi)*z) - I*erf((S.One-I)/2*sqrt(pi)*z))
+
+ assert fresnels(z)._eval_nseries(z, n, None) == z**3*(-z**4)**n*(2**(-2*n-1)*pi**(2*n+1))/((4*n+3)*factorial(2*n+1))
@ness01

ness01 Jan 15, 2012

Contributor

Test an actual series expansion, as in fresnels(z).series(z,n=5) or so.

@raoulb

raoulb Jan 15, 2012

Contributor

Done

@ness01 ness01 commented on an outdated diff Jan 15, 2012

sympy/simplify/hyperexpand.py
@@ -246,6 +247,14 @@ def fm(a, z):
Matrix([[1, 0, 0]]),
Matrix([[-S.Half, S.Half, 0], [0, -S.Half, S.Half], [0, 2*z, 0]]))
+ # FresnelS
+ #add([S(3)/4], [S(3)/2,S(7)/4], fresnels(root(-16*z/pi**2, 4)) * 6/(pi*(root(-16*z/pi**2, 4))**3) )
+ add([S(3)/4], [S(3)/2,S(7)/4], fresnels(2/sqrt(pi)*root(-z,4)) * 6/(pi*8*(-z)**(S(3)/4)/pi**(S(3)/2) ) )
+
+ # FresnelC
+ #add([S(1)/4], [S(1)/2,S(5)/4], fresnelc(root(-16*z/pi**2, 4)) / (pi*(root(-16*z/pi**2, 4))) )
+ add([S(1)/4], [S(1)/2,S(5)/4], fresnelc(2/sqrt(pi)*root(-z,4)) / (2/sqrt(pi)*root(-z,4)) )
@ness01

ness01 Jan 15, 2012

Contributor

Please change these as I explained in another comment; and compute bases (as also explained there).

Contributor

raoulb commented Jan 15, 2012

Now the only things missing are related to the hyperexpand and meijerg stuff. I'll try to
do these soon.

@ness01 : Big Thanks for reviewing my branches!

@ness01 ness01 commented on an outdated diff Jan 15, 2012

sympy/functions/special/error_functions.py
+ .. [1] http://en.wikipedia.org/wiki/Fresnel_integral
+ .. [2] http://dlmf.nist.gov/7
+ .. [3] http://mathworld.wolfram.com/FresnelIntegrals.html
+ .. [4] http://functions.wolfram.com/GammaBetaErf/FresnelS
+ """
+
+ _trigfunc = C.sin
+ _sign = -S.One
+
+ @staticmethod
+ def taylor_term(n, x, *previous_terms):
+ if n < 0:
+ return S.Zero
+ else:
+ x = sympify(x)
+ return x**3*(-x**4)**n*(S(2)**(-2*n-1)*pi**(2*n+1))/((4*n+3)*C.factorial(2*n+1))
@ness01

ness01 Jan 15, 2012

Contributor

I think it is costumary to declare this function @cacheit. Moreover, you should really use previous_terms ... if you want very many terms, evaluating the factorials and powers is probably slowest.

Contributor

ness01 commented Jan 15, 2012

Sure. It's the least I can do.

On 15.01.2012 17:23, raoulb wrote:

Now the only things missing are related to the hyperexpand and meijerg stuff. I'll try to
do these soon.

@ness01 : Big Thanks for reviewing my branches!


Reply to this email directly or view it on GitHub:
#938 (comment)

Contributor

raoulb commented Jan 23, 2012

Rebased on master

Contributor

ness01 commented Jan 23, 2012

Just as a reminder what I think has to be done before this can be pushed:

  • Fix the hyperexpand table entries (bases, coefficients); test that expressions for fresnel(cs) in terms of hypergeometric functions come out "nice". Also test contiguous functions (i.e. demonstrate that the basis is "nice").
  • As a culmination of this, make sure integral(sin(x**2)) etc work as expected.
  • Add a numerical test for fresnel(sc).as_real_imag.
  • Add tests for integration of fresnel(sc) (i.e. test your meijerg table entries).
Contributor

raoulb commented Apr 6, 2012

As a culmination of this, make sure integral(sin(x**2)) etc work as expected.

Done, but some nasty gamma remain. They should cancel to 1.

Add tests for integration of fresnel(sc) (i.e. test your meijerg table entries).

The command

integrate(fresnels(x), x)

gives me a strange error, not sure yet whats wrong here.

Contributor

ness01 commented Apr 7, 2012

I will review this tonight.

On 06.04.2012 16:44, raoulb wrote:

As a culmination of this, make sure integral(sin(x**2)) etc work as expected.

Done, but some nasty gamma remain. They should cancel to 1.

Add tests for integration of fresnel(sc) (i.e. test your meijerg table entries).

The command

integrate(fresnels(x), x)

gives me a strange error, not sure yet whats wrong here.


Reply to this email directly or view it on GitHub:
#938 (comment)

@ness01 ness01 and 2 others commented on an outdated diff Apr 7, 2012

sympy/functions/special/error_functions.py
+ if deep:
+ hints['complex'] = False
+ return (self.expand(deep, **hints), S.Zero)
+ else:
+ return (self, S.Zero)
+ if deep:
+ re, im = self.args[0].expand(deep, **hints).as_real_imag()
+ else:
+ re, im = self.args[0].as_real_imag()
+ return (re, im)
+
+ def as_real_imag(self, deep=True, **hints):
+ x, y = self._as_real_imag(deep=deep, **hints)
+ sq = -y**2/x**2
+ re = S.Half*(self.func(x+x*sqrt(sq))+self.func(x-x*sqrt(sq)))
+ im = x/(2*y) * sqrt(sq) * (self.func(x-x*sqrt(sq)) - self.func(x+x*sqrt(sq)))
@ness01

ness01 Apr 7, 2012

Contributor

Can you add a comment where this comes from?

@raoulb

raoulb Apr 7, 2012

Contributor

These are from here:

    # Fresnel S
    # http://functions.wolfram.com/06.32.19.0003.01
    # http://functions.wolfram.com/06.32.19.0006.01
    # Fresnel C
    # http://functions.wolfram.com/06.33.19.0003.01
    # http://functions.wolfram.com/06.33.19.0006.01 

Should I put this into the code?
Do we want links like these inside the code?

@ness01

ness01 Apr 7, 2012

Contributor

I'm not sure. @asmeurer ?

@asmeurer

asmeurer Apr 7, 2012

Owner

I think it's fine. The Wolfram Functions site is a solid resource.

@ness01 ness01 commented on an outdated diff Apr 7, 2012

sympy/functions/special/error_functions.py
+
+ _trigfunc = C.sin
+ _sign = -S.One
+
+ @staticmethod
+ @cacheit
+ def taylor_term(n, x, *previous_terms):
+ if n < 0:
+ return S.Zero
+ else:
+ x = sympify(x)
+ if len(previous_terms) > 1:
+ p = previous_terms[-1]
+ return (-pi**2*x**4*(4*n - 1)/(8*n*(2*n + 1)*(4*n + 3))) * p
+ else:
+ return x**3*(-x**4)**n*(S(2)**(-2*n-1)*pi**(2*n+1))/((4*n+3)*C.factorial(2*n+1))
@ness01

ness01 Apr 7, 2012

Contributor

PEP8 "requires" a number of spaces here. I suppose use your best judgement.

@ness01 ness01 commented on an outdated diff Apr 7, 2012

sympy/functions/special/tests/test_error_functions.py
+ assert fresnels(2+3*I).as_real_imag() == (fresnels(2 + 3*I)/2 + fresnels(2 - 3*I)/2, I*(fresnels(2 - 3*I) - fresnels(2 + 3*I))/2)
+
+ assert fresnelc(0) == 0
+ assert fresnelc(oo) == S.Half
+ assert fresnelc(-oo) == -S.Half
+
+ assert fresnelc(z) == fresnelc(z)
+ assert fresnelc(-z) == -fresnelc(z)
+ assert fresnelc(I*z) == I*fresnelc(z)
+ assert fresnelc(-I*z) == -I*fresnelc(z)
+
+ assert conjugate(fresnelc(z)) == fresnelc(conjugate(z))
+
+ assert fresnelc(z).diff(z) == cos(pi*z**2/2)
+
+ assert fresnelc(z).rewrite(erf) == (S.One-I)/4 * (erf((S.One+I)/2*sqrt(pi)*z) + I*erf((S.One-I)/2*sqrt(pi)*z))
@ness01

ness01 Apr 7, 2012

Contributor

again, spaces

@ness01 ness01 commented on an outdated diff Apr 7, 2012

sympy/integrals/tests/test_meijerint.py
@@ -586,3 +586,10 @@ def test_3153():
assert not expr.has(hyper)
# XXX the expression is a mess, but actually upon differentiation and
# putting in numerical values seems to work...
+
+
+def test_fresnel():
+ from sympy import fresnels, fresnelc
+
+ assert integrate(sin(pi*x**2/2),x) == 3*fresnels(x)*gamma(S(3)/4)/(4*gamma(S(7)/4))
@ness01

ness01 Apr 7, 2012

Contributor

apply expand_func to get rid of the gammas

Contributor

ness01 commented Apr 7, 2012

Regarding the strange integration problem: this is because of the prefactors you put in. They are more complicated my code can handle, and moreover they are not right here. You will see that there is an expression like t**(9/4)*(t**2)**(3/4)*(-t)**(-3/4). Note the unnatural way in which the exponents are written, etc. This is because on the wolfram function site, the variable is an ordinary complex number, and all the exponentiation is ordinary complex exponentiation, with branch cuts and everything. We are using polar numbers, so you can most likely just replace -1 by polar_lift(-1) throughout and then simplify "naively".

I understand that it may not at all be clear what to do, neither about the meijerint nor the hyperexpand table. If you want me to, I can write these entries for you (it's not trivial and will probably take one or two hours) and perhaps an extended documentation at the same time. Of course, if you want to "figure it out yourself", even better :-).

Member

jrioux commented May 1, 2012

SymPy Bot Summary: 🔴 There were test failures.

@raoulb: Please fix the test failures.

Test results html report: http://reviews.sympy.org/report/agZzeW1weTNyDAsSBFRhc2sY-OEWDA

Interpreter: /usr/bin/python (2.7.2-final-0)
Architecture: Linux (64-bit)
Cache: yes
Test command: setup.py test
master hash: 8426fe8
branch hash: 3609cc697de3b6baac5bc810c62c96295653f3af

Automatic review by SymPy Bot.

Member

jrioux commented May 1, 2012

SymPy Bot Summary: 🔴 There were test failures.

@raoulb: Please fix the test failures.

Test results html report: http://reviews.sympy.org/report/agZzeW1weTNyDAsSBFRhc2sYp6sWDA

Interpreter: /usr/bin/python (2.7.0-final-0)
Architecture: Linux (32-bit)
Cache: yes
Test command: setup.py test
master hash: 8426fe8
branch hash: 3609cc697de3b6baac5bc810c62c96295653f3af

Automatic review by SymPy Bot.

Contributor

raoulb commented May 7, 2012

SymPy Bot Summary: There were test failures.

@raoulb: Please fix the test failures.

Test results html report: http://reviews.sympy.org/report/agZzeW1weTNyDAsSBFRhc2sYzpUYDA

Interpreter: /usr/bin/python (2.7.3-candidate-2)
Architecture: Linux (32-bit)
Cache: yes
Test command: setup.py test
master hash: 25904b8
branch hash: 00895342003618a7e9a97d89b18d4e17705d77e1

Automatic review by SymPy Bot.

Contributor

raoulb commented May 10, 2012

This PR should be ready for a new final review.

There are still some (or rather many) nice integrals we can not solve.
For example S(z)^2, C(z)^2, S(z^2), C(z^2), S(z)*C(z) ...
All have solutions in terms of Fresenl and trigonometric functions.

Contributor

ness01 commented May 11, 2012

This needs to be rebased.

Contributor

raoulb commented May 11, 2012

Rebased

@ness01 ness01 commented on an outdated diff May 11, 2012

sympy/core/expr.py
@@ -554,7 +554,7 @@ def _eval_is_positive(self):
return False
try:
# check to see that we can get a value
- n2 = self._eval_evalf(1)
+ n2 = self.n(1)
@ness01

ness01 May 11, 2012

Contributor

@smichr @asmeurer

I'm somewhat unsure about this. The problem manifests as follows: gamma(3*exp_polar(I*pi)/2)._eval_evalf throws some mpmath internal exception. The problem is that _eval_evalf evaluates the argument to 1 binary digit accuracy, which apparently is -1, and then evaluates gamma at -1, causing the exception. I believe calling n instead alleviates the problem, since then a special routine for evaluating the gamma function is called, which knows about the exception, and tries evaluating at higher precision. [Which works, of course.]

I wonder if something else should be fixed.

@ness01 ness01 commented on the diff May 11, 2012

sympy/functions/__init__.py
@@ -23,9 +22,8 @@
asinh, acosh, atanh, acoth)
from sympy.functions.elementary.integers import floor, ceiling
from sympy.functions.elementary.piecewise import Piecewise, piecewise_fold
-
-from sympy.functions.special.error_functions import (erf, Ei, expint, E1,
- Si, Ci, Shi, Chi)
+from sympy.functions.special.error_functions import (erf, Ei, expint,
@ness01

ness01 May 11, 2012

Contributor

Can you add the sympy. back in? I know it does not really matter, but I don't see why we should break the pattern either.

@raoulb

raoulb May 11, 2012

Contributor

Should be ok with latest push

@ness01 ness01 commented on the diff May 11, 2012

sympy/functions/special/error_functions.py
+ -fresnels(z)
+
+ >>> fresnels(I*z)
+ -I*fresnels(z)
+
+ The Fresnel S integral obeys the mirror symmetry:
+
+ >>> from sympy import conjugate
+ >>> conjugate(fresnels(z))
+ fresnels(conjugate(z))
+
+ Differentiation with respect to z is supported:
+
+ >>> from sympy import diff
+ >>> diff(fresnels(z), z)
+ sin(pi*z**2/2)
@ness01

ness01 May 11, 2012

Contributor

Maybe add an integration example as well? Not too important.

@raoulb

raoulb May 11, 2012

Contributor

Done, I just used the defining integral.
We could also show the Laplace transform.

@ness01 ness01 and 1 other commented on an outdated diff May 11, 2012

sympy/functions/special/tests/test_error_functions.py
@@ -1,6 +1,6 @@
-from sympy import (symbols, expand, erf, nan, oo, Float, conjugate, sqrt, exp, pi, O, I, Ei,
- exp_polar, polar_lift, Symbol, I, exp, uppergamma, expint, log, loggamma, limit,
- meijerg, gamma, S, Shi, Chi, Si, Ci, E1, sin, cos, sinh, cosh)
+from sympy import (symbols, expand, expand_func, erf, nan, oo, Float, conjugate, sqrt, sin, cos, pi, re, im, Abs, O,
+ factorial, exp_polar, polar_lift, Symbol, I, integrate, exp, uppergamma, expint, log, loggamma, limit,
+ hyper, meijerg, gamma, S, Shi, Chi, Si, Ci, E1, Ei, sin, cos, sinh, cosh, fresnels, fresnelc)
@ness01

ness01 May 11, 2012

Contributor

These lines need to be wrapped.

@raoulb

raoulb May 11, 2012

Contributor

Done

@ness01 ness01 and 1 other commented on an outdated diff May 11, 2012

sympy/integrals/tests/test_transforms.py
def test_inverse_laplace_transform():
from sympy import (expand, sinh, cosh, besselj, besseli, exp_polar,
- unpolarify, simplify)
+ unpolarify, simplify, fresnels, fresnelc)
@ness01

ness01 May 11, 2012

Contributor

This does not seem necessary? [presumably you meant to check inverse transforms and realised they did not work ^^]

@raoulb

raoulb May 11, 2012

Contributor

No, just pasted into the wrong function and forgot about it.

@ness01 ness01 and 1 other commented on an outdated diff May 11, 2012

sympy/simplify/hyperexpand.py
@@ -162,6 +163,13 @@ def addb(ap, bq, B, C, M):
# This one is redundant.
add([-S.Half], [S.Half], exp(z) - sqrt(pi*z)*(-I)*erf(I*sqrt(z)))
+ # Added to get nice results for Laplace transform of Fresnel functions
+ # http://functions.wolfram.com/07.22.03.6437.01
+ add([1], [S(3)/4, S(5)/4],
+ sqrt(pi) * (cos(2*sqrt(polar_lift(-1)*z))*fresnelc(2*root(polar_lift(-1)*z,4)/sqrt(pi)) +
+ sin(2*sqrt(polar_lift(-1)*z))*fresnels(2*root(polar_lift(-1)*z,4)/sqrt(pi)))
+ / (2*root(polar_lift(-1)*z,4)))
@ness01

ness01 May 11, 2012

Contributor

I hate to say this, but you should convert this entry to addb format ...

@raoulb

raoulb May 11, 2012

Contributor

I tried to do, but it was not so obvious to me what to put into the basis. I have only three slots in the basis
but a few more function-terms which could be put there. Maybe we should discuss this on irc. But I won't be
there tonight.

@ness01

ness01 May 11, 2012

Contributor

Let w = exp(I*pi/4). I believe

sqrt(pi)/(2*w) z**(-1/4) * (cosh(2*sqrt(z))*fresnelc(2*w*z**(1/4)/sqr(pi)) + I*sinh(2*sqrt(z))*fresnels(2*w/sqrt(pi)*z**(1/4)))

sqrt(pi)/(2*w) z**(1/4) * (sinh(2*sqrt(z))*fresnelc(2*w*z**(1/4)/sqr(pi)) + I*cosh(2*sqrt(z))*fresnels(2*w/sqrt(pi)*z**(1/4)))

1

(the constant function 1) should make a good basis. Then the first row of M would be something like [-1/4, 1, 1/4], the last row would be [0, 0, 0] and the middle row (as well as checking my first row) I leave to you :-).

@raoulb

raoulb May 12, 2012

Contributor

Thanks, that was helpful! I didn't notice that grouping the terms this way leads to a good basis.

@ness01 ness01 commented on the diff May 11, 2012

sympy/utilities/lambdify.py
@@ -50,7 +50,7 @@
"Shi":"shi",
"Chi":"chi",
"Si":"si",
- "Ci":"ci",
+ "Ci":"ci"
@ness01

ness01 May 11, 2012

Contributor

this looks like an unrelated and unnecessary change …?

@raoulb

raoulb May 11, 2012

Contributor

Do I really have to undo this now?

@ness01

ness01 May 11, 2012

Contributor

On 11.05.2012 14:18, raoulb wrote:

@@ -50,7 +50,7 @@
"Shi":"shi",
"Chi":"chi",
"Si":"si",

  • "Ci":"ci",
  • "Ci":"ci"

Do I really have to undo this now?


Reply to this email directly or view it on GitHub:
https://github.com/sympy/sympy/pull/938/files#r808348

I guess not

Contributor

ness01 commented May 12, 2012

(I'm waiting for the test runner but to come round once more)

raoulb added some commits Jan 6, 2012

@raoulb raoulb First implementation of Fresnel integrals fb2b533
@raoulb raoulb Complex expansion for Fresnel integrals 7d75819
@raoulb raoulb Asymptotic series expansions for Fresnel integrals b5d6444
@raoulb raoulb Leading term expansion for Fresnel integrals 9dc7c0f
@raoulb raoulb Numerical evaluation for Fresnel integrals using mpmath ab67ab7
@raoulb raoulb Series expansions for Fresnel integrals e5f83d9
@raoulb raoulb Documentation and examples for Fresnel integrals bc3fb2d
@raoulb raoulb Better asymptotic series expansions for Fresnel integrals e4953d0
@raoulb raoulb Examples in an IPython notebook 75c3ed3
@raoulb raoulb Fixed numerical evaluation of Fresnel integrals. 5f909a4
@raoulb raoulb Updated examples in an IPython notebook 8c5a8e2
@raoulb raoulb Renamed Fresnel integrals to 'fresnels' and 'fresnelc' 379873a
@raoulb raoulb Add testcases for the Fresnel integrals 5080d5a
@raoulb raoulb Adding Fresnel integrals to the hyperexpand Tables a4e7985
@raoulb raoulb Declared the Fresnel integrals to be unbranched. 53c03b1
@raoulb raoulb Apply comments from pull request review cd79acc
@raoulb raoulb Added Meijer G representation of Fresnel Integrals e3b854a
@raoulb raoulb Caching 'taylor_term' functions 19685ac
@raoulb raoulb Improved Taylor terms of Fresnel integrals 86741fe
@raoulb raoulb Test integrals defining Fresnel functions e002dd8
@raoulb raoulb Addressing latest comments on code 0cf95a1
@raoulb raoulb Rules to rewrite fresnel{s,c} functions in terms of 1F2 hypergeometri…
…c functions
b9d4933
@raoulb raoulb Fixed hypergeometric formula for Fresnel functions 694d99d
@raoulb raoulb Fix test_args test failures d26024f
@raoulb raoulb Fix integration of Fresnel functions a44422f
@raoulb raoulb Add tests for Fresnel integration 92f48ca
Contributor

raoulb commented May 12, 2012

With this change, I hope the Fresnel stuff is finished now.

Contributor

ness01 commented May 12, 2012

Actually, no. Something you did makes test_cg hang. I have no idea what, but you'll have to bisect this. [Note that the test runs in less then a few seconds in master.]

Contributor

raoulb commented May 12, 2012

Ah, so I'm responsible that it hangs :-/

Contributor

raoulb commented May 12, 2012

BTW: See also my new Fresnel examples notebook

Contributor

ness01 commented May 12, 2012

Well it could just as easily be one of the patches I sent you, but it is confined to your branch, yes.

Contributor

raoulb commented May 12, 2012

This one hangs:

cg_simp(2 * Sum(CG(1,alpha,0,0,1,alpha), (alpha,-1,1))) == 6

Contributor

raoulb commented May 12, 2012

Actually, I think I found it to be your patch for the gamma poles.
At least without this change it runs fine.

Contributor

ness01 commented May 12, 2012

On 12.05.2012 16:16, raoulb wrote:

Actually, I think I found it to be your patch for the gamma poles.
At least without this change it runs fine.

I'll look into it again, then …

Contributor

ness01 commented May 12, 2012

Can you try this patch instead of the problematic commit:

diff --git a/sympy/core/function.py b/sympy/core/function.py
index 01c942d..22eb252 100644
--- a/sympy/core/function.py
+++ b/sympy/core/function.py
@@ -366,8 +366,12 @@ def _eval_evalf(self, prec):
                 return

         # Convert all args to mpf or mpc
+        # Convert the arguments to *higher* precision than requested for the
+        # final result.
+        # XXX + 5 is a guess, it is similar to what is used in evalf.py. Should
+        #     we be more intelligent about it?
         try:
-            args = [arg._to_mpmath(prec) for arg in self.args]
+            args = [arg._to_mpmath(prec + 5) for arg in self.args]
         except ValueError:
             return

raoulb added some commits May 8, 2012

@raoulb raoulb Tests for the gamma pole patch 2f579d2
@raoulb raoulb A new rule in hyperexpand to get nice results for the Laplace transfo…
…rm of Fresnel functions

Examples without this formula:

  laplace_transform(fresnels(t), t, s, noconds=True)

   ⎛         ⎛         │    4 ⎞        ⎛  2⎞        ⎛  2⎞⎞
   ⎜     ┌─  ⎜   1     │  -s  ⎟        ⎜ s ⎟        ⎜ s ⎟⎟
  -⎜2⋅s⋅ ├─  ⎜         │ ─────⎟ - π⋅sin⎜───⎟ - π⋅cos⎜───⎟⎟
   ⎜    1╵ 2 ⎜3/4, 5/4 │     2⎟        ⎝2⋅π⎠        ⎝2⋅π⎠⎟
   ⎝         ⎝         │ 16⋅π ⎠                          ⎠
  ────────────────────────────────────────────────────────
                           2⋅π⋅s

  laplace_transform(fresnelc(t), t, s, noconds=True)

            ⎛         │    4 ⎞           ⎛  2⎞           ⎛  2⎞
     3  ┌─  ⎜   1     │  -s  ⎟      2    ⎜ s ⎟      2    ⎜ s ⎟
  2⋅s ⋅ ├─  ⎜         │ ─────⎟ - 3⋅π ⋅sin⎜───⎟ + 3⋅π ⋅cos⎜───⎟
       1╵ 2 ⎜5/4, 7/4 │     2⎟           ⎝2⋅π⎠           ⎝2⋅π⎠
            ⎝         │ 16⋅π ⎠
  ────────────────────────────────────────────────────────────
                                2
                             6⋅π ⋅s

And with this formula we can expand the 1F2 above:

  laplace_transform(fresnels(t), t, s, noconds=True)

       ⎛  2⎞                  ⎛  2⎞      ⎛  2⎞                  ⎛  2⎞
       ⎜ s ⎟         ⎛s⎞      ⎜ s ⎟      ⎜ s ⎟         ⎛s⎞      ⎜ s ⎟
    sin⎜───⎟⋅fresnels⎜─⎟   sin⎜───⎟   cos⎜───⎟⋅fresnelc⎜─⎟   cos⎜───⎟
       ⎝2⋅π⎠         ⎝π⎠      ⎝2⋅π⎠      ⎝2⋅π⎠         ⎝π⎠      ⎝2⋅π⎠
  - ──────────────────── + ──────── - ──────────────────── + ────────
             s               2⋅s               s               2⋅s

  laplace_transform(fresnelc(t), t, s, noconds=True)

     ⎛  2⎞                  ⎛  2⎞      ⎛  2⎞                  ⎛  2⎞
     ⎜ s ⎟         ⎛s⎞      ⎜ s ⎟      ⎜ s ⎟         ⎛s⎞      ⎜ s ⎟
  sin⎜───⎟⋅fresnelc⎜─⎟   sin⎜───⎟   cos⎜───⎟⋅fresnels⎜─⎟   cos⎜───⎟
     ⎝2⋅π⎠         ⎝π⎠      ⎝2⋅π⎠      ⎝2⋅π⎠         ⎝π⎠      ⎝2⋅π⎠
  ──────────────────── - ──────── - ──────────────────── + ────────
           s               2⋅s               s               2⋅s

This result is nice and correct.
14ddbc1
@raoulb raoulb Adapt tests to new Laplace transformations of Fresnels 2a47c0b
@raoulb raoulb Fix a remaining test failure 748d08b
@raoulb raoulb Rewrite Fresnel functions in terms of Meijer G e370757
@raoulb raoulb Tests for Fresnel rewrites 248365e
@raoulb raoulb A different patch for the gamma pole issue (by ness) f3d3bc6
@raoulb raoulb Fresnel function latex printing b3ba93f
@raoulb raoulb Add integral example to Fresnel docstring 62bd77c
@raoulb raoulb Addressed points from review discussion 008560c
@raoulb raoulb Manually tuned rule for Laplace transform of Fresnels 56c6af0
Contributor

raoulb commented May 12, 2012

SymPy Bot Summary: All tests have passed.

Test results html report: http://reviews.sympy.org/report/agZzeW1weTNyDAsSBFRhc2sYh7AXDA

Interpreter: /usr/bin/python (2.7.3-candidate-2)
Architecture: Linux (32-bit)
Cache: yes
Test command: setup.py test
master hash: 4c213da
branch hash: 56c6af0

Automatic review by SymPy Bot.

Contributor

raoulb commented May 12, 2012

Seems to work ok.

Contributor

ness01 commented May 12, 2012

SymPy Bot Summary: ✳️ All tests have passed.

Test results html report: http://reviews.sympy.org/report/agZzeW1weTNyDAsSBFRhc2sYyO4XDA

Interpreter: /usr/bin/python (2.7.3-candidate-2)
Architecture: Linux (64-bit)
Cache: yes
Test command: setup.py test
master hash: 4c213da
branch hash: 56c6af0

Automatic review by SymPy Bot.

@ness01 ness01 added a commit that referenced this pull request May 12, 2012

@ness01 ness01 Merge pull request #938 from raoulb/fresnel
Fresnel Integral Functions
30cb220

@ness01 ness01 merged commit 30cb220 into sympy:master May 12, 2012

Contributor

ness01 commented May 12, 2012

This is in. Good work!

Contributor

raoulb commented May 12, 2012

Thanks :-)
And also thanks for your support.

@raoulb raoulb deleted the raoulb:fresnel branch Mar 26, 2015

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment