Linear Differential Equations as a Data-Structure Bruno Salvy - PowerPoint PPT Presentation

Linear Differential Equations 
 as a Data-Structure Bruno Salvy Inria & ENS de Lyon FoCM, July 14, 2017

Computer Algebra Effective mathematics: what can we compute exactly? And complexity: how fast? (also, how big is the result?) Systems with several million users 50+ years of algorithmic progress in computational mathematics! 1/26

Sources of Linear Differential Equations Generating functions 
 in combinatorics M. Bousquet-Mélou Periods Classical elementary 
 and special functions A. Bostan (small order) P. Lairez 2/26

LDEs as a Data-Structure Numerical evaluation Polynomial equations Local and asymptotic Linear expansions Diagonals Differential Equations Proofs of identities Definite sums and integrals Closed forms Conversions M. Singer Solutions called differentially finite (abbrev. D-finite) 3/26

A. Using Linear Differential Equations Exactly

A. Using Linear Differential Equations Exactly I. Numerical Values

̃ ̃ ̃ Fast Computation with Linear Recurrences (70’s and 80’s) 1. Multiplication of integers is fast (Fast Fourier Transform): 
 millions of digits ≪ 1sec. 2. n! in complexity O(n) by divide-and-conquer Notation: O(n) means n ! := n ⇥ · · · ⇥ d n/ 2 e ⇥ ( d n/ 2 e � 1) ⇥ · · · ⇥ 1 O(n log k n) for some k | {z } | {z } size O ( n log n ) size O ( n log n ) 3. Linear recurrence: convert into 1st order recurrence on vectors and apply the same idea. n 1 satisfies a 2nd order rec, computed via Ex: X e n := k ! ✓ e n ◆ ✓ n + 1 ◆ ✓ e n − 1 ◆ ✓ 1 ◆ k = 0 = 1 − 1 = 1 n ! A !( n ) 0 0 e n − 1 n e n − 2 . n | {z } A ( n ) Conclusion: Nth element in O(N) ops. 4/26

Numerical evaluation of solutions of LDEs N ∞ X X Principle : a n x n a n x n f ( x ) = + n = 0 n = N + 1 | {z } | {z } fast evaluation good bounds f solution of a LDE with coeffs in ℚ (x) 1. linear recurrence in N for the first sum (easy); 2. tight bounds on the tail (technical); 3. extend to ℂ by analytic continuation. arctan(1+i) Computation on integers. No roundoff errors. ˜ Conclusion: value anywhere with digits in ops. O ( N ) N Sage code available M. Mezzarobba 5/26 [Chudnovsky-Chudnovsky87;van der Hoeven99;Mezzarobba-S.10;Mezzarobba16]

A. Using Linear Differential Equations Exactly II. Local and Asymptotic Expansions

Dynamic Dictionary of 
 Mathematical Functions • User need • Recent algorithmic progress • Maths on the web http://ddmf.msr-inria.inria.fr/ [Benoit-Chyzak-Darrasse-Gerhold-Mezzarobba-S.2010] 6/26

A. Using Linear Differential Equations Exactly III. Proofs of Identities

Proof technique > series(sin(x)^2+cos(x)^2-1,x,4); O(x 4 ) f satisfies a LDE 
 ⟺ 
 f,f’,f’’,… live in a 
 Why is this a proof? finite-dim. vector space 1. sin and cos satisfy a 2nd order LDE: y’’+y=0; 2. their squares and their sum satisfy a 3rd order LDE; 3. the constant -1 satisfies y’=0; 4. thus sin 2 +cos 2 -1 satisfies a LDE of order at most 4; 5. the Cauchy-Lipschitz theorem concludes. Proofs of non-linear identities by linear algebra! 7/26


 
 Mehler’s identity for Hermite polynomials ⇣ ⌘ 4u ( xy − u ( x 2 + y 2 )) exp ∞ H n ( x ) H n ( y ) u n 1 − 4u 2 X n ! = √ 1 − 4u 2 n = 0 1. Definition of Hermite polynomials: 
 recurrence of order 2; 2. Product by linear algebra: H n+k (x)H n+k (y)/(n+k)!, k ∈ ℕ
 generated over (x,n) by 
 Q H n ( x ) H n ( y ) , H n + 1 ( x ) H n ( y ) , H n ( x ) H n + 1 ( y ) , H n + 1 ( x ) H n + 1 ( y ) n ! n ! n ! n ! → recurrence of order at most 4; 3. Translate into differential equation. 8/26

