Runge-Kutta Theory and Constraint Programming Julien Alexandre dit Sandretto Alexandre ...

45 Pages·2016·2.23 MB·English

Runge-Kutta Theory and Constraint Programming
Julien Alexandre dit Sandretto
Alexandre Chapoutot
DepartmentU2IS
ENSTAParisTech
SCAN2016-Uppsala
Contents
Numerical integration
Runge-Kutta with interval coefficients
Constraint approach to define new schemes
Experimentations
Cost function to define optimal schemes
Experimentations
Conclusion
JulienAlexandreditSandretto-Runge-KuttaTheoryandConstraintProgramming September28,2016-2
Numericalintegration
Numerical integration
Initial value problem
y˙ = f(y) with y(0) = y
0
Time discretization t = 0 < t < < t = T
0 1 n
···
Compute a sequence of values: y ,y ,...,y such that
0 1 n
i 0,n , y y(t ;y ) .
i i 0
∀ ∈ { } ≈
h
s.t. y y(t +h;y ) = y + f(y(s))ds
n+1 ≈ n n n 0
(cid:82)
h
Tool: integration scheme to approx f(y(s))ds
0
with an error (hp+1) (cid:82)
O
(cid:73) h is the integration step-size
(cid:73) p is the order of the method (Taylor series or Runge-Kutta)
JulienAlexandreditSandretto-Runge-KuttaTheoryandConstraintProgramming September28,2016-3
Numericalintegration
Numerical integration
Initial value problem
y˙ = f(y) with y(0) = y
0
Time discretization t = 0 < t < < t = T
0 1 n
···
Compute a sequence of values: y ,y ,...,y such that
0 1 n
i 0,n , y y(t ;y ) .
i i 0
∀ ∈ { } ≈
h
s.t. y y(t +h;y ) = y + f(y(s))ds
n+1 ≈ n n n 0
(cid:82)
h
Tool: integration scheme to approx f(y(s))ds
0
with an error (hp+1) (cid:82)
O
(cid:73) h is the integration step-size
(cid:73) p is the order of the method (Taylor series or Runge-Kutta)
JulienAlexandreditSandretto-Runge-KuttaTheoryandConstraintProgramming September28,2016-3
Numericalintegration
Numerical integration
Initial value problem
y˙ = f(y) with y(0) = y
0
Time discretization t = 0 < t < < t = T
0 1 n
···
Compute a sequence of values: y ,y ,...,y such that
0 1 n
i 0,n , y y(t ;y ) .
i i 0
∀ ∈ { } ≈
h
s.t. y y(t +h;y ) = y + f(y(s))ds
n+1 ≈ n n n 0
(cid:82)
h
Tool: integration scheme to approx f(y(s))ds
0
with an error (hp+1) (cid:82)
O
(cid:73) h is the integration step-size
(cid:73) p is the order of the method (Taylor series or Runge-Kutta)
JulienAlexandreditSandretto-Runge-KuttaTheoryandConstraintProgramming September28,2016-3
Numericalintegration
Runge-Kutta schemes
s-stage Runge-Kutta described by a Butcher tableau:
c a a a
1 11 12 1s
. . . ··· . j
. . . .
. . . .
c a a a
s s1 s2 ··· ss i
b b b
1 2 s
···
The integration scheme is given by:
s s
k = f t +c h , y +h a k y = y +h b k
i n i n n ij j n+1 n i i
 
