1ttyoff: true;
2/* This file contains functions and option settings for
3variational optimization using the calculus of variations and
4the maximum principle.  For a description of its usage see the text
5file  OPTVAR USAGE. */
6
7/* Set options to automatically print cpu time in milliseconds, force
8attempted equation solution even when there are more variables than
9unknowns, when an equation involves logs or exponentials, or when a
10coefficient matrix is singular: */
11time:grindswitch:solveradcan:singsolve:true$
12/* establish  D  as an alias for the differentiation function: */
13alias(d,diff)$
14/* The name of this function has changed */
15ic(soln,xa,ya,dya):=ic2(soln,xa,ya,dya)$
16
17eval_when([translate,batch,demo,load,loadfile],
18dv(a)::=buildq([a],define_variable(a,'a,any)))$
19
20dv(aux)$ dv(c)$ dv(dd)$ dv(dydt)$ dv(k)$
21
22ham(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
46el(f, ylist, tlist) := block( /* This function computes the
47Euler-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
85convert(odes, ylist, t) := block([answ],
86/* This function converts output of EL or HAM to form required by
87DESOLVE  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)) $
92ttyoff: false $
93