login  home  contents  what's new  discussion  bug reports help  links  subscribe  changes  refresh  edit

fricas
++
++ 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)
 (1)
Type: List(BasicOperator?)
fricas
x:=map(operator,mks2 x)
 (2)
Type: List(BasicOperator?)
fricas
y:=map(operator,mks2 y)
 (3)
Type: List(BasicOperator?)
fricas
l:=mks2 l
 (4)
Type: List(Symbol)
fricas
m:=mks2 m
 (5)
Type: List(Symbol)
fricas
-- Make functions of t
qt:=[f(t) for f in q]
 (6)
Type: List(Expression(Integer))
fricas
xt:=[f(t) for f in x]
 (7)
Type: List(Expression(Integer))
fricas
yt:=[f(t) for f in y]
 (8)
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])
 (9)
Type: List(Equation(Expression(Integer)))
fricas
eqy:=mkeq(yt,[-l.1*cos qt.1, -l.1*cos qt.1 - l.2*cos qt.2])
 (10)
Type: List(Equation(Expression(Integer)))
fricas
-- All in a matrix
meqxy:= matrix [eqx,eqy]
 (11)
Type: Matrix(Equation(Expression(Integer)))
fricas
-- Potential energy
V:=operator 'V
 (12)
Type: BasicOperator?
fricas
eqV:=V(concat[xt,yt])=m.1*g*yt.1 + m.2*g*yt.2
 (13)
Type: Equation(Expression(Integer))
fricas
-- Kinetic energy
T:=operator 'T
 (14)
Type: BasicOperator?
fricas
v:=D(matrix [xt,yt],t) -- velocities, col1=v1, col2=v2
 (15)
Type: SquareMatrix?(2,Expression(Integer))
fricas
vv:=diagonal(transpose(v)*v)
 (16)
Type: DirectProduct?(2,Expression(Integer))
fricas
eqT:=T(concat listOfLists transpose v)=(1/2)*(m.1*vv.1 + m.2*vv.2)
 (17)
Type: Equation(Expression(Integer))
fricas
-- DF(eqT,D(xt.1,t,1)) -> m1*x1'
-- Lagrangian L=T-V
L:=operator 'L
 (18)
Type: BasicOperator?
fricas
Largs:=concat [concat[xt,yt],concat listOfLists transpose v]
 (19)
Type: List(Expression(Integer))
fricas
eqL:=L(args)=rhs(eqT)-rhs(eqV)
 (20)
Type: Equation(Expression(Integer))
fricas
-- Side conditions
sc1:=map(normalize, eqx.1^2 + eqy.1^2)
 (21)
Type: Equation(Expression(Integer))
fricas
sc2:=map(normalize, (eqx.2-eqx.1)^2 + (eqy.2-eqy.1)^2)
 (22)
Type: Equation(Expression(Integer))
fricas
L2:=operator 'L2
 (23)
Type: BasicOperator?
fricas
eqL2:=L2(Largs)=rhs(eqT)-rhs(eqV)+A*(rhs(sc1)-lhs(sc1))+B*(rhs(sc2)-lhs(sc2))
 (24)
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))
 (25)
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]
 (26)
Type: List(Equation(Expression(Integer)))
fricas
eqL2x.1
 (27)
Type: Equation(Expression(Integer))
fricas
eqL2x.2
 (28)
Type: Equation(Expression(Integer))
fricas
eqL2y.1
 (29)
Type: Equation(Expression(Integer))
fricas
eqL2y.2
 (30)
Type: Equation(Expression(Integer))

 Subject:   Be Bold !! ( 14 subscribers )