share/contrib/lindstedt.mac

   1 /*
   2    This is a first pass at a Lindstedt code.  It can solve problems
   3    with initial conditions entered, which can be arbitrary constants,
   4    (just not %k1 and %k2) where the initial conditions on the perturbation
   5    equations are z[i]=0,z'[i]=0 for i >0. Examples are:
   6    Lindstedt('diff(x,t,2)+x-(e*x^3)/6,e,2,[1,0]);
   7
   8    Problems occur when initial conditions are not given, as the constants
   9    in the perturbation equations are the same as the zero order equation
  10    solution.  Also, problems occur when the initial conditions for the
  11    perturbation equations are not z[i]=0,z'[i]=0 for i >0, such as the
  12    Van der Pol equation.
  13
  14    Remaining work:  Fix initial condition and limit cycle determination.
  15
  16
  17 Copyright (C) 2001  Dan Stanger
  18
  19 This library is free software; you can redistribute it and/or
  20 modify it under the terms of the GNU Lesser General Public
  21 License as published by the Free Software Foundation; either
  22 version 2.1 of the License, or (at your option) any later version.
  23
  24 This library is distributed in the hope that it will be useful,
  25 but WITHOUT ANY WARRANTY; without even the implied warranty of
  26 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  27 Lesser General Public License for more details.
  28
  29 You should have received a copy of the GNU Lesser General Public
  30 License along with this library; if not, write to the Free Software
  31 Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA
  32
  33 */
  34
  35 second(l):=first(rest(l));
  36
  37 /* This routine finds the dependent and independent variables in a
  38    equation by finding derivatives and looking at the variables in them. */
  39 LindstedtFindDiffVars(exp):=block([dvars:[],ivars:[],d:part('diff(x,y),0),p],
  40    scanmap(
  41       lambda([x],
  42          if not mapatom(x) then
  43             (if inpart(x,0) = d then
  44                block([dvar:inpart(x,1), ivar:inpart(x,2)],
  45                   if not member(dvar, dvars) then dvars:cons(dvar,dvars),
  46                   if not member(ivar, ivars) then ivars:cons(ivar,ivars))),
  47          x),
  48       exp),
  49    [dvars,ivars]);
  50
  51 /* This routine will process the equation eq as follows.
  52    First derivatives will be replaced by zp, second derivatives with z2p,
  53    and the dependent variable with z.  The independent variable will be
  54    replaced by T. */
  55 LindstedtProcessEquation(eq, dvar, ivar, z, zp, z2p, T):=block(
  56    [%zp%:?gensym(),%z2p:?gensym()],
  57    eq:subst(%zp%,'diff(dvar,ivar),eq),
  58    eq:subst(%z2p%, 'diff(dvar,ivar,2),eq),
  59    sublis([dvar=z,ivar=T,%zp%=zp,%z2p%=z2p],eq));
  60
  61 LindstedtProcessIC2(sol, ic, dvar, ivar):=
  62    ic2(sol,ivar=0,dvar=first(ic),'diff(dvar,ivar)=second(ic));
  63
  64 /* This routine finds and solves for secular terms.  For now,
  65    it only looks for powers of the ivar */
  66 LindstedtRemoveSecularTerms(sol,ivar):=block(
  67    /* get the high power of ivar in the solution and loop down */
  68    [h:hipow(sol,ivar), param],
  69    for i:h step -1 thru 1 do block([secularTerm,vars],
  70       /* find the power of ivar */
  71       secularTerm:coeff(sol,ivar,i),
  72       /* get the variables in the term */
  73       vars:listofvars(secularTerm),
  74       /* hopefully find a parameter to solve for */
  75       vars:delete(ivar,vars),
  76       /* print("vars",vars), */
  77       if vars = [] then return(sol),
  78       param:solve(secularTerm, vars),
  79       /* now param is the solution list, which is a array element.
  80          to assign to it I am using the following hack,  but i am
  81          sure there is a better way of doing this */
  82       apply(ev,subst(":","=",param))),
  83       /* This was the previous method I used to assign to the array
  84       map(lambda([eq],block([l:lhs(eq),r:rhs(eq),a,s],
  85                  a:part(l,0),s:part(l,1),
  86                  arraymake(a,[s])::r)),
  87          param)), */
  88    return(param)
  89 );
  90
  91 Lindstedt(eq,pvar, torder, [ic]):=
  92    block([dvars, ivar, l:LindstedtFindDiffVars(eq),solutionList],
  93       dvars:first(l), ivar:second(l),
  94       if length(ivar) # 1 then error("Number of independent variables not 1"),
  95       ivar:first(ivar),
  96       if length(dvars) = 0 then error("Number of dependent variables is 0"),
  97       if length(dvars) = 1 then
  98          /* since there is 1 dependent variable, name it z0, and form a
  99             solution in the form of sum(z[i](lamda*t)).  to avoid
 100             problems, lamda is used instead of lambda. */
 101          block([v,sol,T,dvar:first(dvars),
 102             lamda,z,zSum,lamdaSum,zPrimeSum,zPrimePrimeSum],
 103             local(lamda,z),
 104             array([lamda,z],torder+1),
 105             lamda[0]:1,
 106             lamdaSum:sum(lamda[i]*pvar^i,i,0,torder),
 107             zSum:sum(z[i]*pvar^i,i,0,torder),
 108             depends(z,T),
 109             zPrimeSum:sum(diff(z[i],T)*pvar^i,i,0,torder),
 110             zPrimePrimeSum:sum(diff(z[i],T,2)*pvar^i,i,0,torder),
 111             zPrimeSum:lamdaSum*zPrimeSum,
 112             zPrimePrimeSum:lamdaSum*lamdaSum*zPrimePrimeSum,
 113             eq:LindstedtProcessEquation( eq, dvar, ivar, zSum, zPrimeSum,
 114                zPrimePrimeSum, T),
 115             eq:expand(eq),
 116             /* extract the equation for z[0]. */
 117             v:coeff(eq,pvar,0),
 118             /* print("z[0] ode",v), */
 119             sol:ode2(v, z[0], T), /* solve the ode */
 120             /* if initial conditions are supplied, then apply ic2
 121                to the solution and initial condition list.  Then
 122                remove z[0] from the solution altogether.  A sanity
 123                check for the removal of z[0] could be done here. */
 124             if length(ic) > 0 then
 125                 z[0]:rhs(LindstedtProcessIC2(sol,first(ic), z[0],T))
 126                 else z[0]:rhs(sol),
 127             /* now remove z[0] from the equation and expand again */
 128             eq:ev(ratsubst(z[0], 'z[0], eq),nouns),
 129             eq:expand(trigreduce(eq)),
 130             /* start loop for higher order terms */
 131             /* print("Looping for higher order terms"), */
 132             for i:1 thru torder do block([param],
 133                /* determine the ode */
 134                v:coeff(eq, pvar, i),
 135                sol:ode2(v, z[i], T),
 136                sol:expand(rhs(ic2(sol, T=0, z[i]=0, 'diff(z[i], T) = 0))),
 137                param:LindstedtRemoveSecularTerms(sol,T),
 138                /* replace the following with lratsubst when its implemented */
 139                map(lambda([eq],block([l:lhs(eq),r:rhs(eq)],
 140                    sol:ratsubst(r,l,sol))),
 141                   param),
 142                z[i]:sol,
 143                /* eq:ev(ratsubst(z[i],'z[i],eq),NOUNS), */
 144                eq:ev(eq,nouns),
 145                eq:expand(trigreduce(eq))
 146             ),
 147             solutionList:[[[sum(z[i]*pvar^i,i,0,torder)],
 148                 T=ivar*sum(lamda[i]*pvar^i,i,0,torder)]]
 149          ));