| 1 | r""" |
| 2 | Exponential Integrals |
| 3 | |
| 4 | AUTHORS: |
| 5 | |
| 6 | - Benjamin Jones (2011-06-12) |
| 7 | |
| 8 | This module provides easy access to many exponential integral |
| 9 | special functions. It utilizes Maxima's `special functions package`_ and |
| 10 | the `mpmath library`_. |
| 11 | |
| 12 | REFERENCES: |
| 13 | |
| 14 | - [AS]_ Abramowitz and Stegun: *Handbook of Mathematical Functions* |
| 15 | - Wikipedia Entry: http://en.wikipedia.org/wiki/Exponential_integral |
| 16 | - Online Encyclopedia of Special Function: http://algo.inria.fr/esf/index.html |
| 17 | - NIST Digital Library of Mathematical Functions: http://dlmf.nist.gov/ |
| 18 | - Maxima `special functions package`_ |
| 19 | - `mpmath library`_ |
| 20 | |
| 21 | .. [AS] 'Handbook of Mathematical Functions', Milton Abramowitz and Irene |
| 22 | A. Stegun, National Bureau of Standards Applied Mathematics Series, 55. |
| 23 | See also http://www.math.sfu.ca/~cbm/aands/. |
| 24 | .. _`special functions package`: http://maxima.sourceforge.net/docs/manual/en/maxima_15.html |
| 25 | .. _`mpmath library`: http://code.google.com/p/mpmath/ |
| 26 | |
| 27 | AUTHORS: |
| 28 | |
| 29 | - Benjamin Jones |
| 30 | |
| 31 | Implementations of the classes ``Function_exp_integral_*``. |
| 32 | |
| 33 | - David Joyner and William Stein |
| 34 | |
| 35 | Authors of the code which was moved from special.py and trans.py. |
| 36 | Implementation of :meth:`exp_int` (from sage/functions/special.py). |
| 37 | Implementation of :meth:`exponential_integral_1` (from |
| 38 | sage/functions/transcendental.py). |
| 39 | |
| 40 | """ |
| 41 | |
| 42 | #***************************************************************************** |
| 43 | # Copyright (C) 2011 Benjamin Jones <benjaminfjones@gmail.com> |
| 44 | # |
| 45 | # Distributed under the terms of the GNU General Public License (GPL) |
| 46 | # |
| 47 | # This code is distributed in the hope that it will be useful, |
| 48 | # but WITHOUT ANY WARRANTY; without even the implied warranty of |
| 49 | # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
| 50 | # General Public License for more details. |
| 51 | # |
| 52 | # The full text of the GPL is available at: |
| 53 | # |
| 54 | # http://www.gnu.org/licenses/ |
| 55 | #***************************************************************************** |
| 56 | |
| 57 | import sage.interfaces.all |
| 58 | from sage.misc.sage_eval import sage_eval |
| 59 | from sage.symbolic.function import BuiltinFunction, is_inexact |
| 60 | from sage.calculus.calculus import maxima |
| 61 | from sage.symbolic.expression import Expression |
| 62 | from sage.structure.parent import Parent |
| 63 | from sage.structure.coerce import parent |
| 64 | from sage.libs.mpmath import utils as mpmath_utils |
| 65 | mpmath_utils_call = mpmath_utils.call # eliminate some overhead in _evalf_ |
| 66 | |
| 67 | from sage.rings.rational_field import RationalField |
| 68 | from sage.rings.real_mpfr import RealField |
| 69 | from sage.rings.complex_field import ComplexField |
| 70 | from sage.rings.all import ZZ, QQ, RR, RDF |
| 71 | from sage.functions.log import exp, log |
| 72 | from sage.functions.trig import sin, cos |
| 73 | from sage.functions.hyperbolic import sinh, cosh |
| 74 | |
| 75 | |
| 76 | class Function_exp_integral_e(BuiltinFunction): |
| 77 | r""" |
| 78 | The generalized complex exponential integral `E_n(z)` defined by |
| 79 | |
| 80 | .. math:: |
| 81 | |
| 82 | \operatorname{E_n}(z) = \int_1^{\infty} \frac{e^{-z t}}{t^n} \; dt |
| 83 | |
| 84 | for complex numbers `n` and `z`, see [AS]_ 5.1.4. |
| 85 | |
| 86 | The special case where `n = 1` is denoted in Sage by |
| 87 | ``exp_integral_e1``. |
| 88 | |
| 89 | EXAMPLES: |
| 90 | |
| 91 | Numerical evaluation is handled using mpmath:: |
| 92 | |
| 93 | sage: N(exp_integral_e(1,1)) |
| 94 | 0.219383934395520 |
| 95 | sage: exp_integral_e(1, RealField(100)(1)) |
| 96 | 0.21938393439552027367716377546 |
| 97 | |
| 98 | We can compare this to PARI's evaluation of |
| 99 | :meth:`exponential_integral_1`:: |
| 100 | |
| 101 | sage: N(exponential_integral_1(1)) |
| 102 | 0.219383934395520 |
| 103 | |
| 104 | We can verify one case of [AS]_ 5.1.45, i.e. |
| 105 | `E_n(z) = z^{n-1}\Gamma(1-n,z)`:: |
| 106 | |
| 107 | sage: N(exp_integral_e(2, 3+I)) |
| 108 | 0.00354575823814662 - 0.00973200528288687*I |
| 109 | sage: N((3+I)*gamma(-1, 3+I)) |
| 110 | 0.00354575823814662 - 0.00973200528288687*I |
| 111 | |
| 112 | Maxima returns the following improper integral as a multiple of |
| 113 | ``exp_integral_e(1,1)``:: |
| 114 | |
| 115 | sage: uu = integral(e^(-x)*log(x+1),x,0,oo) |
| 116 | sage: uu |
| 117 | e*exp_integral_e(1, 1) |
| 118 | sage: uu.n(digits=30) |
| 119 | 0.596347362323194074341078499369 |
| 120 | |
| 121 | Symbolic derivatives and integrals are handled by Sage and Maxima:: |
| 122 | |
| 123 | sage: x = var('x') |
| 124 | sage: f = exp_integral_e(2,x) |
| 125 | sage: f.diff(x) |
| 126 | -exp_integral_e(1, x) |
| 127 | |
| 128 | sage: f.integrate(x) |
| 129 | -exp_integral_e(3, x) |
| 130 | |
| 131 | sage: f = exp_integral_e(-1,x) |
| 132 | sage: f.integrate(x) |
| 133 | Ei(-x) - gamma(-1, x) |
| 134 | |
| 135 | Some special values of ``exp_integral_e`` can be simplified. |
| 136 | [AS]_ 5.1.23:: |
| 137 | |
| 138 | sage: exp_integral_e(0,x) |
| 139 | e^(-x)/x |
| 140 | |
| 141 | [AS]_ 5.1.24:: |
| 142 | |
| 143 | sage: exp_integral_e(6,0) |
| 144 | 1/5 |
| 145 | sage: nn = var('nn') |
| 146 | sage: assume(nn > 1) |
| 147 | sage: f = exp_integral_e(nn,0) |
| 148 | sage: f.simplify() |
| 149 | 1/(nn - 1) |
| 150 | |
| 151 | |
| 152 | ALGORITHM: |
| 153 | |
| 154 | Numerical evaluation is handled using mpmath, but symbolics are handled |
| 155 | by Sage and Maxima. |
| 156 | |
| 157 | """ |
| 158 | def __init__(self): |
| 159 | """ |
| 160 | See the docstring for :meth:`Function_exp_integral_e`. |
| 161 | |
| 162 | EXAMPLES:: |
| 163 | |
| 164 | sage: exp_integral_e(1,0) |
| 165 | exp_integral_e(1, 0) |
| 166 | |
| 167 | """ |
| 168 | BuiltinFunction.__init__(self, "exp_integral_e", nargs=2, |
| 169 | latex_name=r'exp_integral_e', |
| 170 | conversions=dict(maxima='expintegral_e')) |
| 171 | |
| 172 | def _eval_(self, n, z): |
| 173 | """ |
| 174 | EXAMPLES:: |
| 175 | |
| 176 | sage: exp_integral_e(1.0, x) |
| 177 | exp_integral_e(1.00000000000000, x) |
| 178 | sage: exp_integral_e(x, 1.0) |
| 179 | exp_integral_e(x, 1.00000000000000) |
| 180 | sage: exp_integral_e(1.0, 1.0) |
| 181 | 0.219383934395520 |
| 182 | |
| 183 | """ |
| 184 | if not isinstance(n, Expression) and not isinstance(z, Expression) and \ |
| 185 | (is_inexact(n) or is_inexact(z)): |
| 186 | coercion_model = sage.structure.element.get_coercion_model() |
| 187 | n, z = coercion_model.canonical_coercion(n, z) |
| 188 | return self._evalf_(n, z, parent(n)) |
| 189 | |
| 190 | z_zero = False |
| 191 | # special case: z == 0 and n > 1 |
| 192 | if isinstance(z, Expression): |
| 193 | if z.is_trivial_zero(): |
| 194 | z_zero = True # for later |
| 195 | if n > 1: |
| 196 | return 1/(n-1) |
| 197 | else: |
| 198 | if not z: |
| 199 | z_zero = True |
| 200 | if n > 1: |
| 201 | return 1/(n-1) |
| 202 | |
| 203 | # special case: n == 0 |
| 204 | if isinstance(n, Expression): |
| 205 | if n.is_trivial_zero(): |
| 206 | if z_zero: |
| 207 | return None |
| 208 | else: |
| 209 | return exp(-z)/z |
| 210 | else: |
| 211 | if not n: |
| 212 | if z_zero: |
| 213 | return None |
| 214 | else: |
| 215 | return exp(-z)/z |
| 216 | |
| 217 | return None # leaves the expression unevaluated |
| 218 | |
| 219 | def _evalf_(self, n, z, parent=None): |
| 220 | """ |
| 221 | EXAMPLES:: |
| 222 | |
| 223 | sage: N(exp_integral_e(1, 1+I)) |
| 224 | 0.000281624451981418 - 0.179324535039359*I |
| 225 | sage: exp_integral_e(1, RealField(100)(1)) |
| 226 | 0.21938393439552027367716377546 |
| 227 | |
| 228 | """ |
| 229 | import mpmath |
| 230 | return mpmath_utils.call(mpmath.expint, n, z, parent=parent) |
| 231 | |
| 232 | def _derivative_(self, n, z, diff_param=None): |
| 233 | """ |
| 234 | If `n` is an integer strictly larger than 0, then the derivative of |
| 235 | `E_n(z)` with respect to `z` is |
| 236 | `-E_{n-1}(z)`. See [AS]_ 5.1.26. |
| 237 | |
| 238 | EXAMPLES:: |
| 239 | |
| 240 | sage: x = var('x') |
| 241 | sage: f = exp_integral_e(2,x) |
| 242 | sage: f.diff(x) |
| 243 | -exp_integral_e(1, x) |
| 244 | |
| 245 | sage: f = exp_integral_e(2,sqrt(x)) |
| 246 | sage: f.diff(x) |
| 247 | -1/2*exp_integral_e(1, sqrt(x))/sqrt(x) |
| 248 | |
| 249 | """ |
| 250 | if n in ZZ and n > 0: |
| 251 | return -1*exp_integral_e(n-1,z) |
| 252 | else: |
| 253 | raise NotImplementedError("The derivative of this function is only implemented for n = 1, 2, 3, ...") |
| 254 | |
| 255 | exp_integral_e = Function_exp_integral_e() |
| 256 | |
| 257 | |
| 258 | class Function_exp_integral_e1(BuiltinFunction): |
| 259 | r""" |
| 260 | The generalized complex exponential integral `E_1(z)` defined by |
| 261 | |
| 262 | .. math:: |
| 263 | |
| 264 | \operatorname{E_1}(z) = \int_z^\infty \frac{e^{-t}}{t}\; dt |
| 265 | |
| 266 | see [AS]_ 5.1.4. |
| 267 | |
| 268 | EXAMPLES: |
| 269 | |
| 270 | Numerical evaluation is handled using mpmath:: |
| 271 | |
| 272 | sage: N(exp_integral_e1(1)) |
| 273 | 0.219383934395520 |
| 274 | sage: exp_integral_e1(RealField(100)(1)) |
| 275 | 0.21938393439552027367716377546 |
| 276 | |
| 277 | We can compare this to PARI's evaluation of |
| 278 | :meth:`exponential_integral_1`:: |
| 279 | |
| 280 | sage: N(exp_integral_e1(2.0)) |
| 281 | 0.0489005107080611 |
| 282 | sage: N(exponential_integral_1(2.0)) |
| 283 | 0.0489005107080611 |
| 284 | |
| 285 | Symbolic derivatives and integrals are handled by Sage and Maxima:: |
| 286 | |
| 287 | sage: x = var('x') |
| 288 | sage: f = exp_integral_e1(x) |
| 289 | sage: f.diff(x) |
| 290 | -e^(-x)/x |
| 291 | |
| 292 | sage: f.integrate(x) |
| 293 | -exp_integral_e(2, x) |
| 294 | |
| 295 | ALGORITHM: |
| 296 | |
| 297 | Numerical evaluation is handled using mpmath, but symbolics are handled |
| 298 | by Sage and Maxima. |
| 299 | |
| 300 | """ |
| 301 | def __init__(self): |
| 302 | """ |
| 303 | See the docstring for :class:`Function_exp_integral_e1`. |
| 304 | |
| 305 | EXAMPLES:: |
| 306 | |
| 307 | sage: exp_integral_e1(1) |
| 308 | exp_integral_e1(1) |
| 309 | |
| 310 | """ |
| 311 | BuiltinFunction.__init__(self, "exp_integral_e1", nargs=1, |
| 312 | latex_name=r'exp_integral_e1', |
| 313 | conversions=dict(maxima='expintegral_e1')) |
| 314 | |
| 315 | def _eval_(self, z): |
| 316 | """ |
| 317 | EXAMPLES:: |
| 318 | |
| 319 | sage: exp_integral_e1(x) |
| 320 | exp_integral_e1(x) |
| 321 | sage: exp_integral_e1(1.0) |
| 322 | 0.219383934395520 |
| 323 | |
| 324 | """ |
| 325 | if not isinstance(z, Expression) and is_inexact(z): |
| 326 | return self._evalf_(z, parent(z)) |
| 327 | |
| 328 | return None # leaves the expression unevaluated |
| 329 | |
| 330 | def _evalf_(self, z, parent=None): |
| 331 | """ |
| 332 | EXAMPLES:: |
| 333 | |
| 334 | sage: N(exp_integral_e1(1+I)) |
| 335 | 0.000281624451981418 - 0.179324535039359*I |
| 336 | sage: exp_integral_e1(RealField(200)(0.5)) |
| 337 | 0.55977359477616081174679593931508523522684689031635351524829 |
| 338 | |
| 339 | """ |
| 340 | import mpmath |
| 341 | return mpmath_utils_call(mpmath.e1, z, parent=parent) |
| 342 | |
| 343 | def _derivative_(self, z, diff_param=None): |
| 344 | """ |
| 345 | The derivative of `E_1(z)` is `-e^{-z}/z`. See [AS], 5.1.26. |
| 346 | |
| 347 | EXAMPLES:: |
| 348 | |
| 349 | sage: x = var('x') |
| 350 | sage: f = exp_integral_e1(x) |
| 351 | sage: f.diff(x) |
| 352 | -e^(-x)/x |
| 353 | |
| 354 | sage: f = exp_integral_e1(x^2) |
| 355 | sage: f.diff(x) |
| 356 | -2*e^(-x^2)/x |
| 357 | |
| 358 | """ |
| 359 | return -exp(-z)/z |
| 360 | |
| 361 | exp_integral_e1 = Function_exp_integral_e1() |
| 362 | |
| 363 | |
| 364 | class Function_log_integral(BuiltinFunction): |
| 365 | r""" |
| 366 | The logarithmic integral `\operatorname{li}(z)` defined by |
| 367 | |
| 368 | .. math:: |
| 369 | |
| 370 | \operatorname{li}(x) = \int_0^z \frac{dt}{\ln(t)} = \operatorname{Ei}(\ln(x)) |
| 371 | |
| 372 | for x > 1 and by analytic continuation for complex arguments z (see [AS]_ 5.1.3). |
| 373 | |
| 374 | EXAMPLES: |
| 375 | |
| 376 | Numerical evaluation for real and complex arguments is handled using mpmath:: |
| 377 | |
| 378 | sage: N(log_integral(3)) |
| 379 | 2.16358859466719 |
| 380 | sage: N(log_integral(3), digits=30) |
| 381 | 2.16358859466719197287692236735 |
| 382 | sage: log_integral(ComplexField(100)(3+I)) |
| 383 | 2.2879892769816826157078450911 + 0.87232935488528370139883806779*I |
| 384 | |
| 385 | Symbolic derivatives and integrals are handled by Sage and Maxima:: |
| 386 | |
| 387 | sage: x = var('x') |
| 388 | sage: f = log_integral(x) |
| 389 | sage: f.diff(x) |
| 390 | 1/log(x) |
| 391 | |
| 392 | sage: f.integrate(x) |
| 393 | x*log_integral(x) - expintegral_ei(2*log(x)) |
| 394 | |
| 395 | Here is a test from the mpmath documentation. There are |
| 396 | 1,925,320,391,606,803,968,923 many prime numbers less than 1e23. The |
| 397 | value of ``log_integral(1e23)`` is very close to this:: |
| 398 | |
| 399 | sage: log_integral(1e23) |
| 400 | 1.92532039161405e21 |
| 401 | |
| 402 | ALGORITHM: |
| 403 | |
| 404 | Numerical evaluation is handled using mpmath, but symbolics are handled |
| 405 | by Sage and Maxima. |
| 406 | |
| 407 | REFERENCES: |
| 408 | |
| 409 | - http://en.wikipedia.org/wiki/Logarithmic_integral_function |
| 410 | - mpmath documentation: `logarithmic-integral`_ |
| 411 | |
| 412 | .. _`logarithmic-integral`: http://mpmath.googlecode.com/svn/trunk/doc/build/functions/expintegrals.html#logarithmic-integral |
| 413 | |
| 414 | |
| 415 | """ |
| 416 | def __init__(self): |
| 417 | """ |
| 418 | See the docstring for ``Function_log_integral``. |
| 419 | |
| 420 | EXAMPLES:: |
| 421 | |
| 422 | sage: log_integral(3) |
| 423 | log_integral(3) |
| 424 | |
| 425 | """ |
| 426 | BuiltinFunction.__init__(self, "log_integral", nargs=1, |
| 427 | latex_name=r'log_integral', |
| 428 | conversions=dict(maxima='expintegral_li')) |
| 429 | |
| 430 | def _eval_(self, z): |
| 431 | """ |
| 432 | EXAMPLES:: |
| 433 | |
| 434 | sage: z = var('z') |
| 435 | sage: log_integral(z) |
| 436 | log_integral(z) |
| 437 | sage: log_integral(3.0) |
| 438 | 2.16358859466719 |
| 439 | |
| 440 | """ |
| 441 | if not isinstance(z, Expression) and is_inexact(z): |
| 442 | return self._evalf_(z, parent(z)) |
| 443 | |
| 444 | return None # leaves the expression unevaluated |
| 445 | |
| 446 | def _evalf_(self, z, parent=None): |
| 447 | """ |
| 448 | EXAMPLES:: |
| 449 | |
| 450 | sage: N(log_integral(1e6)) |
| 451 | 78627.5491594622 |
| 452 | sage: log_integral(RealField(200)(1e6)) |
| 453 | 78627.549159462181919862910747947261161321874382421767074759 |
| 454 | |
| 455 | """ |
| 456 | import mpmath |
| 457 | return mpmath_utils_call(mpmath.li, z, parent=parent) |
| 458 | |
| 459 | def _derivative_(self, z, diff_param=None): |
| 460 | """ |
| 461 | The derivative of `\operatorname{li}(z) is `1/log(z)`. |
| 462 | |
| 463 | EXAMPLES:: |
| 464 | |
| 465 | sage: x = var('x') |
| 466 | sage: f = log_integral(x) |
| 467 | sage: f.diff(x) |
| 468 | 1/log(x) |
| 469 | |
| 470 | sage: f = log_integral(x^2) |
| 471 | sage: f.diff(x) |
| 472 | 2*x/log(x^2) |
| 473 | |
| 474 | """ |
| 475 | return 1/log(z) |
| 476 | |
| 477 | li = log_integral = Function_log_integral() |
| 478 | |
| 479 | |
| 480 | class Function_sin_integral(BuiltinFunction): |
| 481 | r""" |
| 482 | The trigonometric integral `\operatorname{Si}(z)` defined by |
| 483 | |
| 484 | .. math:: |
| 485 | |
| 486 | \operatorname{Si}(z) = \int_0^z \frac{\sin(t)}{t}\; dt, |
| 487 | |
| 488 | see [AS]_ 5.2.1. |
| 489 | |
| 490 | EXAMPLES: |
| 491 | |
| 492 | Numerical evaluation for real and complex arguments is handled using mpmath:: |
| 493 | |
| 494 | sage: sin_integral(0) |
| 495 | 0 |
| 496 | sage: sin_integral(0.0) |
| 497 | 0.000000000000000 |
| 498 | sage: sin_integral(3.0) |
| 499 | 1.84865252799947 |
| 500 | sage: N(sin_integral(3), digits=30) |
| 501 | 1.84865252799946825639773025111 |
| 502 | sage: sin_integral(ComplexField(100)(3+I)) |
| 503 | 2.0277151656451253616038525998 + 0.015210926166954211913653130271*I |
| 504 | |
| 505 | The limit of `\operatorname{Si}(z)` as `z \to \infty` is `\pi/2`:: |
| 506 | |
| 507 | sage: N(sin_integral(1e23)) |
| 508 | 1.57079632679490 |
| 509 | sage: N(pi/2) |
| 510 | 1.57079632679490 |
| 511 | |
| 512 | At 200 bits of precision `\operatorname{Si}(10^{23})` agrees with `\pi/2` up to |
| 513 | `10^{-24}`:: |
| 514 | |
| 515 | sage: sin_integral(RealField(200)(1e23)) |
| 516 | 1.5707963267948966192313288218697837425815368604836679189519 |
| 517 | sage: N(pi/2, prec=200) |
| 518 | 1.5707963267948966192313216916397514420985846996875529104875 |
| 519 | |
| 520 | The exponential sine integral is analytic everywhere:: |
| 521 | |
| 522 | sage: sin_integral(-1.0) |
| 523 | -0.946083070367183 |
| 524 | sage: sin_integral(-2.0) |
| 525 | -1.60541297680269 |
| 526 | sage: sin_integral(-1e23) |
| 527 | -1.57079632679490 |
| 528 | |
| 529 | Symbolic derivatives and integrals are handled by Sage and Maxima:: |
| 530 | |
| 531 | sage: x = var('x') |
| 532 | sage: f = sin_integral(x) |
| 533 | sage: f.diff(x) |
| 534 | sin(x)/x |
| 535 | |
| 536 | sage: f.integrate(x) |
| 537 | x*sin_integral(x) + cos(x) |
| 538 | |
| 539 | sage: integrate(sin(x)/x, x) |
| 540 | 1/2*I*Ei(-I*x) - 1/2*I*Ei(I*x) |
| 541 | |
| 542 | Compare values of the functions `\operatorname{Si}(x)` and |
| 543 | `f(x) = (1/2)i \cdot \operatorname{Ei}(-ix) - (1/2)i \cdot |
| 544 | \operatorname{Ei}(ix) - \pi/2`, which are both anti-derivatives of |
| 545 | `\sin(x)/x`, at some random positive real numbers:: |
| 546 | |
| 547 | sage: f(x) = 1/2*I*Ei(-I*x) - 1/2*I*Ei(I*x) - pi/2 |
| 548 | sage: g(x) = sin_integral(x) |
| 549 | sage: R = [ abs(RDF.random_element()) for i in range(100) ] |
| 550 | sage: all(abs(f(x) - g(x)) < 1e-10 for x in R) |
| 551 | True |
| 552 | |
| 553 | The Nielsen spiral is the parametric plot of (Si(t), Ci(t)):: |
| 554 | |
| 555 | sage: x=var('x') |
| 556 | sage: f(x) = sin_integral(x) |
| 557 | sage: g(x) = cos_integral(x) |
| 558 | sage: P = parametric_plot([f, g], (x, 0.5 ,20)) |
| 559 | sage: show(P, frame=True, axes=False) |
| 560 | |
| 561 | ALGORITHM: |
| 562 | |
| 563 | Numerical evaluation is handled using mpmath, but symbolics are handled |
| 564 | by Sage and Maxima. |
| 565 | |
| 566 | REFERENCES: |
| 567 | |
| 568 | - http://en.wikipedia.org/wiki/Trigonometric_integral |
| 569 | - mpmath documentation: `si`_ |
| 570 | |
| 571 | .. _`si`: http://mpmath.googlecode.com/svn/trunk/doc/build/functions/expintegrals.html#si |
| 572 | |
| 573 | """ |
| 574 | def __init__(self): |
| 575 | """ |
| 576 | See the docstring for ``Function_sin_integral``. |
| 577 | |
| 578 | EXAMPLES:: |
| 579 | |
| 580 | sage: sin_integral(1) |
| 581 | sin_integral(1) |
| 582 | |
| 583 | """ |
| 584 | BuiltinFunction.__init__(self, "sin_integral", nargs=1, |
| 585 | latex_name=r'\operatorname{Si}', |
| 586 | conversions=dict(maxima='expintegral_si')) |
| 587 | |
| 588 | def _eval_(self, z): |
| 589 | """ |
| 590 | EXAMPLES:: |
| 591 | |
| 592 | sage: z = var('z') |
| 593 | sage: sin_integral(z) |
| 594 | sin_integral(z) |
| 595 | sage: sin_integral(3.0) |
| 596 | 1.84865252799947 |
| 597 | sage: sin_integral(0) |
| 598 | 0 |
| 599 | |
| 600 | """ |
| 601 | if not isinstance(z, Expression) and is_inexact(z): |
| 602 | return self._evalf_(z, parent(z)) |
| 603 | |
| 604 | # special case: z = 0 |
| 605 | if isinstance(z, Expression): |
| 606 | if z.is_trivial_zero(): |
| 607 | return z |
| 608 | else: |
| 609 | if not z: |
| 610 | return z |
| 611 | |
| 612 | return None # leaves the expression unevaluated |
| 613 | |
| 614 | def _evalf_(self, z, parent=None): |
| 615 | """ |
| 616 | EXAMPLES: |
| 617 | |
| 618 | The limit `\operatorname{Si}(z)` as `z \to \infty` is `\pi/2`:: |
| 619 | |
| 620 | sage: N(sin_integral(1e23) - pi/2) |
| 621 | 0.000000000000000 |
| 622 | |
| 623 | At 200 bits of precision `\operatorname{Si}(10^{23})` agrees with `\pi/2` up to |
| 624 | `10^{-24}`:: |
| 625 | |
| 626 | sage: sin_integral(RealField(200)(1e23)) |
| 627 | 1.5707963267948966192313288218697837425815368604836679189519 |
| 628 | sage: N(pi/2, prec=200) |
| 629 | 1.5707963267948966192313216916397514420985846996875529104875 |
| 630 | |
| 631 | The exponential sine integral is analytic everywhere, even on the |
| 632 | negative real axis:: |
| 633 | |
| 634 | sage: sin_integral(-1.0) |
| 635 | -0.946083070367183 |
| 636 | sage: sin_integral(-2.0) |
| 637 | -1.60541297680269 |
| 638 | sage: sin_integral(-1e23) |
| 639 | -1.57079632679490 |
| 640 | |
| 641 | """ |
| 642 | import mpmath |
| 643 | return mpmath_utils_call(mpmath.si, z, parent=parent) |
| 644 | |
| 645 | def _derivative_(self, z, diff_param=None): |
| 646 | """ |
| 647 | The derivative of `\operatorname{Si}(z)` is `\sin(z)/z` if `z` is not zero. The derivative |
| 648 | at `z = 0` is `1` (but this exception is not currently implimented). |
| 649 | |
| 650 | EXAMPLES:: |
| 651 | |
| 652 | sage: x = var('x') |
| 653 | sage: f = sin_integral(x) |
| 654 | sage: f.diff(x) |
| 655 | sin(x)/x |
| 656 | |
| 657 | sage: f = sin_integral(x^2) |
| 658 | sage: f.diff(x) |
| 659 | 2*sin(x^2)/x |
| 660 | |
| 661 | """ |
| 662 | return sin(z)/z |
| 663 | |
| 664 | si = sin_integral = Function_sin_integral() |
| 665 | |
| 666 | |
| 667 | class Function_cos_integral(BuiltinFunction): |
| 668 | r""" |
| 669 | The trigonometric integral `\operatorname{Ci}(z)` defined by |
| 670 | |
| 671 | .. math:: |
| 672 | |
| 673 | \operatorname{Ci}(z) = \gamma + \log(z) + \int_0^z \frac{\cos(t)-1}{t}\; dt, |
| 674 | |
| 675 | where `\gamma` is the Euler gamma constant (``euler_gamma`` in Sage), |
| 676 | see [AS]_ 5.2.1. |
| 677 | |
| 678 | EXAMPLES: |
| 679 | |
| 680 | Numerical evaluation for real and complex arguments is handled using mpmath:: |
| 681 | |
| 682 | sage: cos_integral(3.0) |
| 683 | 0.119629786008000 |
| 684 | |
| 685 | Compare ``cos_integral(3.0)`` to the definition of the value using |
| 686 | numerical integration:: |
| 687 | |
| 688 | sage: N(euler_gamma + log(3.0) + integrate((cos(x)-1)/x, x, 0, 3.0) - cos_integral(3.0)) < 1e-14 |
| 689 | True |
| 690 | |
| 691 | Arbitrary precision and complex arguments are handled:: |
| 692 | |
| 693 | sage: N(cos_integral(3), digits=30) |
| 694 | 0.119629786008000327626472281177 |
| 695 | sage: cos_integral(ComplexField(100)(3+I)) |
| 696 | 0.078134230477495714401983633057 - 0.37814733904787920181190368789*I |
| 697 | |
| 698 | The limit `\operatorname{Ci}(z)` as `z \to \infty` is zero:: |
| 699 | |
| 700 | sage: N(cos_integral(1e23)) |
| 701 | -3.24053937643003e-24 |
| 702 | |
| 703 | Symbolic derivatives and integrals are handled by Sage and Maxima:: |
| 704 | |
| 705 | sage: x = var('x') |
| 706 | sage: f = cos_integral(x) |
| 707 | sage: f.diff(x) |
| 708 | cos(x)/x |
| 709 | |
| 710 | sage: f.integrate(x) |
| 711 | x*cos_integral(x) - sin(x) |
| 712 | |
| 713 | The Nielsen spiral is the parametric plot of (Si(t), Ci(t)):: |
| 714 | |
| 715 | sage: t=var('t') |
| 716 | sage: f(t) = sin_integral(t) |
| 717 | sage: g(t) = cos_integral(t) |
| 718 | sage: P = parametric_plot([f, g], (t, 0.5 ,20)) |
| 719 | sage: show(P, frame=True, axes=False) |
| 720 | |
| 721 | ALGORITHM: |
| 722 | |
| 723 | Numerical evaluation is handled using mpmath, but symbolics are handled |
| 724 | by Sage and Maxima. |
| 725 | |
| 726 | REFERENCES: |
| 727 | |
| 728 | - http://en.wikipedia.org/wiki/Trigonometric_integral |
| 729 | - mpmath documentation: `ci`_ |
| 730 | |
| 731 | .. _`ci`: http://mpmath.googlecode.com/svn/trunk/doc/build/functions/expintegrals.html#ci |
| 732 | |
| 733 | """ |
| 734 | def __init__(self): |
| 735 | """ |
| 736 | See the docstring for :class:`Function_cos_integral`. |
| 737 | |
| 738 | EXAMPLES:: |
| 739 | |
| 740 | sage: cos_integral(1) |
| 741 | cos_integral(1) |
| 742 | |
| 743 | """ |
| 744 | BuiltinFunction.__init__(self, "cos_integral", nargs=1, |
| 745 | latex_name=r'\operatorname{Ci}', |
| 746 | conversions=dict(maxima='expintegral_ci')) |
| 747 | |
| 748 | def _eval_(self, z): |
| 749 | """ |
| 750 | EXAMPLES:: |
| 751 | |
| 752 | sage: z = var('z') |
| 753 | sage: cos_integral(z) |
| 754 | cos_integral(z) |
| 755 | sage: cos_integral(3.0) |
| 756 | 0.119629786008000 |
| 757 | sage: cos_integral(0) |
| 758 | cos_integral(0) |
| 759 | sage: N(cos_integral(0)) |
| 760 | -infinity |
| 761 | |
| 762 | """ |
| 763 | if not isinstance(z, Expression) and is_inexact(z): |
| 764 | return self._evalf_(z, parent(z)) |
| 765 | |
| 766 | return None # leaves the expression unevaluated |
| 767 | |
| 768 | def _evalf_(self, z, parent=None): |
| 769 | """ |
| 770 | EXAMPLES:: |
| 771 | |
| 772 | sage: N(cos_integral(1e23)) < 1e-20 |
| 773 | True |
| 774 | sage: N(cos_integral(1e-10), digits=30) |
| 775 | -22.4486352650389235918737540487 |
| 776 | sage: cos_integral(ComplexField(100)(I)) |
| 777 | 0.83786694098020824089467857943 + 1.5707963267948966192313216916*I |
| 778 | |
| 779 | """ |
| 780 | import mpmath |
| 781 | return mpmath_utils_call(mpmath.ci, z, parent=parent) |
| 782 | |
| 783 | def _derivative_(self, z, diff_param=None): |
| 784 | """ |
| 785 | The derivative of `\operatorname{Ci}(z)` is `\cos(z)/z` if `z` is not zero. |
| 786 | |
| 787 | EXAMPLES:: |
| 788 | |
| 789 | sage: x = var('x') |
| 790 | sage: f = cos_integral(x) |
| 791 | sage: f.diff(x) |
| 792 | cos(x)/x |
| 793 | |
| 794 | sage: f = cos_integral(x^2) |
| 795 | sage: f.diff(x) |
| 796 | 2*cos(x^2)/x |
| 797 | |
| 798 | """ |
| 799 | return cos(z)/z |
| 800 | |
| 801 | ci = cos_integral = Function_cos_integral() |
| 802 | |
| 803 | |
| 804 | class Function_sinh_integral(BuiltinFunction): |
| 805 | r""" |
| 806 | The trigonometric integral `\operatorname{Shi}(z)` defined by |
| 807 | |
| 808 | .. math:: |
| 809 | |
| 810 | \operatorname{Shi}(z) = \int_0^z \frac{\sinh(t)}{t}\; dt, |
| 811 | |
| 812 | see [AS]_ 5.2.3. |
| 813 | |
| 814 | EXAMPLES: |
| 815 | |
| 816 | Numerical evaluation for real and complex arguments is handled using mpmath:: |
| 817 | |
| 818 | sage: sinh_integral(3.0) |
| 819 | 4.97344047585981 |
| 820 | sage: sinh_integral(1.0) |
| 821 | 1.05725087537573 |
| 822 | sage: sinh_integral(-1.0) |
| 823 | -1.05725087537573 |
| 824 | |
| 825 | Compare ``sinh_integral(3.0)`` to the definition of the value using |
| 826 | numerical integration:: |
| 827 | |
| 828 | sage: N(integrate((sinh(x))/x, x, 0, 3.0) - sinh_integral(3.0)) < 1e-14 |
| 829 | True |
| 830 | |
| 831 | Arbitrary precision and complex arguments are handled:: |
| 832 | |
| 833 | sage: N(sinh_integral(3), digits=30) |
| 834 | 4.97344047585980679771041838252 |
| 835 | sage: sinh_integral(ComplexField(100)(3+I)) |
| 836 | 3.9134623660329374406788354078 + 3.0427678212908839256360163759*I |
| 837 | |
| 838 | The limit `\operatorname{Shi}(z)` as `z \to \infty` is `\infty`:: |
| 839 | |
| 840 | sage: N(sinh_integral(Infinity)) |
| 841 | +infinity |
| 842 | |
| 843 | Symbolic derivatives and integrals are handled by Sage and Maxima:: |
| 844 | |
| 845 | sage: x = var('x') |
| 846 | sage: f = sinh_integral(x) |
| 847 | sage: f.diff(x) |
| 848 | sinh(x)/x |
| 849 | |
| 850 | sage: f.integrate(x) |
| 851 | x*sinh_integral(x) - cosh(x) |
| 852 | |
| 853 | ALGORITHM: |
| 854 | |
| 855 | Numerical evaluation is handled using mpmath, but symbolics are handled |
| 856 | by Sage and Maxima. |
| 857 | |
| 858 | REFERENCES: |
| 859 | |
| 860 | - http://en.wikipedia.org/wiki/Trigonometric_integral |
| 861 | - mpmath documentation: `shi`_ |
| 862 | |
| 863 | .. _`shi`: http://mpmath.googlecode.com/svn/trunk/doc/build/functions/expintegrals.html#shi |
| 864 | |
| 865 | """ |
| 866 | def __init__(self): |
| 867 | """ |
| 868 | See the docstring for ``Function_sinh_integral``. |
| 869 | |
| 870 | EXAMPLES:: |
| 871 | |
| 872 | sage: sinh_integral(1) |
| 873 | sinh_integral(1) |
| 874 | |
| 875 | """ |
| 876 | BuiltinFunction.__init__(self, "sinh_integral", nargs=1, |
| 877 | latex_name=r'\operatorname{Shi}', |
| 878 | conversions=dict(maxima='expintegral_shi')) |
| 879 | |
| 880 | def _eval_(self, z): |
| 881 | """ |
| 882 | EXAMPLES:: |
| 883 | |
| 884 | sage: z = var('z') |
| 885 | sage: sinh_integral(z) |
| 886 | sinh_integral(z) |
| 887 | sage: sinh_integral(3.0) |
| 888 | 4.97344047585981 |
| 889 | sage: sinh_integral(0) |
| 890 | 0 |
| 891 | |
| 892 | """ |
| 893 | if not isinstance(z, Expression) and is_inexact(z): |
| 894 | return self._evalf_(z, parent(z)) |
| 895 | |
| 896 | # special case: z = 0 |
| 897 | if isinstance(z, Expression): |
| 898 | if z.is_trivial_zero(): |
| 899 | return z |
| 900 | else: |
| 901 | if not z: |
| 902 | return z |
| 903 | |
| 904 | return None # leaves the expression unevaluated |
| 905 | |
| 906 | def _evalf_(self, z, parent=None): |
| 907 | """ |
| 908 | EXAMPLES:: |
| 909 | |
| 910 | sage: N(sinh_integral(1e-10), digits=30) |
| 911 | 1.00000000000000003643219731550e-10 |
| 912 | sage: sinh_integral(ComplexField(100)(I)) |
| 913 | 0.94608307036718301494135331382*I |
| 914 | |
| 915 | """ |
| 916 | import mpmath |
| 917 | return mpmath_utils_call(mpmath.shi, z, parent=parent) |
| 918 | |
| 919 | def _derivative_(self, z, diff_param=None): |
| 920 | """ |
| 921 | The derivative of `\operatorname{Shi}(z)` is `\sinh(z)/z`. |
| 922 | |
| 923 | EXAMPLES:: |
| 924 | |
| 925 | sage: x = var('x') |
| 926 | sage: f = sinh_integral(x) |
| 927 | sage: f.diff(x) |
| 928 | sinh(x)/x |
| 929 | |
| 930 | sage: f = sinh_integral(ln(x)) |
| 931 | sage: f.diff(x) |
| 932 | sinh(log(x))/(x*log(x)) |
| 933 | |
| 934 | """ |
| 935 | return sinh(z)/z |
| 936 | |
| 937 | shi = sinh_integral = Function_sinh_integral() |
| 938 | |
| 939 | |
| 940 | class Function_cosh_integral(BuiltinFunction): |
| 941 | r""" |
| 942 | The trigonometric integral `\operatorname{Chi}(z)` defined by |
| 943 | |
| 944 | .. math:: |
| 945 | |
| 946 | \operatorname{Chi}(z) = \gamma + \log(z) + \int_0^z \frac{\cosh(t)-1}{t}\; dt, |
| 947 | |
| 948 | see [AS]_ 5.2.4. |
| 949 | |
| 950 | EXAMPLES: |
| 951 | |
| 952 | Numerical evaluation for real and complex arguments is handled using mpmath:: |
| 953 | |
| 954 | sage: cosh_integral(1.0) |
| 955 | 0.837866940980208 |
| 956 | |
| 957 | Here is an example from the mpmath documentation:: |
| 958 | |
| 959 | sage: f(x) = cosh_integral(x) |
| 960 | sage: find_root(f, 0.1, 1.0) |
| 961 | 0.5238225713894826 |
| 962 | |
| 963 | Compare ``cosh_integral(3.0)`` to the definition of the value using |
| 964 | numerical integration:: |
| 965 | |
| 966 | sage: N(euler_gamma + log(3.0) + integrate((cosh(x)-1)/x, x, 0, 3.0) - |
| 967 | ... cosh_integral(3.0)) < 1e-14 |
| 968 | True |
| 969 | |
| 970 | Arbitrary precision and complex arguments are handled:: |
| 971 | |
| 972 | sage: N(cosh_integral(3), digits=30) |
| 973 | 4.96039209476560976029791763669 |
| 974 | sage: cosh_integral(ComplexField(100)(3+I)) |
| 975 | 3.9096723099686417127843516794 + 3.0547519627014217273323873274*I |
| 976 | |
| 977 | The limit of `\operatorname{Chi}(z)` as `z \to \infty` is `\infty`:: |
| 978 | |
| 979 | sage: N(cosh_integral(Infinity)) |
| 980 | +infinity |
| 981 | |
| 982 | Symbolic derivatives and integrals are handled by Sage and Maxima:: |
| 983 | |
| 984 | sage: x = var('x') |
| 985 | sage: f = cosh_integral(x) |
| 986 | sage: f.diff(x) |
| 987 | cosh(x)/x |
| 988 | |
| 989 | sage: f.integrate(x) |
| 990 | x*cosh_integral(x) - sinh(x) |
| 991 | |
| 992 | ALGORITHM: |
| 993 | |
| 994 | Numerical evaluation is handled using mpmath, but symbolics are handled |
| 995 | by Sage and Maxima. |
| 996 | |
| 997 | REFERENCES: |
| 998 | |
| 999 | - http://en.wikipedia.org/wiki/Trigonometric_integral |
| 1000 | - mpmath documentation: `chi`_ |
| 1001 | |
| 1002 | .. _`chi`: http://mpmath.googlecode.com/svn/trunk/doc/build/functions/expintegrals.html#chi |
| 1003 | |
| 1004 | """ |
| 1005 | def __init__(self): |
| 1006 | """ |
| 1007 | See the docstring for ``Function_cosh_integral``. |
| 1008 | |
| 1009 | EXAMPLES:: |
| 1010 | |
| 1011 | sage: cosh_integral(1) |
| 1012 | cosh_integral(1) |
| 1013 | |
| 1014 | """ |
| 1015 | BuiltinFunction.__init__(self, "cosh_integral", nargs=1, |
| 1016 | latex_name=r'\operatorname{Chi}', |
| 1017 | conversions=dict(maxima='expintegral_chi')) |
| 1018 | |
| 1019 | def _eval_(self, z): |
| 1020 | """ |
| 1021 | EXAMPLES:: |
| 1022 | |
| 1023 | sage: z = var('z') |
| 1024 | sage: cosh_integral(z) |
| 1025 | cosh_integral(z) |
| 1026 | sage: cosh_integral(3.0) |
| 1027 | 4.96039209476561 |
| 1028 | |
| 1029 | """ |
| 1030 | if not isinstance(z, Expression) and is_inexact(z): |
| 1031 | return self._evalf_(z, parent(z)) |
| 1032 | |
| 1033 | return None |
| 1034 | |
| 1035 | def _evalf_(self, z, parent=None): |
| 1036 | """ |
| 1037 | EXAMPLES:: |
| 1038 | |
| 1039 | sage: N(cosh_integral(1e-10), digits=30) |
| 1040 | -22.4486352650389235918737540487 |
| 1041 | sage: cosh_integral(ComplexField(100)(I)) |
| 1042 | 0.33740392290096813466264620389 + 1.5707963267948966192313216916*I |
| 1043 | |
| 1044 | """ |
| 1045 | import mpmath |
| 1046 | return mpmath_utils_call(mpmath.chi, z, parent=parent) |
| 1047 | |
| 1048 | def _derivative_(self, z, diff_param=None): |
| 1049 | """ |
| 1050 | The derivative of `\operatorname{Chi}(z)` is `\cosh(z)/z`. |
| 1051 | |
| 1052 | EXAMPLES:: |
| 1053 | |
| 1054 | sage: x = var('x') |
| 1055 | sage: f = cosh_integral(x) |
| 1056 | sage: f.diff(x) |
| 1057 | cosh(x)/x |
| 1058 | |
| 1059 | sage: f = cosh_integral(ln(x)) |
| 1060 | sage: f.diff(x) |
| 1061 | cosh(log(x))/(x*log(x)) |
| 1062 | |
| 1063 | """ |
| 1064 | return cosh(z)/z |
| 1065 | |
| 1066 | chi = cosh_integral = Function_cosh_integral() |
| 1067 | |
| 1068 | |
| 1069 | ############################################# |
| 1070 | ## Code below here was moved from either: |
| 1071 | ## - sage/functions/transcendental.py |
| 1072 | ## - sage/functions/special.py |
| 1073 | ## This occured as part of Trac #11143. |
| 1074 | ############################################# |
| 1075 | |
| 1076 | # moved here from sage/functions/special.py |
| 1077 | def exp_int(t): |
| 1078 | r""" |
| 1079 | The exponential integral `\int_t^\infty e^{-x}/x \; dx` (`t` |
| 1080 | belongs to `\mathbb{R}`). This function is deprecated - please use |
| 1081 | ``Ei`` or ``exponential_integral_1`` as needed instead. |
| 1082 | |
| 1083 | EXAMPLES:: |
| 1084 | |
| 1085 | sage: exp_int(6) |
| 1086 | doctest:...: DeprecationWarning: The method exp_int() is deprecated. Use -Ei(-x) or exponential_integral_1(x) as needed instead. |
| 1087 | 0.000360082452162659 |
| 1088 | """ |
| 1089 | from sage.misc.misc import deprecation |
| 1090 | deprecation("The method exp_int() is deprecated. Use -Ei(-x) or exponential_integral_1(x) as needed instead.") |
| 1091 | try: |
| 1092 | return t.eint1() |
| 1093 | except AttributeError: |
| 1094 | from sage.libs.pari.all import pari |
| 1095 | try: |
| 1096 | return pari(t).eint1() |
| 1097 | except: |
| 1098 | raise NotImplementedError |
| 1099 | |
| 1100 | # moved here from sage/functions/transcendental.py |
| 1101 | # |
| 1102 | # This class has a name which is not specific enough |
| 1103 | # see Function_exp_integral_e above, for example, which |
| 1104 | # is the "generalized" exponential integral function. We |
| 1105 | # are leaving the name the same for backwards compatibility |
| 1106 | # purposes. |
| 1107 | class Function_exp_integral(BuiltinFunction): |
| 1108 | r""" |
| 1109 | The generalized complex exponential integral Ei(z) defined by |
| 1110 | |
| 1111 | .. math:: |
| 1112 | |
| 1113 | \operatorname{Ei}(x) = \int_{-\infty}^x \frac{e^t}{t}\; dt |
| 1114 | |
| 1115 | for x > 0 and for complex arguments by analytic continuation, |
| 1116 | see [AS]_ 5.1.2. |
| 1117 | |
| 1118 | EXAMPLES:: |
| 1119 | |
| 1120 | sage: Ei(10) |
| 1121 | Ei(10) |
| 1122 | sage: Ei(I) |
| 1123 | Ei(I) |
| 1124 | sage: Ei(3+I) |
| 1125 | Ei(I + 3) |
| 1126 | sage: Ei(1.3) |
| 1127 | 2.72139888023202 |
| 1128 | |
| 1129 | The branch cut for this function is along the negative real axis:: |
| 1130 | |
| 1131 | sage: Ei(-3 + 0.1*I) |
| 1132 | -0.0129379427181693 + 3.13993830250942*I |
| 1133 | sage: Ei(-3 - 0.1*I) |
| 1134 | -0.0129379427181693 - 3.13993830250942*I |
| 1135 | |
| 1136 | ALGORITHM: Uses mpmath. |
| 1137 | |
| 1138 | """ |
| 1139 | def __init__(self): |
| 1140 | """ |
| 1141 | Return the value of the complex exponential integral Ei(z) at a |
| 1142 | complex number z. |
| 1143 | |
| 1144 | EXAMPLES:: |
| 1145 | |
| 1146 | sage: Ei(10) |
| 1147 | Ei(10) |
| 1148 | sage: Ei(I) |
| 1149 | Ei(I) |
| 1150 | sage: Ei(3+I) |
| 1151 | Ei(I + 3) |
| 1152 | sage: Ei(1.3) |
| 1153 | 2.72139888023202 |
| 1154 | |
| 1155 | The branch cut for this function is along the negative real axis:: |
| 1156 | |
| 1157 | sage: Ei(-3 + 0.1*I) |
| 1158 | -0.0129379427181693 + 3.13993830250942*I |
| 1159 | sage: Ei(-3 - 0.1*I) |
| 1160 | -0.0129379427181693 - 3.13993830250942*I |
| 1161 | |
| 1162 | ALGORITHM: Uses mpmath. |
| 1163 | """ |
| 1164 | BuiltinFunction.__init__(self, "Ei") |
| 1165 | |
| 1166 | def _eval_(self, x ): |
| 1167 | """ |
| 1168 | EXAMPLES:: |
| 1169 | |
| 1170 | sage: Ei(10) |
| 1171 | Ei(10) |
| 1172 | sage: Ei(I) |
| 1173 | Ei(I) |
| 1174 | sage: Ei(1.3) |
| 1175 | 2.72139888023202 |
| 1176 | sage: Ei(10r) |
| 1177 | Ei(10) |
| 1178 | sage: Ei(1.3r) |
| 1179 | 2.7213988802320235 |
| 1180 | """ |
| 1181 | if not isinstance(x, Expression) and is_inexact(x): |
| 1182 | return self._evalf_(x, parent(x)) |
| 1183 | return None |
| 1184 | |
| 1185 | def _evalf_(self, x, parent=None): |
| 1186 | """ |
| 1187 | EXAMPLES:: |
| 1188 | |
| 1189 | sage: Ei(10).n() |
| 1190 | 2492.22897624188 |
| 1191 | sage: Ei(20).n() |
| 1192 | 2.56156526640566e7 |
| 1193 | sage: Ei(I).n() |
| 1194 | 0.337403922900968 + 2.51687939716208*I |
| 1195 | sage: Ei(3+I).n() |
| 1196 | 7.82313467600158 + 6.09751978399231*I |
| 1197 | """ |
| 1198 | import mpmath |
| 1199 | return mpmath_utils_call(mpmath.ei, x, parent=parent) |
| 1200 | |
| 1201 | def __call__(self, x, prec=None, coerce=True, hold=False ): |
| 1202 | """ |
| 1203 | Note that the ``prec`` argument is deprecated. The precision for |
| 1204 | the result is deduced from the precision of the input. Convert |
| 1205 | the input to a higher precision explicitly if a result with higher |
| 1206 | precision is desired. |
| 1207 | |
| 1208 | EXAMPLES:: |
| 1209 | |
| 1210 | sage: t = Ei(RealField(100)(2.5)); t |
| 1211 | 7.0737658945786007119235519625 |
| 1212 | sage: t.prec() |
| 1213 | 100 |
| 1214 | |
| 1215 | sage: Ei(1.1, prec=300) |
| 1216 | doctest:...: DeprecationWarning: The prec keyword argument is deprecated. Explicitly set the precision of the input, for example Ei(RealField(300)(1)), or use the prec argument to .n() for exact inputs, e.g., Ei(1).n(300), instead. |
| 1217 | 2.16737827956340306615064476647912607220394065907142504328679588538509331805598360907980986 |
| 1218 | """ |
| 1219 | if prec is not None: |
| 1220 | from sage.misc.misc import deprecation |
| 1221 | deprecation("The prec keyword argument is deprecated. Explicitly set the precision of the input, for example Ei(RealField(300)(1)), or use the prec argument to .n() for exact inputs, e.g., Ei(1).n(300), instead.") |
| 1222 | |
| 1223 | import mpmath |
| 1224 | return mpmath_utils_call(mpmath.ei, x, prec=prec) |
| 1225 | |
| 1226 | return BuiltinFunction.__call__(self, x, coerce=coerce, hold=hold) |
| 1227 | |
| 1228 | def _derivative_(self, x, diff_param=None): |
| 1229 | """ |
| 1230 | EXAMPLES:: |
| 1231 | |
| 1232 | sage: Ei(x).diff(x) |
| 1233 | e^x/x |
| 1234 | sage: Ei(x).diff(x).subs(x=1) |
| 1235 | e |
| 1236 | sage: Ei(x^2).diff(x) |
| 1237 | 2*e^(x^2)/x |
| 1238 | sage: f = function('f') |
| 1239 | sage: Ei(f(x)).diff(x) |
| 1240 | e^f(x)*D[0](f)(x)/f(x) |
| 1241 | """ |
| 1242 | return exp(x)/x |
| 1243 | |
| 1244 | Ei = exp_integral_ei = Function_exp_integral() |
| 1245 | |
| 1246 | |
| 1247 | # moved here from sage/functions/transcendental.py |
| 1248 | def exponential_integral_1(x, n=0): |
| 1249 | r""" |
| 1250 | Returns the exponential integral `E_1(x)`. If the optional |
| 1251 | argument `n` is given, computes list of the first |
| 1252 | `n` values of the exponential integral |
| 1253 | `E_1(x m)`. |
| 1254 | |
| 1255 | The exponential integral `E_1(x)` is |
| 1256 | |
| 1257 | .. math:: |
| 1258 | |
| 1259 | E_1(x) = \int_{x}^{\infty} e^{-t}/t dt |
| 1260 | |
| 1261 | |
| 1262 | |
| 1263 | INPUT: |
| 1264 | |
| 1265 | |
| 1266 | - ``x`` - a positive real number |
| 1267 | |
| 1268 | - ``n`` - (default: 0) a nonnegative integer; if |
| 1269 | nonzero, then return a list of values E_1(x\*m) for m = |
| 1270 | 1,2,3,...,n. This is useful, e.g., when computing derivatives of |
| 1271 | L-functions. |
| 1272 | |
| 1273 | |
| 1274 | OUTPUT: |
| 1275 | |
| 1276 | |
| 1277 | - ``float`` - if n is 0 (the default) or |
| 1278 | |
| 1279 | - ``list`` - list of floats if n 0 |
| 1280 | |
| 1281 | |
| 1282 | EXAMPLES:: |
| 1283 | |
| 1284 | sage: exponential_integral_1(2) |
| 1285 | 0.04890051070806112 |
| 1286 | sage: w = exponential_integral_1(2,4); w |
| 1287 | [0.04890051070806112, 0.0037793524098489067, 0.00036008245216265873, 3.7665622843924751e-05] # 32-bit |
| 1288 | [0.04890051070806112, 0.0037793524098489063, 0.00036008245216265873, 3.7665622843924534e-05] # 64-bit |
| 1289 | |
| 1290 | IMPLEMENTATION: We use the PARI C-library functions eint1 and |
| 1291 | veceint1. |
| 1292 | |
| 1293 | REFERENCE: |
| 1294 | |
| 1295 | - See page 262, Prop 5.6.12, of Cohen's book "A Course in |
| 1296 | Computational Algebraic Number Theory". |
| 1297 | |
| 1298 | REMARKS: When called with the optional argument n, the PARI |
| 1299 | C-library is fast for values of n up to some bound, then very very |
| 1300 | slow. For example, if x=5, then the computation takes less than a |
| 1301 | second for n=800000, and takes "forever" for n=900000. |
| 1302 | """ |
| 1303 | if isinstance(x, Expression): |
| 1304 | raise NotImplementedError("Use the symbolic exponential integral " + |
| 1305 | "function: exp_integral_e1.") |
| 1306 | from sage.libs.pari.all import pari |
| 1307 | if n <= 0: |
| 1308 | return float(pari(x).eint1()) |
| 1309 | else: |
| 1310 | return [float(z) for z in pari(x).eint1(n)] |