Guess & Prove Continued Fractions 1. Taylor expansion produces first terms (easy): x arctan x = 1 3 x 2 1 + 4 15 x 2 1 + 35 x 2 9 1 + 1 + · · · n 2 2. Guess a formula (easy): a n = 4n 2 − 1 3. Prove that the CF with these a n converges to arctan: ( x 2 + 1)( P n /Q n ) 0 − 1 H n := Q 2 show that � � = O ( x n ) n where P n /Q n is the nth convergent. gfun[ContFrac] Algo ≈ compute a LRE for H n and simplify it. No human intervention needed. 9/26

It Works! • This method has been applied to all explicit 
 C-fractions in Cuyt et alii, starting from either: a Riccati equation: 
 y 0 = A ( z ) + B ( z ) y + C ( z ) y 2 a q -Riccati equation: 
 y ( qz ) = A ( z ) + B ( z ) y ( z ) + C ( z ) y ( z ) y ( qz ) a difference Riccati equation: 
 y ( s + 1) = A ( s ) + B ( s ) y ( s ) + C ( s ) y ( s ) y ( s + 1) • It works in all cases, including Gauss’s CF, Heine’s q- analogue and Brouncker’s CF for Gamma. • In all cases, H n satisfies a recurrence of small order. In progress: 1. explain why this method works so well, 
 2. classify the formulas it yields. [Maulat-S. 15, 17 ] 10/26

B. Conversions (LDE → LDE)

From equations to operators S n ↔ (n ↦ n+1) D x ↔ d/dx 
 x ↔ mult by x n ↔ mult by n product ↔ composition product ↔ composition D x x=xD x +1 S n n=(n+1)S n Taylor morphism: D x ↦ (n+1)S n ; x ↦ S n-1 
 produces linear recurrence from LDE Ex. (erf): D 2 x + 2 xD x 7! ( n + 1) S n ( n + 1) S n + 2 S − 1 n ( n + 1) S n = ( n + 1)( n + 2) S 2 n + 2 n 11/26

