This Maple worksheet accompanies the paper: 

   Paulo D. F. Gouveia, Delfim F. M. Torres, 

   "Computing ODE Symmetries as Abnormal Variational Symmetries", 

   Nonlinear Analysis (2008), doi:10.1016/j.na.2008.10.009. 

Computing ODE symmetries as abnormal variational symmetries 

Paulo D. F. Gouveia*
pgouveia@ipb.pt 

Delfim F. M. Torres**
delfim@ua.pt 

*Bragan?a Polytechnic Institute
5301-854 Bragan?a, Portugal 

**University of Aveiro
3810-193 Aveiro, Portugal 

Abstract 

We give a new computational method to obtain symmetries of ordinary differential equations. The proposed approach appears as an extension of a recent algorithm to compute variational symmetries of optimal control problems [P.D.F. Gouveia, D.F.M. Torres, Automatic computation of conservation laws in the calculus of variations and optimal control, Comput. Methods Appl. Math. 5 (4) (2005) 387-409], and is based on the resolution of a first order linear PDE that arises as a necessary and sufficient condition of invariance for abnormal optimal control problems. A computer algebra procedure is developed, which permits one to obtain ODE symmetries by the proposed method. Examples are given, and results compared with those obtained by previous available methods. 

Mathematics Subject Classification 2000: 34-04; 49-04; 34C14; 49K15. 

 

Keywords: Symmetries, variational symmetries, dynamic symmetries, ODEs, computer algebra systems, optimal control, abnormality. 

Introduction 

Sophus Lie was the first to introduce the use of symmetries into the study of differential equations, Emmy Noether the first to recognize the important role of symmetries in the calculus of variations. Nowadays, all the computer algebra systems which deal with differential equations provide several tools to help the user with the analysis of Lie symmetries. Recently, the authors developed a computer algebra package for the automatic computation of Noether variational symmetries in the calculus of variations and optimal control [5], now available as part of the Maple Application Center at http://www.maplesoft.com/applications/app_center_view.aspx?AID=1983. 

The omnipresent tools for Lie symmetries provide a great help for the search of solutions of ODEs, their classification, order reduction, proof of integrability, or in the construction of first integrals. From the mathematical point of view, a ODE symmetry is described by a group of transformations that keeps the ordinary differential equation invariant. Depending on the type of transformations one is considering, different symmetries are obtained. An important class of symmetries is obtained considering a one-parameter family of transformations which form a local Lie group. Those transformations are often represented by a set of functions known as the infinitesimal generators. From the practical point of view, the determination of the infinitesimal generators that define a symmetry for a given ODE is, in general, a complex task [6, 11]. To address the problem, we follow a different approach. 

We propose a new method for computing symmetries of ODEs by using a Noetherian perspective. Making use of our previous algorithm [5], that has shown up good results for the computation of Noether variational symmetries of problems of the calculus of variations and optimal control, we look to an ODE as being the control system of an optimal control problem. Then, we obtain symmetries for the ODE by computing the abnormal variational symmetries of the associated optimal control problem. 

The paper is organized as follows. In ?2 the necessary concepts associated with variational symmetries in optimal control are reviewed. The new method for computing symmetries of ODEs is explained in ?3. The method is illustrated in ?4, where we compute symmetries for three distinct ODEs and compare the results with the ones obtained by the standard procedures available in Maple. We end the paper with Section 5 of conclusions and final comments.The definitions of the new Maple procedure that implements our method are given in Appendix. 

The Maple package 

The procedure odeSymm, described in the paper and illustrated in the following section, together with some necessary technical routines, have been implemented for the computer algebra system Maple 11. 

odeSymm 

Computes the infinitesimal generators which define the symmetries of the ODE, or system of ODEs, specified in the input. As explained in the paper (section 3), this procedure involves the resolution of a system of partial differential equations. We have used the Maple solver pdsolve, using, as preferential method, the separation of the variables by sum. 

Output: 

• One or more lists of symmetry generators for a given ODE, or system of ODEs ([ =?, =?, =?, . . . , =?], ? ? ? ).

 

Syntax: 

• odeSymm(ode, x(t), opts)

 

Input: 

• ode - ordinary differential equation, or a set or list of ODEs;

 

