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