Inverse variational backward error analysis for the variational trapezoidal rule
Author
Christian Offen
Title
Inverse variational backward error analysis for the variational trapezoidal rule
Description
Accompanying source code for the paper "Variational integration of learned dynamical systems" (arXiv:2112.12619)
Category
Academic Articles & Supplements
Keywords
backward error analysis, variational integrators
URL
http://www.notebookarchive.org/2022-01-608jr4t/
DOI
https://notebookarchive.org/2022-01-608jr4t
Date Added
2022-01-13
Date Last Modified
2022-01-13
File Size
0.79 megabytes
Supplements
Rights
CC BY 4.0
Download
Open in Wolfram Cloud
The notebook is part of accompanying source code for Sina Ober-Blöbaum and Christian Offen, “Variational integration of learned dynamical systems,” https://arxiv.org/abs/2112.12619.
Further source code and experiments can be found on https://github.com/Christian-Offen/LagrangianShadowIntegration.
Further source code and experiments can be found on https://github.com/Christian-Offen/LagrangianShadowIntegration.
Calculation of modified and inverse modified Lagrangian
Calculation of modified and inverse modified Lagrangian
Trapezoidal rule
Input: consistent discretisation of L[y[t],y’[t]] depending on (y(t-h/2),y(t+h/2)), 1 spacial dimension
truncation order o
In[]:=
o=4;
discrete Lagrangian
In[]:=
LdiscTrz=1/2L[y[t-h/2],(y[t+h/2]-y[t-h/2])/h]+1/2L[y[t+h/2],(y[t+h/2]-y[t-h/2])/h]
Out[]=
1
2
h
2
-y-+t+y+t
h
2
h
2
h
1
2
h
2
-y-+t+y+t
h
2
h
2
h
Computation of modified Lagrangian
Expansion of discrete Lagrangian
In[]:=
Ldisc=Series[LdiscTrz,{h,0,o}]
Out[]=
L[y[t],[t]]+[t][y[t],[t]]+3[t][y[t],[t]]+3[t][y[t],[t]]+3[t][y[t],[t]]+5[t][y[t],[t]]+15[t][y[t],[t]]+30[t][t][y[t],[t]]+45[t][y[t],[t]]+60[t][t][y[t],[t]]+30[t][t][y[t],[t]]+90[t][t][y[t],[t]]+15[t][y[t],[t]]+
′
y
1
24
(3)
y
(0,1)
L
′
y
′′
y
(1,0)
L
′
y
2
′
y
(2,0)
L
′
y
2
h
1
5760
(5)
y
(0,1)
L
′
y
2
(3)
y
(0,2)
L
′
y
(4)
y
(1,0)
L
′
y
′′
y
(3)
y
(1,1)
L
′
y
2
′′
y
(2,0)
L
′
y
′
y
(3)
y
(2,0)
L
′
y
2
′
y
(3)
y
(2,1)
L
′
y
2
′
y
′′
y
(3,0)
L
′
y
4
′
y
(4,0)
L
′
y
4
h
5
O[h]
Meshed Lagrangian
In[]:=
Lmesh0=Sum[(2^(1-2i)-1)h^(2i)BernoulliB[2i]/Factorial[2i]D[Ldisc,{t,2i}],{i,0,o/2}];
In[]:=
Lmesh=Simplify[Series[Lmesh0,{h,0,o}]]
Out[]=
L[y[t],[t]]+-[t][y[t],[t]]+2[t][y[t],[t]]-[t][y[t],[t]]+2[t][y[t],[t]]+16[t][y[t],[t]]+7[t][y[t],[t]]+4[t]3[y[t],[t]]+7[t][y[t],[t]]+32[t][t][y[t],[t]]+2[t][t][y[t],[t]]+[t]-[y[t],[t]]+[t][y[t],[t]]-12[t]2[y[t],[t]]+[t]8[y[t],[t]]-[t][y[t],[t]]+8[t][t][y[t],[t]]-2[t]3[y[t],[t]]+2[t][y[t],[t]]-8[t][y[t],[t]]-[t][y[t],[t]]+[t][y[t],[t]]+
′
y
1
24
2
′′
y
(0,2)
L
′
y
′′
y
(1,0)
L
′
y
′
y
(1,1)
L
′
y
2
′
y
(2,0)
L
′
y
2
h
1
5760
2
(3)
y
(0,2)
L
′
y
4
′′
y
(0,4)
L
′
y
3
′′
y
(1,2)
L
′
y
′
y
(1,3)
L
′
y
(3)
y
2
′′
y
(0,3)
L
′
y
′
y
′′
y
(1,2)
L
′
y
′
y
(2,0)
L
′
y
′
y
(2,1)
L
′
y
2
′′
y
(2,0)
L
′
y
′
y
(2,1)
L
′
y
′
y
(2,2)
L
′
y
′′
y
(4)
y
(0,2)
L
′
y
2
′
y
(3,0)
L
′
y
′
y
(3,1)
L
′
y
(4)
y
(1,0)
L
′
y
′
y
(1,1)
L
′
y
4
′
y
(4,0)
L
′
y
4
h
5
O[h]
Reduction of Lmesh to first order
Reduction of Lmesh to first order
In[]:=
Needs["VariationalMethods`"]
Variational derivative of disceret Lagrangian
In[]:=
elLdisc=VariationalD[Normal[Ldisc],y[t],t]
Out[]=
Computation of y’’ in series
In[]:=
y2Ans=Sum[a[k]h^k,{k,0,o}]
Out[]=
a[0]+ha[1]+a[2]+a[3]+a[4]
2
h
3
h
4
h
In[]:=
coeff=(List@@Series[elLdisc/.y''[t]y2Ans,{h,0,o}])[[3]]
Out[]=
In[]:=
sola=Solve[#10&/@coeff,Table[a[k],{k,0,o}]]
Out[]=
a[0] (1,0) L ′ y ′ y (1,1) L ′ y (0,2) L ′ y ⋯1⋯ 24 4 (0,2) L ′ y ⋯1⋯ 5760 7 (0,2) L ′ y | |||||
|
In[]:=
y2=Series[y2Ans/.sola[[1]],{h,0,o}];
output of expression for y'' as a series in h
In[]:=
y2=y2//Simplify
Out[]=
(1,0) L ′ y ′ y (1,1) L ′ y (0,2) L ′ y 2 (4) y 4 (0,2) L ′ y ⋯21⋯ 3 ′ y 3 ⋯1⋯ (3,1) L ′ y 2 h 24 4 (0,2) L ′ y ( ⋯1⋯ 4 h 5760 7 ⋯1⋯ L ⋯1⋯ 5 O[h] | |||||
|
rule for substitution of higher derivatives
In[]:=
replY={Derivative[k_][y][t]/;k≥2->D[y2,{t,k-2}]};
substitution of higher derivatives in Lmesh
In[]:=
Lmod=Lmesh//.replY;
Output of modified Lagrangian
In[]:=
Lmod=FullSimplify[Lmod]
Out[]=
Sanity check: Check if solution to EL(Lmod)==0 solves EL(Ldisc)==0.
In[]:=
err=VariationalD[Lmod,y[t],t]//.replY
Out[]=
5
O[h]
Computation of inverse modified Lagrangian Linv such that (Linv)mod == L
Taylor series ansatz
In[]:=
LinvAns=Sum[h^kSubscript[Ł,k][a0,a1],{k,0,o}]
Out[]=
Ł
0
Ł
1
2
h
Ł
2
3
h
Ł
3
4
h
Ł
4
Substitution of ansatz into Lmod
In[]:=
LmodinvAns=Lmod/.L->Function[{a0,a1},Evaluate[LinvAns]]
Out[]=
In[]:=
eqs=#10&/@PadRight[CoefficientList[LmodinvAns,h][[2;;-1]],o]
Out[]=
In[]:=
varsLinv=Table[Subscript[Ł,k][y[t],y'[t]],{k,0,o}][[2;;-1]]
Out[]=
{[y[t],[t]],[y[t],[t]],[y[t],[t]],[y[t],[t]]}
Ł
1
′
y
Ł
2
′
y
Ł
3
′
y
Ł
4
′
y
Coefficients of ansatz
In[]:=
solLinv=Solve[eqs,varsLinv][[1]]/.{y[t]a0,y'[t]a1}
Out[]=
Substitution rule
In[]:=
replLk=Table[->Function[{a0,a1},Evaluate[solLinv[[k]][[2]]]],{k,1,o}]
Ł
k
Out[]=
Substitution of determined coefficients into ansatz
In[]:=
Linv=LinvAns//.replLk;
In[]:=
Linv=FullSimplify[O[h]^(o+1)+Linv]
Out[]=
Ł
0
1
24
2
[a0,a1]-a1[a0,a1]
(1,0)
Ł
0
(1,1)
Ł
0
(0,2)
Ł
0
2
a1
(2,0)
Ł
0
2
h
1
5760[a0,a1]
5
(0,2)
Ł
0
2
(0,3)
Ł
0
4
[a0,a1]-a1[a0,a1]
(1,0)
Ł
0
(1,1)
Ł
0
(0,2)
Ł
0
(0,3)
Ł
0
2
[a0,a1]-a1[a0,a1]
(1,0)
Ł
0
(1,1)
Ł
0
(1,0)
Ł
0
(1,1)
Ł
0
(1,2)
Ł
0
(0,2)
Ł
0
(2,0)
Ł
0
(0,2)
Ł
0
(2,1)
Ł
0
(0,2)
Ł
0
(0,4)
Ł
0
4
[a0,a1]-a1[a0,a1]
(1,0)
Ł
0
(1,1)
Ł
0
(0,2)
Ł
0
3
(1,0)
Ł
0
(1,2)
Ł
0
(1,3)
Ł
0
2
(1,0)
Ł
0
2
a1
2
(1,2)
Ł
0
(1,1)
Ł
0
(1,2)
Ł
0
(1,3)
Ł
0
(0,2)
Ł
0
(2,0)
Ł
0
(2,1)
Ł
0
(2,2)
Ł
0
(1,0)
Ł
0
2
(1,1)
Ł
0
(1,2)
Ł
0
(1,3)
Ł
0
(1,1)
Ł
0
2
a1
2
(1,2)
Ł
0
(0,2)
Ł
0
(2,0)
Ł
0
(2,1)
Ł
0
2
a1
(2,2)
Ł
0
(0,2)
Ł
0
(1,2)
Ł
0
(2,0)
Ł
0
(2,1)
Ł
0
(0,2)
Ł
0
(3,0)
Ł
0
(3,1)
Ł
0
2
a1
3
(1,1)
Ł
0
(1,2)
Ł
0
(1,3)
Ł
0
2
(1,1)
Ł
0
2
a1
2
(1,2)
Ł
0
(0,2)
Ł
0
(2,0)
Ł
0
(2,1)
Ł
0
(2,2)
Ł
0
(0,2)
Ł
0
(1,1)
Ł
0
(1,2)
Ł
0
(2,0)
Ł
0
(2,1)
Ł
0
(0,2)
Ł
0
(3,0)
Ł
0
(3,1)
Ł
0
2
(0,2)
Ł
0
2
[a0,a1]-a1[a0,a1]
(2,0)
Ł
0
(2,1)
Ł
0
2
a1
(0,2)
Ł
0
(4,0)
Ł
0
4
h
5
O[h]
Validity check:
Lmod/.{y[t]a0,y'[t]a1}/.LFunction[{a0,a1},Evaluate[Normal[Linv]]]
Out[]=
Ł
0
5
O[h]
Nice outputs
Nice outputs
In[]:=
funcs={Plus,Times,Power,Sin,Equal,MatrixForm,List,Rational,SeriesData};SetAttributes[HideArgs,HoldAll];HideArgs[expr_]:=expr/.{x_[__]/;And@@(UnsameQ[x,#]&/@funcs)x}
Lmod
In[]:=
HideArgs[Lmod]
Out[]=
L++2+-16-322--++7-4-3+7+6+33+7+24+16-3-12+9+48+4(4-3)+21+32-6-8(2+)+23+2+3+24-8+4(4-3)+7+16-3-8(2+)+22++
1
24
2
-
(1,0)
L
′
y
(1,1)
L
(0,2)
L
2
()
′
y
(2,0)
L
2
h
1
5760
5
()
(0,2)
L
2
()
(0,3)
L
4
-
(1,0)
L
′
y
(1,1)
L
′
y
(0,2)
L
(0,3)
L
2
-
(1,0)
L
′
y
(1,1)
L
(1,0)
L
′
y
(1,1)
L
(1,2)
L
(0,2)
L
(2,0)
L
′
y
(0,2)
L
(2,1)
L
(0,2)
L
(0,4)
L
4
-
(1,0)
L
′
y
(1,1)
L
(0,2)
L
3
()
(1,0)
L
(1,2)
L
′
y
(1,3)
L
2
()
(1,0)
L
(0,2)
L
(2,0)
L
′
y
(1,1)
L
(1,2)
L
′
y
(1,3)
L
(0,2)
L
(2,1)
L
′
y
2
()
(1,2)
L
(0,2)
L
(2,2)
L
′
y
(1,0)
L
(0,2)
L
(1,1)
L
(2,0)
L
′
y
2
()
(1,1)
L
(1,2)
L
(0,2)
L
(1,1)
L
(2,1)
L
(0,2)
L
(1,2)
L
(2,0)
L
(0,2)
L
(3,0)
L
2
()
′
y
2
()
(1,1)
L
(1,3)
L
(1,1)
L
2
()
(1,2)
L
(0,2)
L
(2,2)
L
(0,2)
L
(1,2)
L
(2,1)
L
(0,2)
L
(3,1)
L
2
()
′
y
(0,2)
L
(2,0)
L
2
()
(1,1)
L
(0,2)
L
(2,0)
L
′
y
3
()
(1,1)
L
(1,2)
L
(0,2)
L
2
()
(1,1)
L
(2,1)
L
2
()
(0,2)
L
(2,0)
L
(2,1)
L
(0,2)
L
(1,1)
L
(1,2)
L
(2,0)
L
(0,2)
L
(3,0)
L
2
()
′
y
3
()
(1,1)
L
(1,3)
L
2
()
(1,1)
L
2
()
(1,2)
L
(0,2)
L
(2,2)
L
(0,2)
L
(1,1)
L
(1,2)
L
(2,1)
L
(0,2)
L
(3,1)
L
2
()
(0,2)
L
2
()
(2,1)
L
(0,2)
L
(4,0)
L
4
h
5
O[h]
Lmodinv
In[]:=
HideArgs[(Linv/.L/.a0y[t]/.a1y'[t])]
Ł
0
Out[]=
L+--2+-4-82--++3+4-3-+11+9-+32-2++9-++2-11+3-2+2++2-13-+6-+3-+11+32-2++213--6-+27+6+
1
24
2
-
(1,0)
L
′
y
(1,1)
L
(0,2)
L
2
()
′
y
(2,0)
L
2
h
1
5760
5
()
(0,2)
L
2
()
(0,3)
L
4
-
(1,0)
L
′
y
(1,1)
L
′
y
(0,2)
L
(0,3)
L
2
-
(1,0)
L
′
y
(1,1)
L
(1,0)
L
′
y
(1,1)
L
(1,2)
L
(0,2)
L
(2,0)
L
′
y
(0,2)
L
(2,1)
L
(0,2)
L
(0,4)
L
4
-
(1,0)
L
′
y
(1,1)
L
(0,2)
L
3
()
(1,0)
L
(1,2)
L
′
y
(1,3)
L
2
()
(1,0)
L
2
()
′
y
2
()
(1,2)
L
′
y
(1,1)
L
(1,2)
L
′
y
(1,3)
L
(0,2)
L
(2,0)
L
′
y
(2,1)
L
′
y
(2,2)
L
′
y
(1,0)
L
′
y
2
()
(1,1)
L
(1,2)
L
′
y
(1,3)
L
(1,1)
L
2
()
′
y
2
()
(1,2)
L
(0,2)
L
(2,0)
L
′
y
(2,1)
L
2
()
′
y
(2,2)
L
′
y
(0,2)
L
(1,2)
L
(2,0)
L
′
y
(2,1)
L
(0,2)
L
(3,0)
L
′
y
(3,1)
L
2
()
′
y
′
y
3
()
(1,1)
L
(1,2)
L
′
y
(1,3)
L
2
()
(1,1)
L
2
()
′
y
2
()
(1,2)
L
(0,2)
L
(2,0)
L
′
y
(2,1)
L
′
y
(2,2)
L
′
y
(0,2)
L
(1,1)
L
(1,2)
L
(2,0)
L
′
y
(2,1)
L
(0,2)
L
(3,0)
L
′
y
(3,1)
L
2
()
(0,2)
L
2
-
(2,0)
L
′
y
(2,1)
L
2
()
′
y
(0,2)
L
(4,0)
L
4
h
5
O[h]
Cite this as: Christian Offen, "Inverse variational backward error analysis for the variational trapezoidal rule" from the Notebook Archive (2022), https://notebookarchive.org/2022-01-608jr4t
Download
