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 */
19asymptotic():=(
20block(
21/* input the problem from the keyboard */
22print("The integrand is of the form: f(t) exp(x phi(t))"),
23f:read("enter f(t)"),
24phi:read("enter phi(t)"),
25a:read("enter the lower limit of integration"),
26b:read("enter the upper limit of integration"),
27print("The integrand is",f*%e^(x*phi)),
28print("integrated from",a,"to",b),
29c:read("enter value of t at which phi =",phi," is maximum"),
30/* set up limits of integration for later use */
31if 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 */
35f1:ev(f,t=u+c),
36phi1:ev(phi,t=u+c),
37trunc:read("enter truncation order"),
38/* determine lowest non-constant term in taylor series
39   for phi about u=0 */
40phi2:taylor(phi1,u,0,2),
41phic:ev(phi,t=c),
42keypower:lopow(phi2-phic,u),
43if 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")),
48phi3:taylor(phi1,u,0,trunc1),
49phikey:coeff(phi3,u,keypower)*u^keypower,
50phirest:phi3-phikey-phic,
51/* set flag so integration routine assumes x is positive */
52assume_pos:true,
53integrand:taylor(f1*%e^(x*phirest),u,0,trunc1),
54/* set u = s/x^(1/keypower) */
55integrand2:ev(integrand,u=s/x^(1/keypower)),
56/* set x = 1/xx in order to truncate */
57integrand3:taylor(ev(integrand2,x=1/xx),xx,0,trunc),
58/* return to x again */
59integrand4:ev(integrand3,xx=1/x),
60approx:%e^(x*phic)/x^(1/keypower)
61       *integrate(integrand4*%e^(ev(phikey,u=s)),s,lowerlim,upperlim),
62approx2:factor(approx)))$
63
64