BeginPackage["VariableGradient`"]
VariableGradient::usage = "VariableGradient[A] determines the time derivative
of the Liapounov Function for an nth order system described by the matrix A.
Any functions in A must be specified with the variables x[1],x[2],...
The output is a numbered list of constraints that should be forced to zero by
direct substitution. To aid in the process use the FindConstraint function."
FindConstraint::usage= "FindConstraint[number,variable] solves the numbered
constraint equation for the desired variable."
Curl::usage = "Curl determines the additional curl contraints on Grad[V] so
that the V function can be found by any path integration. It computes
D[GradV[i],x[j]] - D[GradV[j],x[i]] for all possible combinations of i and j.
It's output consists of these constraints and the form of Grad[V]."
FindV::usage="FindV[ Path List ] integrates the Grad[V] function along the
desired straight-line path. If the path list is {1,2,3}, it integrates along
x[1], with x[2] = x[3] = 0; then along x[2], with x[1] = x[1] and x[3] = 0, etc.
The output is the resulting Liapounov Function V."
ResetVG::usage = "ResetVG eliminates all remaining variables, so that the process
can be restarted."
VariableGradient[A_] := Block[ {X,DelV,DelVC,Vd},
X = Array[Global`x,Length[A]];
DelVC = Array[alpha,{Length[A],Length[A]}];
(* alpha[Length[A],Length[A]] = Length[A]-1;*)
DelV = Partition[ DelVC . X,1];
delV = Collect[DelV,X];
{Vd} = Transpose[delV] . A . X;
Constraint = Partition[Flatten[Table[{(i+j-2),Coefficient[Expand[Vd], Global`x[i] Global`x[j] ]},
{i,1,Length[A]},{j,(i + 1), Length[A]} ] ],2];
Vdot = Collect[Expand[Vd],X];
VDOT = TableForm[ {{"Constraints:"},{TableForm[Constraint]},{"V dot:"},{Vdot}}]
]
FindConstraint[n_,var_] := Solve[Constraint[[ n , 2 ]] == 0, var]
Curl := Block[{},
CR = Flatten[Table[D[delV[[i]],Global`x[j] ] - D[delV[[j]],Global`x[i] ], (* Modified x[i] to Global`x[i] and similiarly for x[j] *)
{i,1,Length[delV]},{j,(i + 1),Length[delV]}]];
DELV = TableForm[ {{"Constraints:"},{TableForm[CR]},{"Grad V:"},{TableForm[delV]}}]
]
FindV[Path_List] := Block[{X,Vx,Xa},
X = Array[Global`x,Length[delV]];
Xa = Array[xa,Length[delV] ];
Vx = Table[Integrate[delV[[ Path[[i]] ]]/.Table[ Which[ j < i, X[[ Path[[ j ]] ]] -> Xa[[ Path[[ j ]] ]], j > i , X[[ Path[[ j ]] ]] -> 0, i == j, X[[ Path[[j]] ]] -> X[[ Path[[j]] ]] ], {j,1,Length[delV]} ],{X[[ Path[[i]] ]],0,Xa[[ Path[[ i ]] ]]}], {i,1,Length[delV]} ];
V = Sum[Vx[[i]], {i,1,Length[delV]}]/.Table[Xa[[j]] -> X[[j]],{j,1,Length[delV]} ]
]
ResetVG := Clear[alpha,delV,Constraint,Vdot,CR,V]
EndPackage[]