Inverse variational backward error analysis for the variational midpoint rule
Author
Christian Offen
Title
Inverse variational backward error analysis for the variational midpoint 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-608giky/
DOI
https://notebookarchive.org/2022-01-608giky
Date Added
2022-01-13
Date Last Modified
2022-01-13
File Size
0.6 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
Midpoint 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[]:=
LdiscMP=L[(y[t+h/2]+y[t-h/2])/2,(y[t+h/2]-y[t-h/2])/h]
Out[]=
Ly-+t+y+t,
1
2
h
2
h
2
-y-+t+y+t
h
2
h
2
h
Computation of modified Lagrangian
Expansion of discrete Lagrangian
In[]:=
Ldisc=Series[LdiscMP,{h,0,o}]
Out[]=
L[y[t],[t]]+[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]]+
′
y
1
24
(3)
y
(0,1)
L
′
y
′′
y
(1,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
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]]-[t][y[t],[t]]+16[t][y[t],[t]]+7[t][y[t],[t]]-8[t][y[t],[t]]+8[t][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]]+6[t]6[y[t],[t]]+4[t][y[t],[t]]+7[t][y[t],[t]]+4[t]2[t][y[t],[t]]+3[t][y[t],[t]]+7[t][y[t],[t]]+7[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
(4)
y
(1,0)
L
′
y
′
y
(4)
y
(1,1)
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
2
′
y
(2,2)
L
′
y
′′
y
(4)
y
(0,2)
L
′
y
2
′
y
(3,0)
L
′
y
3
′
y
(3,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’’ as series in h
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 ⋯28⋯ 2 ′ y 3 ⋯1⋯ (3,0) L ′ y 3 ′ y 3 (0,2) L ′ y (3,1) L ′ y 24 4 (0,2) L ′ y 1 5760 7 (0,2) L ′ y (6) y 7 (0,2) L ′ y (5) y 5 (0,2) L ′ y (0,3) L ′ y (1,0) L ′ y ⋯246⋯ 4 ′ y 6 (0,2) L ′ y (5,0) L ′ y 5 ′ y 6 (0,2) L ′ y (5,1) 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 1 24 4 (0,2) L ′ y (4) y 4 (0,2) L ′ y (0,4) L ′ y 3 (1,0) L ′ y ′ y (0,4) L ′ y 2 (1,0) L ′ y (1,1) L ′ y ⋯20⋯ 3 ′ y 3 (0,2) L ′ y (3,1) L ′ y 2 h 1 5760 7 (0,2) L ′ y (6) y 7 (0,2) L ′ y 2 (3) y 4 (0,2) L ′ y 2 (0,3) L ′ y (1,0) L ′ y ⋯244⋯ 5 ′ y 6 (0,2) L ′ y (5,1) L ′ y 4 h 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]];
Equations to determine coefficients Ł of ansatz
In[]:=
eqs=#10&/@PadRight[CoefficientList[LmodinvAns,h][[2;;-1]],o]
Out[]=
Coefficients of ansatz
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
Solve for coefficients
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;
Output of inverse modified Lagrangian
In[]:=
Linv=FullSimplify[Linv]
Out[]=
Ł
0
1
5760[a0,a1]
5
(0,2)
Ł
0
2
h
2
h
2
(0,3)
Ł
0
4
[a0,a1]-a1[a0,a1]
(1,0)
Ł
0
(1,1)
Ł
0
2
h
(0,2)
Ł
0
3
[a0,a1]-a1[a0,a1]
(1,0)
Ł
0
(1,1)
Ł
0
(0,4)
Ł
0
(1,0)
Ł
0
(1,1)
Ł
0
(0,3)
Ł
0
(1,2)
Ł
0
2
h
2
(0,2)
Ł
0
2
[a0,a1]-a1[a0,a1]
(1,0)
Ł
0
(1,1)
Ł
0
(1,0)
Ł
0
(1,2)
Ł
0
(1,3)
Ł
0
2
(1,2)
Ł
0
(1,1)
Ł
0
(1,2)
Ł
0
(1,3)
Ł
0
(0,3)
Ł
0
(2,0)
Ł
0
(2,1)
Ł
0
2
h
3
(0,2)
Ł
0
(1,0)
Ł
0
(1,1)
Ł
0
(1,0)
Ł
0
(1,1)
Ł
0
(1,2)
Ł
0
(2,0)
Ł
0
(2,1)
Ł
0
2
a1
(1,0)
Ł
0
(1,1)
Ł
0
(2,2)
Ł
0
4
(0,2)
Ł
0
2
(1,0)
Ł
0
(1,0)
Ł
0
(1,1)
Ł
0
2
h
(3,0)
Ł
0
(3,1)
Ł
0
2
a1
2
(1,1)
Ł
0
2
h
2
[a0,a1]-a1[a0,a1]
(2,0)
Ł
0
(2,1)
Ł
0
2
h
(1,1)
Ł
0
(3,0)
Ł
0
(3,1)
Ł
0
5
(0,2)
Ł
0
2
a1
(2,0)
Ł
0
4
a1
2
h
(4,0)
Ł
0
Validity check:
In[]:=
Lmod/.{y[t]a0,y'[t]a1}/.LFunction[{a0,a1},Evaluate[Linv]]
Out[]=
Ł
0
5
O[h]
Cite this as: Christian Offen, "Inverse variational backward error analysis for the variational midpoint rule" from the Notebook Archive (2022), https://notebookarchive.org/2022-01-608giky
Download
