Ask Your Question

Function defined with an integral

asked 2023-11-16 14:26:12 +0100

Sotto gravatar image

updated 2023-11-16 18:00:42 +0100

FrédéricC gravatar image

Dear all, I'm puzzled by this answer from Sage:

sage: def f(x):
....:     return integral(exp(t^2)/sqrt(t),t,1,x^2)
sage: f(4).numerical_approx()
1.18428615262607e109 - 1.28184667602042*I

Why the answer is a complex number? How can I define f as a real function? Thanks.

edit retag flag offensive close merge delete


This just means that approximation works in complex numbers here. Note that the imaginary part is negligible as compared to the real part, and so there is nothing wrong. Since you know that the result is real, you can use f(4).numerical_approx().real() to get the appropriation in real numbers.

Max Alekseyev gravatar imageMax Alekseyev ( 2023-11-16 19:41:28 +0100 )edit

This can be a workaround:

def f(x):
    return real(integral(exp(t^2)/sqrt(t),t,1,x^2))

tolga gravatar imagetolga ( 2023-11-17 04:36:24 +0100 )edit

Also try this:

def f(x):
    return numerical_integral(exp(t^2)/sqrt(t),1,x^2)


As you can see from the result (1.1842861526260861e+109, 4.7486208201754434e+95), the error bound is quite large. You can plot your integrand to see its behavior in the range of your interest.

tolga gravatar imagetolga ( 2023-11-17 04:43:07 +0100 )edit

[ Comment reposted as an answer... ]

Emmanuel Charpentier gravatar imageEmmanuel Charpentier ( 2023-11-17 09:53:36 +0100 )edit

Thank you. So nothing wrong with this function… But why I receive an error message if I try to differentiate it?!? I understand that this is not the original question though… :-)

Sotto gravatar imageSotto ( 2023-11-18 14:15:24 +0100 )edit

1 Answer

Sort by » oldest newest most voted

answered 2023-11-17 10:07:53 +0100

Emmanuel Charpentier gravatar image

updated 2023-11-19 15:28:11 +0100

Complement to Max's and Tolga's comments :

sage: var("t")
sage: with assuming(x>1): (e^(t^2)/sqrt(t)).integrate(t, 1, x^2)
-1/2*(-1)^(3/4)*gamma(1/4, -1) + 1/2*(-x^4)^(3/4)*gamma(1/4, -x^4)/x^3

$$ -\frac{1}{2} \, \left(-1\right)^{\frac{3}{4}} \Gamma\left(\frac{1}{4}, -1\right) + \frac{\left(-x^{4}\right)^{\frac{3}{4}} \Gamma\left(\frac{1}{4}, -x^{4}\right)}{2 \, x^{3}} $$

As you can see, this expression involves different inherently non-real complex factors (fractional powers of negative numbers, $\Gamma$ with negative arguments), whose sums happens to be real. The computation of the numerical value is inexact, so a small inexactitude is enough to result in a numerically non-zero imaginary part.

Showing that this imaginary part is zero seems out of reach of current CASses...

EDIT : To illistrate Max's comment about derivation :

sage: with assuming(x>1): foo=(e^(t^2)/sqrt(t)).integrate(t, 1, x^2) ; foo
-1/2*(-1)^(3/4)*gamma(1/4, -1) + 1/2*(-x^4)^(3/4)*gamma(1/4, -x^4)/x^3
sage: foo.diff(x)
-2*(-x^4)^(3/4)*D[1](gamma)(1/4, -x^4) - 3/2*gamma(1/4, -x^4)/(-x^4)^(1/4) - 3/2*(-x^4)^(3/4)*gamma(1/4, -x^4)/x^4

$$ -2 \, \left(-x^{4}\right)^{\frac{3}{4}} \mathrm{D}_{1}\left(\Gamma\right)\left(\frac{1}{4}, -x^{4}\right) - \frac{3 \, \Gamma\left(\frac{1}{4}, -x^{4}\right)}{2 \, \left(-x^{4}\right)^{\frac{1}{4}}} - \frac{3 \, \left(-x^{4}\right)^{\frac{3}{4}} \Gamma\left(\frac{1}{4}, -x^{4}\right)}{2 \, x^{4}} $$

One notes that Sage does not explicites a derivative of $\Gamma$ wrt one of its arguments...

Suspicious :

sage: giac.diff(*map(giac, (foo, x)))._sage_()
-3/2*gamma(1/4, -x^4)/(-x^4)^(1/4) - 3/2*(-x^4)^(3/4)*gamma(1/4, -x^4)/x^4 + 2*e^(x^4)

$$ -\frac{3 \, \Gamma\left(\frac{1}{4}, -x^{4}\right)}{2 \, \left(-x^{4}\right)^{\frac{1}{4}}} - \frac{3 \, \left(-x^{4}\right)^{\frac{3}{4}} \Gamma\left(\frac{1}{4}, -x^{4}\right)}{2 \, x^{4}} + 2 \, e^{\left(x^{4}\right)} $$

Highly suspicious :

sage: fricas.D(*map(fricas, (foo, x)))._sage_()

and, BTW :

sage: foo._sympy_().diff(x._sympy_())._sage_()
sage: foo._mathematica_().D(x).FullSimplify().sage()

But these expressions tun out to be equal :

sage: (fricas.D(*map(fricas, (foo, x)))._sage_()-giac.diff(*map(giac, (foo, x)))._sage_())
3/2*gamma(1/4, -x^4)/(-x^4)^(1/4) + 3/2*(-x^4)^(3/4)*gamma(1/4, -x^4)/x^4
sage: (fricas.D(*map(fricas, (foo, x)))._sage_()-giac.diff(*map(giac, (foo, x)))._sage_()).factor()

and comparing the Fricas result to a (rough) numerical estimate of the derivative of the function shows a very good agreement : their quotient is 1 up to $10^{-6}$ on the range (1, 4) :

sage: plot(lambda u:((Foo(u+1e-8)-Foo(u))/1e-8).real()/(2*e^(u^4))-1, (1, 4))
Launched png viewer for Graphics object consisting of 1 graphics primitive

image description

Mathematica's, Fricas' and Giac's derivations seem to be exact. The Sage's $\Gamma$ implementation may be perfected...


edit flag offensive delete link more

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools


Asked: 2023-11-16 14:26:12 +0100

Seen: 75 times

Last updated: Nov 19