Skip to content

Instantly share code, notes, and snippets.

@mforets

mforets/integrals.md

Last active April 19, 2018 15:45
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save mforets/e82a3de96a17bc80a4a22ad787bdd104 to your computer and use it in GitHub Desktop.
Save mforets/e82a3de96a17bc80a4a22ad787bdd104 to your computer and use it in GitHub Desktop.
Download ZIP
a bunch of integrals: exploring (non-commercial) symbolic integrators available through SageMath
Raw
integrals.md

t/17469

sage: integrate(exp(-x)*sinh(sqrt(x)), x, 0, oo, algorithm='maxima').simplify() # ok
1/4*(sqrt(pi)*(erf(1) - 1) + sqrt(pi) + 2*e^(-1) - 2)*e^(1/4) - 1/4*(sqrt(pi)*(erf(1) - 1) - sqrt(pi) + 2*e^(-1) - 2)*e^(1/4)
sage: _.canonicalize_radical()
1/2*sqrt(pi)*e^(1/4)
sage: integrate(exp(-x)*sinh(sqrt(x)), x, 0, oo, algorithm='giac').simplify() # ok, fast
1/2*sqrt(pi)*e^(1/4)
sage: integrate(exp(-x)*sinh(sqrt(x)), x, 0, oo, algorithm='sympy') # ok, fast, reduced
1/2*sqrt(pi)*e^(1/4)

t/15504

sage: integrate(log(1+x)/x, x, 0, 1, algorithm='giac') # don't know
integrate(ln(x + 1)/x, x, 0, 1)
sage: integrate(log(1+x)/x, x, 0, 1, algorithm='maxima')  # ok
-1/6*pi^2 + I*pi*log(2) + dilog(2)
sage: integrate(log(1+x)/x, x, 0, 1, algorithm='sympy')  # ok, reduced
1/12*pi^2

t/15747

sage: ys=0.06+2*(x-0.0275)
sage: xs=ys/1.516
sage: i=1/(xs/(1+xs)-x/(1+x))
sage: numerical_integral(i, 0, 0.275)
(12.325862793647284, 4.191806193893343e-11)
sage: integral(i, x, 0, 0.275, algorithm='maxima')  # bug, misuse (*)
(segfault)

(*) integrand contains floating points -- use numerical_integral instead, or:

sage: integral(i._convert(QQ), x, 0, 0.275, algorithm='maxima')
-15715975729016872675835131445395682818116292620105357/5122797704485899255665839525816694137143334669254523*log(1780848786854484) + 120.05607923574551
sage: _.n()
12.3258627936473
sage: integral(i._convert(QQ), x, 0, QQ(0.275), algorithm='maxima')  # exact if endpoints are transformed too
15715975729016872675835131445395682818116292620105357/5122797704485899255665839525816694137143334669254523*log(1967481739716834977/40) - 15715975729016872675835131445395682818
116292620105357/5122797704485899255665839525816694137143334669254523*log(1780848786854484) + 637452858164353703048222377063334101/297169890447534156861399350451376090
sage: _.n()
12.3258627936473
sage: (...) integral(iQQ, x, 0, QQ(0.275), algorithm='maxima') # transform by hand gives reduced
86958139/28344976*log(1381/10) - 86958139/28344976*log(5) + 28551/13310
sage: _.n()
12.3258627936473
sage: integral(i, x, 0, 0.275, algorithm='giac')  # ok
12.3258627936
sage: integral(i, x, 0, 0.275, algorithm='sympy')  # ok
12.3258627936473

t/18822

sage: integral(sqrt(x-1)*sqrt(1/x-1/4), x, algorithm='giac')  # unevaluated (root-of stuff)
...
sage: integral(sqrt(x-1)*sqrt(1/x-1/4), x, algorithm='maxima') # unevaluated
integrate(sqrt(x - 1)*sqrt(1/x - 1/4), x)
sage: integrate(sqrt(-1/4*x^2 + 5/4*x - 1)/x^2, x, 1, 2, algorithm='giac')  # ok
-1/4*pi + 1/20*(20*sqrt(2)*asin(1/3) - 50*sqrt(2)*atan(-5/2*sqrt(2) + 3) + 3*sqrt(2) - 30*asin(1/3) + 75*atan(-5/2*sqrt(2) + 3) - 2)/(2*sqrt(2) - 3) + 5/4*atan(1/2) + 3/1
0
sage: N(-1/4*pi + 1/20*(20*sqrt(2)*asin(1/3) - 50*sqrt(2)*atan(-5/2*sqrt(2) + 3) + 3*sqrt(2) - 30*asin(1/3) + 75*atan(-5/2*sqrt(2) + 3) - 2)/(2*sqrt(2) - 3) + 5/4*atan(1/2) + 3/1
0)  # why does _.n() does not work? (TypeError: cannot evaluate symbolic expression numerically)
0.225112673391975
sage: integrate(sqrt(-1/4*x^2 + 5/4*x - 1)/x^2, x, 1, 2, algorithm='maxima')  # ok better
-sqrt(-1/4*x^2 + 5/4*x - 1)/x + 1/2*arcsin(-2/3*x + 5/3) + 5/8*arcsin(5/3*x/abs(x) - 8/3/abs(x))
sage: _.n()
0.225112673391975