• x(t) - any indeterminate function of one variable, or a list of them, representing the unknowns of the ODE problem;

 

• opts - (optional) specify options for the odeSymm command, where opts is one or more of the following:

 

• allconst - When this argument is given, the output presents all the constants given by the Maple command pdsolve. By default, that is, without option allconst, we eliminate redundant constants; this is done by our Maple procedure reduzConst, that is a technical routine. Essentially, the procedure transforms in one constant each sum of constants.

 

• mindep - When one wants to restrict to the minimum the dependencies of the infinitesimal generators: and . By default, that is, in the absence of options mindep and alldep, the following dependencies are considered: and ;

 

• alldep - All possible dependencies for the infinitesimal generators: and ;

 

• split - When this argument is given, the procedure invoke the split command to divide the resultant set of infinitesimal generators into uncoupled subsets, by fixing the values for all the constants given by the Maple command pdsolve. The procedure split is a technical routine.

 

• showdep - Shows, in the obtained solution, all the dependencies of the generators; otherwise only the name of the generators is shown;

 

• showt - Shows, in the obtained solution, the dependence on the time variable (independent variable); otherwise, the time variable is omitted as a function parameter;

 

• showgen - Shows, in the obtained solution, besides the infinitesimal generators and , the augmented set of variational generators, , and ;

 

• hint=<value> - Indicate a method of solution of the PDE system, where <value> is one of `+`, `*`, or any other expression allowed by the command pdsolve, being also possible to use hint=nohint for the case one wants to use the standard method of resolution of Maple; by default, the system is solved by separating the variables by sum (hint= `+`).

 

 

Definition: 

 

>	odeSymm := proc(ODEs::{`=`, set(`=`), list(`=`)}, depvars::{function,list(function)})

 

>	local n, tt, xx, pp, k, vX, vPSI, syseqd, sol, lstGerad, valGerad, phi, vphi, lpsi, vpsi,Hi, t, Sr, x0, r, aux, mapx, sys, xieta, sol2;

 

>	description "Symmetries in the ODE problems";

 

>	unprotect(Psi); unassign('T'); unassign('X'); unassign('Psi'); unassign('psi');

 

>	Hi:=subs(select(type,[args[3..-1]],`=`),hint); if Hi='hint' then Hi:=`+`; fi;

 

>	n:=nops(depvars);

 

>	if n=1 then x0:=[depvars] else x0:=depvars fi;

 

>	t:=op(1,x0[1]);

 

>	r:=[]:  for aux in x0 do    for k from 1 by 1 while evalb(subs(diff(aux,t$k)=z_z_z,ODEs)<>ODEs) do od;    r:=[r[], k-1]:  od:

 

>	Sr:=sum(r['i'],'i' =1..n);

 

>	mapx:=[seq(x0[i]=_x[1+(sum(r['k'],'k'=1 ..i-1))], i=1..n)];

 

>	mapx:=[mapx[],seq(seq(diff(x0[i],t$j)=_x[j+1+sum(r['k'],'k'=1..i-1)],j=1..r[i]-1),i=1..n)];

 

>	mapx:=[mapx[],seq(diff(x0[i],t$r[i])=_xx[i],i=1..n)];

 

>	mapx:=[seq(mapx[nops(mapx)+1-i],i=1..nops(mapx))];

 

>	sys:=subs(mapx,ODEs);

 

>	if n=1 then solve(sys,{_xx[1]}) else solve({sys[]},{seq(_xx[i],i=1..n)}) fi;

 

>	phi:=subs(%, [seq(_xx[i],i=1..n)]);

 

>	vphi:=Vector([seq([seq(_x[j],j=2+sum(r['k'],'k'=1..i-1)..sum(r['k'],'k'=1..i)),        phi[i]][],i=1..n)]);

 

>	x0:= [seq(_x[i], i = 1 .. Sr)];

 

>	if Sr>1 then lpsi:=[seq(psi[i],i=1..Sr)] else lpsi:=[psi] fi:

 

>	vpsi:=Vector[row](lpsi);

 

>	if member('alldep',[args[3..-1]]) then tt:=t,op(x0),op(lpsi); xx:=tt; pp:=tt;

 

>	elif member('mindep',[args[3..-1]]) then tt:=t; xx:=op(x0); pp:=op(lpsi);

 

