share/diffequations/phypgm.mac

   1 setup_autoload("lodsav",lodesave)$
   2
   3 svhypr(wb,wc):=block([p],swsv:t,print("the type is hypergeometric"),
   4    for i:1 thru 3 do (sl[i]:ssolve(cheq[i],l), la[i,0]:spt[i],
   5           la[i,1]:part(sl[i],1,2), la[i,2]:part(sl[i],2,2)),
   6    print("the solution may be written by Riemann's P-functions as follows"),
   7    pfuncthg(n,dum), replcompform(dum), infosum(dum),
   8    lmin:min(l[1],l[2],l[3],l[4]), lmax:max(l[1],l[2],l[3],l[4]),
   9    if lmin=0 then (j:seekodd(dum),
  10       if j>0 then (ldf3(1),ys:nopmt(j),return(ys))),
  11    ldf3(2),ys:haspmts(dum),return(ys))$
  12
  13 pfuncthg(selv,v):=block([p], p:matrix([la[1,0],la[2,0],la[3,0]],
  14    [la[1,1],la[2,1],la[3,1]],[la[1,2],la[2,2],la[3,2]]),
  15    if selv=y then print("y=",v,"P",p,"(x)") else print("y=P",p,"(x)"))$
  16
  17 infosum(dum):=block([ws],
  18    for i:1 thru 3 do dla[i]:la[i,1]-la[i,2],
  19    ws:dla[1]+dla[2]+dla[3], sum[1]:ratsimp(ws),
  20    ws:dla[1]+dla[2]-dla[3], sum[2]:ratsimp(ws),
  21    ws:dla[1]-dla[2]+dla[3], sum[3]:ratsimp(ws),
  22    ws:dla[1]-dla[2]-dla[3], sum[4]:ratsimp(ws),
  23    for i:1 thru 4 do (var[i]:listofvars(sum[i]),l[i]:length(var[i])))$
  24
  25 excrow(dum):=block([wla], for i:1 thru 3 do exc1c(i))$
  26
  27 ssolve(e,x):=block([l,ws,fws],ws:solve(e,x),l:length(ws),
  28    if l=1 then (fws:first(ws),ws:[fws,fws]),return(ws))$
  29
  30 oddintp(pol,v):=block([c,wtst,wd],c:coeff(pol,v,0),wtst:ratsimp((pol-c)/2),
  31    wd:denom(wtst),if wd=1 and oddp(c) then return(t) else return(f))$
  32
  33 exc1c(i):=block([wla], wla:la[i,1],la[i,1]:la[i,2],
  34                         la[i,2]:wla,     dla[i]:-dla[i])$
  35
  36 seekodd(dum):=block([i,j],i:1,j:0,
  37  loop, if l[i]=0 and oddp(sum[i]) then (j:j+1,jodd[j]:i),
  38        if l[i]=1 then (ans:read("is",sum[i],"an odd integer?"),
  39                        if ans=y then (j:j+1,jodd[j]:i)),
  40    if i<4 then (i:i+1,go(loop)), return(j))$
  41
  42 replcompform(dum):=block([],for i:1 thru 3 do
  43    (nv:listofvars(la[i,1]), if length(nv)>1 then
  44    (ans:read("do you replace in",la[i,1],"?",str4),
  45    if ans=y then (ob:read("replace"),oa:read("by"),
  46    for j:1 thru 2 do la[i,j]:subst(oa,ob,la[i,j])))))$
  47
  48 ldf3(i):=block([], if ldsw[3,i]#y then (ldsw[3,i]:y,
  49    if i=1 then load(pnopm) else
  50    if i=2 then load(phghp) else
  51    if i=3 then load(psleg) else
  52    if i=4 then load(phyp25) else
  53    if i=5 then load(psolh) ))$
  54
  55 lodesave([phypgm,fasl],svhypr,infosum,excrow,ssolve,pfuncthg,
  56         oddintp,exc1c,seekodd,replcompform,ldf3);
  57
  58 /* the sum has no parameter */
  59 nopmt(j):=block([], iodd:minabssum(j),excroots(iodd),
  60    if integerp(sum[iodd]) then (if sum[iodd]>0 then
  61       (excrow(dum),sum[iodd]:-sum[iodd]))
  62    else (ans:read("is",sum[iodd],"positive?",str4), if ans=y then
  63       (excrow(dum),sum[iodd]:-sum[iodd])),
  64    wsum:sum[iodd],xa:x-spt[1],xb:x-spt[2],
  65    print(str1), pfuncthg(n,v), wd:xa*xb,
  66    if sum[iodd]=-1 then ys:elmrephg1(iodd)
  67                    else ys:elmrephg2((1-wsum)/2), return(ys))$
  68
  69 minabssum(j):=block([i,j1,j2,minj],i:1,j1:jodd[i], minj:j1,
  70    loop,if i<j then (i:i+1,j2:jodd[i],
  71    if abs(sum[j1]) > abs(sum[j2]) then minj:j2, go(loop)),
  72    return(minj))$
  73
  74 excroots(iodd):=block([],
  75     if iodd=2 or iodd=4 then exc1c(3),
  76     if iodd=3 or iodd=4 then exc1c(2))$
  77
  78 elmrephg1(iodd):=block([w1,w2,w3,w4],
  79    sdla1:ratsimp(-dla[1]-1), sdmu1:ratsimp(-dla[2]-1),
  80    if la[1,1]=la[2,1] then w1:wd^la[1,1]
  81                       else w1:xa^la[1,1]*xb^la[2,1],
  82    if dla[1]=dla[2] then w2:wd^sdla1 else w2:xa^sdla1*xb^sdmu1,
  83    w3:k1+k2*'integrate(w2,x), w4:w1*w3, print("y=",w4), ys:w4,
  84    reqans1(dum), if ans=y then
  85     (w5:k1+k2*integrate(w2,x), ys:expand(w1*w5),
  86      y1:coeff(ys,k1,1), y2:coeff(ys,k2,1)), return(ys))$
  87
  88 elmrephg2(m):=block([w1,w2,w3,w4,w5],
  89   sdlam :ratsimp(dla[1]+m),  sdmum :ratsimp(dla[2]+m),
  90   sdlam1:ratsimp(sdlam-1),sdmum1:ratsimp(sdmum-1),
  91   if la[1,2]=la[2,2] then w1:wd^la[1,2]
  92                      else w1:xa^la[1,2]*xb^la[2,2],
  93   if dla[1]=dla[2]   then (w2:wd^sdlam1, w3:wd^(-sdlam)) else
  94    (w2:xa^sdlam1*xb^sdmum1,w3:xa^(-sdlam)*xb^(-sdmum)),
  95   w4:k1+k2*'integrate(w3,x), w5:w1*'diff(w2*w4,x,m-1), print("y=",w5),
  96   ys:w5, reqans1(dum), if ans=y then
  97    (w6:k1+k2*'integrate(w3,x), w7:w1* diff(w2*w6,x,m-1), ys:expand(w7),
  98     y1:coeff(ys,k1,1), y2:coeff(ys,k2,1)), return(ys)  )$
  99
 100 lodesave([pnopm,fasl],nopmt,minabssum,excroots,elmrephg1,elmrephg2);
 101
 102 /* the sum has parameters */
 103 haspmts(dum):=block([], y0:1, plist:listofpars(eq),
 104   for i:1 thru 3 do checkcontratn(i),
 105   sortratn(dum), exccolumns(dum), pfuncthg(n,dum),
 106   if not(cdla[1]=cdla[2]) then hgstdd(dum),
 107   make01p(1), make01p(2), pfuncthg(y,y0), xa:x-spt[1],xb:x-spt[2],
 108   if la[1,2]=0 and la[2,2]=0 then (excrow(dum), pfuncthg(y,y0)),
 109   if la[1,1]=0 and la[2,1]=0 and dla[1]=dla[2] then (ans:n,
 110       if not integerp(dla[1]) then
 111          ans:read("Is",dla[1],"an integer?",str4),
 112       if integerp(dla[1]) or ans=y then
 113          (v:-la[3,1],ldf3(3),gsvlgdre(v),return(ys))),
 114   if dla[1]=dla[2] and la[1,3]=2 and la[2,3]=2 then
 115       (ldf3(4),ys:caseof22(dum)),
 116   if ys#f then return(ys) else (hgstdd(dum),
 117        ldf3(5),gsvhg(dum),return(ys)))$
 118
 119 pfdivide(y0,i,lamu):=block([],    wy:ratsimp(y0*(x-la[i,0])^lamu),
 120    la[i,1]:ratsimp(la[i,1]-lamu), la[i,2]:ratsimp(la[i,2]-lamu),
 121    la[3,1]:ratsimp(la[3,1]+lamu), la[3,2]:ratsimp(la[3,2]+lamu),
 122    v:-la[3,1], return(wy))$
 123
 124 listofpars(exp):=block([pl],pl:listofvars(exp), pl:delete(x,pl),
 125    pl:delete(y,pl), lplist:length(pl), return(pl))$
 126
 127 contratn(i):=block([], dexp:dla[i], j:1,
 128    wexp:ratsimp(dexp), if wexp=-1/2 then return(f),
 129    dnop:ratsimp(dexp), if integerp(dnop) then return(dnop),
 130  loop, dnop:ratsimp(dexp-1/j), if psdmint(dnop) then return(1/j),
 131    j:j+1, if j<6 then go(loop),
 132    dnop:ratsimp(dexp-2/5), if psdmint(dnop) then return(2/5),
 133    return(f))$
 134
 135 psdmint(dnop):=block([],wdno:ratsimp(dnop),
 136    if integerp(wdno) then (if wdno<=0 then return(t) else return(f)),
 137    if lplist=0 then return(f),
 138    if lplist=1 then
 139     (if denom(wdno)#1 then return(f),
 140      ans:read("Is",wdno,"a minus integer?",str4),
 141      if ans=y then return(t) else return(f)))$
 142
 143 sortratn(dum):=block([], nrat:0,
 144   for i:1 thru 3 do (if cdla[i]#f then
 145    (nrat:nrat+1, la[i,3]:denom(cdla[i]), if la[i,3]=1 then
 146    (if cdla[i]=0 then la[i,3]:0 else (minla[i]:min(abs(la[i,1]),abs(la[i,2])),
 147                                       la[i,3]:la[i,3]+minla[i])))
 148     else (if la[i,1]*la[i,2]=0 then la[i,3]:8 else la[i,3]:9)),
 149   clsort(dum))$
 150
 151 exccolumns(dum):=block([],
 152   if la[1,4]#1 then (if la[2,4]=1 then exc2c(1,2) else exc2c(1,3)),
 153   if la[2,4]=3 then exc2c(2,3),
 154   for i:1 thru 3 do dla[i]:ratsimp(la[i,1]-la[i,2]))$
 155
 156 exc2c(l,m):=block([],for i:0 thru 4 do exchla(i,l,m))$
 157
 158 exchla(i,l,m):=block([],ws:la[l,i],la[l,i]:la[m,i],la[m,i]:ws)$
 159
 160 wlength(exp):=block([], p1:first(plist),
 161   if not freeof(p1,exp) then (wdeg:hipow(exp,p1),return(wdeg))
 162   else return(0))$
 163
 164 mklntr(x1,x2,x3):=block([a,b,c,d],
 165  if x3=inf then (a:1/(x2-x1),b:-a*x1, c:0,d:1  ) else
 166  if x2=inf then (a:1,        b:-a*x1, c:1,d:-x3) else
 167  if x1=inf then (a:0,        b:x2-x3 ,c:1,d:-x3),
 168  vtr:t=ratsimp((a*x+b)/(c*x+d)) )$
 169
 170 checkcontratn(i):=block([wla],
 171    cdla[i]:contratn(i), if cdla[i]#f then return(t),
 172    exc1c(i), cdla[i]:contratn(i))$
 173
 174 hgstdd(dum):=block([],
 175    if la[1,0]=0 and la[2,0]=1 and la[3,0]=inf then return(f),
 176    mklntr(la[1,0],la[2,0],la[3,0]),
 177    print(str2,vtr), la[1,0]:0,la[2,0]:1,la[3,0]:inf,
 178    print(str1),pfuncthg(n,dum))$
 179
 180 make01p(i):=block([],
 181   if la[i,1]=la[i,2] and la[i,2]#0 then y0:pfdivide(y0,i,la[i,2]) else
 182      (l1:wlength(la[i,1]),l2:wlength(la[i,2]),
 183       if l1>l2 then excrow(dum),
 184       if la[i,1]#0 then y0:pfdivide(y0,i,la[i,2]),
 185       if la[i,2]=0 then exc1c(1)))$
 186
 187 clsort(dum):=block([], for i:1 thru 3 do la[i,4]:3,
 188    for i:1 thru 3 do (i1:nmod(i+1,3),i2:nmod(i+2,3),
 189     if la[i ,3]<=la[i1,3] and la[i,3]<=la[i2,3] then (la[i,4]:1,
 190     if la[i1,3]<=la[i2,3] then la[i1,4]:2 else la[i2,4]:2)))$
 191
 192 nmod(n,k):=block([],if n>k then return(n-k) else return(n))$
 193
 194 lodesave([phghp,fasl],haspmts,pfdivide,listofpars,contratn,
 195    checkcontratn,mklntr,psdmint,sortratn,exccolumns,
 196    exc2c,exchla,wlength,hgstdd,make01p,clsort,nmod);
 197
 198 gsvlgdre(v):=block([],remvalue(L),
 199 print("The solution is representable by the solution of Legendre's eq:
 200        (x^2-1)*y''+2*x*y'-v*(v+1)*y=0"),
 201   if la[1,2]=0 and la[2,2]=0 then svlgdre(v) else
 202  (if la[1,2]=la[2,2] then
 203      (if integerp(la[1,2]) then
 204         (if la[1,2]>0 then (y0:y0*expand(xa*xb),
 205          la[3,1]:la[3,1]+la[1,2], la[3,2]:la[3,2]+la[1,2], v:-la[3,2],
 206          yw:y[L](v,x), ys:y0^la[1,2]*'diff(yw,x,la[1,2]),
 207          print("y=",ys), return(ys)))
 208       else
 209         (ans:read("Is",la[1,2],str5),
 210          if ans=p then (w0:expand(xa*xb),y0:y0*w0^la[1,2] ,
 211             la[3,1]:la[3,1]+la[1,2], la[3,2]:la[3,2]+la[1,2],v:-la[3,2]) else
 212          if ans=m then (la[1,2]:-la[1,2], la[3,1]:la[3,1]-la[1,2],
 213                           la[3,2]:la[3,2]-la[1,2], v:-la[3,2]),
 214          remvalue(L), yw:y[L](v,x), ys:y0*'diff(yw,x,la[1,2]),
 215          print("y=",ys,",where y[L](v,x) is the solution of Legendre's eq."),
 216          return(ys)) ) ))$
 217
 218 svlgdre(v):=block([],
 219    if spt[1]*spt[2]=-1 and (spt[1]=1 or spt[1]=-1) then
 220       (ys:y0*y[L](v,x), print("y=",ys), return(ys))
 221  else (lfrtr:(-2*x+spt[1]+spt[2])/(spt[1]-spt[2]),lfrtr:ratsimp(lfrtr),
 222        vtr:t=lfrtr, ys:y0*y[L](v,t),
 223        print("y=",ys,"where t=",lfrtr),return(ys)))$
 224
 225 lodesave([psleg,fasl],svlgdre,gsvlgdre);
 226
 227 caseof22(dum):=block([], ans:no,
 228   mcdla:ratsimp(-dal[1]+1/2),
 229   if not integerp(mcdla) then
 230      ans:read("Is",mcdla,"a positive integer?",str4),
 231   if (integerp(mcdla) and mcdla>=0) or ans=y then (
 232     w1:sqrt(xa), w2:sqrt(xb),
 233     wy1:(w1+w2)^dla[3], wy2:(w1-w2)^dla[3],
 234     ys:y0*'diff(k1*wy1+k2*wy2,x,mcdla),print("y=",ys),return(ys))
 235   else return(f))$
 236
 237 lodesave([phyp25,fasl],caseof22);
 238
 239 gsvhg(dum):=block([E,K], ys:f, remvalue(E,K),
 240   if la[1,0]=0 and la[2,0]=1 and la[1,1]=0 and la[2,1]=0 then
 241    (print("The solution is representable by Hypergeometric function."),
 242     ha:la[3,1],hb:la[3,2],hg:-la[1,2]+1, remvalue(G), ys:y0*y[G](ha,hb,hg,x),
 243     if ha=-1/2 and hb=1/2 and hg=1 then
 244    (y1:E(sqrt(x)),y2:E(sqrt(1-x))-K(sqrt(1-x)),ys:k1*y1+k2*y2,print("y=",ys),
 245     print("where E and K are elliptic functions of 1st and 2nd kind."))
 246     else print("y=",ys), return(ys)))$
 247
 248 lodesave([psolh,fasl],gsvhg);