lines 72 and 75 have unnecessary horizontal scrollbars
Higher-Order Equation of Motion
Our goal is to solve the equation of motion for the high-order acceleration f, a given higher-order Lagrangian L.
Related worksheet: https://cloud.sagemath.com/projects/b04b5777-e269-4c8f-a4b8-b21dbe1c93c6/files/Two%20Equations%20of%20Lagrange.sagews
Calculations are done using SageManifolds.
%typeset_mode True
from sage.manifolds.utilities import exterior_derivative as d def ev(N): return (lambda x: N.contract(x))
Variables
M = Manifold(5,'M') coord.<t, x, v, a, s> = M.chart()
Vectors (partial derivatives)
[Dt,Dx,Dv,Da,Ds] = coord.frame()
Forms
[dt,dx,dv,da,ds] = coord.coframe()
General Lagrangian
L = M.scalar_field(function('L')(*list(coord))); L.display()
Kinematics
f is an unknown function that is the analogue of acceleration for a higher-order equation of motion.
f = M.scalar_field(function('f')(*list(coord))); f.display()
N = Dt + v*Dx + a*Dv + s*Da + f*Ds; N.display()
The Equation of Lagrange can be defined with the aid of the following auxillary fields
r=Ds(L); r.display()
q=Da(L)-N(r); q.display()
p=Dv(L)-N(q); p.display()
eq1 = (N(p)-Dx(L)); eq1.display()
Number of terms in this expression:
len(eq1.expr().expand())
We want to solve this equation (eq1 = 0) for f in terms of L.
The equation can also be written in the following form:
eq2 = N(Dv(L)) - N(N(Da(L))) + N(N(N(Ds(L)))) - Dx(L)
N.display()
eq1 == eq2
But we would prefer to write it as a polynomial in N(N(f)), N(f), and powers of f. If we write f instead of N(Ds(L)), N(Da(L)) and N(Dv(L)) above then we get new coefficients and terms in powers of f.
N(N(f)).display()
Terms 2$^\textit{nd}orderinfarisefromtheLiebnizruleappliedtothedifferentialoperatorN^2$ on the term Ds f in N(Ds(L)). For example
t2 = eq1.expr().coefficient(diff(f.expr(),t,t))*diff(f.expr(),t,t); M.scalar_field(t2).display()
The necessary coefficient can be read from the above.
c2 = (Ds(Ds(L))*N(N(f))); c2.display()
After removing the term c2 above, the remaining terms are at most first-order in f.
eq1a = eq1-c2; eq1a.display()
Similarly the first-order terms arise from N applied to the term Ds f in N(Da(L)).
N(f).display()
For example:
t1 = eq1a.expr().coefficient(diff(f.expr(),t))*diff(f.expr(),t); M.scalar_field(t1).display()
x1 = eq1a.expr().coefficient(diff(f.expr(),x))*diff(f.expr(),x); M.scalar_field(x1).display()
v1 = eq1a.expr().coefficient(diff(f.expr(),v))*diff(f.expr(),v); M.scalar_field(v1).display()
a1 = eq1a.expr().coefficient(diff(f.expr(),a))*diff(f.expr(),a); M.scalar_field(a1).display()
s1 = eq1a.expr().coefficient(diff(f.expr(),s))*diff(f.expr(),s); M.scalar_field(s1).display()
Liebniz rule applies N to these coefficients.
c1 = 3*N(Ds(Ds(L)))*N(f); c1.display()
bool(c1.expr()==t1+x1+v1+a1+s1)
Removing the 1$^\textit{st}$ order terms:
eq1b = eq1a-c1; eq1b.display()
Only terms algebraic in f remain. The number of terms is:
len(eq1b.expr().expand())
The term of highest degree comes from N3 applied to Ds(L):
t0 = eq1b.expr().coefficient(f.expr()^3)*f.expr()^3; M.scalar_field(t0.expand()).display()
cf3 = Ds(Ds(Ds(Ds(L))))*f^3; cf3.display()
Removing the 3$^{rd}$ degree term:
eq1c = eq1b-cf3; eq1c.display()
leaves terms of at most 2$^{nd}degree.These2^{nd}degreetermscomeforN^2$ applied to Da(L),
t0 = eq1c.expr().coefficient(f.expr()^2)*f.expr()^2; M.scalar_field(t0.expand()).display()
cf2 = N(Ds(Da(L)))*f; cf2.display()
cf2 = Ds(Ds(Da(L)))*f^2; cf2.display()
We cannot solve this for f
algebraically since f
appears as derivatives.
Lagrangian linear in s.
ll=list(coord);ll.remove(s);ll
Ll = M.scalar_field(function('L0')(*ll)) + s * M.scalar_field(function('L1')(*ll)); Ll.display()
solve(eq1.expr().substitute_function(L.expr().operator(),Ll.expr().function(*list(coord))),f.expr())[0]
Second Equation of Lagrange
(Dt(L)-(N(L)-(p*a+q*s+r*f)-(v*N(p)+a*N(q)+s*N(r)))).display()
Multiplying by v
v*(N(p)-Dx(L)) == (Dt(L)-(N(L)-(p*a+q*s+r*f)-(v*N(p)+a*N(q)+s*N(r))))
Checking the calculations from the paper
L = M.scalar_field(function('L')(*list(coord))) p = M.scalar_field(function('p')(*list(coord))) q = M.scalar_field(function('q')(*list(coord))) r = M.scalar_field(function('r')(*list(coord))) t=M.scalar_field(t) x=M.scalar_field(x) v=M.scalar_field(v) a=M.scalar_field(a) s=M.scalar_field(s)
Action differential Form
alpha = L*dt + p*(dx-v*dt) + q*(dv-a*dt) + r*(da-s*dt) alpha.display()
alpha == L*dt+p*dx+q*dv+r*da-(p*v+q*a+r*s)*dt
d(alpha).display()
d(alpha) == d(L).wedge(dt) + d(p).wedge(dx) + d(q).wedge(dv) + d(r).wedge(da) - d(p*v + q*a + r*s).wedge(dt)
ev(alpha)(N)==L
Omega = -(p*dx+q*dv+r*da).wedge(d(t)); Omega.display()
alpha == L*dt + ev(N)(Omega)
Equation of Motion (E=0)
E = ev(N)(d(alpha)) E.display()
Rewriting it in various ways.
E == N(L)*dt - N(t) * d(L) + N(p)*dx - N(x)*d(p) + N(q)*dv - N(v)*d(q) + N(r)*da - N(a)*d(r) - N(p*v+q*a+r*s)*dt + N(t)*d(p*v+q*a+r*s)
E == N(L)*dt - d(L) + N(p)*dx - v*d(p) + N(q)*dv - a*d(q) + N(r)*da - s*d(r) - N(p*v+q*a+r*s)*dt + d(p*v+q*a+r*s)
E == N(L-p*v-q*a-r*s)*dt - d(L - p*v - q*a- r*s) - v*d(p) - a*d(q) - s*d(r) + N(p)*dx + N(q)*dv + N(r)*da
d(p*v) == v*d(p) + p*d(v)
d(q*a) == a*d(q) + q*d(a)
d(r*s) == s*d(r) + r*d(s)
E == N(L - p*v - q*a - r*s)*dt - d(L) + p*dv + q*da + r*ds + N(p)*dx + N(q)*dv + N(r)*da
N(p*v) == p*N(v) + v*N(p)
E == N(L)*dt - (p*a+q*s+r*f)*dt -(v*N(p) + a*N(q) + s*N(r) )*dt - d(L) + p*dv + q*da + r*ds + N(p)*dx + N(q)*dv + N(r)*da
E == ( N(L) - (p*a + q*s + r*f) - (v*N(p) + a*N(q) + s*N(r)))*dt + N(p)*dx + (N(q)+p)*dv + (N(r)+q)*da + r*ds - d(L)
d(L) == Dt(L)*dt + Dx(L)*dx + Dv(L)*dv + Da(L)*da + Ds(L)*ds
r=Ds(L); r.display()
q=Da(L)-N(r); q.display()
p=Dv(L)-N(q); p.display()
For example: The Schiff and Poirier Lagrangian
hbar = var('hbar',latex_name='\hbar') m = var('m') V = M.scalar_field(function('V')(var('x'))) Lp = 1/2*m*v^2 - V - hbar^2/4/m*(s/v^3-5/2*a^2/v^4); Lp.display()
Lp1 = Ds(Lp); Lp1.display()
Lp0 = Lp - s * Lp1; Lp0.display()
Lp == Lp0 + s*Lp1
eq1p = eq1.expr().substitute_function(L.expr().operator(),Lp.expr().function(*list(coord))); eq1p
solve(eq1p,f.expr())