DAETS: a Differential-Algebraic Equation Code in C++ for High Index - PowerPoint PPT Presentation

Overview of DAETS Numerical method Numerical results Event location Summary DAETS: a Differential-Algebraic Equation Code in C++ for High Index and High Accuracy Ned Nedialkov 1 and John Pryce 2 1 Dept. of Computing and Software, McMaster U., Canada nedialk@mcmaster.ca 2 Dept. of Information Systems, Cranfield U., UK j.d.pryce@ntlworld.com DART-IV 4th International Workshop on Differential Algebra and Related Topics October 27-30, 2010, Beijing, China

Overview of DAETS Numerical method Numerical results Event location Summary Outline Overview of DAETS Numerical method Numerical results Event location Summary 2/31

Overview of DAETS Numerical method Numerical results Event location Summary References N. Nedialkov and J. Pryce ◮ Solving differential-algebraic equations by Taylor series (I): computing Taylor coefficients BIT, 45(3), pp. 561–591 , 2005 ◮ (II): computing the System Jacobian BIT, 47(1), pp. 121–135, 2007 ◮ (III): the DAETS code JNAIAM 3 (2008), pp. 61–68 ◮ DAETS User Guide Tech. Report, Dept. Computing & Software, McMaster U. 2008–2009 The DAETS package is available at http://www.cas.mcmaster.ca/~nedialk/daets 3/31

