fricas
(1) -> ++
++ Functional derivatives (Frechet derivative)
++ Sample: Lagrangian of a double pendulum
++ Ugly output of derivatives of (univariate) functions :(
++
++ Macros
macro mksn(s,n) == [s[i] for i in 1..n]
Type: Void
fricas
macro mks2 s == mksn(s,2)
Type: Void
fricas
macro mkeq(s,rhs) == [s.i=rhs.i for i in 1..#s]
Type: Void
fricas
-- Frechet derivative (functional derivative, req. for Lagrangian)
DF(f,g) ==
s:=new()$Symbol -- generate new unique symbol
fs:=subst(f,g=s)
rs:=D(fs,s)
r:=subst(rs,s=g)
return r
Type: Void
fricas
-- test: DF(xt.1^n+sin xt.1,xt.1)
-- Setup
q:=map(operator,mks2 q)
Type: List(BasicOperator
?)
fricas
x:=map(operator,mks2 x)
Type: List(BasicOperator
?)
fricas
y:=map(operator,mks2 y)
Type: List(BasicOperator
?)
fricas
l:=mks2 l
Type: List(Symbol)
fricas
m:=mks2 m
Type: List(Symbol)
fricas
-- Make functions of t
qt:=[f(t) for f in q]
Type: List(Expression(Integer))
fricas
xt:=[f(t) for f in x]
Type: List(Expression(Integer))
fricas
yt:=[f(t) for f in y]
Type: List(Expression(Integer))
fricas
-- Geometry, equations for x1,y1,x2,y2
eqx:=mkeq(xt,[l.1*sin qt.1, l.1*sin qt.1 + l.2*sin qt.2])
Type: List(Equation(Expression(Integer)))
fricas
eqy:=mkeq(yt,[-l.1*cos qt.1, -l.1*cos qt.1 - l.2*cos qt.2])
Type: List(Equation(Expression(Integer)))
fricas
-- All in a matrix
meqxy:= matrix [eqx,eqy]
Type: Matrix(Equation(Expression(Integer)))
fricas
-- Potential energy
V:=operator 'V
fricas
eqV:=V(concat[xt,yt])=m.1*g*yt.1 + m.2*g*yt.2
Type: Equation(Expression(Integer))
fricas
-- Kinetic energy
T:=operator 'T
fricas
v:=D(matrix [xt,yt],t) -- velocities, col1=v1, col2=v2
Type: SquareMatrix
?(2,
Expression(Integer))
fricas
vv:=diagonal(transpose(v)*v)
Type: DirectProduct
?(2,
Expression(Integer))
fricas
eqT:=T(concat listOfLists transpose v)=(1/2)*(m.1*vv.1 + m.2*vv.2)
Type: Equation(Expression(Integer))
fricas
-- DF(eqT,D(xt.1,t,1)) -> m1*x1'
-- Lagrangian L=T-V
L:=operator 'L
fricas
Largs:=concat [concat[xt,yt],concat listOfLists transpose v]
Type: List(Expression(Integer))
fricas
eqL:=L(args)=rhs(eqT)-rhs(eqV)
Type: Equation(Expression(Integer))
fricas
-- Side conditions
sc1:=map(normalize, eqx.1^2 + eqy.1^2)
Type: Equation(Expression(Integer))
fricas
sc2:=map(normalize, (eqx.2-eqx.1)^2 + (eqy.2-eqy.1)^2)
Type: Equation(Expression(Integer))
fricas
L2:=operator 'L2
fricas
eqL2:=L2(Largs)=rhs(eqT)-rhs(eqV)+A*(rhs(sc1)-lhs(sc1))+B*(rhs(sc2)-lhs(sc2))
Type: Equation(Expression(Integer))
fricas
-- Lagrange equations with SC
eqL2x:=[rhs(D(DF(eqL2,D(xt.i,t)),t)-DF(eqL2,xt.i))=0 for i in 1..2]
fricas
Compiling function DF with type (Equation(Expression(Integer)),
Expression(Integer)) -> Equation(Expression(Integer))
Type: List(Equation(Expression(Integer)))
fricas
eqL2y:=[rhs(D(DF(eqL2,D(yt.i,t)),t)-DF(eqL2,yt.i))=0 for i in 1..2]
Type: List(Equation(Expression(Integer)))
fricas
eqL2x.1
Type: Equation(Expression(Integer))
fricas
eqL2x.2
Type: Equation(Expression(Integer))
fricas
eqL2y.1
Type: Equation(Expression(Integer))
fricas
eqL2y.2
Type: Equation(Expression(Integer))