share/contrib/rand/composit.mac

   1 /* Filename composit.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
  19 /* method of composite expansions */
  20 composite():=(
  21 /* input problem from keyboard */
  22 print("The d.e. is: ey''+a(x)y'+b(x)y=0"),
  23 print("with b.c. y(0)=y0 and y(1)=y1"),
  24 a:read("enter a(x) > 0 on [0,1]"),
  25 b:read("enter b(x)"),
  26 y0:read("enter y0"),
  27 y1:read("enter y1"),
  28 print("The d.e. is: ey''+(",a,")y'+(",b,")y=0"),
  29 print("with b.c. y(0)=",y0,"and y(1)=",y1),
  30 trunc:read("enter truncation order"),
  31 /* set up basic form of solution */
  32 y:sum(f[n](x)*e^n,n,0,trunc)+
  33   %e^(-g(x)/e)*sum(h[n](x)*e^n,n,0,trunc),
  34 /* substitute into d.e. */
  35 de1:e*diff(y,x,2)+a*diff(y,x)+b*y,
  36 /* expand and collect terms */
  37 de2:expand(e*de1),
  38 gstuff:coeff(de2,%e^(-g(x)/e)),
  39 other:expand(de2-gstuff*%e^(-g(x)/e)),
  40 gstuff2:taylor(gstuff,e,0,trunc+1),
  41 other2:taylor(other,e,0,trunc+1),
  42 for i:0 thru trunc+1 do
  43    geq[i]:coeff(gstuff2,e,i),
  44 for i:1 thru trunc+1 do
  45    otheq[i]:coeff(other2,e,i),
  46 /* find g(x) */
  47 /* assume a(x) coeff of y' is >0 st bl occurs at x=0 */
  48 geq:geq[0]/h[0](x)/diff(g(x),x),
  49 gsol:ode2(geq,g(x),x),
  50 resultlist:[g(x)=ev(rhs(gsol),%c=0)],
  51 /* find h[i](x)'s */
  52 for i:0 thru trunc do(
  53    heq[i]:ev(geq[i+1],resultlist,diff),
  54    hsol[i]:ode2(heq[i],h[i](x),x),
  55    hsol[i]:ev(hsol[i],%c=hh[i]),
  56    resultlist:append(resultlist,[hsol[i]])),
  57 /* find f[i](x)'s */
  58 /* ff[i] = constant of integration */
  59 for i:0 thru trunc do(
  60    feq[i]:ev(otheq[i+1],resultlist,diff),
  61    fsol[i]:ode2(feq[i],f[i](x),x),
  62    fsol[i]:ev(fsol[i],%c=ff[i]),
  63    resultlist:append(resultlist,[fsol[i]])),
  64 /* boundary conditions */
  65 /* y(0)=y0 and y(1)=y1 */
  66 bc1:ev(y=y0,resultlist,x=0),
  67 bc2:ev(y=y1,resultlist,x=1),
  68 /* kill exponentially small %e terms since taylor won't */
  69 bc2:subst(0,%e,bc2),
  70 for i:0 thru trunc do
  71    (bc1[i]:ratcoef(bc1,e,i),
  72     bc2[i]:ratcoef(bc2,e,i)),
  73 /* solve for unknown consts */
  74 bceqs:makelist(bc1[i],i,0,trunc),
  75 bceqs:append(bceqs,makelist(bc2[i],i,0,trunc)),
  76 bcunk:makelist(hh[i],i,0,trunc),
  77 bcunk:append(bcunk,makelist(ff[i],i,0,trunc)),
  78 const:solve(bceqs,bcunk),
  79 resultlist:ratsimp(ev(resultlist,const)),
  80 ysol:ev(y,resultlist))$
  81