Chebyshev expansions arctan Taylor Chebyshev Z 1 c k = 2 f ( x ) T k ( x ) 1 − x 2 dx √ π − 1 ✓ ◆ T 1 ( x ) T 3 ( x ) T 5 ( x ) z − 1 3z 3 + 1 √ 5z 5 + · · · 2 ( 2 + 1 ) 2 + 3 ) 2 + 2 + 3 ) 3 + · · · √ √ √ − ( 2 2 + 3 ) 3 ( 2 5 ( 2 12/26

Ore fractions Generalize commutative case: R=Q -1 P with P & Q operators. B -1 A=D -1 C when bA=dC with bB=dD=LCLM(B,D). Algorithms for sum and product: B -1 A+D -1 C=LCLM(B,D) -1 (bA+dC), with bB=dD=LCLM(B,D) B -1 AD -1 C=(aB) -1 dC, with aA=dD=LCLM(A,D). [Ore1933] 13/26

Application: Chebyshev expansions Chebyshev Taylor 
 2xT n (x)=T n+1 (x)+T n-1 (x) x n+1 =x · x n ↔ x ↦ X:=S -1 ↔ x ↦ X:=(S n +S n-1 )/2 
 (x n )’=nx n-1 ↔ d/dx ↦ D:=(n+1)S 2(1-x 2 )T n ’(x)=-nT n+1 (x)+nT n-1 (x) 
 ↔ d/dx ↦ D:=(1-X 2 ) -1 n(S n -S n-1 )/2. Prop. If y is a solution of L(x,d/dx), then its Chebyshev coefficients annihilate the numerator of L(X,D). > deqarctan:=(x^2+1)*diff(y(x),x)-1: > diffeqToGFSRec(deqarctan,y(x),u(n),functions=ChebyshevT(n,x)); nu ( n ) + 6( n + 2) u ( n + 2) + ( n + 4) u ( n + 4) Applications to Validated Numerical Approximation [Benoit-S.09;Benoit12;BenoitJoldesMezzarobba17] M. Joldes 14/26

C. Computing Linear Differential Equations (Efficiently)

C. Computing Linear Differential Equations (Efficiently) I. Algebraic Series and Questions of Size

Algebraic Series can be Computed Fast P irreducible P ( X, Y ( X )) = 0 Wanted : the first N Taylor coefficients of Y. Note: P x ( X, Y ( X )) + P y ( X, Y ( X )) · Y 0 ( X ) = 0 F sol LDE 
 ⇒ F(Y(X)) sol LDE Y 0 ( X ) = ( − P x P � 1 mod P )( X, Y ( X )) y (same argument) a polynomial Y ( X ) , Y 0 ( X ) , Y 00 ( X ) , . . . in Vect Q ( X ) (1 , Y, Y 2 , . . . ) finite dimension (deg P) → a LDE by linear algebra 15/26 [Abel1827;Cockle1861;Harley1862;Tannery1875]

Order-Degree Curve D=deg P The cost of minimality degree degree O(D^3) O(D^3) minimal LDE nice recurrence minimal recurrence degree O(D^2) O(D^2) O(D^2) O(D^2) O(D^2) O(D) O(D) O(D) O(D) O(D) O(D) O(D^2) O(D^2) order order O(D^2) O(D^2) O(D^3) order differential equations corresponding recurrences 16/26 [Bostan-Chyzak-Lecerf-S.-Schost07;Chen-Kauers12;Chen-Jaroschek-Kauers-Singer13]

C. Computing Linear Differential Equations (Efficiently) II. Creative Telescoping

Examples n n k ◆ 2 ✓ n + k ◆ 2 ◆ 3 ✓ n ✓ n ◆✓ n + k ◆ ✓ k X X X = k k k k j k = 0 k = 0 j = 0 ✓ j + k ◆✓ r ◆✓ n ◆✓ s + n − j − k ◆ ✓ n + r ◆✓ ◆ s − r X ( − 1 ) j + k = ( − 1 ) l k + l n + l j k m − j m − n − l j , k Z + ∞ xJ 1 ( ax ) I 1 ( ax ) Y 0 ( x ) K 0 ( x ) dx = − ln( 1 − a 4 ) 2 π a 2 0 ⇣ ⌘ I ( 1 + 2xy + 4y 2 ) exp 4x 2 y 2 dy = H n ( x ) 1 1 + 4y 2 3 b n / 2 c ! 2 π i y n + 1 ( 1 + 4y 2 ) 2 n n q k 2 ( − 1 ) k q ( 5k 2 − k ) / 2 X X = ( q ; q ) k ( q ; q ) n − k ( q ; q ) n − k ( q ; q ) n + k k = 0 k = − n Note: at least one 1. Prove them automatically Aims: free variable 2. Find the rhs given the lhs First : find a LDE (or LRE) 17/26

Creative telescoping Z X I ( x ) = f ( x , t ) dt =? U ( n ) = u ( n , k ) =? or k Input : equations (differential for f or recurrence for u ). Output : equations for the sum or the integral. ✓ n ◆ Ex.: U n := X k k ✓ n + 1 ◆ ✓ n ◆ ✓ n + 1 ◆ ✓ n + 1 ◆ ✓ ◆ ✓ n ◆ n X X U n +1 − 2 U n = − 2 = + − − k k k k + 1 k + 1 k k k | {z } | {z } telescopes telescopes Aim: find A(n,S n ) and B(n,k,S n ,S k ) such that Def : ∆ k: =S k -1. (A(n,S n )+ Δ k B(n,k,S n ,S k )) ⋅ u(n,k)=0, certificate then the sum telescopes, leading to A(n,S n ) ⋅ U(n)=0. Integrals : differentiate under the ∫ sign and integrate by parts. 18/26