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