share/diffequations/pivtr.mac

   1 setup_autoload("lodsav",lodesave)$
   2
   3 ivtrns(wb,wc,wr):=block([],tn:ratsimp(tp/10),invsw:n,
   4    if tp<=60 then (neq:dptarg(wb,wc,tn),
   5                      if neq#f then return(neq)),
   6    if tp= 10 then (ldf5(1),trtrig(wb,wc,wr)) else
   7                     (ldf5(0),neq:ivexctrig(wb,wc,wr)),
   8    if tp<100 then return(neq) else return(f))$
   9
  10 dptarg(wb,wc,tn):=block([],wl:listofargs(za*q1+zb*q2+zc,tn),
  11    lg:length(wl),if lg=0 then return(f),wlf:first(wl),
  12  ll, if lg>1 then (wl:rest(wl),wlr:first(wl),wlf:gcd(wlf,wlr),lg:lg-1),
  13    if lg>1 then go(ll),if wlf#x and wlf#2*x then (np:hipow(wlf,x),
  14    if np>1 then (ldf0(14),if powr(wb,wc,np)=t then
  15                 (ldf5(11),trpow(pown),return(neq)) else return(f)),
  16    if np=1 then (ldf5(9),tivln(wb,wc,wr,wlf),return(neq)) else
  17   if np=-1 then (ldf5(7),tivinv(wb,wc,wr),return(neq)) else
  18                  return(f)) else return(f))$
  19
  20 listofargs(exp,n):=block([],wl:[],elmf:matrix([5,sin,cos,tan,cot,0],
  21    [2,%e,0,0,0,0],[2,log,0,0,0,0],[5,sinh,cosh,tanh,coth,0],
  22    [6,sn,cn,dn,pe,th],[2,sqrt,0,0,0,0]),
  23    for i:2 thru elmf[n,1] do wl:append(wl,farg(exp,elmf[n,i])),
  24    return(wl))$
  25
  26 farg(exp,f):=block([],wexp:args(exp),
  27  ll,if atom(wexp) then lg:0 else lg:length(wexp),
  28    if lg>0 then (rexp:rest(wexp),wexp:first(wexp)) else return([]),
  29    if freeof(f,wexp) then (wexp:rexp,go(ll)),
  30    if atom(wexp)     then return(rexp), wexp:args(wexp),
  31    if freeof(f,wexp) then return(wexp) else go(ll))$
  32
  33 ldf5(i):=block([],if ldsw[5,i]#y then (ldsw[5,i]:y,
  34    if i=0 then load(pivtr2) else
  35    if i=1 then load(ptrtrg) else
  36    if i=2 then load(ptrexp) else
  37    if i=3 then load(ptrlog) else
  38    if i=4 then load(ptrhpb) else
  39    if i=5 then load(ptrelp) else
  40    if i=6 then load(ptrsqr) else
  41    if i=7 then load(ptrinv) else
  42    if i=8 then load(ptrasi) else
  43    if i=9 then load(ptrln) else
  44   if i=10 then load(peqp7) else
  45   if i=11 then load(ptrpow) else
  46   if i=12 then load(ptrsym) ))$
  47
  48 lodesave([pivtr,fasl],ivtrns,dptarg,listofargs,farg,ldf5);
  49
  50 ivexctrig(wb,wc,wr):=block([wexp,wb1,wc1],
  51    if tp=15 or tp=16 then (ldf5(1),tivasin(wc,wr)) else
  52    if tp= 20 then (ldf5(2),tivexp(wb,wc,wr))  else
  53    if tp= 30 then (ldf5(3),tivlog(wb,wc))     else
  54    if tp= 40 then (ldf5(4),trhpb(wb,wc,wr))   else
  55    if tp= 50 then (ldf5(5),tivelp(wb,wc)),
  56    if tp= 60 then (ldf5(6),tivsqrt(wb,wc,wr)) else
  57    if tp=70 or tp=80 then (ldf5(11),trpow(pown)) else
  58    if tp= 81 then (ldf5(12),trsymx(pt,qt,1))  else
  59    if tp= 90 then (ldf5(7),tivinv(wb,wc,wr)) else
  60    if tp= 91 or tp=92 then (ldf5(9),tivln(wb,wc,wr,axb)),
  61    return(neq))$
  62
  63 lodesave([pivtr2,fasl],ivexctrig);
  64
  65 trtrig(wb,wc,wr):=block([],fp1:wb-sin(x)/cos(x),fp2:wb-cos(x)/sin(x),
  66    if not freeof(tan,cot,wb*q1+wc) or freeof(sin,cos,fp1)
  67       or  freeof(sin,cos,fp2) then (ldf5(10),eqpt7(wb,wc,wr)),
  68    if 10<tp and tp<15 then return(tp),swtrg:f,
  69    wb1:prptrg(wb),wc1:prptrg(wc),
  70    dbc1:tivsc(1,wb1,wc1),dbc2:tivsc(2,wb1,wc1),
  71    neq:slctrg(dbc1,dbc2),if neq#f then swtrg:t)$
  72
  73 prptrg(exp):=block([w1,w2,w3],trigexpand:true,w1:ratsimp(exp),
  74    w2:ratsubst(sin(x)/cos(x),tan(x),w1),
  75    w3:ratsubst(cos(x)/sin(x),cot(x),w2),trigexpand:false,
  76    return(w3))$
  77
  78 slctrg(w1,w2):=block([w3,w4,wt],wt:t^2-1,
  79    w3:lcm(wt,w1,t),w4:lcm(wt,w2,t),
  80    dg1:hipow(expand(w3),t),dg2:hipow(expand(w4),t),
  81    if dg1>99 and dg2>99 then return(f),
  82    if dg1<=dg2 then (vtr:t=sin(x),i:1) else (vtr:t=cos(x),i:2),
  83    bt:coeff(cnd[i],q,1), ct:coeff(cnd[i],q,0),
  84    bt:ratsimp((t-bt)/wt),ct:ratsimp(-ct/wt),
  85    neq:dyt2+bt*dyt1+ct*y=0)$
  86
  87 tivsc(i,wb,wc):=block([wb1,wc1,wt],wt:1-t^2, if i=1 then
  88    (cs:cos(x),sc:sin(x),ws:+1) else (cs:sin(x),sc:cos(x),ws:-1),
  89    wb1:ratsubst(wt,cs^2,ws*wb*cs),wb1:subst(t,sc,wb1),
  90    wc1:ratsubst(wt,cs^2,wc),      wc1:subst(t,sc,wc1),
  91    cnd[i]:wb1*q+wc1,if freeof(sin,cos,x,cnd[i]) then
  92    return(calca(wb1,wc1,t)) else return(t^99))$
  93
  94 lodesave([ptrtrg,fasl],trtrig,prptrg,slctrg,tivsc);
  95
  96 eqpt7(wb,wc,wr):=block([],
  97   if freeof(sin,cos,fp1) then (fq:-1/2-sin(x)^2/(4*cos(x)^2)
  98     -wb*sin(x)/(2*cos(x))+wc,fq:subst(1-cos(x)^2,sin(x)^2,fq),
  99     fq:expand(fq),if freeof(sin,cos,fq) then
 100     (fp:fp1,fr:wr/sqrt(cos(x)),tp:13,return(t))),
 101   for i:1 thru 2 do
 102     (if i=1 then (tc:tan(x),ws:+1) else (tc:cot(x),ws:-1),
 103      db:hipow(wb,tc),dc:hipow(wc,tc),
 104      if db=1 and dc<=1 and coeff(wb,tc,1)=-2*ws then
 105       (fp:wb+2*ws*tc,if (dc=1 and coeff(wc,tc,1)=-fp*ws) or dc=0 then
 106                      (tp:10+i,fq:coeff(wc,tc,0)+1),return(t))))$
 107
 108 lodesave([peqp7,fasl],eqpt7);
 109
 110 trhpb(wb,wc,wr):=block([],
 111    wb1:prphpb(wb),wc1:prphpb(wc),wt[1]:t^2+1,wt[2]:t^2-1,
 112    dbc1:tivhpb(1,wb1,wc1),dbc2:tivhpb(2,wb1,wc1),
 113    neq:slchpb(dbc1,dbc2),if neq#f then swhpb:t,return(neq))$
 114
 115 prphpb(exp):=block([w1,w2,w3],w1:ratsimp(exp),
 116    w2:ratsubst(sinh(x)/cosh(x),tanh(x),w1),
 117    w3:ratsubst(cosh(x)/sinh(x),coth(x),w2),return(w3))$
 118
 119 slchpb(w1,w2):=block([w3,w4],
 120    w3:lcm(wt[1],w1,t),w4:lcm(wt[2],w2,t),
 121    dg1:hipow(expand(w3),t),dg2:hipow(expand(w4),t),
 122    if dg1>99 and dg2>99 then return(f),
 123    if dg1<dg2 then (vtr:t=sinh(x),i:1) else (vtr:t=cosh(x),i:2),
 124    bt:coeff(cnd[i],q,1), ct:coeff(cnd[i],q,0),
 125    bt:ratsimp((bt+t)/wt[i]),ct:ratsimp(ct/wt[i]),
 126    neq:dyt2+bt*dyt1+ct*y=0)$
 127
 128 tivhpb(i,wb,wc):=block([wb1,wc1], if i=1 then
 129   (cs:cosh(x),sc:sinh(x),ws:+1) else (cs:sinh(x),sc:cosh(x),ws:-1),
 130    wb1:ratsubst(wt[i],cs^2,wb*cs), wb1:subst(t,sc,wb1),
 131    wc1:ratsubst(wt[i],cs^2,wc),    wc1:subst(t,sc,wc1),
 132    cnd[i]:wb1*q+wc1,if freeof(sinh,cosh,x,cnd[i]) then
 133    return(calca(wb1,wc1,t)) else return(t^99))$
 134
 135 lodesave([ptrhpb,fasl],trhpb,prphpb,slchpb,tivhpb);
 136
 137 tivexp(wb,wc,wr):=block([wb1,wc1,wr1,bt,ct,rt],exptsubst:true,
 138    wb1:subst(t,%e^x,wb),wc1:subst(t,%e^x,wc),wr1:subst(t,%e^x,wr),
 139    bt:ratsimp((wb1+1)/t),ct:ratsimp(wc1/t^2),rt:ratsimp(wr1/t^2),
 140    exptsubst:false,neq:dyt2+bt*dyt1+ct*y=rt,
 141    if freeof(%e,x,neq) then (vtr:t=%e^x,return(t)) else neq:f)$
 142
 143 lodesave([ptrexp,fasl],tivexp);
 144
 145 tivlog(wb,wc,wr):=block([wb1,wc1,wr1,bt,ct,rt],
 146    wb1:subst(%e^t,x,x*wb),
 147    wc1:subst(%e^t,x,x^2*wc),wr1:subst(%e^t,x,x^2*wr),
 148    bt:ratsimp(wb1-1),ct:ratsimp(wc1),rt:ratsimp(wr1),
 149    neq:dyt2+bt*dyt1+ct*y=rt,
 150    if freeof(x,neq) then (vtr:t=log(x),return(t)),
 151    wb1:subst(t/x+1,log(x),wb/log(x)+1/(x*(log(x)^2))),
 152    wc1:subst(t/x+1,log(x),wc/log(x)^2),
 153    wr1:subst(t/x+1,log(x),wr/log(x)^2),
 154    bt:ratsimp(wb1),ct:ratsimp(wc1),rt:ratsimp(wr1),
 155    neq:dyt2+bt*dyt1+ct*y=rt,
 156    if freeof(x,neq) then (vtr:t=x*log(x)-x,return(t))
 157                           else return(f) )$
 158
 159 lodesave([prtlog,fasl],tivlog);
 160
 161 tivasin(wc,wr):=block([],
 162    if tp=15 then neq:trasinh(wc,wr) else
 163    if tp=16 then neq:trasin(wc,wr))$
 164
 165 trasinh(wc,wr):=block([],alpha:sqrt(wa2),
 166    vtr:alpha*x=sinh(alpha*t),neq:dyt2+zc*y=0)$
 167
 168 trasin(wc,wr):=block([],vtr:x=sin(x),neq:dyt2-zc*y=0)$
 169
 170 lodesave([ptrasin,fasl],tivasin,trasinh,trasin);
 171
 172 tivinv(wb,wc,wr):=block([wb1,wc1,wr1,bt,ct,rt],
 173    if invsw=y then return(f),
 174    wb1:subst(1/t,x,wb),wc1:subst(1/t,x,wc),wr1:subst(1/t,x,wr),
 175    bt:ratsimp((2*t-wb1)/t^2),ct:ratsimp(wc1/t^4),rt:ratsimp(wr1/t^4),
 176    neq:dyt2+bt*dyt1+ct*y=rt, vtr:t=1/x, invsw:y)$
 177
 178 lodesave([ptrinv,fasl],tivinv);
 179
 180 tivsqrt(wb,wc,wr):=block([wb1,wc1,wr1,bt,ct,rt], vtr:t=sqrt(x),
 181    wb1:subst(t^2,x,wb),wc1:subst(t^2,x,wc),wr1:subst(t^2,x,wr),
 182    bt:ratsimp(2*t*wb1-1/t),ct:ratsimp(wc1*4*t^2),rt:ratsimp(wr1*4*t^2),
 183    neq:dyt2+bt*dyt1+ct*y=rt, neq:subst(t,abs(t),neq))$
 184
 185 lodesave([ptrsqr,fasl],tivsqrt);
 186
 187 /* tivln do t=pa*(x-alpha) */
 188 tivln(wb,wc,wr,x1):=block([wb1,wc1,wr1,bt,ct,rt,x2,ws],
 189    x2:expand(x1),pa:coeff(x2,x,1),ws:solve(x1,x),ws:first(ws),
 190    alpha:part(ws,2), t1:t/pa+alpha, vtr:t=x1,
 191    wb1:subst(t1,x,wb),wc1:subst(t1,x,wc),wr1:subst(t1,x,wr),
 192    bt:ratsimp(wb1/pa),ct:ratsimp(wc1/pa^2),rt:ratsimp(wr1/pa^2),
 193    neq:dyt2+bt*dyt1+ct*y=rt,
 194    if tp=92 then neq:ratsubst(ob,oa,neq))$
 195
 196 lodesave([ptrln,fasl],tivln);
 197
 198 trpow(r):=block([b1,c1],vtr:t=x^r,b1:ratsimp((bt+r-1)/(r*t)),
 199    c1:ratsimp(ct/(r^2*t^2)),d1:ratsimp(rt/(r^2*t^2)),
 200    neq:dyt2+b1*dyt1+c1*y=d1)$
 201
 202 lodesave([ptrpow,fasl],trpow);
 203
 204 trsymx(pt,qt,i):=block([], if i=1 then
 205    (vtr:(x+1/x)/2,neq:(t^2-1)*dyt2+t*dyt1+qt*y=0) else
 206    (vtr:t=log(x),qt:subst(cosh(2*t),t,qt),neq:dyt2+qt*y=0))$
 207
 208 lodesave([ptrsym,fasl],trsymx);