Weird behaviour in integrate code

[ancienthart]

Hi everyone, I've got the following behaviour in sage 4.8.

First I define the following function to calculate the fractional integral of a function:

def fracintegral(func,n, x=x,a=0):
    var('t')
    assume(t > a)
    assume(x > a)
    return integrate((x-t)^n*func(t),t,a,x)/gamma(n)

When I test it on polynomials and logarithms, it seems to work fine.

fracintegral(x,1/2)

 

4/15*x^(5/2)/sqrt(pi)

But when I try it on exponentials or trigonometric functions, I get the following weird error.

fracintegral(sin(x),1/2)

 

 

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "_sage_input_22.py", line 10, in <module>
    exec compile(u'open("___code___.py","w").write("# -*- coding: utf-8 -*-\\n" + _support_.preparse_worksheet_cell(base64.b64decode("ZnJhY2ludGVncmFsKHNpbih4KSwxLzIp"),globals())+"\\n"); execfile(os.path.abspath("___code___.py"))
  File "", line 1, in <module>
    
  File "/tmp/tmpeCkAoa/___code___.py", line 3, in <module>
    exec compile(u'fracintegral(sin(x),_sage_const_1 /_sage_const_2 )
  File "", line 1, in <module>
    
  File "/home/joal/bin/sage-4.8-linux-64bit-ubuntu_10.04.3_lts-x86_64-Linux/local/lib/python2.6/site-packages/sage/misc/latex.py", line 2243, in pretty_print
    view(object)
  File "/home/joal/bin/sage-4.8-linux-64bit-ubuntu_10.04.3_lts-x86_64-Linux/local/lib/python2.6/site-packages/sage/misc/latex.py", line 1969, in view
    s = _latex_file_(objects, title=title, sep=sep, tiny=tiny, debug=debug, **latex_options)
  File "/home/joal/bin/sage-4.8-linux-64bit-ubuntu_10.04.3_lts-x86_64-Linux/local/lib/python2.6/site-packages/sage/misc/latex.py", line 1624, in _latex_file_
    L = latex(x)
  File "/home/joal/bin/sage-4.8-linux-64bit-ubuntu_10.04.3_lts-x86_64-Linux/local/lib/python2.6/site-packages/sage/misc/latex.py", line 875, in __call__
    return LatexExpr(x._latex_())
  File "expression.pyx", line 667, in sage.symbolic.expression.Expression._latex_ (sage/symbolic/expression.cpp:4121)
  File "ring.pyx", line 596, in sage.symbolic.ring.SymbolicRing._latex_element_ (sage/symbolic/ring.cpp:6461)
  File "pynac.pyx", line 433, in sage.symbolic.pynac.py_latex_function (sage/symbolic/pynac.cpp:4702)
  File "pynac.pyx", line 407, in sage.symbolic.pynac.py_latex_function_pystring (sage/symbolic/pynac.cpp:4328)
TypeError: _limit_latex_() takes exactly 4 arguments (5 given)

But if I type in

integrate((x-t)^(1/2)*sin(t),t)

it seems to work, but with a really weird result (expintegral_e?)

-1/2*sqrt(t - x)*(((expintegral_e(-1/2, -I*t + I*x) +
expintegral_e(-1/2, I*t - I*x))*t - (expintegral_e(-1/2, -I*t + I*x) +
expintegral_e(-1/2, I*t - I*x))*x)*sin(x) - ((-I*expintegral_e(-1/2,
-I*t + I*x) + I*expintegral_e(-1/2, I*t - I*x))*x +
(I*expintegral_e(-1/2, -I*t + I*x) - I*expintegral_e(-1/2, I*t -
I*x))*t)*cos(x))

But the really weird bit, is that if I then enter my original command again, I now get a somewhat sensible answer.

-1/2*(x^(3/2)*sin(x)*expintegral_e(-1/2, -I*x) +
x^(3/2)*sin(x)*expintegral_e(-1/2, I*x) +
I*x^(3/2)*cos(x)*expintegral_e(-1/2, -I*x) -
I*x^(3/2)*cos(x)*expintegral_e(-1/2, I*x) -
2*limit(-1/2*sqrt(-t)*t*sin(x)*expintegral_e(-1/2, -I*t) -
1/2*sqrt(-t)*t*sin(x)*expintegral_e(-1/2, I*t) +
1/2*I*sqrt(-t)*t*cos(x)*expintegral_e(-1/2, -I*t) -
1/2*I*sqrt(-t)*t*cos(x)*expintegral_e(-1/2, I*t), t, 0, minus))/sqrt(pi)

Any ideas?

 

[kcrisman]

fracintegral(x,1/2)

 

4/15*x^(5/2)/sqrt(pi)

Though you should get a deprecation error.

 

But when I try it on exponentials or trigonometric functions, I get the following weird error.

fracintegral(sin(x),1/2)

 

I can't reproduce this.

 

But if I type in

integrate((x-t)^(1/2)*sin(t),t)

it seems to work, but with a really weird result (expintegral_e?)

-1/2*sqrt(t - x)*(((expintegral_e(-1/2, -I*t + I*x) +
expintegral_e(-1/2, I*t - I*x))*t - (expintegral_e(-1/2, -I*t + I*x) +
expintegral_e(-1/2, I*t - I*x))*x)*sin(x) - ((-I*expintegral_e(-1/2,
-I*t + I*x) + I*expintegral_e(-1/2, I*t - I*x))*x +
(I*expintegral_e(-1/2, -I*t + I*x) - I*expintegral_e(-1/2, I*t -
I*x))*t)*cos(x))

See http://trac.sagemath.org/sage_trac/ticket/11143.  This is

The generalized complex exponential integral `E_n(z)`

We just need to finish a few things on this, but you are welcome to use it now.  

[ancienthart]

On Wednesday, 21 March 2012 22:42:16 UTC+10, kcrisman wrote:

fracintegral(x,1/2)

 

4/15*x^(5/2)/sqrt(pi)

Though you should get a deprecation error.

Yeah, I do. I don't know the correct way to collect a variable argument from inside a function call, and apply that to a named argument call. :(
(t-x)^n*func(x=t) works fine, but don't know how to change "x=" to "y=" without using a really ugly apply call.
 

 

But when I try it on exponentials or trigonometric functions, I get the following weird error.

fracintegral(sin(x),1/2)

 

I can't reproduce this.

I've found that it seems to only show up for me if I start sage and run the defined function first. It is guaranteed to give me this error first time and will keep on giving me errors until I run the indefinite integral form of the same name function. Then it will sometimes give me correct-ish answers from then on. I've attached my sage worksheet. I'm on a 64-bit machine. Possible issue?
 

 

But if I type in

integrate((x-t)^(1/2)*sin(t),t)

it seems to work, but with a really weird result (expintegral_e?)

-1/2*sqrt(t - x)*(((expintegral_e(-1/2, -I*t + I*x) +
expintegral_e(-1/2, I*t - I*x))*t - (expintegral_e(-1/2, -I*t + I*x) +
expintegral_e(-1/2, I*t - I*x))*x)*sin(x) - ((-I*expintegral_e(-1/2,
-I*t + I*x) + I*expintegral_e(-1/2, I*t - I*x))*x +
(I*expintegral_e(-1/2, -I*t + I*x) - I*expintegral_e(-1/2, I*t -
I*x))*t)*cos(x))

See http://trac.sagemath.org/sage_trac/ticket/11143.  This is

The generalized complex exponential integral `E_n(z)`

We just need to finish a few things on this, but you are welcome to use it now.  

Ohhhh, this looks goooooood. :D
Is it possible to do this over a binary installation, or will I need to install and build from source? I'm asking because as well as my Linux install at home, I'm running a binary virtual image on Windows at work.

Thanks for this kcrisman,
Joal Heagney

[kcrisman]

On Thursday, March 22, 2012 5:46:58 AM UTC-4, ancienthart wrote:

On Wednesday, 21 March 2012 22:42:16 UTC+10, kcrisman wrote:

fracintegral(x,1/2)

 

4/15*x^(5/2)/sqrt(pi)

Though you should get a deprecation error.

Yeah, I do. I don't know the correct way to collect a variable argument from inside a function call, and apply that to a named argument call. :(
(t-x)^n*func(x=t) works fine, but don't know how to change "x=" to "y=" without using a really ugly apply call.
 

I'm not sure either offhand, and don't have time to think about it :(  I think this is an excellent type of question for ask.sagemath.org, especially since the new sage-support interface has been annoying for a lot of people on sage-support.

 

I've found that it seems to only show up for me if I start sage and run the defined function first. It is guaranteed to give me this error first time and will keep on giving me errors until I run the indefinite integral form of the same name function. Then it will sometimes give me correct-ish answers from then on. I've attached my sage worksheet. I'm on a 64-bit machine. Possible issue?

No idea, though hopefully someone will know.  It looks like you are "pretty printing", or at least so the traceback indicates, so maybe this is part of the problem.  Again, one would have to dive into the details.  Here I hope someone will take the time to do it on this list.  I can't even get RSS for the new Google groups interface, truly bizarre.

The generalized complex exponential integral `E_n(z)`

We just need to finish a few things on this, but you are welcome to use it now.  

Ohhhh, this looks goooooood. :D
Is it possible to do this over a binary installation, or will I need to install and build from source? I'm asking because as well as my Linux install at home, I'm running a binary virtual image on Windows at work.

You should be able to import the patch using hg_sage (see the developer guide for how to apply a patch, though I think you know how to do this) and then rebuild Sage using "sage -b", which will change the relevant files and then you should be in better shape with that.

[ancienthart]

Thanks for this. I'm very aware of how rushed you must be at the moment, so I really appreciate such a quick assist. (Reminds me of when Guido Van Rossum used to haunt the general python newsgroup.)
I am pretty printing, will attempt this again with pretty printing turned off and post the result back here. Will take care of the other points as you have suggested.
Again thank you for your help,

Joal Heagney

[ancienthart]

Hi kcrisman, just some followups:

Have added the following to ask.sagemath.org
http://ask.sagemath.org/question/1265/how-to-dynamically-substitute-a-variable-in-a

Have applied the http://trac.sagemath.org/sage_trac/ticket/11143 patch (along with this patch as well http://trac.sagemath.org/sage_trac/ticket/11513 )
and get this for a 1/2 integration of sin(x)

-1/2*(sqrt(x)*sin(x)*exp_integral_e(1/2, -I*x) +
sqrt(x)*sin(x)*exp_integral_e(1/2, I*x) +
I*sqrt(x)*cos(x)*exp_integral_e(1/2, -I*x) -
I*sqrt(x)*cos(x)*exp_integral_e(1/2, I*x) -
2*limit(1/2*sqrt(-t)*sin(x)*exp_integral_e(1/2, -I*t) +
1/2*sqrt(-t)*sin(x)*exp_integral_e(1/2, I*t) -
1/2*I*sqrt(-t)*cos(x)*exp_integral_e(1/2, -I*t) +
1/2*I*sqrt(-t)*cos(x)*exp_integral_e(1/2, I*t), t, 0, minus))/sqrt(pi)

This seems to be promising, except that minus needs to be replaced by -Infinity, and then I get an error message relating to the limit function.

ValueError: call the limit function like this, e.g. limit(expr, x=2).

I've also confirmed that the _limit_latex() error IS related to pretty printing.

[ancienthart]

Or rather, the limit syntax seems to be completely out of whack.

[kcrisman]

This seems to be promising, except that minus needs to be replaced by -Infinity, and then I get an error message relating to the limit function.

ValueError: call the limit function like this, e.g. limit(expr, x=2).

I've also confirmed that the _limit_latex() error IS related to pretty printing.

Hmm, you may want to post this at #11143 as well.   I don't get updates currently on this forum so sorry I haven't seen this.  That will help us there make this better.