In[]:=
ϕ[x_]:=
&,ℒ[n_,a_,b_]:=Module{q,M,V,θ,Θ,z},q=Table[Slot[i],{i,n}];=Table[Slot[n+i],{i,n}];Θ=Table[,{k,n}];z={Cos@#,Sin@#}&/@Θ;M=Accumulate[MapThread[Dot,{ϕ/@q,z}]];V=Accumulate[MapThread[Dot,{ϕ'/@q,z}]];ConstructFunction,SimplifyPlus@@MapThread[Dot,{V,V}]-Plus@@({0,1}.M)/.Cos[]:>a[[k]],Sin[]:>b[[k]]EulerLagrange[n_,a_,b_]:=Table[Dt[(Derivative@@Table[KroneckerDelta[n+k,j],{j,2n}])[ℒ[n,a,b]]@@Table[Slot[i][t],{i,2n}],t,Constants->Join[Table[{a[[k]],b[[k]]},{k,n}]]]==(Derivative@@Table[KroneckerDelta[k,j],{j,2n}])[ℒ[n,a,b]]@@Table[Slot[i][t],{i,2n}]/.Slot[i:_]/;i>n:>Slot[i-n]',{k,n}]
#1 | -#2 |
#2 | #1 |
1-
2
x
1+
2
x
2x
1+
2
x
q
θ
k
q
1
2
θ
(k:_)
θ
(k:_)
In[]:=
Module[{n=3,tmax=60,q},Block[{s=ParametricNDSolve[Join[Construct[Function,EulerLagrange[n,Table[,{k,n}],Table[,{k,n}]]]@@Table[,{k,n}],Table[[0]==0,{k,n}],Table['[0]==If[k==3,.4,0],{k,n}]],Table[,{k,n}],{t,tmax},Join[Table[,{k,n}],Table[,{k,n}]],Method->{"EquationSimplification"->"Residual"}],params=Join[ConstantArray[0,n],ConstantArray[-1,n]]},Plot[Evaluate[#[t]&/@(Table[@@params,{k,n}]/.s)],{t,0,tmax}]]]
a
k
b
k
q
k
q
k
q
k
q
k
a
k
b
k
q
k
Out[]=