>	else tt:=t; xx:=t,op(x0); pp:=t,op(lpsi); fi:

 

>	if Sr>1 then vX:=Vector([seq(X[i](xx), i=1..Sr)]);          else vX:=Vector([X(xx)]); fi;

 

>	if Sr>1 then vPSI:=Vector[row]([seq(PSI[i](pp), i=1..Sr)]);          else vPSI:=Vector[row]([PSI(pp)]); fi;

 

>	syseqd:={ vpsi.( map(diff,vphi,t)*T(tt)+Matrix([seq(map(diff,vphi,i),i=x0)]).vX    +vphi*diff(T(tt),t)-map(diff,vX,t) )+vPSI.vphi,

 

>	convert(-vPSI+(vpsi.vphi)*Vector[row]([seq(diff(T(tt),i),i=x0)])    -vpsi.Matrix([seq(map(diff,vX,i),i=x0)]), 'list')[],

 

>	convert((vpsi.vphi)*Vector[row]([seq(diff(T(tt),i),i=lpsi)])    -vpsi.Matrix([seq(map(diff,vX,i),i=lpsi)]), 'list')[]} minus {0}:

 

>	lstGerad:=[T(tt), convert(vX,'list')[], convert(vPSI,'list')[]];

 

>	if Hi='nohint' then sol:=pdsolve(syseqd, lstGerad);

 

>	else sol:=pdsolve(syseqd, lstGerad, HINT=Hi); fi;

 

>	if not member('allconst',[args[3..-1]]) then sol:=reduzConst(sol); fi:

 

>	valGerad:=subs(sol,lstGerad); sol:=[(lstGerad[i]=valGerad[i])$i=1..nops(lstGerad)];

 

>	sol:=collect(expand(simplify(sol)),[t,op(x0),op(lpsi)]);

 

>	if not member('showdep',[args[3..-1]]) then    xieta:=[xi,seq(eta[i],i=1..n)]; sol:=subs(map(i->i=op(0,i),lstGerad),sol);  else xieta:=[xi(tt),seq(eta[i](xx),i=1..n)]; fi;

 

>	if n=1 then xieta:=subs(eta[1]=eta,xieta) fi;

 

>	sol:=subs('PSI'='Psi', sol);

 

>	sol:=subs(map(i->rhs(i)=lhs(i),mapx),sol); xieta:=subs(map(i->rhs(i)=lhs(i),mapx),xieta);

 

>	sol2:=[xieta[1]=rhs(sol[1]),seq(xieta[i+1]=rhs(sol[2+sum(r['k'],'k'=1 ..i-1)]), i=1..n)];

 

>	if member('split',[args[3..-1]]) then sol2:=[split(sol2)] else sol2:=[sol2] fi;

 

>	if member('showgen',[args[3..-1]]) then sol:=[sol,sol2[]];

 

>	else sol:=sol2 fi;

 

>	if n=1 then x0:=op(0,depvars) else x0:=map(i->op(0,i),depvars)[] fi;

 

>	sol:=subs({map(i->i(t)=i,[x0,op(lpsi)])[]}, sol);

 

>	if member('showt',[args[3..-1]]) then sol:=subs({map(i->i=i(t),[x0,op(lpsi)])[]},sol) fi;

 

>	return sol[];

 

>

end proc:

 

Technical routines 

Essentially, the first three routines, reduzConst, levantamento and convertSums, are used to transform in one constant each sum of constants not repeated in a set of algebraic expressions. The constants in Maple notation are converted to a more usual mathematical notation. 

The others two routines, split and levantamentoCn, are used to divide a set of infinitesimal generatores into uncoupled subsets. 

 

reduzConst 

• Reduces the number of integration constants.

 

• Begins by passing to the format C_n (n=1, 2...), with the aid of function levantamento, all the constants in Maple notation (_Ci, with i=1, 2...) with more than an occurrence in the group of expressions specified in the input (cc), after which, through the function convertSums, each sum of constants _Ci  is turned into one constant Cn. Finally, all individual constants in the format _Ci  are also converted to the format Cn.

 

• Return the input set of expressions with the altered constants.

 

>	reduzConst:=proc(cc::set)

 

