doc/info/Integration.texi.m4

   1 @c -*- mode: texinfo -*-
   2 @menu
   3 * Introduction to Integration::
   4 * Functions and Variables for Integration::
   5 * Introduction to QUADPACK::
   6 * Functions and Variables for QUADPACK::
   7 @end menu
   8
   9 @c -----------------------------------------------------------------------------
  10 @node Introduction to Integration, Functions and Variables for Integration, Integration, Integration
  11 @section Introduction to Integration
  12 @c -----------------------------------------------------------------------------
  13
  14 Maxima has several routines for handling integration.
  15 The @mref{integrate} function makes use of most of them.  There is also the
  16 @mref{antid} package, which handles an unspecified function (and its
  17 derivatives, of course).  For numerical uses,
  18 there is a set of adaptive integrators from QUADPACK, named @mrefcomma{quad_qag}
  19 @mrefcomma{quad_qags} etc., which are described under the heading @code{QUADPACK}.
  20 Hypergeometric functions are being worked on,
  21 see @mref{specint} for details.
  22 Generally speaking, Maxima only handles integrals which are
  23 integrable in terms of the "elementary functions" (rational functions,
  24 trigonometrics, logs, exponentials, radicals, etc.) and a few
  25 extensions (error function, dilogarithm).  It does not handle
  26 integrals in terms of unknown functions such as @code{g(x)} and @code{h(x)}.
  27
  28 @c end concepts Integration
  29
  30 @c -----------------------------------------------------------------------------
  31 @node Functions and Variables for Integration, Introduction to QUADPACK, Introduction to Integration, Integration
  32 @section Functions and Variables for Integration
  33 @c -----------------------------------------------------------------------------
  34
  35 @c NEEDS WORK
  36
  37 @c -----------------------------------------------------------------------------
  38 @anchor{changevar}
  39 @deffn {Function} changevar (@var{expr}, @var{f(x,y)}, @var{y}, @var{x})
  40
  41 Makes the change of variable given by @code{@var{f(x,y)} = 0} in all integrals
  42 occurring in @var{expr} with integration with respect to @var{x}.
  43 The new variable is @var{y}.
  44
  45 The change of variable can also be written @code{@var{f(x)} = @var{g(y)}}.
  46
  47 @c HMM, THIS EXAMPLE YIELDS A CORRECT BUT SLIGHTLY STRANGE RESULT...
  48 @c ===beg===
  49 @c assume(a > 0)$
  50 @c 'integrate (%e**sqrt(a*y), y, 0, 4);
  51 @c changevar (%, y-z^2/a, z, y);
  52 @c ===end===
  53 @example
  54 (%i1) assume(a > 0)$
  55 @group
  56 (%i2) 'integrate (%e**sqrt(a*y), y, 0, 4);
  57                       4
  58                      /
  59                      [    sqrt(a) sqrt(y)
  60 (%o2)                I  %e                dy
  61                      ]
  62                      /
  63                       0
  64 @end group
  65 @group
  66 (%i3) changevar (%, y-z^2/a, z, y);
  67                       0
  68                      /
  69                      [                abs(z)
  70                    2 I            z %e       dz
  71                      ]
  72                      /
  73                       - 2 sqrt(a)
  74 (%o3)            - ----------------------------
  75                                 a
  76 @end group
  77 @end example
  78
  79 An expression containing a noun form, such as the instances of @code{'integrate}
  80 above, may be evaluated by @code{ev} with the @code{nouns} flag.
  81 For example, the expression returned by @code{changevar} above may be evaluated
  82 by @code{ev (%o3, nouns)}.
  83
  84 @code{changevar} may also be used to make changes in the indices of a sum or
  85 product.  However, it must be realized that when a change is made in a
  86 sum or product, this change must be a shift, i.e., @code{i = j+ ...}, not a
  87 higher degree function.  E.g.,
  88
  89 @c ===beg===
  90 @c sum (a[i]*x^(i-2), i, 0, inf);
  91 @c changevar (%, i-2-n, n, i);
  92 @c ===end===
  93 @example
  94 @group
  95 (%i4) sum (a[i]*x^(i-2), i, 0, inf);
  96                          inf
  97                          ====
  98                          \         i - 2
  99 (%o4)                     >    a  x
 100                          /      i
 101                          ====
 102                          i = 0
 103 @end group
 104 @group
 105 (%i5) changevar (%, i-2-n, n, i);
 106                         inf
 107                         ====
 108                         \               n
 109 (%o5)                    >      a      x
 110                         /        n + 2
 111                         ====
 112                         n = - 2
 113 @end group
 114 @end example
 115
 116 @opencatbox{Categories:}
 117 @category{Integral calculus}
 118 @closecatbox
 119 @end deffn
 120
 121 @c THIS ITEM IS A MESS, BUT DON'T BOTHER TO CLEAN IT UP:
 122 @c THE GAUSS-KRONROD FUNCTIONS (QUADPACK) MAKE THIS OBSOLETE
 123
 124 @c -----------------------------------------------------------------------------
 125 @anchor{dblint}
 126 @deffn {Function} dblint (@var{f}, @var{r}, @var{s}, @var{a}, @var{b})
 127
 128 A double-integral routine which was written in
 129 top-level Maxima and then translated and compiled to machine code.
 130 Use @code{load ("dblint")} to access this package.  It uses the Simpson's rule
 131 method in both the x and y directions to calculate
 132
 133 @tex
 134 $$\int_a^b \int_{r\left(x\right)}^{s\left(x\right)} f\left(x,y\right) \, dy \, dx.$$
 135 @end tex
 136 @ifnottex
 137 @example
 138 @group
 139 /b /s(x)
 140 |  |
 141 |  |    f(x,y) dy dx
 142 |  |
 143 /a /r(x)
 144 @end group
 145 @end example
 146 @end ifnottex
 147
 148 The function @var{f} must be a translated or compiled function of two variables,
 149 and @var{r} and @var{s} must each be a translated or compiled function of one
 150 variable, while @var{a} and @var{b} must be floating point numbers.  The routine
 151 has two global variables which determine the number of divisions of the x and y
 152 intervals: @code{dblint_x} and @code{dblint_y}, both of which are initially 10,
 153 and can be changed independently to other integer values (there are
 154 @code{2*dblint_x+1} points computed in the x direction, and @code{2*dblint_y+1}
 155 in the y direction).  The routine subdivides the X axis and then for each value
 156 of X it first computes @code{@var{r}(x)} and @code{@var{s}(x)}; then the Y axis
 157 between @code{@var{r}(x)} and @code{@var{s}(x)} is subdivided and the integral
 158 along the Y axis is performed using Simpson's rule; then the integral along the
 159 X axis is done using Simpson's rule with the function values being the
 160 Y-integrals.  This procedure may be numerically unstable for a great variety of
 161 reasons, but is reasonably fast: avoid using it on highly oscillatory functions
 162 and functions with singularities (poles or branch points in the region).  The Y
 163 integrals depend on how far apart @code{@var{r}(x)} and @code{@var{s}(x)} are,
 164 so if the distance @code{@var{s}(x) - @var{r}(x)} varies rapidly with X, there
 165 may be substantial errors arising from truncation with different step-sizes in
 166 the various Y integrals.  One can increase @code{dblint_x} and @code{dblint_y}
 167 in an effort to improve the coverage of the region, at the expense of
 168 computation time.  The function values are not saved, so if the function is very
 169 time-consuming, you will have to wait for re-computation if you change anything
 170 (sorry).  It is required that the functions @var{f}, @var{r}, and @var{s} be
 171 either translated or compiled prior to calling @code{dblint}.  This will result
 172 in orders of magnitude speed improvement over interpreted code in many cases!
 173
 174 @code{demo ("dblint")} executes a demonstration of @code{dblint} applied to an
 175 example problem.
 176 @c demo (dblint_1) FAILS WITH Could not find `fltdfnk.mc' -- DON'T BOTHER TO MENTION IT. !!!
 177 @c @code{demo (dblint_1)} executes another demonstration.
 178
 179 @opencatbox{Categories:}
 180 @category{Integral calculus}
 181 @closecatbox
 182 @end deffn
 183
 184 @c -----------------------------------------------------------------------------
 185 @anchor{defint}
 186 @deffn {Function} defint (@var{expr}, @var{x}, @var{a}, @var{b})
 187
 188 Attempts to compute a definite integral.  @code{defint} is called by
 189 @code{integrate} when limits of integration are specified, i.e., when
 190 @code{integrate} is called as
 191 @code{integrate (@var{expr}, @var{x}, @var{a}, @var{b})}.
 192 Thus from the user's point of view, it is sufficient to call @code{integrate}.
 193 @c SHOULD WE BOTHER TO DOCUMENT defint ??? NO FUNCTIONALITY HERE THAT IS NOT ALREADY PRESENT IN integrate !!!
 194
 195 @code{defint} returns a symbolic expression, either the computed integral or the
 196 noun form of the integral.  See @mref{quad_qag} and related functions for
 197 numerical approximation of definite integrals.
 198
 199 @opencatbox{Categories:}
 200 @category{Integral calculus}
 201 @closecatbox
 202 @end deffn
 203
 204 @c -----------------------------------------------------------------------------
 205 @anchor{erfflag}
 206 @defvr {Option variable} erfflag
 207 Default value: @code{true}
 208
 209 When @code{erfflag} is @code{false}, prevents @code{risch} from introducing the
 210 @code{erf} function in the answer if there were none in the integrand to
 211 begin with.
 212
 213 @opencatbox{Categories:}
 214 @category{Integral calculus}
 215 @closecatbox
 216 @end defvr
 217
 218 @c NEEDS WORK
 219
 220 @c -----------------------------------------------------------------------------
 221 @anchor{ilt}
 222 @deffn {Function} ilt (@var{expr}, @var{s}, @var{t})
 223
 224 Computes the inverse Laplace transform of @var{expr} with
 225 respect to @var{s} and parameter @var{t}. @var{expr} must be a ratio of
 226 polynomials whose denominator has only linear and quadratic factors;
 227 there is an extension of @code{ilt}, called @mref{pwilt} (Piece-Wise
 228 Inverse Laplace Transform) that handles several other cases where
 229 @code{ilt} fails.
 230
 231 By using the functions @code{laplace} and @code{ilt} together with the
 232 @code{solve} or @code{linsolve} functions the user can solve a single
 233 differential or convolution integral equation or a set of them.
 234
 235 @c ===beg===
 236 @c 'integrate (sinh(a*x)*f(t-x), x, 0, t) + b*f(t) = t**2;
 237 @c laplace (%, t, s);
 238 @c linsolve ([%], ['laplace(f(t), t, s)]);
 239 @c ilt (rhs (first (%)), s, t);
 240 @c input:pos;
 241 @c ===end===
 242 @example
 243 @group
 244 (%i1) 'integrate (sinh(a*x)*f(t-x), x, 0, t) + b*f(t) = t**2;
 245               t
 246              /
 247              [                                    2
 248 (%o1)        I  f(t - x) sinh(a x) dx + b f(t) = t
 249              ]
 250              /
 251               0
 252 @end group
 253 @group
 254 (%i2) laplace (%, t, s);
 255                                a laplace(f(t), t, s)   2
 256 (%o2)  b laplace(f(t), t, s) + --------------------- = --
 257                                        2    2           3
 258                                       s  - a           s
 259 @end group
 260 @group
 261 (%i3) linsolve ([%], ['laplace(f(t), t, s)]);
 262                                         2      2
 263                                      2 s  - 2 a
 264 (%o3)     [laplace(f(t), t, s) = --------------------]
 265                                     5         2     3
 266                                  b s  + (a - a  b) s
 267 @end group
 268 @group
 269 (%i4) ilt (rhs (first (%)), s, t);
 270 Is  a b (a b - 1)  positive, negative, or zero?
 271
 272 pos;
 273                sqrt(a b (a b - 1)) t
 274         2 cosh(---------------------)       2
 275                          b               a t
 276 (%o4) - ----------------------------- + -------
 277               3  2      2               a b - 1
 278              a  b  - 2 a  b + a
 279
 280                                                        2
 281                                              + ------------------
 282                                                 3  2      2
 283                                                a  b  - 2 a  b + a
 284 @end group
 285 @end example
 286
 287 @opencatbox{Categories:}
 288 @category{Laplace transform}
 289 @closecatbox
 290 @end deffn
 291
 292 @c -----------------------------------------------------------------------------
 293 @anchor{intanalysis}
 294 @defvr {Option variable} intanalysis
 295 Default value: @code{true}
 296
 297 When @code{true}, definite integration tries to find poles in the integrand in
 298 the interval of integration.  If there are, then the integral is evaluated
 299 appropriately as a principal value integral.  If intanalysis is @code{false},
 300 this check is not performed and integration is done assuming there are no poles.
 301
 302 See also @mrefdot{ldefint}
 303
 304 Examples:
 305
 306 Maxima can solve the following integrals, when @mref{intanalysis} is set to
 307 @code{false}:
 308
 309 @c ===beg===
 310 @c integrate(1/(sqrt(x+1)+1),x,0,1);
 311 @c integrate(1/(sqrt(x)+1),x,0,1),intanalysis:false;
 312 @c integrate(cos(a)/sqrt((tan(a))^2+1),a,-%pi/2,%pi/2);
 313 @c intanalysis:false$
 314 @c integrate(cos(a)/sqrt((tan(a))^2 +1),a,-%pi/2,%pi/2);
 315 @c ===end===
 316 @example
 317 (%i1) integrate(1/(sqrt(x)+1),x,0,1);
 318                                 1
 319                                /
 320                                [       1
 321 (%o1)                          I  ----------- dx
 322                                ]  sqrt(x) + 1
 323                                /
 324                                 0
 325
 326 (%i2) integrate(1/(sqrt(x)+1),x,0,1),intanalysis:false;
 327 (%o2)                            2 - 2 log(2)
 328
 329 (%i3) integrate(cos(a)/sqrt((tan(a))^2 +1),a,-%pi/2,%pi/2);
 330 The number 1 isn't in the domain of atanh
 331  -- an error. To debug this try: debugmode(true);
 332
 333 (%i4) intanalysis:false$
 334 (%i5) integrate(cos(a)/sqrt((tan(a))^2+1),a,-%pi/2,%pi/2);
 335                                       %pi
 336 (%o5)                                 ---
 337                                        2
 338 @end example
 339
 340 @opencatbox{Categories:}
 341 @category{Integral calculus}
 342 @closecatbox
 343 @end defvr
 344
 345 @c -----------------------------------------------------------------------------
 346 @anchor{integrate}
 347 @deffn  {Function} integrate @
 348 @fname{integrate} (@var{expr}, @var{x}) @
 349 @fname{integrate} (@var{expr}, @var{x}, @var{a}, @var{b})
 350
 351 Attempts to symbolically compute the integral of @var{expr} with respect to
 352 @var{x}.  @code{integrate (@var{expr}, @var{x})} is an indefinite integral,
 353 while @code{integrate (@var{expr}, @var{x}, @var{a}, @var{b})} is a definite
 354 integral, with limits of integration @var{a} and @var{b}.  The limits should
 355 not contain @var{x}, although @code{integrate} does not enforce this
 356 restriction.  @var{a} need not be less than @var{b}.
 357 If @var{b} is equal to @var{a}, @code{integrate} returns zero.
 358
 359 See @mref{quad_qag} and related functions for numerical approximation of
 360 definite integrals.  See @mref{residue} for computation of residues
 361 (complex integration).  See @mref{antid} for an alternative means of computing
 362 indefinite integrals.
 363
 364 The integral (an expression free of @code{integrate}) is returned if
 365 @code{integrate} succeeds.  Otherwise the return value is
 366 the noun form of the integral (the quoted operator @code{'integrate})
 367 or an expression containing one or more noun forms.
 368 The noun form of @code{integrate} is displayed with an integral sign.
 369
 370 In some circumstances it is useful to construct a noun form by hand, by quoting
 371 @code{integrate} with a single quote, e.g.,
 372 @code{'integrate (@var{expr}, @var{x})}.  For example, the integral may depend
 373 on some parameters which are not yet computed.
 374 The noun may be applied to its arguments by @code{ev (@var{i}, nouns)}
 375 where @var{i} is the noun form of interest.
 376
 377 @c BEGIN EXPOSITION ON HEURISTICS
 378 @code{integrate} handles definite integrals separately from indefinite, and
 379 employs a range of heuristics to handle each case.  Special cases of definite
 380 integrals include limits of integration equal to zero or infinity (@mref{inf} or
 381 @mref{minf}), trigonometric functions with limits of integration equal to zero
 382 and @code{%pi} or @code{2 %pi}, rational functions, integrals related to the
 383 definitions of the @mref{beta} and @mref{psi} functions, and some logarithmic
 384 and trigonometric integrals.  Processing rational functions may include
 385 computation of residues.  If an applicable special case is not found, an attempt
 386 will be made to compute the indefinite integral and evaluate it at the limits of
 387 integration.  This may include taking a limit as a limit of integration goes to
 388 infinity or negative infinity; see also @mrefdot{ldefint}
 389
 390 Special cases of indefinite integrals include trigonometric functions,
 391 exponential and logarithmic functions,
 392 and rational functions.
 393 @code{integrate} may also make use of a short table of elementary integrals.
 394
 395 @code{integrate} may carry out a change of variable
 396 if the integrand has the form @code{f(g(x)) * diff(g(x), x)}.
 397 @code{integrate} attempts to find a subexpression @code{g(x)} such that
 398 the derivative of @code{g(x)} divides the integrand.
 399 This search may make use of derivatives defined by the @code{gradef} function.
 400 See also @mref{changevar} and @mrefdot{antid}
 401
 402 If none of the preceding heuristics find the indefinite integral, the Risch
 403 algorithm is executed.  The flag @mref{risch} may be set as an @mrefcomma{evflag}
 404 in a call to @code{ev} or on the command line, e.g.,
 405 @code{ev (integrate (@var{expr}, @var{x}), risch)} or
 406 @code{integrate (@var{expr}, @var{x}), risch}.  If @code{risch} is present,
 407 @code{integrate} calls the @mref{risch} function without attempting heuristics
 408 first.  See also @mrefdot{risch}
 409 @c END EXPOSITION ON HEURISTICS
 410
 411 @code{integrate} works only with functional relations represented explicitly
 412 with the @code{f(x)} notation.  @code{integrate} does not respect implicit
 413 dependencies established by the @mref{depends} function.
 414
 415 @code{integrate} may need to know some property of a parameter in the integrand.
 416 @code{integrate} will first consult the @mref{assume} database,
 417 and, if the variable of interest is not there,
 418 @code{integrate} will ask the user.
 419 Depending on the question,
 420 suitable responses are @code{yes;} or @code{no;},
 421 or @code{pos;}, @code{zero;}, or @code{neg;}.
 422
 423 @code{integrate} is not, by default, declared to be linear.  See @code{declare}
 424 and @code{linear}.
 425
 426 @code{integrate} attempts integration by parts only in a few special cases.
 427
 428 Examples:
 429
 430 @itemize @bullet
 431 @item
 432 Elementary indefinite and definite integrals.
 433
 434 @c ===beg===
 435 @c integrate (sin(x)^3, x);
 436 @c integrate (x/ sqrt (b^2 - x^2), x);
 437 @c integrate (cos(x)^2 * exp(x), x, 0, %pi);
 438 @c integrate (x^2 * exp(-x^2), x, minf, inf);
 439 @c ===end===
 440 @example
 441 @group
 442 (%i1) integrate (sin(x)^3, x);
 443                            3
 444                         cos (x)
 445 (%o1)                   ------- - cos(x)
 446                            3
 447 @end group
 448 @group
 449 (%i2) integrate (x/ sqrt (b^2 - x^2), x);
 450                                  2    2
 451 (%o2)                    - sqrt(b  - x )
 452 @end group
 453 @group
 454 (%i3) integrate (cos(x)^2 * exp(x), x, 0, %pi);
 455                                %pi
 456                            3 %e      3
 457 (%o3)                      ------- - -
 458                               5      5
 459 @end group
 460 @group
 461 (%i4) integrate (x^2 * exp(-x^2), x, minf, inf);
 462                             sqrt(%pi)
 463 (%o4)                       ---------
 464                                 2
 465 @end group
 466 @end example
 467
 468 @item
 469 Use of @code{assume} and interactive query.
 470
 471 @c ===beg===
 472 @c assume (a > 1)$
 473 @c integrate (x**a/(x+1)**(5/2), x, 0, inf);
 474 @c input:no;
 475 @c input:neg;
 476 @c ===end===
 477 @example
 478 (%i1) assume (a > 1)$
 479 @group
 480 (%i2) integrate (x**a/(x+1)**(5/2), x, 0, inf);
 481     2 a + 2
 482 Is  -------  an integer?
 483        5
 484
 485 no;
 486 Is  2 a - 3  positive, negative, or zero?
 487
 488 neg;
 489                                    3
 490 (%o2)                  beta(a + 1, - - a)
 491                                    2
 492 @end group
 493 @end example
 494
 495 @item
 496 Change of variable.  There are two changes of variable in this example:
 497 one using a derivative established by @mrefcomma{gradef} and one using the
 498 derivation @code{diff(r(x))} of an unspecified function @code{r(x)}.
 499
 500 @c ===beg===
 501 @c gradef (q(x), sin(x**2));
 502 @c diff (log (q (r (x))), x);
 503 @c integrate (%, x);
 504 @c ===end===
 505 @example
 506 @group
 507 (%i3) gradef (q(x), sin(x**2));
 508 (%o3)                         q(x)
 509 @end group
 510 @group
 511 (%i4) diff (log (q (r (x))), x);
 512                       d               2
 513                      (-- (r(x))) sin(r (x))
 514                       dx
 515 (%o4)                ----------------------
 516                             q(r(x))
 517 @end group
 518 @group
 519 (%i5) integrate (%, x);
 520 (%o5)                     log(q(r(x)))
 521 @end group
 522 @end example
 523
 524 @item
 525 Return value contains the @code{'integrate} noun form.  In this example, Maxima
 526 can extract one factor of the denominator of a rational function, but cannot
 527 factor the remainder or otherwise find its integral.  @mref{grind} shows the
 528 noun form @code{'integrate} in the result.  See also
 529 @mref{integrate_use_rootsof} for more on integrals of rational functions.
 530
 531 @c ===beg===
 532 @c expand ((x-4) * (x^3+2*x+1));
 533 @c integrate (1/%, x);
 534 @c grind (%);
 535 @c ===end===
 536 @example
 537 @group
 538 (%i1) expand ((x-4) * (x^3+2*x+1));
 539                     4      3      2
 540 (%o1)              x  - 4 x  + 2 x  - 7 x - 4
 541 @end group
 542 @group
 543 (%i2) integrate (1/%, x);
 544                               /  2
 545                               [ x  + 4 x + 18
 546                               I ------------- dx
 547                               ]  3
 548                  log(x - 4)   / x  + 2 x + 1
 549 (%o2)            ---------- - ------------------
 550                      73               73
 551 @end group
 552 @group
 553 (%i3) grind (%);
 554 log(x-4)/73-('integrate((x^2+4*x+18)/(x^3+2*x+1),x))/73$
 555 @end group
 556 @end example
 557
 558 @item
 559 Defining a function in terms of an integral.  The body of a function is not
 560 evaluated when the function is defined.  Thus the body of @code{f_1} in this
 561 example contains the noun form of @code{integrate}.  The quote-quote operator
 562 @code{'@w{}'} causes the integral to be evaluated, and the result becomes the
 563 body of @code{f_2}.
 564
 565 @c ===beg===
 566 @c f_1 (a) := integrate (x^3, x, 1, a);
 567 @c ev (f_1 (7), nouns);
 568 @c /* Note parentheses around integrate(...) here */      f_2 (a) := ''(integrate (x^3, x, 1, a));
 569 @c f_2 (7);
 570 @c ===end===
 571 @example
 572 @group
 573 (%i1) f_1 (a) := integrate (x^3, x, 1, a);
 574                                      3
 575 (%o1)           f_1(a) := integrate(x , x, 1, a)
 576 @end group
 577 @group
 578 (%i2) ev (f_1 (7), nouns);
 579 (%o2)                          600
 580 @end group
 581 @group
 582 (%i3) /* Note parentheses around integrate(...) here */
 583       f_2 (a) := ''(integrate (x^3, x, 1, a));
 584                                    4
 585                                   a    1
 586 (%o3)                   f_2(a) := -- - -
 587                                   4    4
 588 @end group
 589 @group
 590 (%i4) f_2 (7);
 591 (%o4)                          600
 592 @end group
 593 @end example
 594 @end itemize
 595
 596 @opencatbox{Categories:}
 597 @category{Integral calculus}
 598 @closecatbox
 599 @end deffn
 600
 601 @c -----------------------------------------------------------------------------
 602 @anchor{integration_constant}
 603 @defvr {System variable} integration_constant
 604 Default value: @code{%c}
 605
 606 When a constant of integration is introduced by indefinite integration of an
 607 equation, the name of the constant is constructed by concatenating
 608 @code{integration_constant} and @code{integration_constant_counter}.
 609
 610 @code{integration_constant} may be assigned any symbol.
 611
 612 Examples:
 613
 614 @c ===beg===
 615 @c integrate (x^2 = 1, x);
 616 @c integration_constant : 'k;
 617 @c integrate (x^2 = 1, x);
 618 @c ===end===
 619 @example
 620 @group
 621 (%i1) integrate (x^2 = 1, x);
 622                            3
 623                           x
 624 (%o1)                     -- = x + %c1
 625                           3
 626 @end group
 627 @group
 628 (%i2) integration_constant : 'k;
 629 (%o2)                           k
 630 @end group
 631 @group
 632 (%i3) integrate (x^2 = 1, x);
 633                             3
 634                            x
 635 (%o3)                      -- = x + k2
 636                            3
 637 @end group
 638 @end example
 639
 640 @opencatbox{Categories:}
 641 @category{Integral calculus}
 642 @closecatbox
 643 @end defvr
 644
 645 @c -----------------------------------------------------------------------------
 646 @anchor{integration_constant_counter}
 647 @defvr {System variable} integration_constant_counter
 648 Default value: 0
 649
 650 When a constant of integration is introduced by indefinite integration of an
 651 equation, the name of the constant is constructed by concatenating
 652 @code{integration_constant} and @code{integration_constant_counter}.
 653
 654 @code{integration_constant_counter} is incremented before constructing the next
 655 integration constant.
 656
 657 Examples:
 658
 659 @c ===beg===
 660 @c integrate (x^2 = 1, x);
 661 @c integrate (x^2 = 1, x);
 662 @c integrate (x^2 = 1, x);
 663 @c reset (integration_constant_counter);
 664 @c integrate (x^2 = 1, x);
 665 @c ===end===
 666 @example
 667 @group
 668 (%i1) integrate (x^2 = 1, x);
 669                            3
 670                           x
 671 (%o1)                     -- = x + %c1
 672                           3
 673 @end group
 674 @group
 675 (%i2) integrate (x^2 = 1, x);
 676                            3
 677                           x
 678 (%o2)                     -- = x + %c2
 679                           3
 680 @end group
 681 @group
 682 (%i3) integrate (x^2 = 1, x);
 683                            3
 684                           x
 685 (%o3)                     -- = x + %c3
 686                           3
 687 @end group
 688 @group
 689 (%i4) reset (integration_constant_counter);
 690 (%o4)            [integration_constant_counter]
 691 @end group
 692 @group
 693 (%i5) integrate (x^2 = 1, x);
 694                            3
 695                           x
 696 (%o5)                     -- = x + %c1
 697                           3
 698 @end group
 699 @end example
 700
 701 @opencatbox{Categories:}
 702 @category{Integral calculus}
 703 @closecatbox
 704 @end defvr
 705
 706 @c -----------------------------------------------------------------------------
 707 @anchor{integrate_use_rootsof}
 708 @defvr {Option variable} integrate_use_rootsof
 709 Default value: @code{false}
 710
 711 When @code{integrate_use_rootsof} is @code{true} and the denominator of
 712 a rational function cannot be factored, @mref{integrate} returns the integral
 713 in a form which is a sum over the roots (not yet known) of the denominator.
 714
 715 For example, with @code{integrate_use_rootsof} set to @code{false},
 716 @code{integrate} returns an unsolved integral of a rational function in noun
 717 form:
 718
 719 @c ===beg===
 720 @c integrate_use_rootsof: false$
 721 @c integrate (1/(1+x+x^5), x);
 722 @c ===end===
 723 @example
 724 (%i1) integrate_use_rootsof: false$
 725 @group
 726 (%i2) integrate (1/(1+x+x^5), x);
 727         /  2
 728         [ x  - 4 x + 5
 729         I ------------ dx                            2 x + 1
 730         ]  3    2                2            5 atan(-------)
 731         / x  - x  + 1       log(x  + x + 1)          sqrt(3)
 732 (%o2)   ----------------- - --------------- + ---------------
 733                 7                 14             7 sqrt(3)
 734 @end group
 735 @end example
 736
 737 Now we set the flag to be true and the unsolved part of the integral will be
 738 expressed as a summation over the roots of the denominator of the rational
 739 function:
 740
 741 @c ===beg===
 742 @c integrate_use_rootsof: true$
 743 @c integrate (1/(1+x+x^5), x);
 744 @c ===end===
 745 @example
 746 (%i3) integrate_use_rootsof: true$
 747 @group
 748 (%i4) integrate (1/(1+x+x^5), x);
 749       ====        2
 750       \       (%r4  - 4 %r4 + 5) log(x - %r4)
 751        >      -------------------------------
 752       /                    2
 753       ====            3 %r4  - 2 %r4
 754                         3      2
 755       %r4 in rootsof(%r4  - %r4  + 1, %r4)
 756 (%o4) ----------------------------------------------------------
 757                7
 758
 759                                                       2 x + 1
 760                                   2            5 atan(-------)
 761                              log(x  + x + 1)          sqrt(3)
 762                            - --------------- + ---------------
 763                                    14             7 sqrt(3)
 764 @end group
 765 @end example
 766
 767 Alternatively the user may compute the roots of the denominator separately,
 768 and then express the integrand in terms of these roots, e.g.,
 769 @code{1/((x - a)*(x - b)*(x - c))} or @code{1/((x^2 - (a+b)*x + a*b)*(x - c))}
 770 if the denominator is a cubic polynomial.
 771 Sometimes this will help Maxima obtain a more useful result.
 772
 773 @opencatbox{Categories:}
 774 @category{Integral calculus}
 775 @closecatbox
 776 @end defvr
 777
 778 @c -----------------------------------------------------------------------------
 779 @anchor{laplace}
 780 @deffn {Function} laplace (@var{expr}, @var{t}, @var{s})
 781
 782 Attempts to compute the Laplace transform of @var{expr} with respect to the
 783 variable @var{t} and transform parameter @var{s}.  The Laplace
 784 transform of the function @code{f(t)} is the one-sided transform defined by
 785 m4_displaymath(
 786 <<<F(s) = \int_0^{\infty} f(t) e^{-st} dt>>>,
 787 <<<
 788 @example
 789 F(s) = integrate(f(t) * exp(-s*t), t, 0, inf)
 790 @end example
 791 >>>)
 792
 793 where @math{F(s)} is the transform of @math{f(t)}, represented by @var{expr}.
 794
 795 @code{laplace} recognizes in @var{expr} the functions @mrefcomma{delta} @mrefcomma{exp}
 796 @mrefcomma{log} @mrefcomma{sin} @mrefcomma{cos} @mrefcomma{sinh} @mrefcomma{cosh} and @mrefcomma{erf}
 797 as well as @code{derivative}, @mrefcomma{integrate} @mrefcomma{sum} and @mrefdot{ilt} If
 798 @code{laplace} fails to find a transform the function @mref{specint} is called.
 799 @mref{specint} can find the laplace transform for expressions with special
 800 functions like the bessel functions @mrefcomma{bessel_j} @mrefcomma{bessel_i} @dots{}
 801 and can handle the @mref{unit_step} function.  See also @mrefdot{specint}
 802
 803 If @mref{specint} cannot find a solution too, a noun @code{laplace} is returned.
 804
 805 @c REPHRASE THIS
 806 @var{expr} may also be a linear, constant coefficient differential equation in
 807 which case @mref{atvalue} of the dependent variable is used.
 808 @c "used" -- USED HOW ??
 809 The required atvalue may be supplied either before or after the transform is
 810 computed.  Since the initial conditions must be specified at zero, if one has
 811 boundary conditions imposed elsewhere he can impose these on the general
 812 solution and eliminate the constants by solving the general solution
 813 for them and substituting their values back.
 814
 815 @code{laplace} recognizes convolution integrals of the form
 816 m4_displaymath(
 817 <<<\int_0^t f(x) g(t-x) dx>>>,
 818 <<<
 819 @example
 820 integrate (f(x) * g(t - x), x, 0, t)
 821 @end example
 822 >>>)
 823
 824 Other kinds of convolutions are not recognized.
 825
 826 Functional relations must be explicitly represented in @var{expr};
 827 implicit relations, established by @mrefcomma{depends} are not recognized.
 828 That is, if @math{f} depends on @math{x} and @math{y},
 829 @math{f (x, y)} must appear in @var{expr}.
 830
 831 See also @mrefcomma{ilt} the inverse Laplace transform.
 832
 833 Examples:
 834
 835 @c ===beg===
 836 @c laplace (exp (2*t + a) * sin(t) * t, t, s);
 837 @c laplace ('diff (f (x), x), x, s);
 838 @c diff (diff (delta (t), t), t);
 839 @c laplace (%, t, s);
 840 @c assume(a>0)$
 841 @c declare(a, integer)$
 842 @c laplace(gamma_incomplete(a,t),t,s),gamma_expand:true;
 843 @c factor(laplace(gamma_incomplete(1/2,t),t,s));
 844 @c assume(exp(%pi*s)>1)$
 845 @c laplace(sum((-1)^n*unit_step(t-n*%pi)*sin(t),n,0,inf),t,s),
 846 @c    simpsum;
 847 @c ===end===
 848 @example
 849 (%i1) laplace (exp (2*t + a) * sin(t) * t, t, s);
 850                             a
 851                           %e  (2 s - 4)
 852 (%o1)                    ---------------
 853                            2           2
 854                          (s  - 4 s + 5)
 855 (%i2) laplace ('diff (f (x), x), x, s);
 856 (%o2)             s laplace(f(x), x, s) - f(0)
 857 (%i3) diff (diff (delta (t), t), t);
 858                           2
 859                          d
 860 (%o3)                    --- (delta(t))
 861                            2
 862                          dt
 863 (%i4) laplace (%, t, s);
 864                             !
 865                d            !         2
 866 (%o4)        - -- (delta(t))!      + s  - delta(0) s
 867                dt           !
 868                             !t = 0
 869 (%i5) assume(a>0)$
 870 (%i6) laplace(gamma_incomplete(a,t),t,s),gamma_expand:true;
 871                                               - a - 1
 872                          gamma(a)   gamma(a) s
 873 (%o6)                    -------- - -----------------
 874                             s            1     a
 875                                         (- + 1)
 876                                          s
 877 (%i7) factor(laplace(gamma_incomplete(1/2,t),t,s));
 878                                               s + 1
 879                       sqrt(%pi) (sqrt(s) sqrt(-----) - 1)
 880                                                 s
 881 (%o7)                 -----------------------------------
 882                                 3/2      s + 1
 883                                s    sqrt(-----)
 884                                            s
 885 (%i8) assume(exp(%pi*s)>1)$
 886 (%i9) laplace(sum((-1)^n*unit_step(t-n*%pi)*sin(t),n,0,inf),t,s),
 887          simpsum;
 888 @group
 889                          %i                         %i
 890               ------------------------ - ------------------------
 891                               - %pi s                    - %pi s
 892               (s + %i) (1 - %e       )   (s - %i) (1 - %e       )
 893 (%o9)         ---------------------------------------------------
 894                                        2
 895 @end group
 896 (%i9) factor(%);
 897                                       %pi s
 898                                     %e
 899 (%o9)                   -------------------------------
 900                                              %pi s
 901                         (s - %i) (s + %i) (%e      - 1)
 902
 903 @end example
 904
 905 @opencatbox{Categories:}
 906 @category{Laplace transform}
 907 @category{Differential equations}
 908 @closecatbox
 909 @end deffn
 910
 911 @c NEEDS EXAMPLES
 912
 913 @c -----------------------------------------------------------------------------
 914 @anchor{ldefint}
 915 @deffn {Function} ldefint (@var{expr}, @var{x}, @var{a}, @var{b})
 916
 917 Attempts to compute the definite integral of @var{expr} by using @mref{limit}
 918 to evaluate the indefinite integral of @var{expr} with respect to @var{x}
 919 at the upper limit @var{b} and at the lower limit @var{a}.
 920 If it fails to compute the definite integral,
 921 @code{ldefint} returns an expression containing limits as noun forms.
 922
 923 @code{ldefint} is not called from @mrefcomma{integrate} so executing
 924 @code{ldefint (@var{expr}, @var{x}, @var{a}, @var{b})} may yield a different
 925 result than @code{integrate (@var{expr}, @var{x}, @var{a}, @var{b})}.
 926 @code{ldefint} always uses the same method to evaluate the definite integral,
 927 while @code{integrate} may employ various heuristics and may recognize some
 928 special cases.
 929
 930 @opencatbox{Categories:}
 931 @category{Integral calculus}
 932 @closecatbox
 933 @end deffn
 934
 935 @c UMM, IS THERE SOME TEXT MISSING HERE ???
 936 @c WHAT IS THIS ABOUT EXACTLY ??
 937
 938 @c -----------------------------------------------------------------------------
 939 @anchor{pwilt}
 940 @deffn {Function} pwilt (@var{expr}, @var{s}, @var{t})
 941
 942 Computes the inverse Laplace transform of @var{expr} with
 943 respect to @var{s} and parameter @var{t}. Unlike @mrefcomma{ilt}
 944 @code{pwilt} is able to return piece-wise and periodic functions
 945 and can also handle some cases with polynomials of degree greater than 3
 946 in the denominator.
 947
 948 Two examples where @code{ilt} fails:
 949 @c ===beg===
 950 @c pwilt (exp(-s)*s/(s^3-2*s-s+2), s, t);
 951 @c pwilt ((s^2+2)/(s^2-1), s, t);
 952 @c ===end===
 953 @example
 954 (%i1) pwilt (exp(-s)*s/(s^3-2*s-s+2), s, t);
 955                                        t - 1       - 2 (t - 1)
 956                              (t - 1) %e        2 %e
 957 (%o1)         hstep(t - 1) (--------------- - ---------------)
 958                                     3                 9
 959
 960 (%i2) pwilt ((s^2+2)/(s^2-1), s, t);
 961                                          t       - t
 962                                      3 %e    3 %e
 963 (%o2)                    delta(t) + ----- - -------
 964                                        2        2
 965 @end example
 966
 967 @opencatbox{Categories:}
 968 @category{Laplace transform}
 969 @closecatbox
 970 @end deffn
 971
 972 @c -----------------------------------------------------------------------------
 973 @anchor{potential}
 974 @deffn {Function} potential (@var{givengradient})
 975
 976 The calculation makes use of the global variable @code{potentialzeroloc[0]}
 977 which must be @code{nonlist} or of the form
 978
 979 @example
 980 [indeterminatej=expressionj, indeterminatek=expressionk, ...]
 981 @end example
 982
 983 the former being equivalent to the nonlist expression for all right-hand
 984 sides in the latter.  The indicated right-hand sides are used as the
 985 lower limit of integration.  The success of the integrations may
 986 depend upon their values and order.  @code{potentialzeroloc} is initially set
 987 to 0.
 988 @end deffn
 989
 990 @c -----------------------------------------------------------------------------
 991 @anchor{prefer_d}
 992 @defvr {Option variable} prefer_d
 993 Default value: @code{false}
 994
 995 When @code{prefer_d} is @code{true}, @mref{specint} will prefer to
 996 express solutions using @mref{parabolic_cylinder_d} rather than
 997 hypergeometric functions.
 998
 999 In the example below, the solution contains @mref{parabolic_cylinder_d}
1000 when @code{prefer_d} is @code{true}.
1001
1002 @c ===beg===
1003 @c assume(s>0);
1004 @c factor(ex:specint(%e^-(t^2/8)*exp(-s*t),t));
1005 @c specint(ex,t),prefer_d=true;
1006 @c ===end===
1007 @example
1008 @group
1009 (%i1) assume(s>0);
1010 (%o1)                               [s > 0]
1011 @end group
1012 @group
1013 (%i2) factor(specint(ex:%e^-(t^2/8)*exp(-s*t),t));
1014                                          2
1015                                       2 s
1016 (%o2)           - sqrt(2) sqrt(%pi) %e     (erf(sqrt(2) s) - 1)
1017 @end group
1018 @group
1019 (%i3) specint(ex,t),prefer_d=true;
1020                                                           2
1021                                                          s
1022                                                          --
1023                                                  s       8
1024                     parabolic_cylinder_d(- 1, -------) %e
1025                                               sqrt(2)
1026 (%o3)               ---------------------------------------
1027                                     sqrt(2)
1028
1029 @end group
1030 @end example
1031
1032 @opencatbox{Categories:}
1033 @category{Laplace transform}
1034 @category{Special functions}
1035 @closecatbox
1036 @end defvr
1037
1038 @c -----------------------------------------------------------------------------
1039 @anchor{residue}
1040 @deffn {Function} residue (@var{expr}, @var{z}, @var{z_0})
1041
1042 Computes the residue in the complex plane of the expression @var{expr} when the
1043 variable @var{z} assumes the value @var{z_0}.  The residue is the coefficient of
1044 @code{(@var{z} - @var{z_0})^(-1)} in the Laurent series for @var{expr}.
1045
1046 @c ===beg===
1047 @c residue (s/(s**2+a**2), s, a*%i);
1048 @c residue (sin(a*x)/x**4, x, 0);
1049 @c ===end===
1050 @example
1051 @group
1052 (%i1) residue (s/(s**2+a**2), s, a*%i);
1053                                 1
1054 (%o1)                           -
1055                                 2
1056 @end group
1057 @group
1058 (%i2) residue (sin(a*x)/x**4, x, 0);
1059                                  3
1060                                 a
1061 (%o2)                         - --
1062                                 6
1063 @end group
1064 @end example
1065
1066 @opencatbox{Categories:}
1067 @category{Integral calculus}
1068 @category{Complex variables}
1069 @closecatbox
1070 @end deffn
1071
1072 @c -----------------------------------------------------------------------------
1073 @anchor{risch}
1074 @deffn {Function} risch (@var{expr}, @var{x})
1075
1076 Integrates @var{expr} with respect to @var{x} using the
1077 transcendental case of the Risch algorithm.  (The algebraic case of
1078 the Risch algorithm has not been implemented.)  This currently
1079 handles the cases of nested exponentials and logarithms which the main
1080 part of @code{integrate} can't do.  @mref{integrate} will automatically apply
1081 @code{risch} if given these cases.
1082
1083 @code{erfflag}, if @code{false}, prevents @code{risch} from introducing the
1084 @code{erf} function in the answer if there were none in the integrand to begin
1085 with.
1086
1087 @c ===beg===
1088 @c risch (x^2*erf(x), x);
1089 @c diff(%, x), ratsimp;
1090 @c ===end===
1091 @example
1092 @group
1093 (%i1) risch (x^2*erf(x), x);
1094                                                         2
1095              3                      2                - x
1096         %pi x  erf(x) + (sqrt(%pi) x  + sqrt(%pi)) %e
1097 (%o1)   -------------------------------------------------
1098                               3 %pi
1099 @end group
1100 @group
1101 (%i2) diff(%, x), ratsimp;
1102                              2
1103 (%o2)                       x  erf(x)
1104 @end group
1105 @end example
1106
1107 @opencatbox{Categories:}
1108 @category{Integral calculus}
1109 @closecatbox
1110 @end deffn
1111
1112 @anchor{specint}
1113 @deffn {Function} specint (exp(- s*@var{t}) * @var{expr}, @var{t})
1114
1115 Compute the Laplace transform of @var{expr} with respect to the variable @var{t}.
1116 The integrand @var{expr} may contain special functions.   The
1117 parameter @var{s} maybe be named something else; it is determined
1118 automatically, as the examples below show where @var{p} is used in
1119 some places.
1120
1121 The following special functions are handled by @code{specint}: incomplete gamma
1122 function, error functions (but not the error function @code{erfi}, it is easy to
1123 transform @code{erfi} e.g. to the error function @code{erf}), exponential
1124 integrals, bessel functions (including products of bessel functions), hankel
1125 functions, hermite and the laguerre polynomials.
1126
1127 Furthermore, @code{specint} can handle the hypergeometric function
1128 @code{%f[p,q]([],[],z)}, the Whittaker function of the first kind
1129 @code{%m[u,k](z)} and of the second kind @code{%w[u,k](z)}.
1130
1131 The result may be in terms of special functions and can include unsimplified
1132 hypergeometric functions.  If variable @mref{prefer_d} is @code{true}
1133 then the @mref{parabolic_cylinder_d} function may be used in the result
1134 in preference to hypergeometric functions.
1135
1136 When @mref{laplace} fails to find a Laplace transform, @code{specint} is called.
1137 Because @mref{laplace} knows more general rules for Laplace transforms, it is
1138 preferable to use @mref{laplace} and not @code{specint}.
1139
1140 @code{demo("hypgeo")} displays several examples of Laplace transforms computed by
1141 @code{specint}.
1142
1143 Examples:
1144 @c ===beg===
1145 @c assume (p > 0, a > 0)$
1146 @c specint (t^(1/2) * exp(-a*t/4) * exp(-p*t), t);
1147 @c specint (t^(1/2) * bessel_j(1, 2 * a^(1/2) * t^(1/2))
1148 @c               * exp(-p*t), t);
1149 @c ===end===
1150 @example
1151 (%i1) assume (p > 0, a > 0)$
1152 @group
1153 (%i2) specint (t^(1/2) * exp(-a*t/4) * exp(-p*t), t);
1154                            sqrt(%pi)
1155 (%o2)                     ------------
1156                                  a 3/2
1157                           2 (p + -)
1158                                  4
1159 @end group
1160 @group
1161 (%i3) specint (t^(1/2) * bessel_j(1, 2 * a^(1/2) * t^(1/2))
1162               * exp(-p*t), t);
1163                                    - a/p
1164                          sqrt(a) %e
1165 (%o3)                    ---------------
1166                                 2
1167                                p
1168 @end group
1169 @end example
1170
1171 Examples for exponential integrals:
1172
1173 @example
1174 (%i4) assume(s>0,a>0,s-a>0)$
1175 (%i5) ratsimp(specint(%e^(a*t)
1176                       *(log(a)+expintegral_e1(a*t))*%e^(-s*t),t));
1177                              log(s)
1178 (%o5)                        ------
1179                              s - a
1180 (%i6) logarc:true$
1181
1182 (%i7) gamma_expand:true$
1183
1184 radcan(specint((cos(t)*expintegral_si(t)
1185                      -sin(t)*expintegral_ci(t))*%e^(-s*t),t));
1186                              log(s)
1187 (%o8)                        ------
1188                               2
1189                              s  + 1
1190 ratsimp(specint((2*t*log(a)+2/a*sin(a*t)
1191                       -2*t*expintegral_ci(a*t))*%e^(-s*t),t));
1192                                2    2
1193                           log(s  + a )
1194 (%o9)                     ------------
1195                                 2
1196                                s
1197 @end example
1198
1199 Results when using the expansion of @mref{gamma_incomplete} and when changing
1200 the representation to @mref{expintegral_e1}:
1201
1202 @example
1203 (%i10) assume(s>0)$
1204 (%i11) specint(1/sqrt(%pi*t)*unit_step(t-k)*%e^(-s*t),t);
1205                                             1
1206                             gamma_incomplete(-, k s)
1207                                             2
1208 (%o11)                      ------------------------
1209                                sqrt(%pi) sqrt(s)
1210
1211 (%i12) gamma_expand:true$
1212 (%i13) specint(1/sqrt(%pi*t)*unit_step(t-k)*%e^(-s*t),t);
1213                               erfc(sqrt(k) sqrt(s))
1214 (%o13)                        ---------------------
1215                                      sqrt(s)
1216
1217 (%i14) expintrep:expintegral_e1$
1218 (%i15) ratsimp(specint(1/(t+a)^2*%e^(-s*t),t));
1219                               a s
1220                         a s %e    expintegral_e1(a s) - 1
1221 (%o15)                - ---------------------------------
1222                                         a
1223 @end example
1224
1225 @opencatbox{Categories:}
1226 @category{Laplace transform}
1227 @closecatbox
1228 @end deffn
1229
1230 @c NEEDS EXPANSION, CLARIFICATION, AND EXAMPLES
1231
1232 @c -----------------------------------------------------------------------------
1233 @anchor{tldefint}
1234 @deffn {Function} tldefint (@var{expr}, @var{x}, @var{a}, @var{b})
1235
1236 Equivalent to @code{ldefint} with @code{tlimswitch} set to @code{true}.
1237
1238 @opencatbox{Categories:}
1239 @category{Integral calculus}
1240 @closecatbox
1241 @end deffn
1242
1243 @footnotestyle end
1244
1245 @c -----------------------------------------------------------------------------
1246 @node Introduction to QUADPACK, Functions and Variables for QUADPACK, Functions and Variables for Integration, Integration
1247 @section Introduction to QUADPACK
1248 @c -----------------------------------------------------------------------------
1249
1250 @c FOLLOWING TEXT ADAPTED WITH HEAVY MODIFICATION FROM https://www.netlib.org/slatec/src/qpdoc.f
1251
1252 QUADPACK is a collection of functions for the numerical
1253 computation of one-dimensional definite integrals.
1254 It originated from a joint project of
1255 R. Piessens @footnote{Applied Mathematics and Programming Division, K.U. Leuven},
1256 E. de Doncker @footnote{Applied Mathematics and Programming Division, K.U. Leuven},
1257 C. Ueberhuber @footnote{Institut f@"ur Mathematik, T.U. Wien},
1258 and D. Kahaner @footnote{National Bureau of Standards, Washington, D.C., U.S.A}.
1259
1260 The QUADPACK library included in Maxima is an automatic translation (via the
1261 program @code{f2cl}) of the Fortran source code of QUADPACK as it appears in
1262 the SLATEC Common Mathematical Library, Version 4.1 @footnote{@url{https://www.netlib.org/slatec}}.
1263 The SLATEC library is dated July 1993, but the QUADPACK functions
1264 were written some years before.
1265 There is another version of QUADPACK at Netlib @footnote{@url{https://www.netlib.org/quadpack}};
1266 it is not clear how that version differs from the SLATEC version.
1267
1268 The QUADPACK functions included in Maxima are all automatic, in the sense that
1269 these functions attempt to compute a result to a specified accuracy, requiring
1270 an unspecified number of function evaluations.  Maxima's Lisp translation of
1271 QUADPACK also includes some non-automatic functions, but they are not exposed
1272 at the Maxima level.
1273
1274 Further information about QUADPACK can be found in the QUADPACK book
1275 @footnote{R. Piessens, E. de Doncker-Kapenga, C.W. Uberhuber, and D.K. Kahaner.
1276 @i{QUADPACK: A Subroutine Package for Automatic Integration.}
1277 Berlin: Springer-Verlag, 1983, ISBN 0387125531.}.
1278
1279 @c -----------------------------------------------------------------------------
1280 @subsection Overview
1281 @c -----------------------------------------------------------------------------
1282
1283 @table @code
1284 @item @mref{quad_qag}
1285 Integration of a general function over a finite interval.
1286 @mref{quad_qag} implements a simple globally adaptive integrator using the
1287 strategy of Aind (Piessens, 1973).
1288 The caller may choose among 6 pairs of Gauss-Kronrod quadrature
1289 formulae for the rule evaluation component.
1290 The high-degree rules are suitable for strongly oscillating integrands.
1291
1292 @item @mref{quad_qags}
1293 Integration of a general function over a finite interval.
1294 @mref{quad_qags} implements globally adaptive interval subdivision with
1295 extrapolation (de Doncker, 1978) by the Epsilon algorithm (Wynn, 1956).
1296
1297 @item @mref{quad_qagi}
1298 Integration of a general function over an infinite or semi-infinite interval.
1299 The interval is mapped onto a finite interval and
1300 then the same strategy as in @code{quad_qags} is applied.
1301
1302 @item @mref{quad_qawo}
1303
1304 Integration of
1305 m4_math(<<<\cos(\omega x) f(x)>>>,<<<@math{cos(omega x) f(x)}>>>)
1306 or
1307 m4_math(<<<\sin(\omega x) f(x)>>>,<<<@math{sin(omega x) f(x)}>>>)
1308 over a
1309 finite interval, where
1310 m4_math(\omega, @math{omega})
1311 is a constant.
1312 The rule evaluation component is based on the modified Clenshaw-Curtis
1313 technique.  @mref{quad_qawo} applies adaptive subdivision with extrapolation,
1314 similar to @mrefdot{quad_qags}
1315
1316 @item @mref{quad_qawf}
1317 Calculates a Fourier cosine or Fourier sine transform on a semi-infinite
1318 interval.  The same approach as in @mref{quad_qawo} is applied on successive
1319 finite intervals, and convergence acceleration by means of the Epsilon algorithm
1320 (Wynn, 1956) is applied to the series of the integral contributions.
1321
1322 @item @mref{quad_qaws}
1323 Integration of
1324 m4_math(w(x)f(x), @math{w(x) f(x)})
1325 over a finite interval @math{[a, b]}, where
1326 @math{w} is a function of the form
1327 m4_math((x-a)^\alpha (b-x)^\beta v(x), @math{(x - a)^alpha (b -
1328 x)^beta v(x)})
1329 and
1330 @math{v(x)} is 1 or
1331 m4_math(\log(x-a), @math{log(x - a)})
1332 or
1333 m4_math(\log(b-x), @math{log(b - x)})
1334 or
1335 m4_mathcomma(\log(x-a)\log(b-x), @math{log(x - a) log(b - x)})
1336 and
1337 m4_math(\alpha > -1, @math{alpha > -1})
1338 and
1339 m4_mathdot(\beta > -1, @math{beta > -1})
1340
1341 A globally adaptive subdivision strategy is applied, with modified
1342 Clenshaw-Curtis integration on the subintervals which contain @math{a}
1343 or @math{b}.
1344
1345 @item @mref{quad_qawc}
1346 Computes the Cauchy principal value of @math{f(x)/(x - c)} over a finite
1347 interval @math{(a, b)} and specified @math{c}.
1348 The strategy is globally adaptive, and modified
1349 Clenshaw-Curtis integration is used on the subranges
1350 which contain the point @math{x = c}.
1351
1352 @item @mref{quad_qagp}
1353 Basically the same as @mref{quad_qags} but points of singularity or
1354 discontinuity of the integrand must be supplied.  This makes it easier
1355 for the integrator to produce a good solution.
1356 @end table
1357
1358
1359 @opencatbox{Categories:}
1360 @category{Integral calculus}
1361 @category{Numerical methods}
1362 @category{Share packages}
1363 @category{Package quadpack}
1364 @closecatbox
1365
1366 @c -----------------------------------------------------------------------------
1367 @node Functions and Variables for QUADPACK, , Introduction to QUADPACK, Integration
1368 @section Functions and Variables for QUADPACK
1369 @c -----------------------------------------------------------------------------
1370
1371 @c THERE ARE OPTIONAL ARGUMENTS WHICH MAKES LISTING THE VARIANTS A LITTLE TEDIOUS
1372 @c NEED A MORE CONVENIENT (AND NONAMBIGUOUS) NOTATION FOR OPTIONAL ARGUMENTS
1373
1374 @c -----------------------------------------------------------------------------
1375 @anchor{quad_qag}
1376 @deffn  {Function} quad_qag @
1377 @fname{quad_qag} (@var{f(x)}, @var{x}, @var{a}, @var{b}, @var{key}, [@var{epsrel}, @var{epsabs}, @var{limit}]) @
1378 @fname{quad_qag} (@var{f}, @var{x}, @var{a}, @var{b}, @var{key}, [@var{epsrel}, @var{epsabs}, @var{limit}])
1379
1380 Integration of a general function over a finite interval.  @code{quad_qag}
1381 implements a simple globally adaptive integrator using the strategy of Aind
1382 (Piessens, 1973).  The caller may choose among 6 pairs of Gauss-Kronrod
1383 quadrature formulae for the rule evaluation component.  The high-degree rules
1384 are suitable for strongly oscillating integrands.
1385
1386 @code{quad_qag} computes the integral
1387
1388 m4_displaymath(
1389 <<<\int_a^b f(x)\, dx>>>,
1390 @math{integrate (f(x), x, a, b)}>>>)
1391
1392 The function to be integrated is @math{f(x)}, with dependent
1393 variable @math{x}, and the function is to be integrated between the
1394 limits @math{a} and @math{b}.  @var{key} is the integrator to be used
1395 and should be an integer between 1 and 6, inclusive.  The value of
1396 @var{key} selects the order of the Gauss-Kronrod integration rule.
1397 High-order rules are suitable for strongly oscillating integrands.
1398
1399 The integrand may be specified as the name of a Maxima or Lisp function or
1400 operator, a Maxima lambda expression, or a general Maxima expression.
1401
1402 The numerical integration is done adaptively by subdividing the
1403 integration region into sub-intervals until the desired accuracy is
1404 achieved.
1405
1406 The keyword arguments are optional and may be specified in any order.
1407 They all take the form @code{key=val}.  The keyword arguments are:
1408
1409 @table @code
1410 @item epsrel
1411 Desired relative error of approximation.  Default is 1d-8.
1412 @item epsabs
1413 Desired absolute error of approximation.  Default is 0.
1414 @item limit
1415 Size of internal work array.  @var{limit} is the
1416 maximum number of subintervals to use.  Default is 200.
1417 @end table
1418
1419 @code{quad_qag} returns a list of four elements:
1420
1421 @itemize
1422 @item
1423 an approximation to the integral,
1424 @item
1425 the estimated absolute error of the approximation,
1426 @item
1427 the number integrand evaluations,
1428 @item
1429 an error code.
1430 @end itemize
1431
1432 The error code (fourth element of the return value) can have the values:
1433
1434 @table @code
1435 @item 0
1436 if no problems were encountered;
1437 @item 1
1438 if too many sub-intervals were done;
1439 @item 2
1440 if excessive roundoff error is detected;
1441 @item 3
1442 if extremely bad integrand behavior occurs;
1443 @item 6
1444 if the input is invalid.
1445
1446 @end table
1447
1448 @c NEED CROSS REFS HERE -- EITHER CROSS REF A QUADPACK OVERVIEW, OR CROSS REF EACH OF THE quad_* FUNCTIONS
1449
1450 Examples:
1451
1452 @c ===beg===
1453 @c quad_qag (x^(1/2)*log(1/x), x, 0, 1, 3, 'epsrel=5d-8);
1454 @c integrate (x^(1/2)*log(1/x), x, 0, 1);
1455 @c ===end===
1456 @example
1457 @group
1458 (%i1) quad_qag (x^(1/2)*log(1/x), x, 0, 1, 3, 'epsrel=5d-8);
1459 (%o1)    [.4444444444492108, 3.1700968502883E-9, 961, 0]
1460 @end group
1461 @group
1462 (%i2) integrate (x^(1/2)*log(1/x), x, 0, 1);
1463                                 4
1464 (%o2)                           -
1465                                 9
1466 @end group
1467 @end example
1468
1469 @opencatbox{Categories:}
1470 @category{Numerical methods}
1471 @category{Package quadpack}
1472 @closecatbox
1473 @end deffn
1474
1475 @c THERE ARE OPTIONAL ARGUMENTS WHICH MAKES LISTING THE VARIANTS A LITTLE TEDIOUS
1476 @c NEED A MORE CONVENIENT (AND NONAMBIGUOUS) NOTATION FOR OPTIONAL ARGUMENTS
1477
1478 @c -----------------------------------------------------------------------------
1479 @anchor{quad_qags}
1480 @deffn  {Function} quad_qags @
1481 @fname{quad_qags} (@var{f(x)}, @var{x}, @var{a}, @var{b}, [@var{epsrel}, @var{epsabs}, @var{limit}]) @
1482 @fname{quad_qags} (@var{f}, @var{x}, @var{a}, @var{b}, [@var{epsrel}, @var{epsabs}, @var{limit}])
1483
1484 Integration of a general function over a finite interval.
1485 @code{quad_qags} implements globally adaptive interval subdivision with
1486 extrapolation (de Doncker, 1978) by the Epsilon algorithm (Wynn, 1956).
1487
1488 @code{quad_qags} computes the integral
1489
1490 m4_displaymath(
1491 <<<\int_a^b f(x)\, dx>>>,
1492 @math{integrate (f(x), x, a, b)}>>>)
1493
1494 The function to be integrated is @math{f(x)}, with
1495 dependent variable @math{x}, and the function is to be integrated
1496 between the limits @math{a} and @math{b}.
1497
1498 The integrand may be specified as the name of a Maxima or Lisp function or
1499 operator, a Maxima lambda expression, or a general Maxima expression.
1500
1501 The keyword arguments are optional and may be specified in any order.
1502 They all take the form @code{key=val}.  The keyword arguments are:
1503
1504 @table @code
1505 @item epsrel
1506 Desired relative error of approximation.  Default is 1d-8.
1507 @item epsabs
1508 Desired absolute error of approximation.  Default is 0.
1509 @item limit
1510 Size of internal work array.  @var{limit} is the
1511 maximum number of subintervals to use.  Default is 200.
1512 @end table
1513
1514 @code{quad_qags} returns a list of four elements:
1515
1516 @itemize
1517 @item
1518 an approximation to the integral,
1519 @item
1520 the estimated absolute error of the approximation,
1521 @item
1522 the number integrand evaluations,
1523 @item
1524 an error code.
1525 @end itemize
1526
1527 The error code (fourth element of the return value) can have the values:
1528
1529 @table @code
1530 @item 0
1531 no problems were encountered;
1532 @item 1
1533 too many sub-intervals were done;
1534 @item 2
1535 excessive roundoff error is detected;
1536 @item 3
1537 extremely bad integrand behavior occurs;
1538 @item 4
1539 failed to converge
1540 @item 5
1541 integral is probably divergent or slowly convergent
1542 @item 6
1543 if the input is invalid.
1544 @end table
1545
1546 @c NEED CROSS REFS HERE -- EITHER CROSS REF A QUADPACK OVERVIEW, OR CROSS REF EACH OF THE quad_* FUNCTIONS
1547
1548 Examples:
1549
1550 @c ===beg===
1551 @c quad_qags (x^(1/2)*log(1/x), x, 0, 1, 'epsrel=1d-10);
1552 @c ===end===
1553 @example
1554 @group
1555 (%i1) quad_qags (x^(1/2)*log(1/x), x, 0, 1, 'epsrel=1d-10);
1556 (%o1)   [.4444444444444448, 1.11022302462516E-15, 315, 0]
1557 @end group
1558 @end example
1559
1560 Note that @code{quad_qags} is more accurate and efficient than @code{quad_qag} for this integrand.
1561
1562 @opencatbox{Categories:}
1563 @category{Numerical methods}
1564 @category{Package quadpack}
1565 @closecatbox
1566 @end deffn
1567
1568 @c THERE ARE OPTIONAL ARGUMENTS WHICH MAKES LISTING THE VARIANTS A LITTLE TEDIOUS
1569 @c NEED A MORE CONVENIENT (AND NONAMBIGUOUS) NOTATION FOR OPTIONAL ARGUMENTS
1570
1571 @c -----------------------------------------------------------------------------
1572 @anchor{quad_qagi}
1573 @deffn  {Function} quad_qagi @
1574 @fname{quad_qagi} (@var{f(x)}, @var{x}, @var{a}, @var{b}, [@var{epsrel}, @var{epsabs}, @var{limit}]) @
1575 @fname{quad_qagi} (@var{f}, @var{x}, @var{a}, @var{b}, [@var{epsrel}, @var{epsabs}, @var{limit}])
1576
1577 Integration of a general function over an infinite or semi-infinite interval.
1578 The interval is mapped onto a finite interval and
1579 then the same strategy as in @code{quad_qags} is applied.
1580
1581 @code{quad_qagi} evaluates one of the following integrals
1582
1583 m4_displaymath(
1584 <<<\int_a^\infty f(x) \, dx>>>,
1585 <<<@math{integrate (f(x), x, a, inf)}>>>)
1586
1587 m4_displaymath(
1588 <<<\int_\infty^a f(x) \, dx>>>,
1589 <<<@math{integrate (f(x), x, minf, a)}>>>)
1590
1591 m4_displaymath(
1592 <<<\int_{-\infty}^\infty f(x) \, dx>>>,
1593 <<<@math{integrate (f(x), x, minf, inf)}>>>)
1594
1595 using the Quadpack QAGI routine.  The function to be integrated is
1596 @math{f(x)}, with dependent variable @math{x}, and the function is to
1597 be integrated over an infinite range.
1598
1599 The integrand may be specified as the name of a Maxima or Lisp function or
1600 operator, a Maxima lambda expression, or a general Maxima expression.
1601
1602 One of the limits of integration must be infinity.  If not, then
1603 @code{quad_qagi} will just return the noun form.
1604
1605 The keyword arguments are optional and may be specified in any order.
1606 They all take the form @code{key=val}.  The keyword arguments are:
1607
1608 @table @code
1609 @item epsrel
1610 Desired relative error of approximation.  Default is 1d-8.
1611 @item epsabs
1612 Desired absolute error of approximation.  Default is 0.
1613 @item limit
1614 Size of internal work array.  @var{limit} is the
1615 maximum number of subintervals to use.  Default is 200.
1616 @end table
1617
1618 @code{quad_qagi} returns a list of four elements:
1619
1620 @itemize
1621 @item
1622 an approximation to the integral,
1623 @item
1624 the estimated absolute error of the approximation,
1625 @item
1626 the number integrand evaluations,
1627 @item
1628 an error code.
1629 @end itemize
1630
1631 The error code (fourth element of the return value) can have the values:
1632
1633 @table @code
1634 @item 0
1635 no problems were encountered;
1636 @item 1
1637 too many sub-intervals were done;
1638 @item 2
1639 excessive roundoff error is detected;
1640 @item 3
1641 extremely bad integrand behavior occurs;
1642 @item 4
1643 failed to converge
1644 @item 5
1645 integral is probably divergent or slowly convergent
1646 @item 6
1647 if the input is invalid.
1648
1649 @end table
1650
1651 @c NEED CROSS REFS HERE -- EITHER CROSS REF A QUADPACK OVERVIEW, OR CROSS REF EACH OF THE quad_* FUNCTIONS
1652
1653 Examples:
1654
1655 @c ===beg===
1656 @c quad_qagi (x^2*exp(-4*x), x, 0, inf, 'epsrel=1d-8);
1657 @c integrate (x^2*exp(-4*x), x, 0, inf);
1658 @c ===end===
1659 @example
1660 @group
1661 (%i1) quad_qagi (x^2*exp(-4*x), x, 0, inf, 'epsrel=1d-8);
1662 (%o1)        [0.03125, 2.95916102995002E-11, 105, 0]
1663 @end group
1664 @group
1665 (%i2) integrate (x^2*exp(-4*x), x, 0, inf);
1666                                1
1667 (%o2)                          --
1668                                32
1669 @end group
1670 @end example
1671
1672 @opencatbox{Categories:}
1673 @category{Numerical methods}
1674 @category{Package quadpack}
1675 @closecatbox
1676 @end deffn
1677
1678 @c THERE ARE OPTIONAL ARGUMENTS WHICH MAKES LISTING THE VARIANTS A LITTLE TEDIOUS
1679 @c NEED A MORE CONVENIENT (AND NONAMBIGUOUS) NOTATION FOR OPTIONAL ARGUMENTS
1680
1681 @c -----------------------------------------------------------------------------
1682 @anchor{quad_qawc}
1683 @deffn  {Function} quad_qawc @
1684 @fname{quad_qawc} (@var{f(x)}, @var{x}, @var{c}, @var{a}, @var{b}, [@var{epsrel}, @var{epsabs}, @var{limit}]) @
1685 @fname{quad_qawc} (@var{f}, @var{x}, @var{c}, @var{a}, @var{b}, [@var{epsrel}, @var{epsabs}, @var{limit}])
1686
1687 Computes the Cauchy principal value of @math{f(x)/(x - c)} over a finite
1688 interval.  The strategy is globally adaptive, and modified
1689 Clenshaw-Curtis integration is used on the subranges
1690 which contain the point @math{x = c}.
1691
1692 @code{quad_qawc} computes the Cauchy principal value of
1693
1694 m4_displaymath(
1695 <<<\int_{a}^{b}{{{f\left(x\right)}\over{x-c}}\>dx}>>>,
1696 <<<@math{integrate (f(x)/(x - c), x, a, b)}>>>)
1697
1698 using the Quadpack QAWC routine.  The function to be integrated is
1699 @math{f(x)/(x-c)}, with dependent variable @math{x}, and the
1700 function is to be integrated over the interval @math{a} to @math{b}.
1701
1702 The integrand may be specified as the name of a Maxima or Lisp function or
1703 operator, a Maxima lambda expression, or a general Maxima expression.
1704
1705 The keyword arguments are optional and may be specified in any order.
1706 They all take the form @code{key=val}.  The keyword arguments are:
1707
1708 @table @code
1709 @item epsrel
1710 Desired relative error of approximation.  Default is 1d-8.
1711 @item epsabs
1712 Desired absolute error of approximation.  Default is 0.
1713 @item limit
1714 Size of internal work array.  @var{limit} is the
1715 maximum number of subintervals to use.  Default is 200.
1716 @end table
1717
1718 @code{quad_qawc} returns a list of four elements:
1719
1720 @itemize
1721 @item
1722 an approximation to the integral,
1723 @item
1724 the estimated absolute error of the approximation,
1725 @item
1726 the number integrand evaluations,
1727 @item
1728 an error code.
1729 @end itemize
1730
1731 The error code (fourth element of the return value) can have the values:
1732
1733 @table @code
1734 @item 0
1735 no problems were encountered;
1736 @item 1
1737 too many sub-intervals were done;
1738 @item 2
1739 excessive roundoff error is detected;
1740 @item 3
1741 extremely bad integrand behavior occurs;
1742 @item 6
1743 if the input is invalid.
1744
1745 @end table
1746
1747 Examples:
1748
1749 @c ===beg===
1750 @c quad_qawc (2^(-5)*((x-1)^2+4^(-5))^(-1), x, 2, 0, 5, 'epsrel=1d-7);
1751 @c integrate (2^(-alpha)*(((x-1)^2 + 4^(-alpha))*(x-2))^(-1), x, 0, 5);
1752 @c ev (%, alpha=5, numer);
1753 @c ===end===
1754 @example
1755 @group
1756 (%i1) quad_qawc (2^(-5)*((x-1)^2+4^(-5))^(-1), x, 2, 0, 5,
1757                  'epsrel=1d-7);
1758 (%o1)    [- 3.130120337415925, 1.306830140249558E-8, 495, 0]
1759 @end group
1760 @group
1761 (%i2) integrate (2^(-alpha)*(((x-1)^2 + 4^(-alpha))*(x-2))^(-1),
1762       x, 0, 5);
1763 Principal Value
1764                        alpha
1765         alpha       9 4                 9
1766        4      log(------------- + -------------)
1767                       alpha           alpha
1768                   64 4      + 4   64 4      + 4
1769 (%o2) (-----------------------------------------
1770                         alpha
1771                      2 4      + 2
1772
1773        3 alpha                       3 alpha
1774        -------                       -------
1775           2            alpha/2          2          alpha/2
1776     2 4        atan(4 4       )   2 4        atan(4       )   alpha
1777   - --------------------------- - -------------------------)/2
1778               alpha                        alpha
1779            2 4      + 2                 2 4      + 2
1780 @end group
1781 @group
1782 (%i3) ev (%, alpha=5, numer);
1783 (%o3)                    - 3.130120337415917
1784 @end group
1785 @end example
1786
1787 @opencatbox{Categories:}
1788 @category{Numerical methods}
1789 @category{Package quadpack}
1790 @closecatbox
1791 @end deffn
1792
1793 @c THERE ARE OPTIONAL ARGUMENTS WHICH MAKES LISTING THE VARIANTS A LITTLE TEDIOUS
1794 @c NEED A MORE CONVENIENT (AND NONAMBIGUOUS) NOTATION FOR OPTIONAL ARGUMENTS
1795
1796 @c -----------------------------------------------------------------------------
1797 @anchor{quad_qawf}
1798 @deffn  {Function} quad_qawf @
1799 @fname{quad_qawf} (@var{f(x)}, @var{x}, @var{a}, @var{omega}, @var{trig}, [@var{epsabs}, @var{limit}, @var{maxp1}, @var{limlst}]) @
1800 @fname{quad_qawf} (@var{f}, @var{x}, @var{a}, @var{omega}, @var{trig}, [@var{epsabs}, @var{limit}, @var{maxp1}, @var{limlst}])
1801
1802 Calculates a Fourier cosine or Fourier sine transform on a semi-infinite
1803 interval using the Quadpack QAWF function.  The same approach as in
1804 @code{quad_qawo} is applied on successive finite intervals, and convergence
1805 acceleration by means of the Epsilon algorithm (Wynn, 1956) is applied to the
1806 series of the integral contributions.
1807
1808 @code{quad_qawf} computes the integral
1809
1810 m4_displaymath(
1811 <<<\int_a^\infty f(x) \, w(x) \, dx>>>,
1812 <<<@math{integrate (f(x)*w(x), x, a, inf)}>>>)
1813
1814 The weight function @math{w} is selected by @var{trig}:
1815
1816 @table @code
1817 @item cos
1818 m4_math(w(x) = \cos\omega x, @math{w(x) = cos (omega x)})
1819 @item sin
1820 m4_math(w(x) = \sin\omega x, @math{w(x) = sin (omega x)})
1821 @end table
1822
1823 The integrand may be specified as the name of a Maxima or Lisp function or
1824 operator, a Maxima lambda expression, or a general Maxima expression.
1825
1826 The keyword arguments are optional and may be specified in any order.
1827 They all take the form @code{key=val}.  The keyword arguments are:
1828
1829 @table @code
1830 @item epsabs
1831 Desired absolute error of approximation.  Default is 1d-10.
1832 @item limit
1833 Size of internal work array.  (@var{limit} - @var{limlst})/2 is the
1834 maximum number of subintervals to use.  Default is 200.
1835 @item maxp1
1836 Maximum number of Chebyshev moments.  Must be greater than 0.  Default
1837 is 100.
1838 @item limlst
1839 Upper bound on the number of cycles.  Must be greater than or equal to
1840 3.  Default is 10.
1841 @end table
1842
1843 @code{quad_qawf} returns a list of four elements:
1844
1845 @itemize
1846 @item
1847 an approximation to the integral,
1848 @item
1849 the estimated absolute error of the approximation,
1850 @item
1851 the number integrand evaluations,
1852 @item
1853 an error code.
1854 @end itemize
1855
1856 The error code (fourth element of the return value) can have the values:
1857
1858 @table @code
1859 @item 0
1860 no problems were encountered;
1861 @item 1
1862 too many sub-intervals were done;
1863 @item 2
1864 excessive roundoff error is detected;
1865 @item 3
1866 extremely bad integrand behavior occurs;
1867 @item 6
1868 if the input is invalid.
1869
1870 @end table
1871
1872 Examples:
1873
1874 @c ===beg===
1875 @c quad_qawf (exp(-x^2), x, 0, 1, 'cos, 'epsabs=1d-9);
1876 @c integrate (exp(-x^2)*cos(x), x, 0, inf);
1877 @c ev (%, numer);
1878 @c ===end===
1879 @example
1880 @group
1881 (%i1) quad_qawf (exp(-x^2), x, 0, 1, 'cos, 'epsabs=1d-9);
1882 (%o1)   [.6901942235215714, 2.84846300257552E-11, 215, 0]
1883 @end group
1884 @group
1885 (%i2) integrate (exp(-x^2)*cos(x), x, 0, inf);
1886                           - 1/4
1887                         %e      sqrt(%pi)
1888 (%o2)                   -----------------
1889                                 2
1890 @end group
1891 @group
1892 (%i3) ev (%, numer);
1893 (%o3)                   .6901942235215714
1894 @end group
1895 @end example
1896
1897 @opencatbox{Categories:}
1898 @category{Numerical methods}
1899 @category{Package quadpack}
1900 @closecatbox
1901 @end deffn
1902
1903 @c THERE ARE OPTIONAL ARGUMENTS WHICH MAKES LISTING THE VARIANTS A LITTLE TEDIOUS
1904 @c NEED A MORE CONVENIENT (AND NONAMBIGUOUS) NOTATION FOR OPTIONAL ARGUMENTS
1905
1906 @c -----------------------------------------------------------------------------
1907 @anchor{quad_qawo}
1908 @deffn  {Function} quad_qawo @
1909 @fname{quad_qawo} (@var{f(x)}, @var{x}, @var{a}, @var{b}, @var{omega}, @var{trig}, [@var{epsrel}, @var{epsabs}, @var{limit}, @var{maxp1}, @var{limlst}]) @
1910 @fname{quad_qawo} (@var{f}, @var{x}, @var{a}, @var{b}, @var{omega}, @var{trig}, [@var{epsrel}, @var{epsabs}, @var{limit}, @var{maxp1}, @var{limlst}])
1911
1912 Integration of
1913 m4_math(<<<\cos(\omega x) f(x)>>>,<<<@math{cos (omega x) f(x)}>>>)
1914 or
1915 m4_math(<<<\sin(\omega x)>>>, <<<@math{sin (omega x) f(x)}>>>)
1916 over a finite interval,
1917 where
1918 m4_math(\omega, @math{omega})
1919 is a constant.
1920 The rule evaluation component is based on the modified
1921 Clenshaw-Curtis technique.  @code{quad_qawo} applies adaptive subdivision with
1922 extrapolation, similar to @code{quad_qags}.
1923
1924 @code{quad_qawo} computes the integral using the Quadpack QAWO
1925 routine:
1926
1927 m4_displaymath(
1928 <<<\int_a^b f(x) \, w(x) \, dx>>>,
1929 <<<@math{integrate (f(x)*w(x), x, a, b)}>>>)
1930
1931 The weight function @math{w} is selected by @var{trig}:
1932
1933 @table @code
1934 @item cos
1935 m4_math(w(x) = \cos\omega x, @math{w(x) = cos (omega x)})
1936 @item sin
1937 m4_math(w(x) = \sin\omega x, @math{w(x) = sin (omega x)})
1938 @end table
1939
1940 The integrand may be specified as the name of a Maxima or Lisp function or
1941 operator, a Maxima lambda expression, or a general Maxima expression.
1942
1943 The keyword arguments are optional and may be specified in any order.
1944 They all take the form @code{key=val}.  The keyword arguments are:
1945
1946 @table @code
1947 @item epsrel
1948 Desired relative error of approximation.  Default is 1d-8.
1949 @item epsabs
1950 Desired absolute error of approximation.  Default is 0.
1951 @item limit
1952 Size of internal work array.  @var{limit}/2 is the
1953 maximum number of subintervals to use.  Default is 200.
1954 @item maxp1
1955 Maximum number of Chebyshev moments.  Must be greater than 0.  Default
1956 is 100.
1957 @item limlst
1958 Upper bound on the number of cycles.  Must be greater than or equal to
1959 3.  Default is 10.
1960 @end table
1961
1962 @code{quad_qawo} returns a list of four elements:
1963
1964 @itemize
1965 @item
1966 an approximation to the integral,
1967 @item
1968 the estimated absolute error of the approximation,
1969 @item
1970 the number integrand evaluations,
1971 @item
1972 an error code.
1973 @end itemize
1974
1975 The error code (fourth element of the return value) can have the values:
1976
1977 @table @code
1978 @item 0
1979 no problems were encountered;
1980 @item 1
1981 too many sub-intervals were done;
1982 @item 2
1983 excessive roundoff error is detected;
1984 @item 3
1985 extremely bad integrand behavior occurs;
1986 @item 6
1987 if the input is invalid.
1988
1989 @end table
1990
1991 Examples:
1992
1993 @c ===beg===
1994 @c quad_qawo (x^(-1/2)*exp(-2^(-2)*x), x, 1d-8, 20*2^2, 1, cos);
1995 @c rectform (integrate (x^(-1/2)*exp(-2^(-alpha)*x) * cos(x), x, 0, inf));
1996 @c input:pos;
1997 @c ev (%, alpha=2, numer);
1998 @c ===end===
1999 @example
2000 @group
2001 (%i1) quad_qawo (x^(-1/2)*exp(-2^(-2)*x), x, 1d-8, 20*2^2, 1, cos);
2002 (%o1)     [1.376043389877692, 4.72710759424899E-11, 765, 0]
2003 @end group
2004 @group
2005 (%i2) rectform (integrate (x^(-1/2)*exp(-2^(-alpha)*x) * cos(x),
2006       x, 0, inf));
2007                    alpha/2 - 1/2            2 alpha
2008         sqrt(%pi) 2              sqrt(sqrt(2        + 1) + 1)
2009 (%o2)   -----------------------------------------------------
2010                                2 alpha
2011                          sqrt(2        + 1)
2012 @end group
2013 @group
2014 (%i3) ev (%, alpha=2, numer);
2015 (%o3)                     1.376043390090716
2016 @end group
2017 @end example
2018
2019 @opencatbox{Categories:}
2020 @category{Numerical methods}
2021 @category{Package quadpack}
2022 @closecatbox
2023 @end deffn
2024
2025 @c THERE ARE OPTIONAL ARGUMENTS WHICH MAKES LISTING THE VARIANTS A LITTLE TEDIOUS
2026 @c NEED A MORE CONVENIENT (AND NONAMBIGUOUS) NOTATION FOR OPTIONAL ARGUMENTS
2027
2028 @c -----------------------------------------------------------------------------
2029 @anchor{quad_qaws}
2030 @deffn  {Function} quad_qaws @
2031 @fname{quad_qaws} (@var{f(x)}, @var{x}, @var{a}, @var{b}, @var{alpha}, @var{beta}, @var{wfun}, [@var{epsrel}, @var{epsabs}, @var{limit}]) @
2032 @fname{quad_qaws} (@var{f}, @var{x}, @var{a}, @var{b}, @var{alpha}, @var{beta}, @var{wfun}, [@var{epsrel}, @var{epsabs}, @var{limit}])
2033
2034 Integration of @math{w(x) f(x)} over a finite interval, where @math{w(x)} is a
2035 certain algebraic or logarithmic function.  A globally adaptive subdivision
2036 strategy is applied, with modified Clenshaw-Curtis integration on the
2037 subintervals which contain the endpoints of the interval of integration.
2038
2039 @code{quad_qaws} computes the integral using the Quadpack QAWS routine:
2040
2041 m4_displaymath(
2042 <<<\int_a^b f(x) \, w(x) \, dx>>>,
2043 <<<@math{integrate (f(x)*w(x), x, a, b)}>>>)
2044
2045 The weight function @math{w} is selected by @var{wfun}:
2046
2047 @table @code
2048 @item 1
2049 m4_math(w(x) = (x - a)^\alpha (b - x)^\beta, @math{w(x) = (x - a)^alpha (b - x)^beta})
2050 @item 2
2051 m4_math(w(x) = (x - a)^\alpha (b - x)^\beta \log(x - a),
2052 @math{w(x) = (x - a)^alpha (b - x)^beta log(x - a)})
2053 @item 3
2054 m4_math(w(x) = (x - a)^\alpha (b - x)^\beta \log(b - x), @math{w(x) = (x - a)^alpha (b - x)^beta log(b - x)})
2055 @item 4
2056 m4_math(w(x) = (x - a)^\alpha (b - x)^\beta \log(x - a) \log(b - x),
2057 @math{w(x) = (x - a)^alpha (b - x)^beta log(x - a) log(b - x)})
2058 @end table
2059
2060 The integrand may be specified as the name of a Maxima or Lisp function or
2061 operator, a Maxima lambda expression, or a general Maxima expression.
2062
2063 The keyword arguments are optional and may be specified in any order.
2064 They all take the form @code{key=val}.  The keyword arguments are:
2065
2066 @table @code
2067 @item epsrel
2068 Desired relative error of approximation.  Default is 1d-8.
2069 @item epsabs
2070 Desired absolute error of approximation.  Default is 0.
2071 @item limit
2072 Size of internal work array.  @var{limit}is the
2073 maximum number of subintervals to use.  Default is 200.
2074 @end table
2075
2076 @code{quad_qaws} returns a list of four elements:
2077
2078 @itemize
2079 @item
2080 an approximation to the integral,
2081 @item
2082 the estimated absolute error of the approximation,
2083 @item
2084 the number integrand evaluations,
2085 @item
2086 an error code.
2087 @end itemize
2088
2089 The error code (fourth element of the return value) can have the values:
2090
2091 @table @code
2092 @item 0
2093 no problems were encountered;
2094 @item 1
2095 too many sub-intervals were done;
2096 @item 2
2097 excessive roundoff error is detected;
2098 @item 3
2099 extremely bad integrand behavior occurs;
2100 @item 6
2101 if the input is invalid.
2102
2103 @end table
2104
2105 Examples:
2106
2107 @c ===beg===
2108 @c quad_qaws (1/(x+1+2^(-4)), x, -1, 1, -0.5, -0.5, 1, 'epsabs=1d-9);
2109 @c integrate ((1-x*x)^(-1/2)/(x+1+2^(-alpha)), x, -1, 1);
2110 @c input:pos;
2111 @c ev (%, alpha=4, numer);
2112 @c ===end===
2113 @example
2114 @group
2115 (%i1) quad_qaws (1/(x+1+2^(-4)), x, -1, 1, -0.5, -0.5, 1,
2116                  'epsabs=1d-9);
2117 (%o1)     [8.750097361672832, 1.24321522715422E-10, 170, 0]
2118 @end group
2119 @group
2120 (%i2) integrate ((1-x*x)^(-1/2)/(x+1+2^(-alpha)), x, -1, 1);
2121        alpha
2122 Is  4 2      - 1  positive, negative, or zero?
2123
2124 pos;
2125                           alpha         alpha
2126                    2 %pi 2      sqrt(2 2      + 1)
2127 (%o2)              -------------------------------
2128                                alpha
2129                             4 2      + 2
2130 @end group
2131 @group
2132 (%i3) ev (%, alpha=4, numer);
2133 (%o3)                     8.750097361672829
2134 @end group
2135 @end example
2136
2137 @opencatbox{Categories:}
2138 @category{Numerical methods}
2139 @category{Package quadpack}
2140 @closecatbox
2141 @end deffn
2142
2143 @c THERE ARE OPTIONAL ARGUMENTS WHICH MAKES LISTING THE VARIANTS A LITTLE TEDIOUS
2144 @c NEED A MORE CONVENIENT (AND NONAMBIGUOUS) NOTATION FOR OPTIONAL ARGUMENTS
2145
2146 @c -----------------------------------------------------------------------------
2147 @anchor{quad_qagp}
2148 @deffn  {Function} quad_qagp @
2149 @fname{quad_qagp} (@var{f(x)}, @var{x}, @var{a}, @var{b}, @var{points}, [@var{epsrel}, @var{epsabs}, @var{limit}]) @
2150 @fname{quad_qagp} (@var{f}, @var{x}, @var{a}, @var{b}, @var{points}, [@var{epsrel}, @var{epsabs}, @var{limit}])
2151
2152 Integration of a general function over a finite interval.
2153 @code{quad_qagp} implements globally adaptive interval subdivision with
2154 extrapolation (de Doncker, 1978) by the Epsilon algorithm (Wynn, 1956).
2155
2156 @code{quad_qagp} computes the integral
2157
2158 m4_displaymath(
2159 <<<\int_a^b f(x) \, dx>>>,
2160 <<<@math{integrate (f(x), x, a, b)}>>>)
2161
2162 The function to be integrated is @math{f(x)}, with
2163 dependent variable @math{x}, and the function is to be integrated
2164 between the limits @math{a} and @math{b}.
2165
2166 The integrand may be specified as the name of a Maxima or Lisp function or
2167 operator, a Maxima lambda expression, or a general Maxima expression.
2168
2169 To help the integrator, the user must supply a list of points where
2170 the integrand is singular or discontinuous.  The list is provided by
2171 @var{points}.  It may be an empty list.  The elements of the list must
2172 be between @var{a} and @var{b}, exclusive.  An error is thrown if
2173 there are elements out of range.  The list points may be in any order.
2174
2175 The keyword arguments are optional and may be specified in any order.
2176 They all take the form @code{key=val}.  The keyword arguments are:
2177
2178 @table @code
2179 @item epsrel
2180 Desired relative error of approximation.  Default is 1d-8.
2181 @item epsabs
2182 Desired absolute error of approximation.  Default is 0.
2183 @item limit
2184 Size of internal work array.  @var{limit} is the
2185 maximum number of subintervals to use.  Default is 200.
2186 @end table
2187
2188 @code{quad_qagp} returns a list of four elements:
2189
2190 @itemize
2191 @item
2192 an approximation to the integral,
2193 @item
2194 the estimated absolute error of the approximation,
2195 @item
2196 the number integrand evaluations,
2197 @item
2198 an error code.
2199 @end itemize
2200
2201 The error code (fourth element of the return value) can have the values:
2202
2203 @table @code
2204 @item 0
2205 no problems were encountered;
2206 @item 1
2207 too many sub-intervals were done;
2208 @item 2
2209 excessive roundoff error is detected;
2210 @item 3
2211 extremely bad integrand behavior occurs;
2212 @item 4
2213 failed to converge
2214 @item 5
2215 integral is probably divergent or slowly convergent
2216 @item 6
2217 if the input is invalid.
2218 @end table
2219
2220 @c NEED CROSS REFS HERE -- EITHER CROSS REF A QUADPACK OVERVIEW, OR CROSS REF EACH OF THE quad_* FUNCTIONS
2221
2222 Examples:
2223
2224 @c ===beg===
2225 @c quad_qagp(x^3*log(abs((x^2-1)*(x^2-2))),x,0,3,[1,sqrt(2)]);
2226 @c quad_qags(x^3*log(abs((x^2-1)*(x^2-2))), x, 0, 3);
2227 @c ===end===
2228 @example
2229 @group
2230 (%i1) quad_qagp(x^3*log(abs((x^2-1)*(x^2-2))),x,0,3,[1,sqrt(2)]);
2231 (%o1)   [52.74074838347143, 2.6247632689546663e-7, 1029, 0]
2232 @end group
2233 @group
2234 (%i2) quad_qags(x^3*log(abs((x^2-1)*(x^2-2))), x, 0, 3);
2235 (%o2)   [52.74074847951494, 4.088443219529836e-7, 1869, 0]
2236 @end group
2237 @end example
2238
2239 The integrand has singularities at @code{1} and @code{sqrt(2)} so we supply
2240 these points to @code{quad_qagp}.  We also note that @code{quad_qagp} is
2241 more accurate and more efficient that @mrefdot{quad_qags}
2242
2243 @opencatbox{Categories:}
2244 @category{Numerical methods}
2245 @category{Package quadpack}
2246 @closecatbox
2247 @end deffn
2248
2249 @c -----------------------------------------------------------------------------
2250 @anchor{quad_control}
2251 @deffn  {Function} quad_control (@var{parameter}, [@var{value}])
2252
2253 Control error handling for quadpack.  The parameter should be one of
2254 the following symbols:
2255
2256 @table @code
2257 @item current_error
2258 The current error number
2259 @item control
2260 Controls if messages are printed or not.  If it is set to zero or
2261 less, messages are suppressed.
2262 @item max_message
2263 The maximum number of times any message is to be printed.
2264 @end table
2265
2266 If @var{value} is not given, then the current value of the
2267 @var{parameter} is returned.  If @var{value} is given, the value of
2268 @var{parameter} is set to the given value.
2269
2270 @opencatbox{Categories:}
2271 @category{Numerical methods}
2272 @category{Package quadpack}
2273 @closecatbox
2274 @end deffn
2275