share/diffequations/pdvtr.mac

   1 setup_autoload("lodsav",lodesave)$
   2
   3 /* try transformations for dependent variables. */
   4 dvtrns(a,wb,wc,wr,ro,lam1):=block([neq,alpha],
   5    if tp=12 or tp=13 or tp=30 then (ldf6(1),dvt60(dum),return(neq)),
   6    neq:f, wro:ro, i:1, tblspt(a,wb,wc,x,n),
   7    if swro=1 or tp<60 then go(ll2),
   8    if eqpt6(wb,wc)=t then (wro:1,lam1:ratsimp(-pc/pa),inflam:det),
   9    if not integerp(wro) then wro:wro-1/2,
  10    if wro=swro then wro:wro-1, owro:wro,
  11 ll1,   neq:dvtexp(wb,wc,wr,wro,lam1), inflam:n,
  12    if swsv=tdv then neq:dvtexp(wb,wc,wr,1,lam1),
  13    if neq#f then (wexp:part(neq,1),
  14    if coeff(wexp,u,1)=0 then (swro:wro, return(neq))),
  15    if wro>1 then (wro:wro-1,go(ll1)), swro:wro,
  16    if owro>1 then (neq:dvtexp(wb,wc,wr,owro,lam1),return(neq))
  17                else  return(neq),
  18 ll2, if i< nspt then (alpha:spt[i], neq:dvtmxp(a,wb,wc,wr,alpha)),
  19    wexp:part(neq,1), w1:coeff(wexp,dux1,1),w2:coeff(wexp,u),
  20    w3:calca(w1,w2,x),w3:subst(spt[i],x,w3),
  21    if coeff(wexp,u,1)=0 or w3#0 then return(neq) else (i:i+1,go(ll2)),
  22    return(f) )$
  23
  24 /* eq=k54 y"+(a*x+b)*y'+(c*x+d)*y=0 */
  25 eqpt6(wb,wc):=block([], if dega=0 and degb=1 and degc=1 then
  26    (wb:expand(wb), wc:expand(wc), pa:coeff(wb,x,1), pb:coeff(wb,x,0),
  27     pc:coeff(wc,x,1), pd:coeff(wc,x,0), return(t)) else return(f))$
  28
  29 /* dvtexp do y=exp(l*x^r)*u. */
  30 dvtexp(wb,wc,wr,r,lam1):=block([l,de],
  31    if swsv=tdv or inflam=det then l:lam1,
  32    w1:2*r*x^(r-1)*l+wb, w3:wr*%e^(-l*x^r),
  33    w2:r^2*x^(2*r-2)*l^2+(r*(r-1)*x^(r-2)+r*wb*x^(r-1))*l+wc,
  34    w1:ratsimp(w1),w2:ratsimp(w2),w3:ratsimp(w3),
  35    if swsv=tdv or inflam=det then (de:dux2+w1*dux1+w2*u=w3,lam:l)
  36                                  else  de:degenc(w1,w2,w3),
  37    if de=f then return(de) else (vtr:y=%e^(lam*x^r)*u, return(de)))$
  38
  39 /* dvtmxp do y=(x-alpha)^l*u. */
  40 dvtmxp(a,wb,wc,wr,alpha):=block([l,de],
  41    xa:x-alpha, w1:ratsimp((2*l+xa*wb)/xa),
  42    w2:ratsimp((l^2+(wb*xa-1)*l+xa^2*wc)/xa^2),w3:ratsimp(wr/xa^l),
  43    de:degenc(w1,w2,w3), if de=f then return(de),
  44    vtr:y=xa^lam*u, return(de))$
  45
  46 degenc(w1,w2,w3):=block([wc1,wc2,dwc1,dwc2],
  47    wc:num(w2),wq:expand(wc), deg:hipow(wq,x),
  48   ll,co:coeff(wq,x,deg), deg:deg-1,
  49    if freeof(l,co) then if deg>=0 then go(ll) else return(f),
  50    ws:solve(co,l), lg:length(ws),
  51    if lg=1 then (lam1:part(ws[1],2), lam2:lam1)
  52            else (lam1:part(ws[1],2), lam2:part(ws[2],2)),
  53    if lg=1 and lam1=0 then return(f),
  54    wc1:subst(lam1,l,wq), wc2:subst(lam2,l,wq),
  55    dwc1:hipow(wc1,x),        dwc2:hipow(wc2,x),
  56    if dwc1<dwc2 then (lam:lam1,clam:lam2) else
  57    if dwc1>dwc2 then (lam:lam2,clam:lam1) else
  58    if dwc1=dwc2 then (if wc1=0 then (lam:lam1,clam:lam2) else
  59                      (lam:min(lam1,lam2),clam:max(lam1,lam2))),
  60    if lam=0 then lam:clam,wp:ratsubst(lam,l,w1),wr:ratsubst(lam,l,w3),
  61    wq:ratsubst(lam,l,w2),de:dux2+wp*dux1+wq*u=wr, return(de))$
  62
  63 ldf6(i):=block([],if ldsw[6,1]#y then (ldsw[6,i]:y,
  64    if i=1 then load(pdvt60)))$
  65
  66 lodesave([pdvtr,fasl],dvtexp,dvtmxp,degenc,dvtrns,eqpt6,ldf6);
  67
  68 dvt60(dum):=block([],if tp=11 then vtr:y=u/cos(x),
  69    if tp=12 then vtr:y=u/sin(x), if tp=13 then vtr:y=sqrt(cos(x))*u,
  70    if 10<tp and tp<14 then (neq:dux2+fp*ux1+fq*u=fr,return(neq)),
  71    if tp=30 then (w1:ratsimp(2/(x*log(x))+wb),
  72       w2:ratsimp(-1/(x^2*log(x))+wb/(x*log(x))+wc),
  73       w3:ratsimp(wr/log(x)),
  74       if freeof(log(x),w1*q1+w2*q2+w3) or w2=0 then
  75       (vtr:y=log(x)*u,neq:dux2+w1*dux1+w2*u=w3, return(neq))))$
  76
  77 lodesave([pdvt60,fasl],dvt60);
  78
  79 gptm(wb,wc):=block([],rptm:f,
  80    if gpt078(wb,wc)#f then go(ll),
  81    if gpt128(wb,wc)#f then go(ll),
  82    if gpt442(wb,wc)#f then go(ll) else return(f),
  83  ll,rptm:y, return(t))$
  84
  85 gpt078(wb,wc):=block([],w1:wc-wb^2/4-diff(wb,x)/2, w1:ratsimp(w1),
  86    if freeof(x,w1) then (vtr:y=%e^(-integrate(wb/2,x))*u,
  87                          neq:dux2+w1*u=0,return(t))
  88                    else return(f))$
  89
  90 gpt128(wb,wc):=block([],w1:ratsimp(wc*x),w2:ratsimp(wb-2/x),
  91    if w1=w2 then (vtr:y=u/x,neq:dux2+w1*dux1=0,
  92                   return(t)) else return(f))$
  93
  94 gpt442(wb,wc):=block([], if wc=0 then return(f),
  95    w1:ratsimp(-1/wc),w2:ratsimp(-wb/wc),w3:w2-x,
  96    if hipow(w3,x)#0 then return(f), wf:ratsimp(w2/w1),
  97    w1:integrate(wf,x), w3:%e^(-w1), w4:integrate(w3,x),
  98    y1:w2,y2:y1*w4, ys:k1*y1+k2*y2)$
  99
 100 lodesave([pgptm,fasl],gptm,gpt078,gpt128,gpt442);
 101