>	local L0, LL, LLr, aux, ss, sss, termo, indexConst, sol;

 

>	LL:={}; LLr:={}; sol:={};

 

>	for aux in cc do

 

>	LL, LLr:=levantamento(aux, LL, LLr);

 

>

od;

 

>	ss:=convert(LLr,'list');

 

>	sss:={seq(ss[i]=C[i],i=1..nops(ss))};

 

>	L0:=subs(sss,cc);

 

>	indexConst:=nops(ss);

 

>	for aux in L0 do

 

>	termo, indexConst := convertSums(aux,indexConst);

 

>	sol:={termo} union sol;

 

>

od:

 

>	return sol;

 

>

end proc:

 

 

levantamento 

• Identify repeated constants in an expression.

 

• Add to set conj all the constants _Ci present in the expression termo, and to the set conjr the constants with more than an occurrence in that expression.

 

• Return conj and conjr updated.

 

>	levantamento:=proc(termo, conj::set(symbol), conjr::set(symbol))

 

>	local aux, conj2, conjr2;

 

>	if nops(termo)=1 then

 

>	if type(termo,'symbol') and

 

>	StringTools[IsPrefix]("_C",termo) and

 

>	StringTools[IsDigit](substring(termo,3..-1)) then

 

>	if evalb(termo in conj) then return conj, {termo} union conjr;

 

>	else return {termo} union conj, conjr; fi:

 

>	else return conj, conjr;

 

>

fi:

 

>

else

 

>	conj2, conjr2 := conj, conjr;

 

>	for aux in op(termo) do

 

>	conj2, conjr2 :=levantamento(aux, conj2, conjr2);

 

>

od:

 

>	return conj2, conjr2;

 

>

fi:

 

>

end proc:

 

 

convertSums 

• Convert sums of constants in an individual constant.

 

• Represent each sum of constants of the type _Ci , present in the expression cc, by a constant C_n. The remaining constants  _Ci are also converted to the format C_n. The n is initialized with the value of indC, and increased whenever a new constant C_n is created.
  Return cc and indC updated.

 

>	convertSums:=proc(cc, indC::integer)

 

>	local tipo, soma, i, aux, auxcc, flag, indexConst;

 

>	indexConst:=indC;

 

>	tipo:=op(0,cc);

 

>	if type(cc, extended_numeric) then return cc, indexConst;

 

>	elif tipo='symbol' then return cc, indexConst;

 

>	elif tipo=`+` then

 

>	soma:=0; flag:=false;

 

>	for i from 1 to nops(cc) do

 

>	aux, indexConst := convertSums(op(i,cc),indexConst);

 

>	if type(aux,'symbol') and

 

>	StringTools[IsPrefix]("_C",aux) and

 

>	StringTools[IsDigit](substring(aux,3..-1)) then flag:=true;

 

>	else soma:=soma+aux; fi:

 

>

od:

 

>	if flag then indexConst:=indexConst+1; soma:=soma+C[indexConst]; fi:

 

>	return soma, indexConst;

 

>

else

 

>	auxcc:=cc;

 

>	for i from 1 to nops(cc) do

 

>	aux, indexConst := convertSums(op(i,cc),indexConst);

 

>	if type(aux,'symbol') and

 

>	StringTools[IsPrefix]("_C",aux) and

 

>	StringTools[IsDigit](substring(aux,3..-1)) then

 

>	indexConst:=indexConst+1; auxcc:=subsop(i=C[indexConst],auxcc);

 

>	else auxcc:=subsop(i=aux,auxcc); fi:

 

>

od:

 

>	return auxcc, indexConst;

 

>

fi:

 

>

end proc:

 

 

split 

• Divide an algebraic expession into uncoupled subexpressions, or a list of algebraic expressions into uncoupled sublists. Each uncoupled term is the expression of the input paramenter (al), where was fixed values for the constants in format (identified by the technical routine LevantamentoCn): one constant is fixed to the value one, and all the other constants are fixed to zero

 

>	split := proc(al::{algebraic,`=`,list(algebraic),list(`=`)})

 

>	local n, N, cc, LL, LLr, aux;

 

>	description "Divide a set of infinitesimal generatores into uncoupled subsets";

 

>	LL:={}; LLr:={};

 

>	for aux in al do

 