t/21057

sage: integrate(1/(x*(2+x)^(1/3)), x, algorithm='giac')  # complicated rootof stuff
NotImplementedError: Unable to parse Giac output: 3*(rootof([[1,-4,-10,19,53,-61],[1,0,-9,-4,27,-36,-23]])/612*ln(((x+2)^(1/3))^2+2^(1/3)*(x+2)^(1/3)+2^(1/3)*2^(1/3))+roo
tof([[4,1,-23,-26,25,-125],[1,0,-9,-4,27,-36,-23]])/306*atan(((x+2)^(1/3)+1/2*2^(1/3))/sqrt(3)*2/2^(1/3))+2^(1/3)*2^(1/3)/6*ln(abs((x+2)^(1/3)-2^(1/3))))
sage: integrate(1/(x*(2+x)^(1/3)), x, algorithm='maxima')
1/2*sqrt(3)*2^(2/3)*arctan(1/6*sqrt(3)*2^(2/3)*(2^(1/3) + 2*(x + 2)^(1/3))) - 1/4*2^(2/3)*log(2^(2/3) + 2^(1/3)*(x + 2)^(1/3) + (x + 2)^(2/3)) + 1/2*2^(2/3)*log(-2^(1/3)
+ (x + 2)^(1/3))

The following returns both exactly computed and numerical values, feature that should be handled separately:

sage: integrate(1/(x*(2+x)^(1/3)), x, 1, 2, algorithm='giac') # wrong -- in console [0, 0.459755569194]
sage: integrate(1/(x*(2+x)^(1/3)), x, 1, 2, algorithm='sympy')  # ok (?)
1/6*2^(2/3)*(-I*sqrt(3) + 1)*log(-1/4*3^(1/3)*2^(2/3)*(I*sqrt(3) - 1) + 1)*gamma(2/3)/gamma(5/3) + 1/6*2^(2/3)*(I*sqrt(3) + 1)*log(-1/4*3^(1/3)*2^(2/3)*(-I*sqrt(3) - 1) + 1)*gamma(2/3)/gamma(5/3) - 1/6*2^(2/3)*(-I*sqrt(3) + 1)*log(-1/2*2^(1/3)*(I*sqrt(3) - 1) + 1)*gamma(2/3)/gamma(5/3) - 1/6*2^(2/3)*(I*sqrt(3) + 1)*log(-1/2*2^(1/3)*(-I*sqrt(3) - 1) + 1)*gamma(2/3)/gamma(5/3) - 1/3*2^(2/3)*log(1/2*3^(1/3)*2^(2/3)*exp_polar(-I*pi) + 1)*gamma(2/3)/gamma(5/3) + 1/3*2^(2/3)*log(2^(1/3)*exp_polar(-I*pi) + 1)*gamma(2/3)/gamma(5/3)
sage: _.n()  # cannot numerically evaluate exp_polar(-I*pi)
...

t/22138

sage: n = var('n'); assume(n, 'integer')
sage: integral(x^n*sin(x), x, algorithm='maxima') # wrong
1/4*(((-1)^n - 1)*gamma(n + 1, I*x) - ((-1)^n - 1)*gamma(n + 1, -I*x))*(-1)^(-1/2*n)
sage: integral(x^n*sin(x), x, algorithm='giac')  # don't know
...
NotImplementedError: Unable to parse Giac output: integrate(x^n*sin(x),x)
sage: integral(x^n*sin(x), x, algorithm='sympy') # ok
1/2*x^2*x^n*hypergeometric((1/2*n + 1,), (3/2, 1/2*n + 2), -1/4*x^2)*gamma(1/2*n + 1)/gamma(1/2*n + 2)

t/13097

sage: integrate(ln(1+4/5*sin(x)), x, -3.1415, 3.1415, algorithm='sympy')  # timeout
# 
sage: integrate(ln(1+4/5*sin(x)), x, -3.1415, 3.1415, algorithm='maxima') # don't know, long
integrate(log(4/5*sin(x) + 1), x, -3.1415, 3.1415)
sage: integrate(ln(1+4/5*sin(x)), x, -3.1415, 3.1415, algorithm='giac') # ok, fast
-1.40205228301
sage: numerical_integral(log(4/5*sin(x) + 1), -3.1415, 3.1415)
(-1.4020522830091542, 5.270707787524382e-08)

