| 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 | sage: log_integral(0) |
| 385 | 0 |
| 386 | |
| 387 | Symbolic derivatives and integrals are handled by Sage and Maxima:: |
| 388 | |
| 389 | sage: x = var('x') |
| 390 | sage: f = log_integral(x) |
| 391 | sage: f.diff(x) |
| 392 | 1/log(x) |
| 393 | |
| 394 | sage: f.integrate(x) |
| 395 | x*log_integral(x) - Ei(2*log(x)) |
| 396 | |
| 397 | Here is a test from the mpmath documentation. There are |
| 398 | 1,925,320,391,606,803,968,923 many prime numbers less than 1e23. The |
| 399 | value of ``log_integral(1e23)`` is very close to this:: |
| 400 | |
| 401 | sage: log_integral(1e23) |
| 402 | 1.92532039161405e21 |
| 403 | |
| 404 | ALGORITHM: |
| 405 | |
| 406 | Numerical evaluation is handled using mpmath, but symbolics are handled |
| 407 | by Sage and Maxima. |
| 408 | |
| 409 | REFERENCES: |
| 410 | |
| 411 | - http://en.wikipedia.org/wiki/Logarithmic_integral_function |
| 412 | - mpmath documentation: `logarithmic-integral`_ |
| 413 | |
| 414 | .. _`logarithmic-integral`: http://mpmath.googlecode.com/svn/trunk/doc/build/functions/expintegrals.html#logarithmic-integral |
| 415 | |
| 416 | |
| 417 | """ |
| 418 | def __init__(self): |
| 419 | """ |
| 420 | See the docstring for ``Function_log_integral``. |
| 421 | |
| 422 | EXAMPLES:: |
| 423 | |
| 424 | sage: log_integral(3) |
| 425 | log_integral(3) |
| 426 | |
| 427 | """ |
| 428 | BuiltinFunction.__init__(self, "log_integral", nargs=1, |
| 429 | latex_name=r'log_integral', |
| 430 | conversions=dict(maxima='expintegral_li')) |
| 431 | |
| 432 | def _eval_(self, z): |
| 433 | """ |
| 434 | EXAMPLES:: |
| 435 | |
| 436 | sage: z = var('z') |
| 437 | sage: log_integral(z) |
| 438 | log_integral(z) |
| 439 | sage: log_integral(3.0) |
| 440 | 2.16358859466719 |
| 441 | sage: log_integral(0) |
| 442 | 0 |
| 443 | |
| 444 | """ |
| 445 | if isinstance(z, Expression): |
| 446 | if z.is_trivial_zero(): # special case: z = 0 |
| 447 | return z |
| 448 | else: |
| 449 | if is_inexact(z): |
| 450 | return self._evalf_(z, parent(z)) |
| 451 | elif not z: |
| 452 | return z |
| 453 | return None # leaves the expression unevaluated |
| 454 | |
| 455 | def _evalf_(self, z, parent=None): |
| 456 | """ |
| 457 | EXAMPLES:: |
| 458 | |
| 459 | sage: N(log_integral(1e6)) |
| 460 | 78627.5491594622 |
| 461 | sage: log_integral(RealField(200)(1e6)) |
| 462 | 78627.549159462181919862910747947261161321874382421767074759 |
| 463 | |
| 464 | """ |
| 465 | import mpmath |
| 466 | return mpmath_utils_call(mpmath.li, z, parent=parent) |
| 467 | |
| 468 | def _derivative_(self, z, diff_param=None): |
| 469 | """ |
| 470 | The derivative of `\operatorname{li}(z) is `1/log(z)`. |
| 471 | |
| 472 | EXAMPLES:: |
| 473 | |
| 474 | sage: x = var('x') |
| 475 | sage: f = log_integral(x) |
| 476 | sage: f.diff(x) |
| 477 | 1/log(x) |
| 478 | |
| 479 | sage: f = log_integral(x^2) |
| 480 | sage: f.diff(x) |
| 481 | 2*x/log(x^2) |
| 482 | |
| 483 | """ |
| 484 | return 1/log(z) |
| 485 | |
| 486 | li = log_integral = Function_log_integral() |
| 487 | |
| 488 | |
| 489 | class Function_sin_integral(BuiltinFunction): |
| 490 | r""" |
| 491 | The trigonometric integral `\operatorname{Si}(z)` defined by |
| 492 | |
| 493 | .. math:: |
| 494 | |
| 495 | \operatorname{Si}(z) = \int_0^z \frac{\sin(t)}{t}\; dt, |
| 496 | |
| 497 | see [AS]_ 5.2.1. |
| 498 | |
| 499 | EXAMPLES: |
| 500 | |
| 501 | Numerical evaluation for real and complex arguments is handled using mpmath:: |
| 502 | |
| 503 | sage: sin_integral(0) |
| 504 | 0 |
| 505 | sage: sin_integral(0.0) |
| 506 | 0.000000000000000 |
| 507 | sage: sin_integral(3.0) |
| 508 | 1.84865252799947 |
| 509 | sage: N(sin_integral(3), digits=30) |
| 510 | 1.84865252799946825639773025111 |
| 511 | sage: sin_integral(ComplexField(100)(3+I)) |
| 512 | 2.0277151656451253616038525998 + 0.015210926166954211913653130271*I |
| 513 | |
| 514 | The alias `Si` can be used instead of `sin_integral`:: |
| 515 | |
| 516 | sage: Si(3.0) |
| 517 | 1.84865252799947 |
| 518 | |
| 519 | The limit of `\operatorname{Si}(z)` as `z \to \infty` is `\pi/2`:: |
| 520 | |
| 521 | sage: N(sin_integral(1e23)) |
| 522 | 1.57079632679490 |
| 523 | sage: N(pi/2) |
| 524 | 1.57079632679490 |
| 525 | |
| 526 | At 200 bits of precision `\operatorname{Si}(10^{23})` agrees with `\pi/2` up to |
| 527 | `10^{-24}`:: |
| 528 | |
| 529 | sage: sin_integral(RealField(200)(1e23)) |
| 530 | 1.5707963267948966192313288218697837425815368604836679189519 |
| 531 | sage: N(pi/2, prec=200) |
| 532 | 1.5707963267948966192313216916397514420985846996875529104875 |
| 533 | |
| 534 | The exponential sine integral is analytic everywhere:: |
| 535 | |
| 536 | sage: sin_integral(-1.0) |
| 537 | -0.946083070367183 |
| 538 | sage: sin_integral(-2.0) |
| 539 | -1.60541297680269 |
| 540 | sage: sin_integral(-1e23) |
| 541 | -1.57079632679490 |
| 542 | |
| 543 | Symbolic derivatives and integrals are handled by Sage and Maxima:: |
| 544 | |
| 545 | sage: x = var('x') |
| 546 | sage: f = sin_integral(x) |
| 547 | sage: f.diff(x) |
| 548 | sin(x)/x |
| 549 | |
| 550 | sage: f.integrate(x) |
| 551 | x*sin_integral(x) + cos(x) |
| 552 | |
| 553 | sage: integrate(sin(x)/x, x) |
| 554 | 1/2*I*Ei(-I*x) - 1/2*I*Ei(I*x) |
| 555 | |
| 556 | Compare values of the functions `\operatorname{Si}(x)` and |
| 557 | `f(x) = (1/2)i \cdot \operatorname{Ei}(-ix) - (1/2)i \cdot |
| 558 | \operatorname{Ei}(ix) - \pi/2`, which are both anti-derivatives of |
| 559 | `\sin(x)/x`, at some random positive real numbers:: |
| 560 | |
| 561 | sage: f(x) = 1/2*I*Ei(-I*x) - 1/2*I*Ei(I*x) - pi/2 |
| 562 | sage: g(x) = sin_integral(x) |
| 563 | sage: R = [ abs(RDF.random_element()) for i in range(100) ] |
| 564 | sage: all(abs(f(x) - g(x)) < 1e-10 for x in R) |
| 565 | True |
| 566 | |
| 567 | The Nielsen spiral is the parametric plot of (Si(t), Ci(t)):: |
| 568 | |
| 569 | sage: x=var('x') |
| 570 | sage: f(x) = sin_integral(x) |
| 571 | sage: g(x) = cos_integral(x) |
| 572 | sage: P = parametric_plot([f, g], (x, 0.5 ,20)) |
| 573 | sage: show(P, frame=True, axes=False) |
| 574 | |
| 575 | ALGORITHM: |
| 576 | |
| 577 | Numerical evaluation is handled using mpmath, but symbolics are handled |
| 578 | by Sage and Maxima. |
| 579 | |
| 580 | REFERENCES: |
| 581 | |
| 582 | - http://en.wikipedia.org/wiki/Trigonometric_integral |
| 583 | - mpmath documentation: `si`_ |
| 584 | |
| 585 | .. _`si`: http://mpmath.googlecode.com/svn/trunk/doc/build/functions/expintegrals.html#si |
| 586 | |
| 587 | """ |
| 588 | def __init__(self): |
| 589 | """ |
| 590 | See the docstring for ``Function_sin_integral``. |
| 591 | |
| 592 | EXAMPLES:: |
| 593 | |
| 594 | sage: sin_integral(1) |
| 595 | sin_integral(1) |
| 596 | |
| 597 | """ |
| 598 | BuiltinFunction.__init__(self, "sin_integral", nargs=1, |
| 599 | latex_name=r'\operatorname{Si}', |
| 600 | conversions=dict(maxima='expintegral_si')) |
| 601 | |
| 602 | def _eval_(self, z): |
| 603 | """ |
| 604 | EXAMPLES:: |
| 605 | |
| 606 | sage: z = var('z') |
| 607 | sage: sin_integral(z) |
| 608 | sin_integral(z) |
| 609 | sage: sin_integral(3.0) |
| 610 | 1.84865252799947 |
| 611 | sage: sin_integral(0) |
| 612 | 0 |
| 613 | |
| 614 | """ |
| 615 | if not isinstance(z, Expression) and is_inexact(z): |
| 616 | return self._evalf_(z, parent(z)) |
| 617 | |
| 618 | # special case: z = 0 |
| 619 | if isinstance(z, Expression): |
| 620 | if z.is_trivial_zero(): |
| 621 | return z |
| 622 | else: |
| 623 | if not z: |
| 624 | return z |
| 625 | |
| 626 | return None # leaves the expression unevaluated |
| 627 | |
| 628 | def _evalf_(self, z, parent=None): |
| 629 | """ |
| 630 | EXAMPLES: |
| 631 | |
| 632 | The limit `\operatorname{Si}(z)` as `z \to \infty` is `\pi/2`:: |
| 633 | |
| 634 | sage: N(sin_integral(1e23) - pi/2) |
| 635 | 0.000000000000000 |
| 636 | |
| 637 | At 200 bits of precision `\operatorname{Si}(10^{23})` agrees with `\pi/2` up to |
| 638 | `10^{-24}`:: |
| 639 | |
| 640 | sage: sin_integral(RealField(200)(1e23)) |
| 641 | 1.5707963267948966192313288218697837425815368604836679189519 |
| 642 | sage: N(pi/2, prec=200) |
| 643 | 1.5707963267948966192313216916397514420985846996875529104875 |
| 644 | |
| 645 | The exponential sine integral is analytic everywhere, even on the |
| 646 | negative real axis:: |
| 647 | |
| 648 | sage: sin_integral(-1.0) |
| 649 | -0.946083070367183 |
| 650 | sage: sin_integral(-2.0) |
| 651 | -1.60541297680269 |
| 652 | sage: sin_integral(-1e23) |
| 653 | -1.57079632679490 |
| 654 | |
| 655 | """ |
| 656 | import mpmath |
| 657 | return mpmath_utils_call(mpmath.si, z, parent=parent) |
| 658 | |
| 659 | def _derivative_(self, z, diff_param=None): |
| 660 | """ |
| 661 | The derivative of `\operatorname{Si}(z)` is `\sin(z)/z` if `z` is not zero. The derivative |
| 662 | at `z = 0` is `1` (but this exception is not currently implimented). |
| 663 | |
| 664 | EXAMPLES:: |
| 665 | |
| 666 | sage: x = var('x') |
| 667 | sage: f = sin_integral(x) |
| 668 | sage: f.diff(x) |
| 669 | sin(x)/x |
| 670 | |
| 671 | sage: f = sin_integral(x^2) |
| 672 | sage: f.diff(x) |
| 673 | 2*sin(x^2)/x |
| 674 | |
| 675 | """ |
| 676 | return sin(z)/z |
| 677 | |
| 678 | Si = sin_integral = Function_sin_integral() |
| 679 | |
| 680 | |
| 681 | class Function_cos_integral(BuiltinFunction): |
| 682 | r""" |
| 683 | The trigonometric integral `\operatorname{Ci}(z)` defined by |
| 684 | |
| 685 | .. math:: |
| 686 | |
| 687 | \operatorname{Ci}(z) = \gamma + \log(z) + \int_0^z \frac{\cos(t)-1}{t}\; dt, |
| 688 | |
| 689 | where `\gamma` is the Euler gamma constant (``euler_gamma`` in Sage), |
| 690 | see [AS]_ 5.2.1. |
| 691 | |
| 692 | EXAMPLES: |
| 693 | |
| 694 | Numerical evaluation for real and complex arguments is handled using mpmath:: |
| 695 | |
| 696 | sage: cos_integral(3.0) |
| 697 | 0.119629786008000 |
| 698 | |
| 699 | The alias `Ci` can be used instead of `cos_integral`:: |
| 700 | |
| 701 | sage: Ci(3.0) |
| 702 | 0.119629786008000 |
| 703 | |
| 704 | Compare ``cos_integral(3.0)`` to the definition of the value using |
| 705 | numerical integration:: |
| 706 | |
| 707 | sage: N(euler_gamma + log(3.0) + integrate((cos(x)-1)/x, x, 0, 3.0) - cos_integral(3.0)) < 1e-14 |
| 708 | True |
| 709 | |
| 710 | Arbitrary precision and complex arguments are handled:: |
| 711 | |
| 712 | sage: N(cos_integral(3), digits=30) |
| 713 | 0.119629786008000327626472281177 |
| 714 | sage: cos_integral(ComplexField(100)(3+I)) |
| 715 | 0.078134230477495714401983633057 - 0.37814733904787920181190368789*I |
| 716 | |
| 717 | The limit `\operatorname{Ci}(z)` as `z \to \infty` is zero:: |
| 718 | |
| 719 | sage: N(cos_integral(1e23)) |
| 720 | -3.24053937643003e-24 |
| 721 | |
| 722 | Symbolic derivatives and integrals are handled by Sage and Maxima:: |
| 723 | |
| 724 | sage: x = var('x') |
| 725 | sage: f = cos_integral(x) |
| 726 | sage: f.diff(x) |
| 727 | cos(x)/x |
| 728 | |
| 729 | sage: f.integrate(x) |
| 730 | x*cos_integral(x) - sin(x) |
| 731 | |
| 732 | The Nielsen spiral is the parametric plot of (Si(t), Ci(t)):: |
| 733 | |
| 734 | sage: t=var('t') |
| 735 | sage: f(t) = sin_integral(t) |
| 736 | sage: g(t) = cos_integral(t) |
| 737 | sage: P = parametric_plot([f, g], (t, 0.5 ,20)) |
| 738 | sage: show(P, frame=True, axes=False) |
| 739 | |
| 740 | ALGORITHM: |
| 741 | |
| 742 | Numerical evaluation is handled using mpmath, but symbolics are handled |
| 743 | by Sage and Maxima. |
| 744 | |
| 745 | REFERENCES: |
| 746 | |
| 747 | - http://en.wikipedia.org/wiki/Trigonometric_integral |
| 748 | - mpmath documentation: `ci`_ |
| 749 | |
| 750 | .. _`ci`: http://mpmath.googlecode.com/svn/trunk/doc/build/functions/expintegrals.html#ci |
| 751 | |
| 752 | """ |
| 753 | def __init__(self): |
| 754 | """ |
| 755 | See the docstring for :class:`Function_cos_integral`. |
| 756 | |
| 757 | EXAMPLES:: |
| 758 | |
| 759 | sage: cos_integral(1) |
| 760 | cos_integral(1) |
| 761 | |
| 762 | """ |
| 763 | BuiltinFunction.__init__(self, "cos_integral", nargs=1, |
| 764 | latex_name=r'\operatorname{Ci}', |
| 765 | conversions=dict(maxima='expintegral_ci')) |
| 766 | |
| 767 | def _eval_(self, z): |
| 768 | """ |
| 769 | EXAMPLES:: |
| 770 | |
| 771 | sage: z = var('z') |
| 772 | sage: cos_integral(z) |
| 773 | cos_integral(z) |
| 774 | sage: cos_integral(3.0) |
| 775 | 0.119629786008000 |
| 776 | sage: cos_integral(0) |
| 777 | cos_integral(0) |
| 778 | sage: N(cos_integral(0)) |
| 779 | -infinity |
| 780 | |
| 781 | """ |
| 782 | if not isinstance(z, Expression) and is_inexact(z): |
| 783 | return self._evalf_(z, parent(z)) |
| 784 | |
| 785 | return None # leaves the expression unevaluated |
| 786 | |
| 787 | def _evalf_(self, z, parent=None): |
| 788 | """ |
| 789 | EXAMPLES:: |
| 790 | |
| 791 | sage: N(cos_integral(1e23)) < 1e-20 |
| 792 | True |
| 793 | sage: N(cos_integral(1e-10), digits=30) |
| 794 | -22.4486352650389235918737540487 |
| 795 | sage: cos_integral(ComplexField(100)(I)) |
| 796 | 0.83786694098020824089467857943 + 1.5707963267948966192313216916*I |
| 797 | |
| 798 | """ |
| 799 | import mpmath |
| 800 | return mpmath_utils_call(mpmath.ci, z, parent=parent) |
| 801 | |
| 802 | def _derivative_(self, z, diff_param=None): |
| 803 | """ |
| 804 | The derivative of `\operatorname{Ci}(z)` is `\cos(z)/z` if `z` is not zero. |
| 805 | |
| 806 | EXAMPLES:: |
| 807 | |
| 808 | sage: x = var('x') |
| 809 | sage: f = cos_integral(x) |
| 810 | sage: f.diff(x) |
| 811 | cos(x)/x |
| 812 | |
| 813 | sage: f = cos_integral(x^2) |
| 814 | sage: f.diff(x) |
| 815 | 2*cos(x^2)/x |
| 816 | |
| 817 | """ |
| 818 | return cos(z)/z |
| 819 | |
| 820 | Ci = cos_integral = Function_cos_integral() |
| 821 | |
| 822 | |
| 823 | class Function_sinh_integral(BuiltinFunction): |
| 824 | r""" |
| 825 | The trigonometric integral `\operatorname{Shi}(z)` defined by |
| 826 | |
| 827 | .. math:: |
| 828 | |
| 829 | \operatorname{Shi}(z) = \int_0^z \frac{\sinh(t)}{t}\; dt, |
| 830 | |
| 831 | see [AS]_ 5.2.3. |
| 832 | |
| 833 | EXAMPLES: |
| 834 | |
| 835 | Numerical evaluation for real and complex arguments is handled using mpmath:: |
| 836 | |
| 837 | sage: sinh_integral(3.0) |
| 838 | 4.97344047585981 |
| 839 | sage: sinh_integral(1.0) |
| 840 | 1.05725087537573 |
| 841 | sage: sinh_integral(-1.0) |
| 842 | -1.05725087537573 |
| 843 | |
| 844 | The alias `Shi` can be used instead of `sinh_integral`:: |
| 845 | |
| 846 | sage: Shi(3.0) |
| 847 | 4.97344047585981 |
| 848 | |
| 849 | Compare ``sinh_integral(3.0)`` to the definition of the value using |
| 850 | numerical integration:: |
| 851 | |
| 852 | sage: N(integrate((sinh(x))/x, x, 0, 3.0) - sinh_integral(3.0)) < 1e-14 |
| 853 | True |
| 854 | |
| 855 | Arbitrary precision and complex arguments are handled:: |
| 856 | |
| 857 | sage: N(sinh_integral(3), digits=30) |
| 858 | 4.97344047585980679771041838252 |
| 859 | sage: sinh_integral(ComplexField(100)(3+I)) |
| 860 | 3.9134623660329374406788354078 + 3.0427678212908839256360163759*I |
| 861 | |
| 862 | The limit `\operatorname{Shi}(z)` as `z \to \infty` is `\infty`:: |
| 863 | |
| 864 | sage: N(sinh_integral(Infinity)) |
| 865 | +infinity |
| 866 | |
| 867 | Symbolic derivatives and integrals are handled by Sage and Maxima:: |
| 868 | |
| 869 | sage: x = var('x') |
| 870 | sage: f = sinh_integral(x) |
| 871 | sage: f.diff(x) |
| 872 | sinh(x)/x |
| 873 | |
| 874 | sage: f.integrate(x) |
| 875 | x*sinh_integral(x) - cosh(x) |
| 876 | |
| 877 | Note that due to some problems with the way Maxima handles these |
| 878 | expressions, definite integrals can sometimes give unexpected |
| 879 | results (typically when using inexact endpoints) due to |
| 880 | inconsistent branching:: |
| 881 | |
| 882 | sage: integrate(sinh_integral(x), x, 0, 1/2) |
| 883 | -cosh(1/2) + 1/2*sinh_integral(1/2) + 1 |
| 884 | sage: integrate(sinh_integral(x), x, 0, 1/2).n() # correct |
| 885 | 0.125872409703453 |
| 886 | sage: integrate(sinh_integral(x), x, 0, 0.5).n() # incorrect! |
| 887 | 0.125872409703453 + 1.57079632679490*I |
| 888 | |
| 889 | ALGORITHM: |
| 890 | |
| 891 | Numerical evaluation is handled using mpmath, but symbolics are handled |
| 892 | by Sage and Maxima. |
| 893 | |
| 894 | REFERENCES: |
| 895 | |
| 896 | - http://en.wikipedia.org/wiki/Trigonometric_integral |
| 897 | - mpmath documentation: `shi`_ |
| 898 | |
| 899 | .. _`shi`: http://mpmath.googlecode.com/svn/trunk/doc/build/functions/expintegrals.html#shi |
| 900 | |
| 901 | """ |
| 902 | def __init__(self): |
| 903 | """ |
| 904 | See the docstring for ``Function_sinh_integral``. |
| 905 | |
| 906 | EXAMPLES:: |
| 907 | |
| 908 | sage: sinh_integral(1) |
| 909 | sinh_integral(1) |
| 910 | |
| 911 | """ |
| 912 | BuiltinFunction.__init__(self, "sinh_integral", nargs=1, |
| 913 | latex_name=r'\operatorname{Shi}', |
| 914 | conversions=dict(maxima='expintegral_shi')) |
| 915 | |
| 916 | def _eval_(self, z): |
| 917 | """ |
| 918 | EXAMPLES:: |
| 919 | |
| 920 | sage: z = var('z') |
| 921 | sage: sinh_integral(z) |
| 922 | sinh_integral(z) |
| 923 | sage: sinh_integral(3.0) |
| 924 | 4.97344047585981 |
| 925 | sage: sinh_integral(0) |
| 926 | 0 |
| 927 | |
| 928 | """ |
| 929 | if not isinstance(z, Expression) and is_inexact(z): |
| 930 | return self._evalf_(z, parent(z)) |
| 931 | |
| 932 | # special case: z = 0 |
| 933 | if isinstance(z, Expression): |
| 934 | if z.is_trivial_zero(): |
| 935 | return z |
| 936 | else: |
| 937 | if not z: |
| 938 | return z |
| 939 | |
| 940 | return None # leaves the expression unevaluated |
| 941 | |
| 942 | def _evalf_(self, z, parent=None): |
| 943 | """ |
| 944 | EXAMPLES:: |
| 945 | |
| 946 | sage: N(sinh_integral(1e-10), digits=30) |
| 947 | 1.00000000000000003643219731550e-10 |
| 948 | sage: sinh_integral(ComplexField(100)(I)) |
| 949 | 0.94608307036718301494135331382*I |
| 950 | |
| 951 | """ |
| 952 | import mpmath |
| 953 | return mpmath_utils_call(mpmath.shi, z, parent=parent) |
| 954 | |
| 955 | def _derivative_(self, z, diff_param=None): |
| 956 | """ |
| 957 | The derivative of `\operatorname{Shi}(z)` is `\sinh(z)/z`. |
| 958 | |
| 959 | EXAMPLES:: |
| 960 | |
| 961 | sage: x = var('x') |
| 962 | sage: f = sinh_integral(x) |
| 963 | sage: f.diff(x) |
| 964 | sinh(x)/x |
| 965 | |
| 966 | sage: f = sinh_integral(ln(x)) |
| 967 | sage: f.diff(x) |
| 968 | sinh(log(x))/(x*log(x)) |
| 969 | |
| 970 | """ |
| 971 | return sinh(z)/z |
| 972 | |
| 973 | Shi = sinh_integral = Function_sinh_integral() |
| 974 | |
| 975 | |
| 976 | class Function_cosh_integral(BuiltinFunction): |
| 977 | r""" |
| 978 | The trigonometric integral `\operatorname{Chi}(z)` defined by |
| 979 | |
| 980 | .. math:: |
| 981 | |
| 982 | \operatorname{Chi}(z) = \gamma + \log(z) + \int_0^z \frac{\cosh(t)-1}{t}\; dt, |
| 983 | |
| 984 | see [AS]_ 5.2.4. |
| 985 | |
| 986 | EXAMPLES: |
| 987 | |
| 988 | Numerical evaluation for real and complex arguments is handled using mpmath:: |
| 989 | |
| 990 | sage: cosh_integral(1.0) |
| 991 | 0.837866940980208 |
| 992 | |
| 993 | The alias `Chi` can be used instead of `cosh_integral`:: |
| 994 | |
| 995 | sage: Chi(1.0) |
| 996 | 0.837866940980208 |
| 997 | |
| 998 | Here is an example from the mpmath documentation:: |
| 999 | |
| 1000 | sage: f(x) = cosh_integral(x) |
| 1001 | sage: find_root(f, 0.1, 1.0) |
| 1002 | 0.5238225713894826 |
| 1003 | |
| 1004 | Compare ``cosh_integral(3.0)`` to the definition of the value using |
| 1005 | numerical integration:: |
| 1006 | |
| 1007 | sage: N(euler_gamma + log(3.0) + integrate((cosh(x)-1)/x, x, 0, 3.0) - |
| 1008 | ... cosh_integral(3.0)) < 1e-14 |
| 1009 | True |
| 1010 | |
| 1011 | Arbitrary precision and complex arguments are handled:: |
| 1012 | |
| 1013 | sage: N(cosh_integral(3), digits=30) |
| 1014 | 4.96039209476560976029791763669 |
| 1015 | sage: cosh_integral(ComplexField(100)(3+I)) |
| 1016 | 3.9096723099686417127843516794 + 3.0547519627014217273323873274*I |
| 1017 | |
| 1018 | The limit of `\operatorname{Chi}(z)` as `z \to \infty` is `\infty`:: |
| 1019 | |
| 1020 | sage: N(cosh_integral(Infinity)) |
| 1021 | +infinity |
| 1022 | |
| 1023 | Symbolic derivatives and integrals are handled by Sage and Maxima:: |
| 1024 | |
| 1025 | sage: x = var('x') |
| 1026 | sage: f = cosh_integral(x) |
| 1027 | sage: f.diff(x) |
| 1028 | cosh(x)/x |
| 1029 | |
| 1030 | sage: f.integrate(x) |
| 1031 | x*cosh_integral(x) - sinh(x) |
| 1032 | |
| 1033 | ALGORITHM: |
| 1034 | |
| 1035 | Numerical evaluation is handled using mpmath, but symbolics are handled |
| 1036 | by Sage and Maxima. |
| 1037 | |
| 1038 | REFERENCES: |
| 1039 | |
| 1040 | - http://en.wikipedia.org/wiki/Trigonometric_integral |
| 1041 | - mpmath documentation: `chi`_ |
| 1042 | |
| 1043 | .. _`chi`: http://mpmath.googlecode.com/svn/trunk/doc/build/functions/expintegrals.html#chi |
| 1044 | |
| 1045 | """ |
| 1046 | def __init__(self): |
| 1047 | """ |
| 1048 | See the docstring for ``Function_cosh_integral``. |
| 1049 | |
| 1050 | EXAMPLES:: |
| 1051 | |
| 1052 | sage: cosh_integral(1) |
| 1053 | cosh_integral(1) |
| 1054 | |
| 1055 | """ |
| 1056 | BuiltinFunction.__init__(self, "cosh_integral", nargs=1, |
| 1057 | latex_name=r'\operatorname{Chi}', |
| 1058 | conversions=dict(maxima='expintegral_chi')) |
| 1059 | |
| 1060 | def _eval_(self, z): |
| 1061 | """ |
| 1062 | EXAMPLES:: |
| 1063 | |
| 1064 | sage: z = var('z') |
| 1065 | sage: cosh_integral(z) |
| 1066 | cosh_integral(z) |
| 1067 | sage: cosh_integral(3.0) |
| 1068 | 4.96039209476561 |
| 1069 | |
| 1070 | """ |
| 1071 | if not isinstance(z, Expression) and is_inexact(z): |
| 1072 | return self._evalf_(z, parent(z)) |
| 1073 | |
| 1074 | return None |
| 1075 | |
| 1076 | def _evalf_(self, z, parent=None): |
| 1077 | """ |
| 1078 | EXAMPLES:: |
| 1079 | |
| 1080 | sage: N(cosh_integral(1e-10), digits=30) |
| 1081 | -22.4486352650389235918737540487 |
| 1082 | sage: cosh_integral(ComplexField(100)(I)) |
| 1083 | 0.33740392290096813466264620389 + 1.5707963267948966192313216916*I |
| 1084 | |
| 1085 | """ |
| 1086 | import mpmath |
| 1087 | return mpmath_utils_call(mpmath.chi, z, parent=parent) |
| 1088 | |
| 1089 | def _derivative_(self, z, diff_param=None): |
| 1090 | """ |
| 1091 | The derivative of `\operatorname{Chi}(z)` is `\cosh(z)/z`. |
| 1092 | |
| 1093 | EXAMPLES:: |
| 1094 | |
| 1095 | sage: x = var('x') |
| 1096 | sage: f = cosh_integral(x) |
| 1097 | sage: f.diff(x) |
| 1098 | cosh(x)/x |
| 1099 | |
| 1100 | sage: f = cosh_integral(ln(x)) |
| 1101 | sage: f.diff(x) |
| 1102 | cosh(log(x))/(x*log(x)) |
| 1103 | |
| 1104 | """ |
| 1105 | return cosh(z)/z |
| 1106 | |
| 1107 | Chi = cosh_integral = Function_cosh_integral() |
| 1108 | |
| 1109 | |
| 1110 | ################################################################### |
| 1111 | ## Code below here was moved from sage/functions/transcendental.py |
| 1112 | ## This occured as part of Trac #11143. |
| 1113 | ################################################################### |
| 1114 | # |
| 1115 | # This class has a name which is not specific enough |
| 1116 | # see Function_exp_integral_e above, for example, which |
| 1117 | # is the "generalized" exponential integral function. We |
| 1118 | # are leaving the name the same for backwards compatibility |
| 1119 | # purposes. |
| 1120 | class Function_exp_integral(BuiltinFunction): |
| 1121 | r""" |
| 1122 | The generalized complex exponential integral Ei(z) defined by |
| 1123 | |
| 1124 | .. math:: |
| 1125 | |
| 1126 | \operatorname{Ei}(x) = \int_{-\infty}^x \frac{e^t}{t}\; dt |
| 1127 | |
| 1128 | for x > 0 and for complex arguments by analytic continuation, |
| 1129 | see [AS]_ 5.1.2. |
| 1130 | |
| 1131 | EXAMPLES:: |
| 1132 | |
| 1133 | sage: Ei(10) |
| 1134 | Ei(10) |
| 1135 | sage: Ei(I) |
| 1136 | Ei(I) |
| 1137 | sage: Ei(3+I) |
| 1138 | Ei(I + 3) |
| 1139 | sage: Ei(1.3) |
| 1140 | 2.72139888023202 |
| 1141 | |
| 1142 | The branch cut for this function is along the negative real axis:: |
| 1143 | |
| 1144 | sage: Ei(-3 + 0.1*I) |
| 1145 | -0.0129379427181693 + 3.13993830250942*I |
| 1146 | sage: Ei(-3 - 0.1*I) |
| 1147 | -0.0129379427181693 - 3.13993830250942*I |
| 1148 | |
| 1149 | ALGORITHM: Uses mpmath. |
| 1150 | |
| 1151 | """ |
| 1152 | def __init__(self): |
| 1153 | """ |
| 1154 | Return the value of the complex exponential integral Ei(z) at a |
| 1155 | complex number z. |
| 1156 | |
| 1157 | EXAMPLES:: |
| 1158 | |
| 1159 | sage: Ei(10) |
| 1160 | Ei(10) |
| 1161 | sage: Ei(I) |
| 1162 | Ei(I) |
| 1163 | sage: Ei(3+I) |
| 1164 | Ei(I + 3) |
| 1165 | sage: Ei(1.3) |
| 1166 | 2.72139888023202 |
| 1167 | |
| 1168 | The branch cut for this function is along the negative real axis:: |
| 1169 | |
| 1170 | sage: Ei(-3 + 0.1*I) |
| 1171 | -0.0129379427181693 + 3.13993830250942*I |
| 1172 | sage: Ei(-3 - 0.1*I) |
| 1173 | -0.0129379427181693 - 3.13993830250942*I |
| 1174 | |
| 1175 | ALGORITHM: Uses mpmath. |
| 1176 | """ |
| 1177 | BuiltinFunction.__init__(self, "Ei", |
| 1178 | conversions=dict(maxima='expintegral_ei')) |
| 1179 | |
| 1180 | def _eval_(self, x ): |
| 1181 | """ |
| 1182 | EXAMPLES:: |
| 1183 | |
| 1184 | sage: Ei(10) |
| 1185 | Ei(10) |
| 1186 | sage: Ei(I) |
| 1187 | Ei(I) |
| 1188 | sage: Ei(1.3) |
| 1189 | 2.72139888023202 |
| 1190 | sage: Ei(10r) |
| 1191 | Ei(10) |
| 1192 | sage: Ei(1.3r) |
| 1193 | 2.7213988802320235 |
| 1194 | """ |
| 1195 | if not isinstance(x, Expression) and is_inexact(x): |
| 1196 | return self._evalf_(x, parent(x)) |
| 1197 | return None |
| 1198 | |
| 1199 | def _evalf_(self, x, parent=None): |
| 1200 | """ |
| 1201 | EXAMPLES:: |
| 1202 | |
| 1203 | sage: Ei(10).n() |
| 1204 | 2492.22897624188 |
| 1205 | sage: Ei(20).n() |
| 1206 | 2.56156526640566e7 |
| 1207 | sage: Ei(I).n() |
| 1208 | 0.337403922900968 + 2.51687939716208*I |
| 1209 | sage: Ei(3+I).n() |
| 1210 | 7.82313467600158 + 6.09751978399231*I |
| 1211 | """ |
| 1212 | import mpmath |
| 1213 | return mpmath_utils_call(mpmath.ei, x, parent=parent) |
| 1214 | |
| 1215 | def __call__(self, x, prec=None, coerce=True, hold=False ): |
| 1216 | """ |
| 1217 | Note that the ``prec`` argument is deprecated. The precision for |
| 1218 | the result is deduced from the precision of the input. Convert |
| 1219 | the input to a higher precision explicitly if a result with higher |
| 1220 | precision is desired. |
| 1221 | |
| 1222 | EXAMPLES:: |
| 1223 | |
| 1224 | sage: t = Ei(RealField(100)(2.5)); t |
| 1225 | 7.0737658945786007119235519625 |
| 1226 | sage: t.prec() |
| 1227 | 100 |
| 1228 | |
| 1229 | sage: Ei(1.1, prec=300) |
| 1230 | 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. |
| 1231 | 2.16737827956340306615064476647912607220394065907142504328679588538509331805598360907980986 |
| 1232 | """ |
| 1233 | if prec is not None: |
| 1234 | from sage.misc.misc import deprecation |
| 1235 | 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.") |
| 1236 | |
| 1237 | import mpmath |
| 1238 | return mpmath_utils_call(mpmath.ei, x, prec=prec) |
| 1239 | |
| 1240 | return BuiltinFunction.__call__(self, x, coerce=coerce, hold=hold) |
| 1241 | |
| 1242 | def _derivative_(self, x, diff_param=None): |
| 1243 | """ |
| 1244 | EXAMPLES:: |
| 1245 | |
| 1246 | sage: Ei(x).diff(x) |
| 1247 | e^x/x |
| 1248 | sage: Ei(x).diff(x).subs(x=1) |
| 1249 | e |
| 1250 | sage: Ei(x^2).diff(x) |
| 1251 | 2*e^(x^2)/x |
| 1252 | sage: f = function('f') |
| 1253 | sage: Ei(f(x)).diff(x) |
| 1254 | e^f(x)*D[0](f)(x)/f(x) |
| 1255 | """ |
| 1256 | return exp(x)/x |
| 1257 | |
| 1258 | Ei = exp_integral_ei = Function_exp_integral() |
| 1259 | |
| 1260 | |
| 1261 | # moved here from sage/functions/transcendental.py |
| 1262 | def exponential_integral_1(x, n=0): |
| 1263 | r""" |
| 1264 | Returns the exponential integral `E_1(x)`. If the optional |
| 1265 | argument `n` is given, computes list of the first |
| 1266 | `n` values of the exponential integral |
| 1267 | `E_1(x m)`. |
| 1268 | |
| 1269 | The exponential integral `E_1(x)` is |
| 1270 | |
| 1271 | .. math:: |
| 1272 | |
| 1273 | E_1(x) = \int_{x}^{\infty} e^{-t}/t dt |
| 1274 | |
| 1275 | |
| 1276 | |
| 1277 | INPUT: |
| 1278 | |
| 1279 | |
| 1280 | - ``x`` - a positive real number |
| 1281 | |
| 1282 | - ``n`` - (default: 0) a nonnegative integer; if |
| 1283 | nonzero, then return a list of values E_1(x\*m) for m = |
| 1284 | 1,2,3,...,n. This is useful, e.g., when computing derivatives of |
| 1285 | L-functions. |
| 1286 | |
| 1287 | |
| 1288 | OUTPUT: |
| 1289 | |
| 1290 | |
| 1291 | - ``float`` - if n is 0 (the default) or |
| 1292 | |
| 1293 | - ``list`` - list of floats if n 0 |
| 1294 | |
| 1295 | |
| 1296 | EXAMPLES:: |
| 1297 | |
| 1298 | sage: exponential_integral_1(2) |
| 1299 | 0.04890051070806112 |
| 1300 | sage: w = exponential_integral_1(2,4); w |
| 1301 | [0.04890051070806112, 0.0037793524098489067, 0.00036008245216265873, 3.7665622843924751e-05] # 32-bit |
| 1302 | [0.04890051070806112, 0.0037793524098489063, 0.00036008245216265873, 3.7665622843924534e-05] # 64-bit |
| 1303 | sage: exponential_integral_1(0) |
| 1304 | +Infinity |
| 1305 | |
| 1306 | IMPLEMENTATION: We use the PARI C-library functions eint1 and |
| 1307 | veceint1. |
| 1308 | |
| 1309 | REFERENCE: |
| 1310 | |
| 1311 | - See page 262, Prop 5.6.12, of Cohen's book "A Course in |
| 1312 | Computational Algebraic Number Theory". |
| 1313 | |
| 1314 | REMARKS: When called with the optional argument n, the PARI |
| 1315 | C-library is fast for values of n up to some bound, then very very |
| 1316 | slow. For example, if x=5, then the computation takes less than a |
| 1317 | second for n=800000, and takes "forever" for n=900000. |
| 1318 | """ |
| 1319 | if isinstance(x, Expression): |
| 1320 | if x.is_trivial_zero(): |
| 1321 | from sage.rings.infinity import Infinity |
| 1322 | return Infinity |
| 1323 | else: |
| 1324 | raise NotImplementedError("Use the symbolic exponential integral " + |
| 1325 | "function: exp_integral_e1.") |
| 1326 | elif not is_inexact(x): # x is exact and not an expression |
| 1327 | if not x: # test if exact x == 0 quickly |
| 1328 | from sage.rings.infinity import Infinity |
| 1329 | return Infinity |
| 1330 | |
| 1331 | # else x is in not an exact 0 |
| 1332 | from sage.libs.pari.all import pari |
| 1333 | if n <= 0: |
| 1334 | return float(pari(x).eint1()) |
| 1335 | else: |
| 1336 | return [float(z) for z in pari(x).eint1(n)] |