>	LL, LLr:=levantamentoCn(aux, LL, LLr);

 

>

od;

 

>	cc:={}; LL:=[op(LL)]; N:=nops(LL);

 

>	for n to N do

 

>	cc:=cc union {subs((LL[i]=0)$i=1..n-1,LL[n]=1,(LL[i]=0)$i=n+1..N, al)};

 

>

od;

 

>	return op(cc)

 

>

end proc:

 

 

levantamentoCn 

• Identify repeated constants in an expression.

 

• Add to set conj all the constants present in the expression termo, and to the set conjr the constants with more than an occurrence in that expression.

 

• Return conj and conjr updated.

 

>	levantamentoCn:=proc(termo, conj::set(indexed), conjr::set(indexed))

 

>	local aux, conj2, conjr2;

 

>	if nops(termo)=1 then

 

>	if type(termo,'indexed') and

 

>	StringTools[IsPrefix]("C[",convert(termo,string)) and

 

>	StringTools[IsSuffix]("]",convert(termo,string)) and

 

>	StringTools[IsDigit](substring(convert(termo,string),3..-2)) then

 

>	if evalb(termo in conj) then return conj, {termo} union conjr;

 

>	else return {termo} union conj, conjr; fi:

 

>	else return conj, conjr;

 

>

fi:

 

>

else

 

>	conj2, conjr2 := conj, conjr;

 

>	for aux in op(termo) do

 

>	conj2, conjr2 :=levantamentoCn(aux, conj2, conjr2);

 

>

od:

 

>	return conj2, conjr2;

 

>

fi:

 

>

end proc:

 

Auxiliary routine 

(Routine of our package of the Calculus of Variations [4], neccessary to run the Axample 3 of this worksheet) 

Construct the system of Euler-Lagrange equations of a higher-order problem of the Calculus of Variations with several dependent variables. 

Output: 

• A set or a vector with the Euler-Lagrange equations.

 

Syntax: 

• EulerLagrange(L, t, x, x1, x2, ..., xr)

 

Input: 

• L - expression of the Lagrangian;

 

• t  - name of the independent variable;

 

• x  - name, list of names or vector of names of the dependent variables;

 

• xi - (i=1, ..., r) name, list of names or vector of names for the ith derivatives of the dependent variables.

 

 

Definition: 

 

>	EulerLagrange := proc(L::algebraic, t::name, x0::{name,list(name),'Vector[column]'(name)}, x1::{name,list(name),'Vector[column]'(name)})

 

>	local xx, n, m, Lxi, xi, V, EL, i, j, k;

 

>	description "Euler-Lagrange Extremals";

 

>	if nargs<4 then print(`N?mero de argumentos insuficiente!`); return;

 

>	elif not type([args[3..-1]],{'list'(name),'listlist'(name),'list'('Vector[column]'(name))})

 

>	then print(`Erro na lista de nomes das var. depend. ou suas derivadas!`); return;

 

>

end if;

 

>	xx:=convert(x0, 'list')[];

 

>	n:=nops([xx]); m:=nargs-3;

 

>	xi:=[seq(Vector(convert(args[i], 'list')),i=3..m+3)];

 

>

V:=[0$n];

 

>	for i from 1 to m do

 

>	Lxi:=[seq(diff(L,k),k=convert(xi[i+1],'list'))]:

 

>	Lxi:=subs({map(k->k=k(t),[xx])[]},Lxi);

 

>	Lxi:=subs({seq(seq(xi[k+1][j]=diff(xi[1][j](t),t$k),j=1..n),k=1..m)},Lxi);

 

>	V:=V+(-1)^i*map(diff,Lxi,t$i);

 

>

end do:

 

>	EL:=[seq(diff(L,k),k=convert(xi[1],'list'))];

 

>	EL:=subs({map(k->k=k(t),[xx])[]},EL);

 

>	EL:=subs({seq(seq(xi[k+1][j]=diff(xi[1][j](t),t$k),j=1..n),k=1..m)},EL);

 

>

EL:=EL+V;

 

>	if type(x0,'Vector') then return convert(map(i->i=0,EL),'Vector[column]');

 

>	elif type(x0,'list') then return convert(map(i->i=0,EL),'set');

 