t/21440

sage: f1 = (2*cos(2*pi*x) - cos(4*pi*x)) / (5 - 4*cos(2*pi*x))
sage: integrate(f1, x, 0, 1, algorithm='maxima')  # wrong
23/24
sage: integrate(f1, x, 0, 1, algorithm='sympy')  # don't know (long) 
-integrate(-cos(4*pi*x)/(4*cos(2*pi*x) - 5), x, 0, 1) - integrate(2*cos(2*pi*x)/(4*cos(2*pi*x) - 5), x, 0, 1)
sage: integrate(f1, x, 0, 1, algorithm='giac')  # ok, fast
1/4
sage: numerical_integral(f1, 0, 1)
(0.24999999999999997, 4.616007731122122e-15)

t/17511

sage: integrate(abs(cos(x)),x,0,pi, algorithm='maxima') # wrong
-1
sage: integrate(abs(cos(x)),x,0,pi, algorithm='sympy') # don't know
integrate(abs(cos(x)), x, 0, pi)
sage: integrate(abs(cos(x)),x,0,pi, algorithm='giac') # ok
2

t/11164

sage: integrate(sin(x)/x, x, -1e-6, 1e-6, algorithm='maxima') # wrong
-I*Ei(1/1000000*I) + I*Ei(-1/1000000*I)
sage: _.n()
3.14159465358979
sage: integrate(sin(x)/x, x, -1e-6, 1e-6, algorithm='sympy') # ok
1.99999999999989e-6
sage: integrate(sin(x)/x, x, -1e-6, 1e-6, algorithm='giac') # ok
2e-06
sage: numerical_integral(sin(x)/x, -1e-6, 1e-6)  # nan (?)
(1.999999999999889e-06, nan)

t/18599

sage: integrate(sqrt(cot(x)^2),x,pi/4,pi/4*3, algorithm='maxima') # wrong
0
sage: integrate(sqrt(cot(x)^2),x,pi/4,pi/4*3, algorithm='sympy') # don't know
integrate(sqrt(cot(x)^2), x, 1/4*pi, 3/4*pi)
sage: integrate(sqrt(cot(x)^2),x,pi/4,pi/4*3, algorithm='giac') # ok
ln(2)
sage: numerical_integral(sqrt(cot(x)^2),pi/4,pi/4*3)
(0.6931471805599454, 7.702172233337023e-15)

t/13718

sage: integrate((2/3)*x^(5/2)*sqrt(x+1), x, 0, 1, algorithm='sympy') # timeout ~10min
#
sage: integrate((2/3)*x^(5/2)*sqrt(x+1), x,0,1, algorithm='maxima') # ok, fast
-5/192*I*pi + 61/288*sqrt(2) - 5/192*log(sqrt(2) + 1) + 5/192*log(-sqrt(2) + 1)
sage: _.n()
0.253633414928700
sage: integrate((2/3)*x^(5/2)*sqrt(x+1), x, 0, 1, algorithm='giac')  # ok, fast, reduced
61/288*sqrt(2) + 5/96*ln(sqrt(2) - 1)
sage: (61/288*sqrt(2) + 5/96*ln(sqrt(2) - 1)).n()
0.253633414928700
sage: numerical_integral((2/3)*x^(5/2)*sqrt(x+1), 0, 1) 
(0.25363341492869607, 5.446801299750804e-14)

but:

sage: integral((2/3)*x^(5/2)*(x+1)^.5, x,0,1, algorithm='giac') # ok
0.253633414929
sage: integral((2/3)*x^(5/2)*(x+1)^.5, x,0,1, algorithm='maxima') # wrong, misuse (*)
-0.8888888888888888*sqrt(2)
sage: integral((2/3)*x^(5/2)*(x+1)^.5, x,0,1, algorithm='sympy') # timeout

(*) shouldn't use floating point numbers for symbolic integration -- use numerical_integral instead:

sage: integral((2/3)*x^(5/2)*(x+1)^(1/2), x,0,1, algorithm='maxima') # ok
-5/192*I*pi + 61/288*sqrt(2) - 5/192*log(sqrt(2) + 1) + 5/192*log(-sqrt(2) + 1)
sage: _.n()
0.253633414928700
sage: numerical_integral((2/3)*x^(5/2)*(x+1)^.5, 0,1)
(0.25363341492869607, 5.446801299750804e-14)

t/17468

sage: integral(log(abs(2*sin(x))), x, 0, pi/3, algorithm='giac') # don't know
integrate(ln(2*sin(x)), x, 0, 1/3*pi)
sage: integral(log(abs(2*sin(x))), x, 0, pi/3, algorithm='sympy') # bug
...
TypeError: __new__() got an unexpected keyword argument 'manual'
sage: integral(log(abs(2*sin(x))), x, 0, pi/3, algorithm='maxima') # (?)
1/36*I*pi^2 + I*dilog(1/2*I*sqrt(3) + 1/2) + I*dilog(-1/2*I*sqrt(3) - 1/2) - limit(x*log(2*sgn(sin(x))*sin(x)), x, 0, plus)