Overview of DAETS Numerical method Numerical results Event location Summary DAETS solves DAEs by Taylor series expansion ◮ DAETS— Differential Algebraic Equations by Taylor Series ◮ Solves DAE initial value problems, for state variables x j ( t ) , j = 1 , . . . , n , of the general form f i ( t , the x j and derivatives of them ) = 0 , i = 1 , . . . , n ◮ Can be fully implicit ◮ d / d t can appear anywhere in the expressions for f i e.g. one of the equations could be � 1 sin t ) ′ � 2 ( x ′ + t 2 cos x 2 = 0 1 + ( x ′ 2 ) 2 4/31

Overview of DAETS Numerical method Numerical results Event location Summary DAETS solves high-index problems ◮ Index ν measures how “hard” DAE is to solve ◮ For traditional methods, ν ≥ 3 considered hard ◮ DAETS is based on structural analysis of DAE + Taylor series ◮ In principle unaffected by index ◮ Have solved artificial problems up to ν = 47 (Any physical sense? . . . is another matter) 5/31

Overview of DAETS Numerical method Numerical results Event location Summary DAETS is a new kind of DAE solver ◮ Excellent at high-index DAEs ◮ Excellent for getting high accuracy ◮ Returns useful data about structure of problem ◮ Doesn’t compete on speed at moderate accuracies ◮ . . . or on handling very large problems 6/31

Overview of DAETS Numerical method Numerical results Event location Summary Numerical method summary ◮ Start with code for the f i that define the DAE ◮ Use automatic differentiation (AD), FADBAD ++ package, = d r f i to evaluate suitable derivatives f ( r ) d t r at given t = t j i ◮ For each step: Equate these to zero “in batches” to get Taylor ◮ coefficients of (vector) solution x ( t ) at current ( t j , x j ) x j + 1 at Sum Taylor series to get approximation � ◮ t j + 1 = t j + h j x j + 1 on DAE’s constraints to get a Project this � ◮ consistent x j + 1 ◮ Repeat, to step along range in usual way 7/31

Overview of DAETS Numerical method Numerical results Event location Summary x j + 1 � x j + 1 x j 8/31

Overview of DAETS Numerical method Numerical results Event location Summary Numerical method, cont. ◮ Before all this, do Structural Analysis: preprocess the DAE code to find the 2 n integer offsets, one for each variable, one for each equation ◮ These prescribe the “batches” in the overall process of solving for TCs ◮ They imply the Initial Values (IVs) data is not a flat vector unlike with most DAE solvers ◮ Following example illustrates 9/31

Overview of DAETS Numerical method Numerical results Event location Summary The simple pendulum Index-3 system with equations 0 = f = x ′′ + x λ 0 = g = y ′′ + y λ − G G = gravity 0 = h = x 2 + y 2 − L 2 L = length of pendulum State variables x ( t ) , y ( t ) , λ ( t ) x y f g h Item λ Offset 2 2 0 0 0 2 10/31

Overview of DAETS Numerical method Numerical results Event location Summary User needs offsets to understand IVs ◮ Offsets tell what initial values (IVs) should be provided ◮ Offsets x y x x ′ λ 2 2 0 y y ′ mean that IVs comprise values for x , x ′ ; y , y ′ ◮ Except when DAETS sees DAE is non-linear x x ′ x ′′ in leading derivatives (here x ′′ , y ′′ , λ ) y y ′ y ′′ It requires an extra set of derivatives E.g. if first equation were λ 0 = ( x ′′ ) 3 + x λ then IVs must comprise x , x ′ , x ′′ ; y , y ′ , y ′′ ; λ ◮ Reason: to assure local uniqueness of solution 11/31

Overview of DAETS Numerical method Numerical results Event location Summary User needs offsets to understand constraints Offsets tell what constraints the provided IVs must meet for consistency f g h Offsets 0 0 2 mean they must satisfy in the linear case in the non-linear case h , h ′ = 0: f ; g ; h , h ′ , h ′′ = 0, so in our example, add these: 0 = h = x 2 + y 2 − L 2 0 = f = ( x ′′ ) 3 + x λ 0 = h ′ = 2 xx ′ + 2 yy ′ 0 = g = y ′′ + y λ − G 0 = h ′′ = 2 ( xx ′′ + ( x ′ ) 2 + yy ′′ + ( y ′ ) 2 ) 12/31

Overview of DAETS Numerical method Numerical results Event location Summary Finding consistent point ◮ Solution x ( t ) must satisfy algebraic constraints for all t to be consistent ◮ Constraint can be obvious, as (for Pendulum) h = 0, or hidden, as h ′ = 0 ◮ Finding initial consistent point can be hardest part of solving DAE ◮ Not built in to most solvers ◮ Fits naturally into DAETS workflow We formulate it as a nonlinear minimization problem and give it to IPOPT package 13/31

Overview of DAETS Numerical method Numerical results Event location Summary Numerical results ◮ Accuracy ◮ Efficiency ◮ DAETS on a High-index problem ◮ DAETS on a Continuation problem 14/31

Overview of DAETS Numerical method Numerical results Event location Summary Accuracy vs. tolerance ◮ We plot Significant Correct Digits: SCD = − log 10 � componentwise rel. error at end of integration � as a function of tolerance ◮ Problem is Transistor Amplifier, index-1 DAE of size n = 8, from Test Set for Initial Value Problem Solvers , Bari University, Italy ◮ “DAETS”, “RADAU”, “DASSL ” curves compare with reference solution in Test Set documentation ◮ “DAETS-2” curve uses reference solution computed by DAETS with tol = 10 − 16 15/31

Overview of DAETS Numerical method Numerical results Event location Summary Plots of accuracy vs. tolerance Transistor amplifier 14 DAETS DAETS-2 12 RADAU DASSL 10 SCD 8 6 4 2 10 -4 10 -6 10 -8 10 -10 10 -12 10 -14 tol 16/31

Overview of DAETS Numerical method Numerical results Event location Summary Efficiency experiments ◮ Plot CPU time vs. SCD ◮ Problems are (from the ODE/DAE test set) index-3 DAE, n = 10 ◮ Car axis index-1 DAE, n = 8 ◮ Transistor amplifier index-1 DAE, n = 6 ◮ Chemical Akzo Nobel ODE, n = 8 ◮ HIRES ◮ Produce work-precision diagrams for DAETS, DASSL, and RADAU (as described in the test set) ◮ Compare with DASSL and RADAU solvers 17/31

Overview of DAETS Numerical method Numerical results Event location Summary Work-precision diagrams Car axis Transistor amplifier 10 0 10 2 DAETS DASSL RADAU 10 1 CPU time CPU time 10 -1 10 0 10 -2 10 -1 10 -3 10 -2 0 2 4 6 8 10 12 14 2 4 6 8 10 12 14 SCD SCD Chemical Akzo HIRES 10 -1 10 1 10 0 CPU time CPU time 10 -2 10 -1 10 -2 10 -3 10 -3 0 2 4 6 8 10 12 14 16 -2 0 2 4 6 8 10 12 SCD SCD 18/31

Overview of DAETS Numerical method Numerical results Event location Summary “Multi-pendula” — a class of high-index problems ◮ System is a “chain” of P simple pendula with coupling ◮ Pendulum 1 is as normal ◮ Length of pendulum p depends on λ p − 1 , for p = 2 , . . . , P ◮ For P = 2 0 = x ′′ 1 + λ 1 x 1 0 = x ′′ 2 + λ 2 x 2 0 = y ′′ 1 + λ 1 y 1 − G 0 = y ′′ 2 + λ 2 y 2 − G and 0 = x 2 1 + y 2 1 − L 2 0 = x 2 2 + y 2 2 − ( L + c λ 1 ) 2 where c is a constant ◮ Chain of length P has size n = 3 P and index 2 P + 1 ◮ DAETS has solved systems for P up to 23 giving index 47 19/31

Overview of DAETS Numerical method Numerical results Event location Summary P = 7, index 15, x 7 ( t ) plot Seven pendula, tol = 10 -9 , x 7 with two slightly differing ICs 5 4 3 2 1 0 -1 -2 -3 -4 -5 0 10 20 30 40 50 60 20/31

Overview of DAETS Numerical method Numerical results Event location Summary Continuation problems ◮ No need for derivatives to be present Can solve n purely algebraic equations f ( λ, x ) = 0 to find x = ( x 1 , . . . , x n ) as a function of λ ◮ To handle turning points, use arc-length continuation ◮ Define Euclidean arc-length s by ( d λ/ d s ) 2 + ( d x 1 / d s ) 2 + . . . + ( d x n / d s ) 2 = 1 ◮ Find ( λ, x ) as function of s ◮ Gives an index-1 DAE of size n + 1 21/31