>	else return op(EL)=0; end if;

 

>

end proc:

 

Illustrative Examples 

In order to show the functionality and the usefulness of our new procedure odeSymm, we consider three concrete problems found in the literature. All the examples were carried out with Maple version 11 on a 1.4GHz 512MB RAM Pentium Centrino. The running time of procedure odeSymm is indicated, for each example, in seconds. 

Example 1  (Kamke?s ODE 120) 

We begin with a first order ODE found in Kamke?s book [7]: 

>	ode:= t*diff(y(t),t)-y(t)*(t*ln(t^2/y(t))+2)=0;

 

(4.1.1)

 

To obtain symmetries of the equation we use our Maple procedure odeSymm with the additional parameter hint=noint. This means that we will use the default method of resolution of PDEs of the Maple solver pdsolve. If the optional parameter hint is not used (see Examples 2 and 3 below), our procedure odeSymm uses the method of separation of variables. We obtain the following infinitesimal generators (0.72 s): 

>	gerad:= odeSymm(ode, y(t), split, hint=nohint);

 

(4.1.2)

 

One can test the validity of the obtained symmetries with the symtest command of the DEtools Maple package: 

>	map(DEtools[symtest], [gerad], ode, y(t));

 

(4.1.3)

 

The symtest confirm that the infinitesimal generators leave the given ODE invariant, i.e., the generators obtained by our method give indeed a symmetry to Kamke?s ODE 120. It is interesting to remark that, without the knowledge of the computed symmetries, the ODE Maple solver dsolve is not able to integrate the ODE: 

>	dsolve(ode, y(t), class);

 

However, when one gives to the Maple solver the infinitesimal generators found by our method, the ODE is correctly solved:  

>	dsolve(ode, y(t), HINT=[gerad]): simplify(%);

 

(4.1.4)

 

It is also interesting to note that our method is able to find one symmetry that is different from the ones obtained using the standard methods of the literature. The Maple system provides nine different algorithms to compute symmetries of ODEs through the command symgen of the DEtools package. All the available schemes for determining the infinitesimal generators ? option way=all ? are not able to identify our pair of infinitesimals : 

>	DEtools[symgen](ode, y(t), way=all);

 

(4.1.5)

 

Example 2 (Damped Harmonic Oscillator) 

We consider a harmonic oscillator with restoring force −kx, emersed in a liquid in such a way that the motion of the mass m is damped by a force proportional to its velocity. Using Newton?s second law one obtains, as the equation of motion, the following second order differential equation [9, pp. 432?434]: 

>	EL:= m*diff(x(t),t,t)+a*diff(x(t),t)+k*x(t)=0;

 

(4.2.1)

 

The symmetries for this equation are easily obtained with our Maple procedure odeSymm (1.21 s) 

>	gerad:= odeSymm(EL, x(t), split);

 

(4.2.2)

 

One can confirm that these infinitesimals represent valid symmetries for the differential equation: 

>	map(DEtools[symtest], [gerad], EL, x(t));

 

(4.2.3)

 

Note that the output of our odeSymm procedure includes a dynamical symmetry: the derivative of the dependent variable is present in the second pair of obtained infinitesimal generators. 

Example 3 (Kepler?s problem) 

We now consider the Kepler?s problem: a problem of the calculus of variations ? see [14, p. 217]. In this case the Lagrangian depends on two dependent variables and : 

, with . 

We will use the proposed method to determine symmetries for the corresponding Euler-Lagrange differential equation. The Euler-Lagrange equation is trivially obtained using our package of the calculus of variations [4, Example 5.2]: 

>	L:= m/2*(v[1]^2+v[2]^2)+K/sqrt(q[1]^2+q[2]^2);

 

(4.3.1)

 

>	EL:= EulerLagrange(L, t, [q[1],q[2]], [v[1],v[2]]);

 

(4.3.2)

 

In this case, the Euler-Lagrange equation is a system of two second order ODEs. Our odeSymm procedure is able to determine symmetries for systems of differential equations as well (13.32 s): 

>	odeSymm(EL, [q[1](t),q[2](t)], split);

 

(4.3.3)

 

It is worth to mention that this example can not be handled by the algorithms available in Maple. Indeed, the Maple command symgen that looks for a symmetry generator for a given ODE is not able to deal with more than one dependent variable. 

