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