share/diffequations/abel.mac

   1 /* -*- Macsyma -*- */
   2 abel(de,a,x):=block([v,coefs,f0,f1,f2,f3,a,%a,%fff,y,u,z,%ff,%dff,
   3                      %v,%uu,newde],
   4     depends(v,x),
   5     de:lhs(de)-rhs(de),
   6     coefs:makecoeflist3(de),
   7     f0:coefs[1],f1:coefs[2],f2:coefs[3],f3:coefs[4],a:coefs[5],
   8     if a#1 then return(abel2(coefs,y,x)),
   9     if f0=0 and f1=0 and freeof(x,%a:diff(%fff:f3/f2,x)/f2) then
  10       return(transform(de,y,x,u,z,[x=z,u=%fff*y])),
  11     %ff:-f2/(3*f3),%dff:diff(%ff,x),
  12     %v:exp(integrate((f1+%ff*f2),x)),
  13     %uu:(f0-%dff+f1*%ff+2*f2*%ff^2/3)/%v^3/f3,
  14     newde:transform(de,y,x,u,z,[y=u*%v+%ff,z=ratsimp(integrate(f3*%v^2,x))]),
  15      return(ratsimp(newde)))$
  16
  17
  18 makecoeflist3(de):=block([a,f0,f1,f2,f3,y,x],
  19         de:expand(de),
  20         a: coeff(de,'diff(y,x),1),
  21
  22         f0:-coeff(de,y,0),
  23         f1:-coeff(de,y,1),
  24         f2:-coeff(de,y,2),
  25         f3:-coeff(de,y,3),
  26 return([f0,f1,f2,f3,a]))$
  27
  28 /*method for non-linear y'*/
  29 nonlin(de,y,x):=block([temp,temp1,newy,dd,newde,dispflag,%p],
  30         de:expand(de),depends(newy,x),dispflag:false,
  31         dd:derivdegree(de,y,x),
  32         if(dd>1) then return(false) else
  33          if(hipow(de,diff(y,x))<2) then return(false),
  34         newde:subst(%p,diff(y,x),de),
  35         depends(%p,x),
  36         newde:lhs(newde)-rhs(newde),
  37         temp:solve(newde=0,%p),
  38         temp:subst('diff(y,x),%p,apply('ev,[temp])),
  39         temp:map(lambda([v],ode2(v,y,x)),apply('ev,[temp,'diff])),
  40         temp:map(lambda([v],solve(v,y)),temp),
  41         return(apply('ev,[temp,'infeval])))$
  42
  43
  44
  45 berno(de,y,x):=block([nn,mm,nni,mmi,%q,ans],ans:[],
  46           if(derivdegree(de,y,x)>1) then return(false),
  47           de:lhs(de)-rhs(de),
  48           nn:expand(coeff(de,diff(y,x))),
  49           mm:expand(radcan(de-nn*diff(y,x))),
  50           nni:apply('append,maplist('elements,nn)),
  51           mmi:apply('append,maplist('elements,mm)),
  52           if(atom(mm) or part(mm,0)="^") then mm:[mm],
  53 /* commented out of DOE MACSYMA because it can't be translated
  54           MAP(LAMBDA([V],IF NOT FREEOF(%Q,RATSUBST(%Q,
  55              DELETE(NUMFACTOR(DIFF(V,X)),DELETE(DIFF(Y,X),DIFF(V,X))),NN))
  56              THEN ANS:CONS(V,ANS)),MM), */
  57 /* and replaced by the equivalent do loop */
  58           for v in mm do (if not freeof(%q,ratsubst(%q,
  59              delete(numfactor(diff(v,x)),delete(diff(y,x),diff(v,x))),nn))
  60              then ans:cons(v,ans)),
  61           return(ans))$
  62
  63 elements(a):=block([],if atom(a) then return([a]),
  64           if(part(a,0) ="*") then return(maplprod(lambda([v],v),a)),
  65           if(part(a,0)="+") then return(maplsum(lambda([v],v),a)),
  66           return([a]))$
  67
  68 /*Method for nonlinear coeffs*/
  69 nonlin1(de,y,x):=block([%v,a1,a2,newde],
  70           de:expand(rhs(de)-lhs(de)),depends(%v,x),
  71           a1:coeff(de,diff(y,x)),
  72           a2:expand(de-a1*diff(y,x)),
  73           if(hipow(a1,x)=hipow(a2,y) and hipow(a2,x)=hipow(a1,y))
  74             then (newde:apply('ev,[a1,y=%v,x=1])*x*diff(%v,x)+
  75                        (apply('ev,[a2,y=%v,x=1])+apply('ev,[a1,y=%v,x=1])*%v),
  76                        a1:ode2(newde=0,%v,x),
  77                        if freeof(nounify('integrate),a1) then
  78                        return(subst([%v=y/x],a1))  else return(false))
  79             else(return(false)))$