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