share/contrib/rand/asympexp.mac

   1 /* Filename asympexp.mac
   2
   3    ***************************************************************
   4    *                                                             *
   5    *                     <package name>                          *
   6    *                <functionality description>                  *
   7    *                                                             *
   8    *          from: "The Use of Symbolic Computation             *
   9    *                 in Perturbation Analysis"                   *
  10    *                 by Rand in Symbolic Computation             *
  11    *                 in Fluid Mechanics and Heat Transfer        *
  12    *                 ed H.H.Bau (ASME 1988)                      *
  13    *                Programmed by Richard Rand                   *
  14    *      These files are released to the public domain          *
  15    *                                                             *
  16    ***************************************************************
  17 */
  18 /* asymptotic expansion of integrals */
  19 asymptotic():=(
  20 block(
  21 /* input the problem from the keyboard */
  22 print("The integrand is of the form: f(t) exp(x phi(t))"),
  23 f:read("enter f(t)"),
  24 phi:read("enter phi(t)"),
  25 a:read("enter the lower limit of integration"),
  26 b:read("enter the upper limit of integration"),
  27 print("The integrand is",f*%e^(x*phi)),
  28 print("integrated from",a,"to",b),
  29 c:read("enter value of t at which phi =",phi," is maximum"),
  30 /* set up limits of integration for later use */
  31 if c=a then (lowerlim:0, upperlim:inf)
  32    else if c=b then (lowerlim:minf, upperlim:0)
  33        else (lowerlim:minf, upperlim:inf),
  34 /* move origin to t=c with u=t-c */
  35 f1:ev(f,t=u+c),
  36 phi1:ev(phi,t=u+c),
  37 trunc:read("enter truncation order"),
  38 /* determine lowest non-constant term in taylor series
  39    for phi about u=0 */
  40 phi2:taylor(phi1,u,0,2),
  41 phic:ev(phi,t=c),
  42 keypower:lopow(phi2-phic,u),
  43 if keypower=1 then trunc1:2*trunc else
  44    if keypower=2 then trunc1:6*trunc else
  45       (print("phi does not behave like t-",c,
  46              "or (t-",c,")^2 near t=",c),
  47        return("program aborted")),
  48 phi3:taylor(phi1,u,0,trunc1),
  49 phikey:coeff(phi3,u,keypower)*u^keypower,
  50 phirest:phi3-phikey-phic,
  51 /* set flag so integration routine assumes x is positive */
  52 assume_pos:true,
  53 integrand:taylor(f1*%e^(x*phirest),u,0,trunc1),
  54 /* set u = s/x^(1/keypower) */
  55 integrand2:ev(integrand,u=s/x^(1/keypower)),
  56 /* set x = 1/xx in order to truncate */
  57 integrand3:taylor(ev(integrand2,x=1/xx),xx,0,trunc),
  58 /* return to x again */
  59 integrand4:ev(integrand3,xx=1/x),
  60 approx:%e^(x*phic)/x^(1/keypower)
  61        *integrate(integrand4*%e^(ev(phikey,u=s)),s,lowerlim,upperlim),
  62 approx2:factor(approx)))$
  63