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 */
20composite():=(
21/* input problem from keyboard */
22print("The d.e. is: ey''+a(x)y'+b(x)y=0"),
23print("with b.c. y(0)=y0 and y(1)=y1"),
24a:read("enter a(x) > 0 on [0,1]"),
25b:read("enter b(x)"),
26y0:read("enter y0"),
27y1:read("enter y1"),
28print("The d.e. is: ey''+(",a,")y'+(",b,")y=0"),
29print("with b.c. y(0)=",y0,"and y(1)=",y1),
30trunc:read("enter truncation order"),
31/* set up basic form of solution */
32y: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. */
35de1:e*diff(y,x,2)+a*diff(y,x)+b*y,
36/* expand and collect terms */
37de2:expand(e*de1),
38gstuff:coeff(de2,%e^(-g(x)/e)),
39other:expand(de2-gstuff*%e^(-g(x)/e)),
40gstuff2:taylor(gstuff,e,0,trunc+1),
41other2:taylor(other,e,0,trunc+1),
42for i:0 thru trunc+1 do
43   geq[i]:coeff(gstuff2,e,i),
44for 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 */
48geq:geq[0]/h[0](x)/diff(g(x),x),
49gsol:ode2(geq,g(x),x),
50resultlist:[g(x)=ev(rhs(gsol),%c=0)],
51/* find h[i](x)'s */
52for 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 */
59for 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 */
66bc1:ev(y=y0,resultlist,x=0),
67bc2:ev(y=y1,resultlist,x=1),
68/* kill exponentially small %e terms since taylor won't */
69bc2:subst(0,%e,bc2),
70for i:0 thru trunc do
71   (bc1[i]:ratcoef(bc1,e,i),
72    bc2[i]:ratcoef(bc2,e,i)),
73/* solve for unknown consts */
74bceqs:makelist(bc1[i],i,0,trunc),
75bceqs:append(bceqs,makelist(bc2[i],i,0,trunc)),
76bcunk:makelist(hh[i],i,0,trunc),
77bcunk:append(bcunk,makelist(ff[i],i,0,trunc)),
78const:solve(bceqs,bcunk),
79resultlist:ratsimp(ev(resultlist,const)),
80ysol:ev(y,resultlist))$
81
82