j=1 i=1
(cid:88) (cid:88)
(cid:73) Explicit method (ERK): a = 0 for i (cid:54) j
ij
(cid:73) Diagonal Implicit method (DIRK): a = 0 for i < j
ij
(cid:73) Implicit method (IRK) otherwise
Order p if y(t ;y ) y = C hp+1 with C R.
n n−1 n
− · ∈
JulienAlexandreditSandretto-Runge-KuttaTheoryandConstraintProgramming September28,2016-4
28 II. NumericalIntegrators
II.1.1 Runge–Kutta Methods
Inthissection, wetreat non-autonomous systemsoffirst-orderordinarydifferential
equations
y˙ = f(t,y), y(t ) = y . (1.1)
0 0
t
The integration of this equation gives y(t ) = y + 1 f(t,y(t))dt, and replacing
1 0 t
0
the integral by the trapezoidal rule, we obtain
!
h
y = y + f(t ,y )+f(t ,y ) . (1.2)
1 0 0 0 1 1
2
" #
This is the implicit trapezoidal rule, which, in addition to its historical impor-
tance for computations in partial differential equations (Crank–Nicolson) and in
A-stability theory (Dahlquist), played a crucial role even earlier in the discovery of
Runge–Kutta methods. It was the starting point of Runge (1895), who “predicted”
the unknown y -value to the right by an Euler step, and obtained the first of the
1
following formulas (the second being the analogous formula for the midpoint rule)
k = f(t ,y ) k = f(t ,y )
1 0 0 1 0 0
k = f(t +h,y +hk ) k = f(t + h,y + hk ) (1.3)
2 0 0 1 2 0 2 0 2 1
y = y + h k +k y = y +hk .
1 0 2 1 2 1 0 2
These methods have a ni"ce geome#tric interpretation (which is illustrated in the first
two pictures of Fig.1.2 for a famous problem, the Riccati equation): they consist
of polygonal lines, which assume the slopes prescribed by the differential equation
evaluated at previous points.
IdeaofHeun(1900)and Kutta(1901):compute several polygonal lines,each start-
ing at y0 and assuminNgumtehricealivnteagrraitoionus slopes kj on portions of the integration interval,
which are proportional to some given constants a ; at the final point of each poly-
Runge-Kutta vs Taylor ij
gon evaluate a new slope k . The last of these polygons, with constants b , deter-
i i
miTneasyltohre: nounmlyeorniceamlseothluotdiocnomypu(tseedeftohredtihffierrdenpticotrudreero(ftiFlli1g2.01.!2)).Thisidealeadsto
1
theync+la1s=s oyfne+xplicsi=it1Rhiufn[i]g(ey–n)K⇒uttacommeptuhtoadtiso,ni.fero.,mfoornmlyuloane(1p.o4i)ntbe!low with aij = 0
forRiungje.-Kutta: many methods with different order (often 4),
com≤putation(cid:80)from different points:
y y y
expl.trap.rule expl.midp.rule b b b
1 2 3
y1 y1
strong stability properties for various kinds
y 1
1
1 1of problems (A-stable, L-stable1, algebraic
a a
31 32
stability, etc.), and may preserve quadratic
k
3
algebraic invariant (symplectic methods)
k
2 c
k1 k1 k2 k1 k2 3
y y y
0 1 0 1 0 a c2
2 t 1 2 t 1 21 t 1
Fig.1.2.Runge–Kuttamethodsfory˙ = t2 +y2,y = 0.46,h = 1;dotted:exactsolution
JulienAlexandreditSandretto-Runge-KuttaTheoryandConstrain0tProgramming September28,2016-5
Numericalintegration
Race for higher order schemes
Why higher order ?
High order implies low difference between solution and
approximation !
A global competition
(cid:73) Explicit Runge-Kutta order 14 with 35 stages with Maple [1]
(cid:73) Radau order 17 with 9 stages with Mathematica [2]
[1]Feagin,Terry,“High-orderExplicitRunge-KuttaMethodsUsingM-Symmetry”,Neural,Parallel&Scientific
Computations,Vol.20,No.4,December2012,pp.437-458
[2]JMart´ın-Vaquero,“A17th-orderRadauIIAmethodforpackageRADAU”,Applicationsinmechanicalsystems,
Computers&MathematicswithApplications,2010
JulienAlexandreditSandretto-Runge-KuttaTheoryandConstraintProgramming September28,2016-6
Numericalintegration
New scheme: a complex problem
Needs to solve constraints
High order polynomials (till p), number of constraints increases
exponentially (4 for p = 3, 8 for p = 4, 17, 37, 85, 200)
1. s b = 1
1 i
2. s b c = 1/2
3.(cid:80)1s bici2 = 1/3 s s b a c = 1/6
(cid:80)1 i i 1 1 i ij j
Cl(cid:80)assical approach(cid:80) (cid:80)
Solve by using polynomials with known exact zeros such as
Legendre (for Gauss) or Jacobi (for Radau) [see Butcher]
Problems
(cid:73) Discovery of new methods guided by solver and not by
requirements !
(cid:73) Solved numerically: additive approximations !
JulienAlexandreditSandretto-Runge-KuttaTheoryandConstraintProgramming September28,2016-7
Numericalintegration
New scheme: a complex problem
Needs to solve constraints
High order polynomials (till p), number of constraints increases
exponentially (4 for p = 3, 8 for p = 4, 17, 37, 85, 200)
1. s b = 1
1 i
2. s b c = 1/2
3.(cid:80)1s bici2 = 1/3 s s b a c = 1/6
(cid:80)1 i i 1 1 i ij j
Cl(cid:80)assical approach(cid:80) (cid:80)
Solve by using polynomials with known exact zeros such as
Legendre (for Gauss) or Jacobi (for Radau) [see Butcher]
Problems
(cid:73) Discovery of new methods guided by solver and not by
requirements !
(cid:73) Solved numerically: additive approximations !
JulienAlexandreditSandretto-Runge-KuttaTheoryandConstraintProgramming September28,2016-7