t/11493

sage: integrate(x^(1/2)*e^(-x-x^2),x,0,infinity, algorithm='giac') # don't know
integrate(sqrt(x)*e^(-x^2 - x), x, 0, +Infinity)
I(1/4, 1/8)*e^(1/8)*gamma(5/4))*gamma(3/4))/sqrt(pi)
sage: integrate(x^(1/2)*e^(-x-x^2),x,0,infinity, algorithm='maxima') # don't know
integrate(sqrt(x)*e^(-x^2 - x), x, 0, +Infinity)
sage: integrate(x^(1/2)*e^(-x-x^2),x,0,infinity, algorithm='sympy') # ok, fast
-1/12*(4*sqrt(pi)*(bessel_I(7/4, 1/8)*e^(1/8)*gamma(7/4) + 13*bessel_I(3/4, 1/8)*e^(1/8)*gamma(7/4))*gamma(5/4) - 3*sqrt(pi)*(bessel_I(5/4, 1/8)*e^(1/8)*gamma(5/4) + 5*bessel_I(1/4, 1/8)*e^(1/8)*gamma(5/4))*gamma(3/4))/sqrt(pi)
sage: _.n()
0.320157090360147
sage: numerical_integral(x^(1/2)*e^(-x-x^2),0,infinity)
(0.3201570903601315, 8.176582744212624e-10)

t/F2FzA0CYZw

sage: integrate(heaviside(x-1)*(x-1)+1,(x,0,2), algorithm='sympy')  # wrong
2
sage: integrate(heaviside(x-1)*(x-1)+1,(x,0,2), algorithm='maxima') # unevaluated
integrate((x - 1)*heaviside(x - 1) + 1, x, 0, 2)
sage: integrate(heaviside(x-1)*(x-1)+1,(x,0,2), algorithm='giac') # ok
5/2
sage: numerical_integral(heaviside(x-1)*(x-1)+1,0,2)
Exception TypeError: TypeError('unable to simplify to float approximation',) in 'sage.calculus.integration.c_ff' ignored
(2.5, 2.7755575615628914e-14)

t/12145

sage: integrate(x/(exp(x)-1), x, 0, oo, algorithm='sympy')  # unevaluated
integrate(x/(e^x - 1), x, 0, +Infinity)
sage: integrate(x/(exp(x)-1), x, 0, oo, algorithm='maxima')  # value of the limit?
-1/6*pi^2 + limit(-1/2*x^2 + x*log(-e^x + 1) + dilog(e^x), x, +Infinity, minus)
sage: integrate(x/(exp(x)-1), x, 0, oo, algorithm='giac')  # ok
1/6*pi^2
sage: _.n()
1.64493406684823
sage: numerical_integral(x/(exp(x)-1), 0, oo)
(1.6449340668482264, 5.935645283617803e-10)

t/17328

sage: integral(1/(ln(x)**2), x, 2, 3, algorithm='giac')  # timeout
..
sage: integral(1/(ln(x)**2), x, 2, 3, algorithm='sympy')  # ok
-3/log(3) + 2/log(2) + integrate(1/log(x), x, 2, 3)
sage: _.n()
1.27309721644711
sage: integral(1/(ln(x)**2), x, 2, 3, algorithm='maxima')  # ok
gamma(-1, -log(3)) - gamma(-1, -log(2))
sage: _.n()
1.27309721644711
sage: numerical_integral(1/ln(x)**2, 2, 3)
(1.273097216447114, 1.4134218422857824e-14)

system

sage: version()
'SageMath version 8.0.beta3, Release Date: 2017-04-23'
sage: giac('version()')
"giac 1.2.3, (c) B. Parisse and R. De Graeve, Institut Fourier, Universite de Grenoble I"
$ ./sage -maxima
...
Maxima 5.39.0 http://maxima.sourceforge.net
using Lisp ECL 16.1.2
sage: import sympy
sage: sympy.__version__
'1.0'

scores

maxima giac sympy
2.5 13 7

where: correct -> +1, fast or reduced -> +0.5, wrong -> -1, unevaluated -> 0.

@mforets
Copy link
Author

mforets commented May 8, 2017

task: transform this raw markdown file into a (sage -t) testable .py file.

@mforets
Copy link
Author

mforets commented May 9, 2017

apply time tag to tests, cf. #19519 and discussion

@rwst
Copy link

rwst commented Apr 19, 2018

Rubi-4.14 cannot solve 17468 either, having no rules with abs(). Solution of 12145 is -exp(1-x)-x*exp(1-x).

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