share/diffequations/pmain.mac

   1 /* ************** beginning of main ******************************** */
   2
   3 lode2(oeq):=block([lhs,ay,ypx],tobegin(dum),
   4    remvalue(y1,y2), swro:0, ys:f, nhstddeq(oeq), if rr#0 then
   5     (print("first",str1,eq)), kt:0, ys:hlode2(eq),
   6    if rr=0 then return(ys), print("the solution of the homog. eq. is",ys),
   7    ldf0(8), ypx:nonhom(rr,y1,y2), return(ys+ypx))$
   8
   9 tobegin(dum):=block([],dyx2:'diff(y,x,2),dyx1:'diff(y,x),
  10    dyt2:'diff(y,t,2),dyt1:'diff(y,t),dux2:'diff(u,x,2),dux1:'diff(u,x),
  11    str1:"we solve",str2:"we use",str3:"the result is",str4:"type y or n",
  12    str5:"positive,zero,or negative? type p,z,or n")$
  13
  14 nhstddeq(oeq):=block([], weq:expand(oeq), lhs:part(weq,1),
  15    rhs:part(weq,2), lhs:lhs-rhs, deqcoeff(lhs), wlhs:lhs,
  16    lhs:za*dyx2+zb*dyx1+zc*y, rhs:-ratsimp(wlhs-lhs),
  17    weq:lhs=rhs, steq:stddeq(weq), lhs:part(steq,1), eq:lhs=0,
  18    if lhs#wlhs then print(str1,steq))$
  19
  20 hlode2(eq):=block([exp,ro], aeq[kt]:eq, swsv:f,
  21  ll1, leq:part(eq,1), tp:brchde1(leq),
  22    if zc=0   then (ldf0( 9),svode1(rb),return(ys)),
  23    if tp= 50 then (ldf0(10),skelp(rb,rc), if ys#f then return(ys)),
  24    if swpt=y then (ldf0(11),gptm(rb,rc),  if ys#f then return(ys)
  25       else if rptm=t then (if freeof(y,vtr) then idvprep(1)
  26                            else idvprep(2),go(ll1))),
  27    if tp>60  then (lv:listofvars(za), if length(lv)>1 then
  28       (ldf0(16), if stddsp(za)=t then go(ll2)), tblspt(za,rb,rc,x,n),
  29       if nspt>3  then (ldf0(12),if grmaps(dum)=y then go(ll1))),
  30    ha:hipow(za*q1+zb*q2+zc*q3,x), netv:listofvars(ha),lg:length(netv),
  31       if lg>0 then (ldf0(13),if ppowr(rb,rc,lg) then tp:70),
  32    if tp<=70 then go(ll2), vecdeg(za,zb,zc,x), tp:brchde2(rb,rc),
  33    if tp=110 then (ldf0(2),svconf(rb,rc),return(ys)) else
  34    if tp=120 then (ldf0(3),svhypr(rb,rc),return(ys)) else
  35    if tp=125 then (idvprep(2),go(ll1)),
  36    if tp=126 then (idvprep(1),go(ll1)),
  37    if tp=130 then return(ys),   if swsv=tdv or tp=990 then go(ll3),
  38   ll2,ldf0(5),neq:ivtrns(rb,rc,rr),if 10<tp and tp<14 then go(ll3),
  39    if neq#f  then (idvprep(1),go(ll1)),   if tp=0 then go(ll4),
  40   ll3, ldf0(6),neq:dvtrns(za,rb,rc,rr,ro,lam1),
  41    if neq#f  then (idvprep(2),go(ll1)),
  42   ll4, nosol(oeq),return(oeq))$
  43
  44 deqcoeff(exp):=block([],za:coeff(exp,dyx2,1),zb:coeff(exp,dyx1,1),
  45    zc:coeff(exp,y,1),rb:ratsimp(zb/za),rc:ratsimp(zc/za))$
  46
  47 brchde1(exp):=block([],
  48    if not freeof(sin,cos,tan,cot,exp)     then return(10),
  49    if not freeof(%e,exp)                  then return(20),
  50    if not freeof(log,exp)                 then return(30),
  51    if not freeof(sinh,cosh,tanh,coth,exp) then return(40),
  52    if not freeof(sn,cn,dn,pe,th,exp)      then return(50),
  53    if not freeof(sqrt(x),exp)             then return(60),
  54    return(990))$
  55
  56 brchde2(wb,wc):=block([ws,wi],
  57    if nspt=1 then (if rk[1]=1 then return(110)),
  58    if nspt=2 then (if rk[1]<=0 and denom(rk[2])<3 and rk[2]<=1 then
  59      (if spt[1]#0 then (axb:x-spt[1],return(91)) else return(110))),
  60    if nspt=3 and rk[1]=0 and rk[2]=0 then (if rk[3]=0 then
  61      (if nf=1 and maxacf=2 then (ldf0(14),if powr(wb,wc,2) then return(80)),
  62       if da=2 and db=1 and dc=0 then
  63          (ldf0(15),tp:eqpt1(dum), if tp#f then return(tp)),
  64       if da=3 and db=2 and dc=1 then
  65          (ldf0(18),tp:eqpt12(wb,wc),if tp#f then return(tp))),
  66     if da=4 and db=3 and dc<=4 then
  67          (ldf0(17),tp:eqpt11(dum),if tp#f then return(tp)) else
  68     if rk[3]=0 then return(120)),
  69    ldf0(7),tp:brchde3(wb,wc),return(tp))$
  70
  71 tblspt(za,wb,wc,x,ksw):=block([k,n,wa,wf,wfi], k:0, n:hipow(za,x),
  72    if n=0 then go(lf), pspt:1, maxacf:0, wf:facza(za), j:0,
  73    for i:1 thru nf do (pwf[i]:part(wf,i),if pwf[i]=x then j:i),
  74    if nf>1 and j>1 then (ws:pwf[1],pwf[1]:pwf[j],pwf[j]:ws),
  75    for i:1 thru nf do (wfi:solve(pwf[i],x),l1:length(wfi),
  76      lacf[i]:l1, l2:polord(za,pwf[i]), maxacf:max(maxacf,l2),
  77      for j:1 thru l1 do (k:k+1,spt[k]:part(wfi[j],2),
  78       if ksw#n then kindspt(wb,wc,k),
  79       if l1=1 then pspt:pspt*spt[k] else pspt:0)),
  80  lf,k:k+1,nspt:k,spt[k]:inf,if ksw#n then (kindspt(wb,wc,k),ro:rk[k]))$
  81
  82 facza(za):=block([w1,w2,w3,w4], w1:za,w2:diff(w1,x),w3:gcd(w1,w2),
  83    w4:ratsimp(w1/w3), wf:factor(w4), if wf=w4 then nf:1 else nf:length(wf),
  84    if nf=1 then wf:[wf], return(wf))$
  85
  86 vecdeg(za,zb,zc,x):=block([],da:hipow(za,x),db:hipow(zb,x),dc:hipow(zc,x))$
  87
  88 nosol(oeq):=block([], print("we cannot solve"), print(oeq))$
  89
  90 idvprep(i):=block([],print(str2,vtr),print(str3,neq),
  91    if i=1 then eq:subst(x,t,neq) else eq:subst(y,u,neq), vprep(dum))$
  92
  93 vprep(dum):=block([],eq:stddeq(eq),print(str1,eq),kt:kt+1,aeq[kt]:eq,
  94    avtr[kt]:vtr)$
  95
  96 stddeq(oeq):=block([wexp,wr,deq],wexp:part(oeq,1),
  97    wr:part(oeq,2),deqcoeff(wexp),rr:ratsimp(wr/za),
  98    za:calca(rb,rc,x), zb:ratsimp(rb*za), zc:ratsimp(rc*za),
  99    deq:dyx2+rb*dyx1+rc*y=rr, return(deq))$
 100
 101 calca(wb,wc,x):=block([],w1:denom(wb),w2:denom(wc),za:lcm(w1,w2,x))$
 102
 103 lcm(p,q,x):=block([r,s],r:gcd(p,q,x),s:ratsimp(p*q/r),return(s))$
 104
 105 polord(exp,sexp):=block([we,n],we:exp,n:0,ll,we:ratsimp(we/sexp),
 106    if denom(we)#1 then return(n), n:n+1,go(ll))$
 107
 108 kindspt(wb,wc,k):=block([bt,ct,b1,c1,dbt,dct],
 109    if spt[k]=inf then alpha:0 else alpha:spt[k],
 110    bt:ratsubst(1/t+alpha,x,wb), b1:ratsimp((2*t-bt)/t^2),
 111    ct:ratsubst(1/t+alpha,x,wc), c1:ratsimp(ct/t^4),if spt[k]=inf then
 112     (bt:subst(x,t,b1),ct:subst(x,t,c1),dbt:ratdeg(wb,x),dct:ratdeg(wc,x))
 113  else (bt:ratsimp(wb),ct:ratsimp(wc),  dbt:ratdeg(b1,t),dct:ratdeg(c1,t)),
 114    rk[k]:max(dbt,dct/2)+1, if rk[k]=0 and freeof(%i,alpha)
 115    then regsingpt(x-alpha,bt,ct,alpha))$
 116
 117 ratdeg(rf,x):=block([wrf,nrf,dgnrf,drf,dgdrf,dgrf],if rf=0 then
 118    return(-99), wrf:ratsimp(rf), nrf:num(wrf), dgnrf:hipow(nrf,x),
 119    drf:denom(wrf),dgdrf:hipow(drf,x), dgrf:ratsimp(dgnrf-dgdrf),
 120    return(dgrf))$
 121
 122 regsingpt(xa,wb,wc,alpha):=block([f1,g1,flx], remvalue(l),
 123    f1:ratsimp(xa*wb), g1:ratsimp(xa^2*wc), flx:l^2+(f1-1)*l+g1,
 124    cheq[k]:subst(alpha,x,flx),cflx[k]:flx)$
 125
 126 reqans1(dum):=block([],ans:read("continue?",str4))$
 127
 128 ldf0(i):=block([],if ldsw[0,i]#y then (ldsw[0,i]:y,
 129    if i= 2 then load(pconfl) else
 130    if i= 3 then load(phypgm) else
 131    if i= 5 then load(pivtr) else
 132    if i= 6 then load(pdvtr) else
 133    if i= 7 then load(pbrnch) else
 134    if i= 8 then load(pnonh) else
 135    if i= 9 then load(psode1) else
 136    if i=10 then load(pskelp) else
 137    if i=11 then load(pgptm) else
 138    if i=12 then load(premap) else
 139    if i=13 then load(ppow) else
 140    if i=14 then load(pnpow) else
 141    if i=15 then load(peqp1) else
 142    if i=16 then load(pstdsp) else
 143    if i=17 then load(peqp11) else
 144    if i=18 then load(peqp12) ))$
 145
 146 lodesave(filename,[lode_functions]):=
 147   if status(feature,its)=true
 148   then apply('fassave,cons(filename,lode_functions))
 149   else block([packagefile:true],
 150              if status(feature,unix)=true
 151              then apply('save,
 152                         cons(concat(first(filename),".l"),lode_functions))
 153              else error("unknown system"))$
 154
 155 lodesave([pmain,fasl],lode2,tobegin,nhstddeq,hlode2,brchde1,vecdeg,
 156    brchde2,deqcoeff,tblspt,facza,nosol,idvprep,vprep,stddeq,calca,
 157    lcm,polord,kindspt,ratdeg,regsingpt,reqans1,ldf0);
 158
 159 lodesave([lodsav,fasl],lodesave)$