share/calculus/optvar.mac

   1 ttyoff: true;
   2 /* This file contains functions and option settings for
   3 variational optimization using the calculus of variations and
   4 the maximum principle.  For a description of its usage see the text
   5 file  OPTVAR USAGE. */
   6
   7 /* Set options to automatically print cpu time in milliseconds, force
   8 attempted equation solution even when there are more variables than
   9 unknowns, when an equation involves logs or exponentials, or when a
  10 coefficient matrix is singular: */
  11 time:grindswitch:solveradcan:singsolve:true$
  12 /* establish  D  as an alias for the differentiation function: */
  13 alias(d,diff)$
  14 /* The name of this function has changed */
  15 ic(soln,xa,ya,dya):=ic2(soln,xa,ya,dya)$
  16
  17 eval_when([translate,batch,demo,load,loadfile],
  18 dv(a)::=buildq([a],define_variable(a,'a,any)))$
  19
  20 dv(aux)$ dv(c)$ dv(dd)$ dv(dydt)$ dv(k)$
  21
  22 ham(odes) := block(
  23 /* This function computes the Hamiltonian & the auxiliary equations */
  24
  25    [t,nsv,statevars,auxvars,answ,elist,auxde], /*declare local vars*/
  26
  27    if not listp(odes) then odes: [odes], /* ensure list argument */
  28    t: part(odes,1,1,2),  /* get independent var from derivative */
  29    nsv: length(odes),  /* determine number of state variables */
  30    /* Form list of state and auxiliary variables: */
  31    statevars: auxvars: elist: [],
  32    for i thru nsv do (statevars: endcons(part(odes,i,1,1), statevars),
  33       auxvars: endcons(aux[i], auxvars)),
  34    answ: [sum(rhs(odes[i])*aux[i], i, 1, nsv)], /* form Hamiltonian */
  35
  36    /* Form list of auxiliary equations and any trivial solutions: */
  37    for i thru nsv do (
  38       auxde: 'diff(aux[i],t) = -diff(answ[1], statevars[i]),
  39       answ: endcons(auxde, answ),
  40       if rhs(auxde)=0 then answ:endcons(aux[i]=c[i],answ)),
  41
  42    /* Form list of E-labels while displaying computed results: */
  43    for item in answ do elist: endcons(first(apply('ldisp,[item])), elist),
  44    elist) $
  45
  46 el(f, ylist, tlist) := block( /* This function computes the
  47 Euler-Lagrange equations and any trivial first integrals: */
  48
  49    [ly,lt,fsub,energycon,fy,answ,elist],  /* declare local variables */
  50
  51    if not listp(tlist) then tlist: [tlist],
  52    if not listp(ylist) then ylist: [ylist],
  53    /* compute number of independent & independent variables: */
  54    ly: length(ylist), lt: length(tlist), fsub: f,
  55
  56    /* no conservation of energy if more than 1 independent var: */
  57    energycon: is(lt=1),
  58    for i thru ly do /* substitute for derivatives: */
  59       for j thru lt do (dd[i,j]: derivdegree(fsub,ylist[i],tlist[j]),
  60          if dd[i,j] > 1 then energycon: false,
  61          for kk thru dd[i,j] do
  62             fsub:subst('diff(ylist[i],tlist[j],kk)=dydt[i,j,kk], fsub)),
  63    /* no conservation of energy if independent var. in integrand: */
  64    if not freeof(tlist[1],fsub) then energycon: false,
  65    answ: if energycon then [fsub] else [], /* form list of results: */
  66    for i thru ly do (fy: diff(fsub,ylist[i]),
  67       answ: endcons(
  68          sum(sum((-1)^(kk-1)*'diff(diff(fsub,dydt[i,j,kk]),tlist[j],kk),
  69          kk,1,dd[i,j]), j, 1, lt) = fy, answ),
  70       if energycon then answ[1]: answ[1] -
  71          diff(fsub,dydt[i,1,1])*'diff(ylist[i],tlist[1]),
  72       if fy=0 and lt=1 and dd[i,1]=1 then /* momentum integral */
  73          answ: endcons(diff(fsub,dydt[i,1,1])=k[i], answ)),
  74    if energycon then answ[1]: answ[1]=k[0],
  75
  76    for i thru ly do  /* resubstitute original derivatives: */
  77       for j thru lt do
  78          for kk thru dd[i,j] do
  79             answ:subst(dydt[i,j,kk]='diff(ylist[i],tlist[j],kk), answ),
  80
  81    elist:[],  /* form list of E-labels while displaying results: */
  82    for eqn in answ do elist: endcons(first(apply('ldisp,[eqn])), elist),
  83    elist) $
  84
  85 convert(odes, ylist, t) := block([answ],
  86 /* This function converts output of EL or HAM to form required by
  87 DESOLVE  from the file  DESOLN. */
  88    if not listp(ylist) then ylist: [ylist],
  89    answ: apply('ev,[odes,'eval]),  /* if E-labels, replace with values */
  90    for yy in ylist do answ: subst(yy=funmake(yy,[t]), answ),
  91    return(answ)) $
  92 ttyoff: false $