share/diffequations/odeaux.mac

   1 /* -*- Mode: macsyma; Package: CL-MAXIMA -*- */
   2 /* A collection of all the best ODE routines, run in the*/
   3 /*  best (hopefully) order  */
   4 /*************************************************/
   5
   6 /* This intitializes the environment when the package gets loaded
   7    or batched */
   8
   9 define_variable(allmethods2,
  10                 '[ode2,whittaker,solvehyper,solfac,diffsol,desol,series],
  11                 list)$
  12 define_variable(allmethods1,
  13                 '[ode2,diffsol,nonlin,nonlin1,riccati],
  14                 list)$
  15
  16 define_variable(odetutor,false,boolean)$
  17
  18 define_variable(method,'method,any)$
  19
  20 define_variable(method2,'method2,any)$
  21
  22 define_variable(sings,'sings,any)$
  23
  24 define_variable(ode_routinelist,[],list)$
  25
  26 eval_when([batch,loadfile,translate],block(
  27 /* for DOE MACSYMA, the autoload properties are in max$disk:[ode]ode.lsp */
  28 if not status(feature,nil) or status(feature,cl)  then (
  29  ode_routinelist:'[[[abel,lisp,dsk,ode],nonlin,nonlin1],
  30                    [[ricsol,lisp,dsk,ode],ricsol],
  31                    [[diffac,lisp,dsk,ode],solfac],
  32                    [[lapl,lisp,dsk,ode],diffsol],
  33                    [[trans,lisp,dsk,ode],desol, deptran,indtran,
  34                         invariant,normalform,schwartzian,bothvartran,dadjoint],
  35                    [[ricati,lisp,dsk,ode],apollo,riccati],
  36                    [[hyper,lisp,dsk,ode],solvehyper],
  37                    [[manyt,lisp,dsk,ode],manrecur],
  38                    [[series,lisp,dsk,ode],series],
  39                    [[whit,lisp,dsk,ode],whittaker],
  40                    [[schmid,lisp,dsk,ode],schmidt],
  41                    [[intfac,lisp,dsk,ode],int_factor],
  42                    [[sings,lisp,dsk,ode],getsings],
  43                    [[tran,lisp,dsk,ode],transform],
  44                    [[hyp,fasl,dsk,ode],hgfred],
  45                    [[remote,lisp,dsk,ell],run_remote],
  46                    [[kamke,lisp,dsk,ode],kamke]],
  47  for i in ode_routinelist do apply('setup_autoload,i)),
  48  sstatus(feature,"ode")))$
  49
  50
  51 eval_when([batch,loadfile,translate],block([],
  52 /*  for DOE MACSYMA, the autoload properties are in max$disk:[ode]ode.lsp */
  53 if   status(feature,cl)  then (
  54  ode_routinelist:'[["maxima-source:maxima;abel",nonlin,nonlin1],
  55                    ["maxima-source:maxima;ricsol",ricsol],
  56                    ["maxima-source:maxima;diffac",solfac],
  57                    ["maxima-source:maxima;lapl",diffsol],
  58                    ["maxima-source:maxima;trans",desol, deptran,indtran,
  59                         invariant,normalform,schwartzian,bothvartran,dadjoint],
  60                    ["maxima-source:maxima;ricati",apollo,riccati],
  61                    ["maxima-source:maxima;hyper",solvehyper],
  62                    ["maxima-source:maxima;manyt",manrecur],
  63                    ["maxima-source:maxima;series",series],
  64                    ["maxima-source:maxima;whit",whittaker],
  65                    ["maxima-source:maxima;schmid",schmidt],
  66                    ["maxima-source:maxima;intfac",int_factor],
  67                    ["maxima-source:maxima;sings",getsings],
  68                    ["maxima-source:maxima;tran",transform],
  69                    ["maxima-source:maxima;hyp",hgfred],
  70                    ["maxima-source:maxima;remote",run_remote],
  71                    ["maxima-source:maxima;kamke",kamke]],
  72  for i in ode_routinelist do apply('setup_autoload,i)),
  73  sstatus(feature,"ode")))$
  74
  75 /* This intitializes the environment for TRANSLATION,
  76    loads modedeclarations, for the whole package, macro packages, etc. */
  77
  78 /* EVAL_WHEN([TRANSLATE],?PRINC('?"
  79 WE ARE TRANSLATING AWAY!
  80 ",?TYO))$ */
  81
  82
  83 odeaux(eq,oldy,oldx,option):=block([dd,ans,solutionlist,x,y],
  84         eq:subsf(eq,oldy,oldx),
  85     return(subst(oldy,y,subst(oldx,x,solveode(eq,y,x,option)))))$
  86
  87 solveode(eq,y,x,option):=block([dd,ans,solutionlist],
  88     eq:eqsimp(eq),
  89     dd:derivdegree(eq,y,x),
  90
  91     if option=[] or option=['any] then (
  92     if(dd=1) then (
  93     if(ans:win('ode2,eq,y,x))#false then return(ans),
  94     if(testli(rhs(eq)-lhs(eq),y,x))=true then
  95     (if(ans:win('diffsol,eq,y,x))#false then
  96           (method:"LAPLACE XFORM", return(ans))) else
  97     (if(ans:win('nonlin,eq,y,x))#false then
  98           (method:"NONLIN IN YPRIME", return (ans)),
  99      if(ans:win('nonlin1,eq,y,x))#false then
 100           (method:"NONLIN XFORM", return (ans)),
 101     if riccati_test(lhs(eq)-rhs(eq),y,x)=true then
 102     if(ans:win('riccati,eq,y,x))#false then
 103           (method:"RICCATI",return(ans)),
 104     if(ans:win('int_factor,eq,y,x))#false then
 105           (method:"FOUND EULER MULTIPLIER",return(ans))
 106           )) else
 107           ( /*These methods are tailored for order = 2  */
 108     if dd=2 then(
 109     if(ans:win('ode2,eq,y,x))#false then return(ans),
 110     if(ans:win('desol,eq,y,x))#false then (method:'method1,return(ans)),
 111     if(ans:win('diffsol,eq,y,x))#false then
 112        (method:"LAPLACE XFORM",return(ans)),
 113
 114
 115   if(length(sings:getsings(eq,y,x))=3) then
 116
 117     if(ans:win('solvehyper,eq,y,x))#false then(method:'HYPERGEOMETRIC,
 118       return(ans)),
 119   if(sings='[0,inf] ) then
 120     if(ans:win('whittaker,eq,y,x))#false then(method:"CONFLUENT HYPER",
 121       return(ans)),
 122     if(ans:win('solfac,eq,y,x))#false then
 123           (method:"FACTOR OPERATOR",return(ans)),
 124     if(ans:win('series,eq,y,x))#false then (method:"SERIES",
 125       return(ans))) else
 126    if (ans:win('diffsol,eq,y,x))#false then (method:"LAPLACE XFORM",
 127        return(ans)),
 128   method:"NONE",return(false)))
 129   else(
 130     if option=['all] then
 131       if dd=2 then option:allmethods2 else if dd=1 then option:allmethods1
 132          else option:['diffsol],
 133     solutionlist:map(lambda([u],win(u,eq,y,x)),option),
 134     if (length(solutionlist)=1) then
 135       if programmode#true then
 136       solutionlist:part(solutionlist,1),
 137     return(solutionlist)))$
 138
 139 win(fun,eq,y,x):= block([ans],
 140     if odetutor=true then print( "TRYING    ",fun),
 141     if((ans:errcatch(block(method2:fun,apply(method2,[eq,y,x])))=[]))then 'fun,
 142 /* commented out of DOE MACSYMA since the GC statistics are on by default
 143        if ?cadr(status(freecore))=0 then error("out of core...sorry"), */
 144     if  (ans=[false] or ans=[]) then return(false)
 145             else return(part(ans,1)))$
 146
 147 homog(eq,y,x):=block([eq1,p1,q1,hom,h1],
 148     eq1:lhs(eq)-rhs(eq),
 149     p1:p(eq1,y,x),q1:q(eq1,y,x),
 150     hom:diff(y,x,2)+p1*diff(y,x)+q1*y,
 151     if((h1:radcan(eq1-hom))=0) then return(false) else
 152        return(hom=-h1))$
 153
 154 purify(sol,y,x):=block([sol1,s1,s2,%k1,%k2],
 155     sol1:rhs(sol),
 156     s1:ratcoef(sol1,%k1),s2:ratcoef(sol1,%k2),
 157     if(s1=0) then if(s2=0) then return(false) else s1:s2,
 158     if(length(s1)>length(s2)) then s1:s2,
 159     return(s1))$
 160
 161 testli(de,y,x):=block([h],
 162     de:expand(de),
 163     h:hipow(de,diff(y,x)),
 164     if(h#1)then return(false),h:[],
 165     h:cons(ratcoef(de,diff(y,x)),h),
 166     h:cons(ratcoef(de,y,1),h),
 167     h:cons(ratcoef(de,y,0),h),
 168 if(hipow(de,y)#1) then return(false),
 169     h:cons(ratcoef(de,diff(y,x,2)),h),
 170       if not freeof(y,h) then return(false) else
 171       return(true))$
 172
 173
 174
 175 /*Simplifies an ODE removes factors common to all terms*/
 176
 177 eqsimp(de):=block([inflag:true],
 178     de:factor(lhs(de)-rhs(de)),
 179     if(atom(de) or not(inpart(de,0)="*")) then return(de=0),
 180     de:map(lambda([v],if freeof(nounify('diff),v) then 1 else v),de),
 181     return(de=0))$
 182 subsf(eq,oldy,oldx):=block([%%newy,y,x],
 183     depends(y,x),
 184     eq:subst(%%newy,oldy,eq),
 185     eq:subst(x,oldx,eq),
 186     eq:subst(y,%%newy,eq),
 187     return(eq))$
 188
 189 riccati_test(de,y,x):=block([q],de:expand(de),
 190     if hipow(de,'diff(y,x))#1 then return(false),
 191     if ((q:lopow(de,'diff(y,x)))#1 and  q#0) then return(false),
 192     if hipow(de,y)=2 then return(true) else return(false))$