Conclusions 

We have used the CAS Maple to define a new computational procedure that determines, in an automatic way, symmetries of ODEs. The automatic calculation of symmetries is a subject much studied under the theory of differential equations, with many results and applications in many different areas. Our main novelty is the presentation of a new algorithm, alternative to existing ones, which looks to symmetries of ODEs as particular cases of Noether-variational symmetries. As explained in the paper (?3), our algorithm involves the resolution of a first order, homogeneous, and linear PDE, which is the abnormal case of the necessary and sufficient condition of invariance for problems of optimal control studied in connection with Noether?s theorem [5, 12]. Interesting points of the proposed method are: (i) it is based on a new approach to the subject ? in particular, it is different from all the nine alternative algorithms available in Maple; (ii) allows us to get dynamic symmetries for ODEs of any order; (iii) allows to determine symmetries for systems of ODEs, when the analog simgen Maple command of the DEtools package can only obtain solutions for a single ODE. 

References 

• E. S. Cheb-Terrab, L. G. S. Duarte and L. A. C. P. da Mota, Computer algebra solving of second order ODEs using symmetry methods, Comput. Phys. Comm. 108 (1998), no. 1, 90?114.

 

• E. S. Cheb-Terrab and K. von Blow, A computational approach for the analytical solving of partial differential equations, Comput. Phys. Comm. 90 (1995), no. 1, 102?116.

 

• D. S. Dukic, Noether?s theorem for optimum control systems, Internat. J. Control (1) 18 (1973), 667?672.

 

• P. D. F. Gouveia and D. F. M. Torres, Algebraic computation in the calculus of variations: determining symmetries and conservation laws, TEMA Tend. Mat. Apl. Comput. 6 (2005), no. 1, 81?90.

 

• P. D. F. Gouveia and D. F. M. Torres, Automatic computation of conservation laws in the calculus of variations and optimal control, Comput. Methods Appl. Math. 5 (2005), no. 4, 387?409.

 

• W. Hereman, Review of symbolic software for the computation of Lie symmetries of differential equations, Euromath Bull. 1 (1994), no. 2, 45?82.

 

• E. Kamke, Differentialgleichungen. Losungsmethoden und Losungen. Teil I: Gewohnliche Differentialgleichungen. 6. Aufl.; Teil II: Partielle Differentialgleichungen erster Ordnung fur eine gesuchte Funktion. 4. Aufl, Geest & Portig, Leipzig, 1959.

 

• P. K. Kythe, P. Puri and M. R. Schaferkotter, Partial differential equations and Mathematica, CRC, Boca Raton, FL, 1997.

 

• J. D. Logan, Applied mathematics, Wiley, New York, 1987.

 

• L. S. Pontryagin, V. G. Boltyanskii, R. V. Gamkrelidze and E. F. Mishchenko, The mathematical theory of optimal processes, Translated from the Russian by K. N. Trirogoff; edited by L.W. Neustadt, Interscience Publishers John Wiley & Sons, Inc. New York, 1962.

 

• A. Samokhin, Full symmetry algebra for ODEs and control systems, Acta Appl. Math. 72 (2002), no. 1-2, 87?99.

 

• D. F. M. Torres, Conservation laws in optimal control, in Dynamics, bifurcations, and control (Kloster Irsee, 2001), 287?296, Lecture Notes in Control and Inform. Sci., 273, Springer, Berlin, 2002.

 

• D. F. M. Torres, Weak conservation laws for minimizers which are not Pontryagin extremals, Proc. of the 2005 International Conference ?Physics and Control? (PhysCon 2005), August 24-26, 2005, Saint Petersburg, Russia. Edited by A.L. Fradkov and A.N. Churilov, 2005 IEEE, pp. 134?138.

 

• B. van Brunt, The calculus of variations, Springer, New York, 2004.

 

• D. Zwillinger, Handbook of differential equations, Academic Press, Boston, MA, 1989.

 

Legal Notice: The copyright for this application is owned by the author(s). Neither Maplesoft nor the author are responsible for any errors contained within and are not liable for any damages resulting from the use of this material. This application is intended for non-commercial, non-profit use only. Contact the author for permission if you wish to use this application in for-profit activities.