Home

Daets User Guide - Computing and Software

1. S ALESSANDRINI A motivational example for the numerical solution of two point boundary value problems SIAM Review 37 1995 pp 423 427 K BRENAN S CAMPBELL AND L PETZOLD Numerical Solution of Initial Value Problems in Differential Algebraic Equations SIAM Philadelphia sec ond ed 1996 S L CAMPBELL AND E GRIEPENTROG Solvability of general differential alge braic equations SIAM J Sci Comput 16 1995 pp 257 270 Cpp Unit project http sourceforge net projects cppunit Gcov a test coverage program nttp gcc gnu org onlinedocs gcc Gcov html A GRIEWANK Evaluating Derivatives Principles and Techniques of Algorithmic Dif ferentiation Frontiers in applied mathematics SIAM Philadelphia PA 2000 E HAIRER AND G WANNER Solving Ordinary Differential Equations II Stiff and Differential Algebraic Problems Springer Verlag Berlin 1991 INTEL CORP Intel C compiler for Linux web page http support intel com support performancetools c linu R JONKER AND A VOLGENANT A shortest augmenting path algorithm for dense and sparse linear assignment problems Computing 38 1987 pp 325 340 LAPACK PROJECT BLAS Basic Linear Algebra Subprograms http www netlib org blas LAPACK PROJECT LAPACK Linear Algebra PACKage http www netlib org lapack F MAzzIA AND F IAVERNARO Test set for initial value problem solvers Tech Rep 40 Department of Mathematics University of Bari
2. 2 4 DAETS default settings The table below lists the default values of all parameters that control the operation of the integration process or are otherwise relevant to the user 2 4 DAETS default settings tolerance accuracy control order min order max length of TS max order stepsize min stepsize max stepsize Summary 36 1078 mixture of relative and absolute selected automatically by DAETS according to 2 3 or set by the user 1 80 changeable by recompiling the source determined by variable determined by DAETS determined by DAETS according to 2 4 the largest finite double precision number All the user callable features of DAETS have been covered in this chapter and Chap ter l Many but by no means all of them are illustrated in the example programs of the next chapter Chapter 3 Examples of DAETS Code We begin in 83 1 with suggestions for how to handle a DAE whose structure is un familiar to you In 93 2 we show code to solve a small DAE from the Test Set 12 illustrating a technique to simplify the translation from Fortran where arrays are usu ally indexed from 1 into C which indexes them from 0 In 93 3 we solve Van Der Pol s ODE showing how parameters are passed to a DAE function and how one step mode may be coded B 4 shows how DAETS can solve continuation problems In 43 35 we give code to solve a model of putting a golf ball originally approximated in by an ODE We show
3. n one gets n equations that in principle by the Implicit Function Theorem should be solvable for the d th derivative of x j 1 n in terms of lower derivatives thus converting the DAE to an ODE This involves an n x n matrix the System Jacobian J Broadly the numerical method succeeds if this is found to be nonsingular The d th derivatives are the leading derivatives of the DAE If they occur in a jointly linear way in the original undifferentiated f then the DAE is called quasi linear Given the same offsets a quasilinear problem needs fewer IVs than one that is not quasilinear Counter intuitively a leading derivative need not occur anywhere in the DAE see FAQ To recognize quasilinearity ask Can the DAE be written in matrix form Ax b 0 where x is the column vector of leading derivatives and the matrix A and vector b are independent of x Example For the pendulum system L2 the offsets of x y A are 2 2 0 so the leading derivatives are x y and A In this example we have f 10 x g 0 g 0 1 y y G h 00 0 x y L so the answer is Yes and the DAE is quasilinear If say the first equation is replaced by 0 f 2 x or 0 f x y x the system ceases to be quasilinear The rule for specifying IVs can now be stated For the reason behind it see 1 e If the DAE is quasilinear IVs for each x and its derivatives up to but not including the d th must be gi
4. 0 taking G 9 8 and L 10 and displays solution values at the end of the integration The interface to DAETS is defined by the header file DAEsolver h code line B The source code of DAETS is in namespace daets hence the line B Lines HLI are the function that defines the DAE x y Its name is fcn here but this is arbitrary as are the names of the arguments t z f param Input argu ment t represents the independent variable t Argu ment z is an array of a templated type T representing Y the state variable vector z x y The code com putes the corresponding values of the functions f g h Figure 1 1 Simple pendulum in L2 and puts them in the array f which also has type T The type T is interpreted as several actual types at compile time one of them is the basis for the numerical computation with Taylor series and the others support symbolic manipulation of the fcn code during the preliminary structural analysis process Arrays are indexed the C way from zero so z 0 means x f 0 means f etc The operation Diff q means d dt so Diff z 0 2 means 2 i e d x dt Argument param is explained in Its use is illustrated in the Van Der Pol example of 3 3 and the continuation example of BA The main program is in lines BHA Lines T5HI6 construct the two key objects for problem solution Solver an object of the DAEsolver class which understands the DAE and x of the DAEsolution class which may be
5. 85 86 87 88 89 90 91 92 93 94 95 3 2 Chemical Akzo Nobel problem 40 gt fcn evaluates the Akzo Nobel DAE functions template lt typename T gt void fcn T t const T y T f void param fcn0 t y f for int i 0 i lt 5 i f i Diff y il 1 FORWARD FUNCTION DECLARATIONS double compSignificantCorrectDigits DAEsolution amp x void setInitialValues DAEsolution amp x A eee VAIN PROGRAM gt int main int n 6 double tend 180 0 Create solver and solution point and set initial values DAEsolver Solver n DAE_FCN fcn DAEsolution x Solver setInitialValues x Integrate SolverExitFlag flag Solver integrate x tend flag if flag success printSolverExitFlag flag return 1 Print solution and statistics about the integration x printSolution x printStats How did it compare with reference solution print Significant correct digits 5 1f n compSignificantCorrectDigits x return 0 SET INITIAL VALUES void setInitialValues DAEsolution amp x double t 0 x setX 0 0 0 444 setX 0 1 0 setX 1 0 0 00123 setX 1 1 0 setX 2 0 0 0 setX 2 1 0 setX 3 0 0 007 setX 3 1 0 setX 4 0 0 0 setX 4 1 0 setX 5 0 0 setT t gt Merci coninipiiceccaitces CHECK ACCURACY Check solution against reference solution from
6. Ti g t Ly 22 l Ly T3 where g t is an arbitrary given function The offsets are d 2 1 0 and c 2 1 0 Thus x is a leading derivative though it does not occur in the system By adding extra equations 75 t4 one can make an arbitrarily large discrepancy between the leading derivative order and the order that occurs in the equations 5 2 7 Not many standard functions are allowed in fcn Are more planned These have sufficed for the examples we have done so far If you need more let us know 5 2 Frequently asked questions 67 5 2 8 5 2 9 Have you plans for an event location facility Yes and efficient step by step output or plotting without using one step mode and banded sparse linear algebra and other features that modern ODE solvers have Please show us some real applications and say what your needs are so we can set our priorities Can DAETS handle diagnostics silently That is can it be put in a mode where it handles exceptional conditions not by printing a diagnostic message but by returning a diagnostic data structure so that the calling program can take appropriate action This allows the code to be used as a numerical engine at the heart of some application This will be in the next version The exception handling mechanism of C provides most of the needed features 5 2 10 Have you tips for debugging code that uses DAETS One of us JDP is a C beginner He has found that
7. lation of these That is it succeeds in exact arithmetic ill conditioning of the System Jacobian may cause it to fail in finite precision 6 An upper bound for the standard differentiation index va is given by the Taylor index 1 if some d is zero 0 otherwise 5 6 Vp max c 2 For the commonest kinds of DAE vr vg but the difference can be arbitrarily large see FAQ 5 215 5 2 Frequently asked questions 64 7 It was explained in 41 2 2Jthat required IVs which are in fact the values carried along from step to step throughout the integration process comprise the vector di a f d2 a dn a eee eo ee L 5 2 Bn D n 5 7 where a 1 if the DAE is quasilinear 0 otherwise and for consistency X must satisfy the set of equations N 9 qu 7 dr 5 8 In the quasilinear case an x whose d is zero it must have o 0 for all i and thus is a purely algebraic variable does not appear in the vector X Similarly an f with c 0 does not appear in F fi means df dt treating the variables and their derivatives as unknown func tions of t for instance if f x x1 3 then fi x 2 13 2 23 similarly for higher derivatives The reason for the extra components in the non quasilinear case is as follows Consider a particular independent variable value t If a set of values X in 57 is consistent with some solution of the DAE at t then in the linear cas
8. 3 2 Even for n 10 this is considered quite difficult We use the approach of seeking a fixed point of the parameterized problem x Ag x that is a root of x 0 3 3 where f A x x Ag x When 0 equation 3 3 has the trivial solution x 0 and we hope to track a solution all the way to the desired root at 1 DAETS can handle this directly taking as the independent variable and solving for x x A This formulation gives an index 0 system all the offsets c and d are zero and the System Jacobian is J L 28 Ox Ox For general problems of the form 3 3 this approach works as long as fx Of Ox is nonsingular but typically it has a serious weakness one reaches values of where f becomes singular while the n x n 1 matrix f f retains full rank n These are turning points with respect to A powerful alternative is to treat all of x1 n A as on the same footing instead of viewing as special and to use arc length continuation The system 8 3 may be rewritten as n equations in n 1 unknowns f y 0 3 4 where y x A Generically fy is always of full row rank and there is a unique except for direction solution path through any point satisfying 3 4 tangential to the 1 dimensional null space of fy Invent a new independent variable s and add to BT an equation dy ds 1 that is gt y I 3 5 where means d ds This defines s to be Eu
9. Italy 2003 http pitagora dm uniba it testset N S NEDIALKOV Computing Rigorous Bounds on the Solution of an Initial Value Problem for an Ordinary Differential Equation PhD thesis Department of Computer Science University of Toronto Toronto Canada M5S 3G4 February 1999 N S NEDIALKOV AND J D PRYCE Solving differential algebraic equations by Taylor series I Computing Taylor coefficients BIT 45 2005 pp 561 591 73 BIBLIOGRAPHY 74 15 16 I7 18 19 20 21 22 23 24 Solving differential algebraic equations by Taylor series II Computing the Sys tem Jacobian BIT 47 2007 pp 121 135 Solving differential algebraic equations by Taylor series III The DAETS code J Numerical Analysis Industrial and Applied Mathematics 3 2008 pp 61 68 C C PANTELIDES The consistent initialization of differential algebraic systems SIAM J Sci Stat Comput 9 1988 pp 213 231 PathScale compiler suite http www pathscale co J D PRYCE A simple structural analysis method for DAEs BIT 41 2001 pp 364 394 G REISSIG W S MARTINSON AND P I BARTON Differential algebraic equations of index 1 may have an arbitrarily high structural index SIAM J Sci Comput 21 2000 pp 1987 1990 O STAUNING AND C BENDTSEN FADBAD web page FADBAD is available at A WACHTER An Interior Point Algorithm for Large Scale Nonlinear Optimization wit
10. See P T 6land 93 3 for details Result This constructor does much work behind the scenes It performs the struc tural analysis of the DAE and stores the offsets and determines whether the system is quasilinear Thus the constructed object knows the scheme for gener ating the Taylor series including what set of x and their derivatives comprise valid IVs and the equations they must satisfy for consistency Note For convenience the macro DAE_FCN is introduced to simplify a call to this constructor as e g DAEsolver Solver n DAE_FCN fcn and the user is recommended to call it in this way Ill posed DAEs This constructor recognizes when a DAE is ill posed see 45 1 for a definition If such a DAE is encoded a solver object is still created but it is in an illposed state This can be detected by a call to the isI11Posed method see 2 1 2 1 3 Basic classes and methods 10 integrate method void integrate DAEsolution amp x double tend SolverExitFlag amp flag x Input contains the initial t and the initial values of the solution Output contains the reached t and the computed values of the solution at t tend Input the desired final value of t If tend lt t integration proceeds in the negative direction If tend t the code uses the supplied x to compute a consistent initial point if x is not already consistent and returns this in x leaving t unchanged flag Output reports the success of failu
11. let DAETS do the analysis as suggested in 8 1 4 Form the n x n System Jacobian matrix Cn o i H oe Z A RN equivalently J or a C di qn J i eee 0 otherwise 5 4 Here the x and derivatives thereof are treated as unrelated independent vari ables within f For instance if f is 1 111 then af 8xP 0Of1 0x x x5 If J thus defined is not identically singular the SA based method almost cer tainly succeeds To be precise if there is a consistent point as defined below where J is nonsingular a solution to the DAE exists through that point it is lo cally analytic and the DAETS algorithm in exact arithmetic can find its Taylor series to arbitrary order For the pendulum 5 4 gives the system Jacobian afd 0 OFf OX Tr 0 lt J 0 0g 0y 0g 9A 0 1 yl 5 5 dh x Oh dy 0 es 29 This is not identically singular since det J 2 x y Indeed at a consistent point from the third equation of T 2 det J 2L 0 5 The SA method sometimes fails see FAQ B2B However it always succeeds on various common classes of system These include index 0 DAEs implicit ODEs semi explicit index 1 and index 2 DAEs DAEs with Hessenberg structure of arbitrary index of which a particular case is index 3 constrained mechanical systems Mx GTA t g x 0 where G is the Jacobian of g purely algebraic systems f t x 0 continuation problems and the arc length formu
12. T amp y const T amp z void param T amp f T amp dfx T amp dfy T amp dfz Compute both f and grad f where f x y z 0 defines the surface of the green f must be defined so that df dz is positive Input x y z Output E holds x y z s dfx dfy dfz holds grad f This green has the form z f1 2 0 where f1 defines a plane and 2 defines a Gaussian hump centered at x y a b of height h and width squared v double p double param double dgx p 0 dgy plil a p 2 b p 3 h p 4 v p 5 T f1 dgx x dgy y planar part T f2 h exp sqr x a sqr y b 2 v hump part f z f1 2 dfx dgx x a v f2 dfy dgy y b v 2 dfz 1 gt Pin in ti FCN DEFINES EQNS DE MOTION template lt typename T gt void fcn T t const T w T f void param state variables define R w O normal reaction force of green on ball define x w 1 define y w 2 co ordinates of ball on surface of green define z w 3 z points upward const double constants from problem definition G 9 81 gravity in SI units mu 0 04 coefficient of friction between green and ball velocity T vx Diff x 1 T vy Diff y 1 T vz Diff z 1 T v norm2 vx vy vz vector normal to green T fval dfx dfy dfz dfgreen x y z param fval dfx dfy dfz T normdf norm2 dfx dfy dfz f
13. amp mu i x setFirstEntry x updatePoint xp setT t0 Create file name and open a file for writing char buf 20 sprintf buf vdp 1e out mulil FILE fd fopen buf w Restrict the stepsize so the plots appear smooth Solver setHmax hmax i Integrate and write t x x at each step into file vdp mu out while x getT tend i amp amp flag success Solver integrate x tend i flag fprintf fd 3e 6e 6e n x getT x getX 0 0 x getX 0 1 fclose fd if flag success printSolverExitFlag flag return 1 return 0 When passing parameter s to fcn the DAEsolver constructor must be invoked with a third parameter that is the address of a variable line PI that will be passed to fcn Since this parameter is of type void the user must extract the needed value s 3 3 One step mode Passing data to a problem 43 inside fcn line Here inside fcn the data at param is converted to double the type used in this function If integrations with different parameter values are desired these values should be passed by calling the passDataToDAE method line B7 Here we pass in a loop the values 10 and 100 for u WARNING a The DAEsolver constructor and passDataToDAE must be called with the same type of their corresponding parameters b The size n of a problem is fixed at the time a DAEsolver object Solver is con structed You can not use passDataToDAE to change n
14. appear as constants It may help to think of the rod as being ever so slightly elastic Thus at any instant it stretches or compresses infinitesimally to produce that precise force T on the bob which keeps it going in an exact circle Were we to model the dynamics of this elasticity the equations of motion would be an ODE The DAE is the limiting and simpler case where the elasticity goes to zero For a purely mathematical illustration why DAEs are not ODEs see the example DAE in FAQ 5 216 It has no degrees of freedom its unique solution is 7 g t tg g t 13 9 t But if the zero on the left of the first equation is replaced by ea for a nonzero e it is a normal but unpleasant when e is small ODE with three degrees of freedom requiring three IVs to specify a unique solution The pendulum also has a well known description by an ODE d 9 G de T sin 0 where 0 is as in FAQ 5 2012 Many more complex mechanical systems such as robotic arms can be modeled by ODEs typically involving angles of rotation Why use a DAE when an ODE is available One good reason is that as the size of the system increases it becomes increasingly harder to formulate such ODE models whereas the algorithm to formulate the DAE form typically remains straightforward This is true in other application areas such as electrical circuit modeling For simulation packages that let one define a system interactively and that then set up its governing
15. by the Diff function For instance Diff x j 1 2 represents x that is d x dt Notes The code should not contain any branches that is if or switch statements Functions that are defined using these are not covered by the current theory The code may freely declare other variables These will normally be of type int or double or T Any intermediate variable that depends directly or indirectly on the inputs x or t must be of type T 1 3 2 DAEsolver class The constructor and the integrate method are described here They are the only DAEsolver methods used by the program in L2 1 A DAEsolver object has many attributes describing the DAE system the fcn code the offsets etc and holding strategy parameters of the integration accuracy tolerance Taylor series order max imum allowed stepsize etc For the methods related to these see 2 1 1 3 Basic classes and methods 9 Constructor DAEsolver int n T1 fcn1 T2 fcn2 T3 fcn3 T4 fcn4 void param NULL Input the value n as above fcn1 fcn4 Input All four must be the same namely the function fcn defining the DAE as in 1 3 1 with the signature lt typename T gt void fcn T t const T x T f void param The repetition is due to C language limitations See the Note below The types T1 T2 T3 T4 do not concern the user param Input optional pointer to a user defined object of arbitrary type that passes parameters to the definition of the DAE
16. does not change The choice is a matter of convenience 16 2 1 Controlling integration the DAEsolver class 17 2 1 1 Access to structural analysis data The first three methods use the vector class of the C standard library getSigmaMatrix method void getSigmaMatrix vector lt vector lt int gt gt amp s int neginfval 1 throw std logic_error Output A matrix which must be nxn If not an exception std logic_error results neginfval Input Optional the value that encodes oo 1 because int does not have an infinity value and the alternative of using a large negative number looks ugly in printed output If neginfval gt 0 an exception std logic_error results Result Stores in s the signature matrix X 0 of the DAE defined in equation ED getCVector method void getCVector vector lt int gt amp c const c Output A vector which must be of length n If not an exception std logic_error results Result Stores in c the vector of offsets c of the functions f getDVector method void getDVector vector lt int gt amp d const d Output A vector which must be of length n If not an exception std logic_error results Result Stores in d the vector of offsets d of the variables x 2 1 Controlling integration the DAEsolver class 18 getNumDegsOfFreedom method int getNumDegsOfFreedom const Returns The number of degrees of freedom of the DAE g
17. equations automatically this gives the DAE form a decided advantage Appendix A Revision History In this appendix we report changes in the DAETS code since version 1 0 of June 2008 A 1 Version 1 1 May 2009 This version has been extensively tested using the CPPUNIT 4 testing framework and the code has been analyzed with the Gcov 5 test coverage tool This work is reported in BA As a result we have fixed several minor bugs in the code and inaccuracies in the user guide Changes in functionality compared to version 1 0 are as follows 1 The getX method of the DAEsolution class can now return derivatives of an x up to p d that is x for any k 0 p dj where p is the order of the Taylor series see p 27 2 We have added a mechanism for data output inside the integrate function This allows for example one to generate a file with solution data suitable for plotting See p BO and 38 3 The constant in the formula for computing Amin 2 4 is changed from 10 to 16 70 Index Enumerated types ErrorEstType P1 SolutionState VarType 23 SolverExitFlag Functions printSolverExitFlag 13 Macros DAE_FCN El Methods getATol getCPUtime getCVector 17 getDVector 17 getErrorEstType getHmax getHmin getMaxOrder 23 getMinOrder 23 getNumAccSteps getNumDegsOfFreedom 18 getNumDerivatives 33 getNumRejSteps B1 getNumVariables 33 getOrder 23 getRTol getSigmaM
18. for an existing Solver If you attempt this a segmentation fault or undefined behavior is likely to occur An fcn can allow different values of n as in the continuation example of B4 but to use several values of n in one program run you must create a different solver for each n We store the initial values x x in a DAEpoint object lines 23 and 24 Then we create a DAEsolution object and activate the one step mode in line 23 In the for loop before x is reinitialized setFirstEntry must be called line B8 to unlock x for further use Since DAETS does not provide continuous output yet which can be used for producing smooth plots we restrict the largest stepsize the solver takes line 26 The values in hmax are found by trial and error The plots in Figure B 3Jof x versus t and x versus x for u 1 10 and 100 show the behavior for which this ODE is famous From the initial point the path moves towards a limit cycle to which it rapidly converges and which shows abrupt changes evidenced by the corners near 2 0 The larger is u the sharper are these corners The GNUPLOT file used to produce these plots is file vdp gp in daets examples gnuplot This file plots x vs t and x vs x for the Van Der Pol s equation set terminal postscript enh color eps dl 2 Courier 28 mu 1 set title Symbol m 1 font Times 28 set output vdp0x eps set xlabel t font 28 set ylabel y font
19. getType int index int order const throw std out_of_range index order Input such that x getType j 1 k returns the type of the entry of x repre senting ae the kth derivative of the jth variable J Constraint j k must be in the index set J of 2 5 Otherwise an exception std out_of_range results k Returns x getType j 1 k returns the type of the entry of x representing x getX method double getX int index int order const throw std logic_error std out_of_range index order Input values j and k such that x getX j 1 k returns the current value of the entry of x representing 4 the kth derivative of the jth variable Constraint If x is in Initial state j k must be in the index set J of 2 5 Otherwise index 0 n 1 i e j 1 n and k 0 p dj where p is the order of the Taylor series Returns x getX j 1 k returns the current value of the entry of x representing k DAA j An exception std logic_erorr results if e x is in an Initial state and entry j 1 k is not initialized by a call to setX or e x is in an Initial state and j 1 k is not in the index set J For example if x has not been integrated this exception arises when attempting to access an entry with indices outside J An exception std out_of_range results if variable index j 1 n or derivative index k 0 p dj 2 2 Sol
20. left by the relevant offset relative to the k 0 column The process starts at stage k 2 because the largest d is 2 At stage k 2 one must find x and y that satisfy h x y 0 This is one equation for two unknowns so there is 2 1 1 degree of freedom DOF at this stage which is how many of x and y one is allowed to specify as fixed the Fixable row in the table At stage k 1 one must find 2 and y that satisfy h x y 0 Actually h involves x and y also but these have been solved for and are no longer unknowns Again there is one equation for two unknowns so again 2 1 1 DOF which is how many of x and y one is allowed to specify as fixed At stage k 0 one has three equations for three unknowns x y x y also occur in the equations but are no longer unknowns so there are no DOF and one cannot specify any of the unknowns 2 y of this stage as fixed In fact it is clear whether the DAE is quasilinear or not that stage k 0 always has n equations for n unknowns so none of its unknowns can be fixed The general rule is that only stages k lt 0 are relevant to fixing the number of equations at stage k is mj the number of i for which k c gt 0 the number of unknowns at stage k is n the number of j for which k d gt 0 and the number of unknowns that can be fixed at stage k d lt k lt 0 is the difference nj Mp Consider the pendulum exam
21. little later than they might have been causing slight inefficiency Or it is large 5 2 Frequently asked questions 66 and in exact arithmetic J is structurally singular In the latter case in finite precision J will be structurally singular within roundoff a situation that is easily detected numerically 5 2 5 What about Rei ig s example Reifig et al 20 give a family of DAEs for which the structural index can be an arbitrary v gt 1 though the differentiation index is 1 For v 3 the system is t2 3 101 alt a Lo T3 Ta b C Ta T5 T3 t La Ts Ta d 111 Le 11 els i and J gt d 11 a d 1 One easily sees that the SA method succeeds with c 0 0 1 1 2 d 0 1 1 2 2 and gives vr 3 But by inspection one can write the general so lution explicitly in a form involving no derivatives of the driving terms sufficient proof that the index is 1 In fact using two new variables xg 12 x3 and 17 4 x5 to eliminate the common subexpressions t t3 and t4 45 one obtains a system on which SA succeeds It has all c 0 and all d 0 except d d7 1 This raises a different issue from either FAQB2B or FAQ 6 2M but the cure by defining new variables can be seen as an illustration of the former 5 2 6 How can a leading derivative not actually occur in the DAE Consider the classic example of Linda Petzold 0
22. mode is switched off Note The operation of this method may be sensitive to your terminal settings Please notify the authors in case of problems 2 2 Solution path properties the DAEsolution class 30 setOutputFunction method This method sets a user defined function say outfcn for data output inside the integrate method for instance to write solution values to a file for plotting For an example see 3 01 The outfcn function is called at the end of each integration step Such a function should be declared as void outfcn const DAEsolution amp x void out_param void problem_param x Input reference to a solution object out_param Input pointer to a user defined object of arbitrary type that passes output related parameters e g file handles to this function problem_param Input pointer to a user defined object of arbitrary type that passes problem related parameters to this function To make integrate use outfcn call on a DAEsolution object void setOutputFunction OutputFcn outfcn void out_param void problem_param NULL outfcn Input pointer to an output function out_param problem_param Input pointers to user defined objects of arbitrary type see above Here OutputFcn is defined as typedef void DutputFcn const DAEsolution amp void void 2 2 5 Statistics The following statistics are cumulative starting at zero for an object x in the Initial state and added to by each s
23. of a robot arm that is index 5 but seems to be index 3 The Ring Modulator in 12 in case C 0 is an index 2 electrical circuit DAE that seems to be index 1 In our experience this has always been curable by purely algebraic manipulation of the equations to make them more sparse For instance in the robot arm DAE one can make SA find the correct index by naming one common subex pression in the equations as a new algebraic variable In the ring modulator DAE suitably adding to one equation a linear combination of other equations moves the sparsity around to give the desired effect We do not at present know any systematic way to do such things or any theory about when it can be done The signature matrix is sometimes wrong Can this spoil the com puted solution The point here is that the order of derivative of each variable in each equation is found by a pass through the code list treating it symbolically but with out any clever algebraic rearrangement So for instance if equation 1 is 123 tx tr 0 the code does not notice that the x terms cancel and sets 012 to 2 instead of the correct value 1 Such an error can overestimate but not underestimate any gij It is proved in 14 Theorems 5 1 and 5 2 that this cannot deceive the method Just one of two things happens Either the overestimation is small and the method computes the correct Taylor series but with some TCs possibly being computed a
24. related feature is provided by setPrintProgress of 2 2 4 which puts x into a mode where the current t step size etc are displayed at each step with an optional pause to help one read them This can help debug an integration that seems to get stuck without the need to program a one step loop in the main program Apart from the operations in this section DAEsolution objects support the arith metic operations of their base class DAEpoint which just holds an irregular array X Such operations strip out the non DAEpoint attributes For instance let x1 and x2 be two computed solutions at teng of the same DAE We wish to find how different they are by comparing their X arrays To do so form the difference x1 x2 which is a DAEpoint object and then x1 x2 norm which returns the max norm of the difference x1 x2 norm max a ol gt j k EJ in an obvious notation where J is the index set in 2 5 x records statistics such as the number of accepted and rejected steps used so far in the integration It also holds the attribute that determines whether the integration is done in one step mode 2 2 1 Access getT method double getT const throw std logic_error Returns x getT returns the current t stored in x Constraint t must be initialized by a call to setT Otherwise an exception std logic_error results 2 2 Solution path properties the DAEsolution class 27 getType method VarType
25. tolerance and accuracy achieved cannot be 5 2 Frequently asked questions 65 5 2 2 9 2 3 5 2 4 guaranteed The user is strongly recommended to call DAETS with more than one value for tol and to compare the results obtained to estimate their accuracy For instance 1e 7 and 1e 8 if 7 digits of accuracy are wanted If the method underlying DAETS sometimes fails can the code be trusted Yes If a nonsingular system Jacobian J is found at a consistent point in exact arithmetic this constitutes a proof that the DAE is locally solvable and the DAETS method will work see 19 Theorem 4 2 and 14 Theorem 4 2 Then in finite precision things can only go wrong because J is ill conditioned in which case we suspect any DAE method given the same equations and the same precision would have trouble What remedies are there when SA fails The SA method relies on the fact that DAEs from applications are typically sparse i e just a few variables and derivatives occur in each equation If the sparsity is destroyed the signature matrix may not reflect the underlying math ematical structure of the DAE and SA gives the wrong result For a crude example if the pendulum system f O is converted to the equivalent Mf 0 where M is a random nonsingular constant 3 x 3 matrix all three rows of X become 2 2 0 and SA gives the wrong result In interesting cases the failure of SA is more subtle Campbell and Griepentrog 3 give the DAE
26. 0 fval f 1 Diff x 2 R dfx normdf mu R vx v f 2 Diff y 2 R dfy normdf mu R vy v f 3 Diff z 2 R dfz normdf mu R vz v G undef R undef x undef y undef z gt MEE SS IE AUXILIARY FUNCTIONS 2 25232 2 3 gt template lt typename T gt T norm2 const T x const T y const T z return sqrt sqr x sqr y saqr z gt void outputfcn FILE fd const DAEsolution amp x 4 fprintf fd 4e 6e 6e 6e 6e 1 2 6e 6e n 3 5 A golf putting problem 51 x getT x getX 1 0 x getX 2 0 x getX 3 0 x getX 1 1 x getX 2 1 x getX 3 1 MAIN PROGRAM int main int argc char argv Check correct no of inputs on command line if argc 4 cout lt lt n Usage lt lt string argv 0 lt lt u0 vO dt n lt lt di where u0 vO are initial x and y velocities An lt lt dt is time interval for output n n exit 1 gt const int n 4 size of the problem Use atof to convert real number literals to double const double t0 0 0 initial time vx0 atof argv 1 initial x velocity vyO atof argv 2 initial y velocity dt atof argv 3 at intervals of dt output to file double greenparams 0 0 035 plane slightly uphill in y direction 3 1 1 center of hump 0 1 0 25 height and width of hump DAEsolver Solver n DAE_FCN
27. 0060ct21 0THERS blas In Ipopt Fortran_20060ct21 type configure prefix DIRECTORY where DIRECTORY is the directory in which IPOPT is to be installed For example if you wish to install it in subdirectory ipopt of your home directory you can type configure prefix HOME ipopt In Ipopt Fortran_20060ct21 type make install This should install IPOPT in DIRECTORY For more information about the configure options type configure help 4 2 Installation 59 4 2 2 Creating executables Here is the makefile in daets examples used to produce the executables whose results are reported in this user guide file makefile in daets examples This file is used to produce the results in the DAETS user guide We use g on an x86 MacOSX machine SET C COMPILER AND FLAGS CXX g CXXFLAGS g Wall ansi SET DIRECTORIES directory where DAETS is unpacked DAETS HOME daets directory where FADBAD is unpacked FADBAD HOME FADBAD IPOPT include directory IPOPT_INCLUDE HOME ipopt include IPOPT lib directory where libipopt a is located IPOPT_LIB HOME ipopt lib CXXFLAGS I DAETS include I FADBAD I CIPOPT_INCLUDE LDFLAGS L DAETS lib L IPOPT_LIB LIBRARIES LDLIBS ldaets lipopt llapack lblas lg2c 1m IF THE LINKING GIVES UNDEFINED REFERENCES HAVING _gfortran DO LDLIBS lgfortran FOR MAC OSX YOU MAY HAVE TO ADD LDFLAGS bi
28. 1 0 Nedialko S Nedialkov John D Pryce Department of Computing and Department of Information Systems Software Cranfield University McMaster University Shrivenham Campus Hamilton Ontario Swindon L8S 4K1 Canada SN6 8LA UK May 25 2009 Contents i i sage 1 ETE o EE E E ESE 1 A e ai ea 3 12 1 Am example program o A el ES 3 1 2 2 Offsets initial values quasilinearity consistency 6 lasses and methodg socia RR de RE re a CRS 7 ea daa eer ree 7 L 5 DAR solver EA tt e e a hee Gee he ee eo a be D 8 1 3 3 DAEsolutionclasd o poa esrar pakeke bieti rika 10 1 3 4 golyerExitElag UDA A 12 ii CONTENTS iii 35 37 37 39 3 3 One step mode Passing data to a problem 41 3 4 A continuation problem 202 0 2 0000 eee ee 45 3 5 A golf putting proble 48 6 Setting an output functions s 2 s ae d on e saa a oe Roe a oa a 53 TOA wa era WA A P SS es AS 53 4 Installatio 56 41 Obtaining DARE 56 4 1 1 Demo version o ss ca 4 a 644 4 6 vos 4a 4 908 4 a 44 44 56 412 Fullversion s csr astas sparra tua grea raud 56 4 2 In allatio A O e A 2 B a se A Boe Oe ee RB E 57 i 57 59 59 5 Further theory FAQ 61 51 Structural analyse 61 64 70 oe OG Se eee ee Be eee Sr ee os 70 71 Bibliography 73 ersion 1 1 Mas In the electronic version of this document every cro
29. 1 shows those of y and y while no data about is shown This is because the offsets d of the variables see IL2 2 are 2 2 and 0 respectively and the behavior of printSolution is to display up to the d 1 th derivative of the jth variable for each j in case the DAE is quasilinear as it is here see 1 2 2 The program also displays a few lines of splash screen of the IPOPT package 1 2 Quick start guide 5 TABLEAU o0 1 2 c_i ol 0 O il 2x O 0 2 O 0 2 d_jl 2 2 0 Index 3 DOF 2 SOLUTION t 1 000000e 02 x x x0 5 683109e 00 4 563114e 00 x1 5 950172e 00 8 716614e 00 Figure 1 3 Output of the program from Figure 2 Some programming points follow 1 The problem size n 3 line 3 and integration endpoints tO 0 tend 100 line need not be given names For instance one could delete line 3Jand replace line 5 by DAEsolver Solver 3 DAE_FCN fcn 2 Setting the IVs in lines 1YH19 is done by a single statement containing chained calls to setX This could be done by separate statements such as x setT 0 0 followed by x setX 0 0 10 0 etc the choice is a matter of taste 3 To access individual solution components e g for finer control of output than offered by the convenience function printSolution use getT and getX see 2 2 4 The IVs given in lines MAL happen to be consistent the initial position x y 10 0 exactly satisfies the obvious const
30. 13 3 is used to initialize such arrays see also updatePoint P _2 2 1 3 Basic classes and methods The classes constructors and methods required for the simplest use are described here All other methods are described in Chapter B 1 3 1 Defining the DAE system Function fcn defines the f in the system LJ It is a templated function that is instantiated with several different actual types T to support the automatic differenti ation and structural analysis For coding purposes T can be thought of as denoting type double with some restrictions described below The specification of fcn is 1 3 Basic classes and methods 8 template lt typename T gt void fcn T t const T x T f void param t Input denotes the variable t Input Array of length n where x j 1 denotes the variable x for j 1 n Output Array of length n where f 1 must be set to the value of f for VE V rs param Input optional pointer to a user defined object of arbitrary type that passes parameters to the definition of the DAE The code of fcn can use the following operations on variables of type T The arithmetic operations Assignment and the shortcut operators Assigning a double or a literal numeric constant value as in T sum 0 These analytic standard functions Sin cos tan sqr sqrt exp pow log asin acos atan sqr y computes y Differentiation to any non negative integer order done
31. 28 plot vdp 1 0e 00 out u 1 2 notitle with lines lw 3 set output vdpOxxp eps set xlabel y font 28 set ylabel y font 28 plot vdp 1 0e 00 out u 2 3 notitle with lines lw 3 mu 10 set title 4 Symbol m 10 font Times 28 set output vdpix eps set xlabel t font 28 set ylabel y font 28 3 3 One step mode Passing data to a problem 44 o Ra F w D 100 F 50 F 50 100 F 3 150 0 100 Figure 3 1 Plots of x t versus t and x t versus x t for the Van der Pol s equation plot vdp 1 0e 01 out u 1 2 notitle with lines set output vdpixxp eps set xlabel y font 28 set ylabel y font 28 plot vdp 1 0e 01 out u 2 3 notitle with lines mu 100 set title Symbol m 10712 font Times 28 set output vdp2x eps set xlabel t font 28 set ylabel y font 28 plot vdp 1 0e 02 out u 1 2 notitle with lines set output vdp2xxp eps set xlabel y font 28 set ylabel y font 28 plot vdp 1 0e 02 out u 2 3 notitle with lines lw 3 lw 3 lw 3 lw 3 3 4 A continuation problem 45 3 4 A continuation problem This problem comes from Layne Watson 23 We seek a fixed point that is a root of x g x for the nonlinear function g R R defined by OX gi X1 En exp cos i 3 4 T Lie T
32. 3D plot its deficiencies are obvious Where the green is appreciably curved the ball spends most of its time in the air or underground A correct formu lation is as a DAE which is a typical index 3 mechanical system like L2 except for the presence of friction The reaction R acts as a Lagrange multiplier like A in 12 Let x z y z be the position of the ball Taking the mass of the ball as unity one may write the equations of motion as X uRn Rng Gk 3 7 where dot is d dt n is the unit vector tangential to the surface in the direction of the velocity x ns is the unit normal to the surface in the upward direction G is the acceleration of gravity and k is the unit vector vertically upward Thus x Vg x Es 28 Tyga Pa The DAE consists of 3 01 8 7 after substituting the formulae 3 8 into B It is physically clear there are 4 DOF the initial x y and t Formulating as a BVP we remove 2 DOF by specifying x y o yo at t 0 Specifying that the path reach the hole at some chosen time tena would remove the remaining 2 DOF but this is not realistic a golfer does not choose tena but aims to make the ball s speed nearly zero at the moment it reaches the hole so that it does not jump the hole So we treat tena as an extra unknown and remove the extra DOF by the condition X tenq 0 A MATLAB formulation of this as an ODE system has been successfully solved by BVP4C for a number of sha
33. DAETS User Guide Version 1 1 Nedialko S Nedialkov John D Pryce McMaster University Cranfield University Canada United Kingdom Copyright 2008 2009 N S Nedialkov and J D Pryce Preface DAETS Differential Algebraic Equations by Taylor Series is a C package that integrates an Initial Value Problem IVP for a differential algebraic equation DAE system of arbitrary index and order over a range using a Taylor series method It can report detailed information about the structure of the DAE as well as the numerical solution of an IVP either at the end of the range or step by step We believe the first report of solving a DAE by Taylor series in an ad hoc way is in the work of Y F Chang with G F Corliss on the ATOMFT code c 1994 John Pryce JDP learned of this on a visit to Corliss in 1996 during which between them they found the basic principles underlying the structural analysis approach JDP soon developed these into a systematic method and wrote a proof of concept code in MATLAB in 1996 an introductory paper in 1997 and a basic theory paper in 2001 There was little progress until the collaboration with Ned Nedialkov NSN began after JDP acted as external examiner for NSN s PhD thesis 13 on the VNODE vali dated ODE solver in 1999 VNODE s architecture strongly influenced that of DAETS We started work on theory and code of DAETS in 2002 with papers 14 15 6 Ex perience over several years with coding DAE probl
34. DAEsolver h using namespace daets PA A DAE DEFINITION template lt typename T gt void compDAE T s const T x T f void param The param argument conveys the problem size int nplusi intx param int n nplusi 1 T lambda x 0 Define 0 function S that specifies s to be arc length f O 1 0 for int k 0 k lt nplusi k f 0 sqr Diff x k 1 T sum 0 0 for int k 1 k lt n k sum x k The algebraic equations are f 1 to f n for int i 1 i lt n i f i x i lambda exp cos i sum gt gt FORWARD FUNCTION DEFINITION static void outputfcn FILE fd const DAEsolution amp x int nplus1 MAIN PROGRAM int main 4 int n 10 nplus1 n 1 DAEsolver Solver nplusi DAE_FCN compDAE amp nplusi DAEsolution x Solver double tO 0 0 tend int order 15 double tol 1e 7 Initial guess has x and x all O if except x 0 that is d lambda ds is 1 EE which hopefully gets us going in right direction x setX 0 0 0 0 setX 0 1 1 0 for int i 1 i lt n i x setX i 0 0 0 setX i 1 0 0 x setT t0 Solver setOrder order Solver setTol tol x setOneStepMode true FILE fd fopen laynewatson2 out w SolverExitFlag flag success while x getX 0 0 lt 1 amp amp flag success recall x 0 0 is lambda Solver integ
35. T STATE SUCCESS The integration has completed successfuly STATISTICS TIME siii eee 0 016 ORDER sia 15 DEP Sit aoe acd sacha 1 TOLERANCE accepted v v awa eu 1 relative 1 0e 07 rejected 0 absolute 1 0e 07 diia 0 0 STEPSIZES smallest 0 629 largest 0 629 SOLUTION t 8 750393e 01 x x x0 1 000000 e 00 6 442725e 02 x1 1 491914e 00 1 691229e 01 x2 5 066654e 01 7 231676e 02 x3 3 890434e 01 4 516114e 03 x4 9 273171e 01 1 377331e 01 x5 2 419807 e 00 1 464676e 01 x6 2 186966e 00 5 771084e 01 x7 7 729182e 01 3 289186e 01 x8 3 720929e 01 4 789042e 02 3 5 A golf putting problem 48 x9 5 865923e 01 2 006262e 01 x10 1 753840e 00 6 616551e 01 From the program s output to a file the plot in Figure 3 2 was produced Figure 3 2 Layne Watson problem n 10 The paths of x and 219 are plotted against the original parameter Many turning points can be seen The markers show successive integration points Points beyond 1 are part of the event location Tests showed that for reliable results one should impose a stepsize limit by the setHmax method not done in this code Otherwise one may jump from the correct path to another very close one A limit of 0 3 worked for all experiments we tried namely every value of n from 2 to 40 and several larger n With many combinations of order and tolerance we always found the same solution at 1 Many others exist one finds some of them by c
36. a function of the tolerance or may be set by the user but is fixed during one call to integrate The following block diagram shows the basic order of events in solving a problem as illustrated by the example program below Structural Set ning A Initial E Consistent E E Event location finding a t where some function of the solution becomes zero is not currently provided However the integrator can be put into one step mode allowing event location to be coded in the calling program See B A for an example Auxiliary functions to control the integration are described in Chapter B For instance the TS order may be changed as may the tolerance and the type of accuracy control relative absolute or mixed the number of solution components printed etc 1 2 Quick start guide 3 1 2 Quick start guide 1 2 1 An example program A simple program is shown in Figure L2 It can be adapted to solve other simple problems As given it solves a problem with the simple pendulum system a DAE of differential index 3 defined by the equations 0O f zx xA O g9 y yA G 1 2 2h E 1 see Figure LI Here gravity G and the length Z of the pendulum are constants The independent variable is time t The dependent state variables are the coordinates x t y t of the pendulum bob and the Lagrange multiplier A t The program computes the motion on 0 lt t lt 100 with IVs x y 10 0 and 2 y 0 1 att
37. a more realistic formulation as an index 3 DAE though the problem is really a boundary value problem and our main program does not handle this aspect In 36 we illustrate how an output function is set so integrate outputs data after each accepted step All the code shown in this chapter is included in the DAETS distribution 3 1 Handling an unfamiliar DAE If a DAE is new to you you probably do not know the offsets and other results of its structural analysis nor the index set J in 2 5 i e for what j and k must a be initialized Even when this is known it may be unclear what consistent initial values look like To remedy this you might write first a main program that merely creates a Solver and calls the reporting functions in 2 1 2 The output particularly of printDAEpoint Structure will give the needed structural information Note printDAEtableau displays an n x n matrix so it should be used with caution when n is large For the chemakzo program of 43 2 below instead of lines onwards one can initially just write x printDAEinfo x printDAEtableau x printDAEpointStructure return 0 gt which produces 37 3 1 Handling an unfamiliar DAE DAE not quasilinear S126 ie E ee ee ee 6 UNAS AA 1 TABLEAU 0 1 0 0 0 0 11 0 1 0 0 21 0 0 1 0 0 31 0 0 1 0 4 0 0 0 1 0 51 0 0 Ox Variable derivatives x0 0 1 x1 0 1 x2 0 1 38 With the output from this you know which componen
38. a set of h files and a precompiled library libdaets a These files are stored in zipped files as described below 4 1 1 Demo version A free demo version of DAETS is available at http www cas mcmaster ca nedialk daets This version is restricted to solving systems of at most 8 equations It is distributed as zipped files Architecture OS Compiler Zip file Cygwin daets 1 1 demo x86 cygwin gcc zip x86 Linux GNU daets 1 1 demo x86 linux gcc zip Mac OSX daets 1 1 demo x86 osx gcc zip 4 1 2 Full version The full version is available through Canada s innovation portal Flintbox at 56 4 2 Installation 57 http www flintbox com It is distributed in the following zipped files Architecture OS Compiler Zip file Cygwin daets 1 1 x86 cygwin gcc zip x86 Linux GNU daets 1 1 x86 linux gcc zi1p Mac OSX daets 1 1 x86 osx gcc zip A GNU daets 1 1 x86 64 linux gcc zip Pa Li PathScale daets 1 1 x86 64 linux pathCC zip PowerPC Mac OSX GNU daets 1 1 ppc osx gcc zip For other combinations of OS architecture compiler and suggestions on improving DAETS please contact Ned Nedialkov O John Pryce d pryce ntlworld com 4 2 Installation After downloading the corresponding zip file unzip it by e g unzip This will create a directory daets with subdirectories subdirectory contains include h files lib libdaets a examples source code of the examples in this guide and the makefile used to produce the executabl
39. acy test esttype should be set to AbsRel The default value AbsRel should be chosen unless there are good reasons for a different choice 2 1 Controlling integration the DAEsolver class 21 setTol method This method lets the program alter the accuracy tolerance and the type of accuracy control from the default of a mixed test with tolerance 1078 Calling setTo1 with no arguments resets the working values to their defaults void setTol double tol 1e 8 ErrorEstType esttype AbsRel throw std logic_error tol Input the accuracy tolerance default 1078 Constraint tol 107 107 If tol is outside this range an exception std logic_error results esttype Input The type of accuracy control default is a mixture of absolute and relative accuracy control esttype must be one of the three values of the enumeration type typedef enum ErrorEstType Abs Rel AbsRel getTol method double getTol const Returns the current tolerance getRTol method double getRTol const Returns the current 7 in 2 1 getATol method double getATol const Returns the current 7 in 2 1 2 1 Controlling integration the DAEsolver class 22 getErrorEstType method ErrorEstType getErrorEstType const Returns the current type of accuracy control 2 1 4 Controlling the Taylor series order We denote the order integrate uses by p The Taylor series of the jth variable x is comput
40. ass The DAEpoint class is a base class for DAEsolution A DAEpoint object stores an irregular array i e values representing derivatives ay for j k in an index set J see 2 5 2 3 1 Constructors DAEpoint const DAEsolver amp Solver throw std logic_error Solver Input a DAEsolver object Result Constructs a DAEpoint object with the correct shape for Solver If Solver is constructed from an ill posed DAE an exception std logic_error results DAEpoint const DAEpoint amp X X Input a DAEpoint or a DAEsolution object Result Constructs a DAEpoint object with the same shape as X and copies the contents of X into it 2 3 2 Access setX method DAEpoint amp setX int index int order double value throw std out_of_range index order value Input such that x setX j 1 k vu sets the entry of x representing La the kth derivative of the jth variable equal to v This is explained in IL3 5 Constraint j k must be in the index set J of 2 5 Otherwise an exception std out_of_range results Returns A reference to the updated object Returning a reference lets one chain calls as with the setX method of DAEsolution see 1133 2 3 Arithmetic on solutions the DAEpoint class 33 getX method double getX int index int order const throw std out_of_range index order Input such that x getX j 1 k returns the current value of the entry of x representing a the At
41. ation the DAEsolver class 24 By default hmax is set to the largest finite double precision number and can be set to a smaller value by the user The usual reason for doing so is that the solution path is observed to have long smooth stretches during which h builds up alternating with short regions of rapid change that may be missed altogether if h is too large The continuation problem in Bis an example Let 7 max to tenal f where the integration is from to to tena Then hmin 16 T4 T 2 4 where T4 is the next double precision number after T hmin cannot be changed by the user getHmin method double getHmin const Returns current Amin getHmax method double getHmax const Returns current Amax setHmax method void setHmax double hmax hmax Input scalar Result sets max to hmax If the computed by integrate Amin gt Amax an exception std logic_error oc curs In this case Amax should be set to at least Amin 2 2 Solution path properties the DAEsolution class 25 2 1 6 Passing data to a DAE function passDataToDAE method void passDataToDAE void param param Input pointer to data that need to be passed to the function fcn that specifies the DAE through the void param parameter in fcn If param is the NULL pointer a call to passDataToDAE has no effect The user should ensure that data is correctly extracted inside fcn For examples illustrating the use of t
42. atrix getStructuralIndex 18 getTol PI getType B7 getT getX 271 83 integrate islllPosed 18 isQuasilinear 18 norm 85 operator 35 operator e 34 operator e 35 operator 35 operator 81 passDataToDAE 25 printDAEinfo printDAEpointStructure printDAEtableau 19 printSolutionState S1 printSolution 12 printStats setFirstEntry setHmax PZ setOneStepMode 29 setOrder 23 setOutputFunction BO setPrintProgress setTol setT I setX II B3 updatePoint DAEpoint 71 getNumDerivatives 33 getNumVariables 33 getX 33 norm 35 operator 35 operator e operator e 35 operator 35 operator BA setX DAEsolution getCPUtime getNumAccSteps B1 getNumRejSteps B1 getType getT getX BY INDEX printSolutionState printSolution printStats BY setFirstEntry B29 setOneStepMode 29 setT setX updatePoint DAEsolver getATol BI getCVector 17 getDVector getErrorEstType getHmax 24 getHmin 24 getMaxOrder 23 getMinOrder getNumDegsOfFreedon getOrder 23 getRTol Bi getSigmaMatrix 7 getStructuralIndex 18 getTol integrate isIllPosed 18 isQuasilinear 18 passDataToDAE printDAEinfo printDAEpointStructure printDAEtableau 19 setHmax PX setOrder 23 setOutputFunction BO setPrintProgress setTol Bij 72 Bibliography 1 2 10 11 12
43. clidean arc length along the path one can insert weights to scale the components of y if needed This formulation gives an index 1 system of n 1 equations and variables To match array indexing in C we label these from 0 to n with being variable 0 and B 5 being equation 0 The offsets are co 0 C1 Cn 1 and all do d 1 The n 1 x n 1 System Jacobian has 2 y as its first row with f beneath and is nonsingular if f has full rank since y is orthogonal to the rows of fy How does DAETS know which direction to take along the path The code recognizes the resulting DAE as not quasilinear and thus requires values up to stage k 0 as an 000 JO UR O NH e er pa e O OW WW W W W WW W W W N UNO RI NONI NNN NY YN AA A A hp 00 40 0 WNHRrF OTOH AN DA FE wWwnNHrF OCA N DA A O NHO a Fr a 00 410 BR O NHO lt OL OL OL OL OL OL OL OL OL Ot WR V NOn A OW NHO 3 4 A continuation problem 46 initial guess Since all offsets d are 1 these values are Y1 Yn 41 Ylo Unyr not only an initial position but also an initial direction This is just what is needed Code to solve the Layne Watson problem by arc length continuation is given below at a smaller font size because of its length The main program uses two loops to output solution values in one step mode and to locate 1 by a Newton iteration file laynewatson cc in daets examples include lt cmath gt include
44. cout state Input Usually the value returned by a call to the integrate method Input output stream default is cout Result printSolverExitFlag outputs into s a message describing the meaning of state 1 3 5 Fixed and free initial values Clearly there must be limits on how many IVs one can specify as fixed For instance in the pendulum system one cannot give random values to both x and y and fix both of them because the constraint x y L generally will not be satisfied The rules for what one can fix are entirely determined by the offsets One has to explain at this point that the irregular array X of values that define a consistent point see 2 5 divides into disjoint parts that are computed in stages and have a natural indexing by stage numbers k in the range d lt k lt a where d max d and a is 1 if the DAE is quasilinear O otherwise Namely at stage k to have consistent values one must solve the equations f 0 for those i 1 n for which k c gt O qe for those j 1 n for which k d gt 0 1 3 Basic classes and methods 14 We illustrate with the pendulum assuming it has been changed to be non quasilinear The figure below shows the situation Irregular array X consists of x y x y y A k 2 0 1 Unknowns Fixable 2 1 1 2 1 1 3 3 0 The figure also shows the reason for the term offset each list of derivatives is shifted
45. e that solution is unique In the nonlinear case there may be several solutions consistent with these values however including the next level of derivatives in X restores uniqueness This affects the user at the initial point t t The initial conditions you provide are usually only approrimately consistent values which the code corrects by a root finding process By providing hopefully good guesses of these extra derivatives you make it more likely in cases of non uniqueness that the code finds the consistent values hence the DAE solution that you intended An example is the use of DAETS to do arc length continuation example in 34 Here t is arc length from an initial point Xy along a path in some RY There is inherent non uniqueness from Xy you can traverse the path in either direction In this case the extra derivatives DAETS forces you to provide turn out to be a guess of the unit tangent at Xy in the direction you wish to go In most cases you can estimate this a priori 5 2 Frequently asked questions 5 2 1 How well does the accuracy tolerance control the accuracy Our experiments see for instance 16 Section 5 1 suggest that for most prob lems the error in the solution is roughly proportional to the tolerance if the type of control relative absolute or mixed specified by the esttype argument of setTol is unchanged Thus the behavior is essentially the same as for an ODE code However this relation between
46. ed up to the term of order p dj where d is the offset of x The minimum and maximum allowed orders are denoted by Pmin and Pmax respectively Currently Pmin 1 The maximum number of Taylor coefficients that can be generated is set to MaxLength 80 Then the largest allowed value of p for a given DAE is Pmax MaxLength 1 max d 2 2 j When the constructor of DAEsolver is executed it checks that pmax gt Pmin If this is not the case an exception std out_of_range is thrown with a message stating the smallest value for MaxLength that has to be given when DAETS is compiled When integrate is allowed to select the order itself because setOrder below is not called or is called with argument 0 it uses a heuristic formula that has usually performed well p max pmin min 0 5 n tol 1 Pmax 2 3 Here tol is the tolerance given to DAETS and x is the next integer gt x In the absence of the limits pmin and Pmax and with the allowed range of tolerances see setTol method this gives values in the range 2 to 18 Tests on many problems have shown that the CPU time to solve a problem as a function of p for fixed tol usually has a flat minimum with the best order ppest in the range 15 to 20 and slowly increasing as tol is decreased Further it is much safer to take p too large rather than too small A small p like 1 or 2 can increase solution time by orders of magnitude while taking p 2 X Ppest say typ
47. ems and performance tests resulted in many improvements to the user interface and the algorithms The first version 1 0 of the code and user guide were released in June 2008 An improved version 1 1 was released in May 2009 We acknowledge with thanks the support given to JDP by the Leverhulme Trust and the Engineering and Physical Sciences Research Council both of the UK and to NSN by the Canadian Natural Sciences and Engineering Research Council We thank the Fields Institute in Toronto in whose stimulating environment we had the discussions that committed us to this project We acknowledge the help of friends and colleagues George Corliss for having provided the original idea and for having carefully read and constructively criticized our writings over many years including this guide Ole Stauning for providing assistance with his FADBAD automatic differentiation package which was instrumental in the development of both VNODE and DAETS Wanhe Zhang for working on the design and many numerical experiments with an early version of DAETS Andreas W chter for assistance with IPOPT Andrea Walther and Andreas Griewank for helpful discussion on automatic differentiations and their ADOL C package Rob Corless Wayne Enright and Silvana Ilie for insightful discussions on theory especially on methods for numerical error estimation and control Emil Sekerinski for his input on object oriented design Bingzhou Zheng for extensive testing of version
48. es The DAETS code uses the third party components BLAS 0 Lapack 1 Fap BAD I IpoPT 22 and Lar 9 The LAP program is distributed with the source code of DAETS We describe how the rest are installed and show how programs using DAETS are compiled and linked We use the GNU make utility 4 2 1 Third party components LAPACK and BLAS The LAPACK and BLAS libraries are likely installed on your machine If not you can download and install them from http www netlib org lapack and http www netlib org blas respectively 4 2 Installation 58 FADBAD 1 2 Download FADBAD from http www fadbad com fadbad html Extract the FADBAD files by unzip FADBAD 2 1 zip These files will be stored in subdirectory FADBAD of the directory in which you have executed tar IPOPT 1 3 Download the Fortran version Ipopt Fortran_20060ct21 tgz of the IPOPT package from http www coin or org download source Ipopt Fortran Extract the IPOPT files by tar zxvf Ipopt Fortran_20060ct21 tgz This will create a directory Ipopt Fortran_20060ct21 containing the IPOPT files Download mc19ad f ma27ad f and ma47ad f from http www cse clrc ac uk nag hs1 fand save them in Ipopt Fortran_20060ct21 OTHERS HSL For an explanation of how to download these files see Ipopt Fortran_20060ct21 OTHERS HSL INSTALL HSL Download dimach f from http www netlib org blas dimach f and save it in Ipopt Fortran_2
49. f outfile n Figure 3 4 Setting an output function 54 3 7 Summary 10 r T r 50 r r r 8 F 40 6 30 4 7 20 2 10 x OF 7 Pa 0 2 F 10 4 20 6 J 30 8 40 10 1 50 0 5 10 15 20 t 40 30 20 0 gt 0 Di 10 20 30 2 40 f 0 5 0 15 20 0 5 10 15 20 t t 3 T r T 5 r r 2 49 7 2 L lt 149 T a 1 0 5 7 a 0 i 15 i 0 5 10 15 20 0 5 10 15 20 t t Figure 3 5 Plots of x y y A and A versus t for the simple pendulum 55 Chapter 4 Installation Read this chapter if DAETS is not already installed on your system Section I describes how DAETS can be obtained 44 2 describes the installation process and 41 2 3Jexplains how the examples can be executed The lists below show the supported platforms at the time of writing Visit the web sites mentioned for information on the current status of support In our experience the main obstacles to installation come from the third party components mainly IPOPT This is not a criticism of IPOPT which is a powerful and inevitably somewhat complex optimization package See if it is already on your system in which case it may be possible to link to it instead of re installing 4 1 Obtaining DAETS There are demo and full versions of DAETS They are distributed as
50. fcn greenparams create solver analyze DAE Solver printDAEinfo Solver printDAEtableau Solver printDAEpointStructure DAEsolution x Solver create solution object initial position x y 0 0 with a guessed value z 0 initial velocity x y vx0 vy0 with a guessed value z 0 x setT t0 setX 1 0 0 0 Fixed setX 2 0 0 0 Fixed setX 3 0 0 0 setX 1 1 vx0 Fixed setX 2 1 vy0 Fixed setX 3 1 0 0 open output file FILE fd fopen golfputt out w Its first line holds parameters defining the surface for int i 0 i lt 6 i fprintf fd 6e greenparanmsli fprintf fd n declare variables to record change of direction of motion double vx vy double vxold vx0 vyold vy0 v_vold 1 do integration in one step mode upper limit 60 seconds SolverExitFlag flag success x setOneStepMode true double tend t0 dt Solver integrate x t0 flag to output initial consistent point outputfcn fd x while flag success amp amp v_vold gt 0 amp amp tend lt 60 60 sec time limit Solver integrate x tend flag if x getT tend outputfcn fd x tend tend dt vx x getX 1 1 vy x getX 2 1 v_vold vx vxold vy vyold vxold vx vyold vy outputfcn fd x fclose fd output results to screen x printStats x printSolution if flag success printSolverExitFlag flag return 1 return 0 3 5 A golf putting p
51. generated by gnuplot gnuplot laynewatson gp Finally golfputt produces golfputt out and you can generate plots using MATLAB and the program given in examples matlab golfputtOutStudy m If you encounter discrepancies between the numerical results in this user guide and the output of your program executions please contact the authors Chapter 5 Further theory FAQs This chapter is in two parts Read Section B I if you need more detail about the structural analysis theory on which DAETS is based and references to the literature Section B J lists some in our experience frequently asked questions FAQs about the structural analysis approach and about DAEs in general The authors welcome suggestions for useful additions to the FAQs 5 1 Structural analysis Often DAEs are generated in an exploratory way by simulation software that allows basic components mechanical electrical etc to be connected using a graphical interface It may be non obvious what the resulting DAE requires as initial data what its index is and even whether it is ill posed because it contains too little or contradictory information The structural analysis SA performed by DAETS helps to resolve such difficulties Here we give enough theory that a user can on small problems do the analysis by hand For further details see 14 15 79 Our method gives essentially the same result as does e g that of Pantelides 17 but is easier to use Current
52. h Applications in Process Engineering PhD thesis Carnegie Mellon University Pitts burgh PA 2002 L T WATSON A globally convergent algorithm for computing fired points of C maps Appl Math Comput 5 1979 pp 297 311 B ZHENG Testing scientific computing software A case study of a DAE solver MSc thesis in preparation Dept of Computing and Software McMaster University Hamil ton Ontario Canada L8S 4K1 2009
53. h derivative of the jth variable Constraint j k must be in the index set J of 2 5 Otherwise an exception std out_of_range results Returns x getX 1 k returns the current value of the entry of x representing k H j getNumVariables method int getNumVariables const Returns n number of variables a DAEpoint object stores getNumDerivatives method int getNumDerivatives int j const throw std out_of_range j Input variable number Constraint 0 lt j lt n _1 If j lt 0 or j gt n 1 exception std out_of_range is thrown which also prints the range for j Returns d 1 4 number of derivatives a DAEpoint object stores for variable j cf 2 5 Note When getNumDerivatives j returns 0 it implies that variable 2 1 cannot be accessed by setX and getX methods In such a case for example X setX j 0 value would result in std out_of_range exception 2 3 Arithmetic on solutions the DAEpoint class 34 2 3 3 Assignment operator method DAEpoint amp operator const DAEpoint amp X X Input a DAEpoint or a DAEsolution object Result copies the content of X to the object on the left in e g Y X Returns reference to the updated DAEpoint object DAEpoint amp operator double a a Input a double scalar Result X a sets all entries of X equal to a Returns reference to the updated DAEpoint object 2 3 4 Arithmetic In the operators that fo
54. hat show other useful techniques V 0 JO UU R O N H w w A A OW W Y Y wns Nb N DY NY N NN NY NY DN PR Rh Rh H H H H H H H 0 0 JD 0 R OO Ne O 0 ODN DT FF 0 NN Pp O Ob ODN Dot FF Ww N I H O 3 7 Summary file pendulumsimple outfcn cc in daets examples include DAEsolver h using namespace daets template lt typename T gt void fcn T t const T z T f void param const double G 9 8 L 10 0 f 0 Diff z 0 2 zL 0 z 2 f 1 Diff z 1 2 z 1 z 2 G f 2 sqr z 0 sqr z 1 sqr L gt Forward declaration of output function void outfcn const DAEsolution amp x void output_paranm void problem_param int main const int n 3 double tO 0 0 tend 20 0 DAEsolver Solver n DAE_FCN fcn DAEsolution x Solver x setT t0 setX 0 0 9 setX 1 0 0 0 setX 0 1 0 0 setX 1 1 1 0 FILE outfile fopen pend out w x setOutputFunction outfcn outfile set output function SolverExitFlag flag Solver integrate x tend flag if flag success printSolverExitFlag flag fclose outfile return 0 Output function void outfcn const DAEsolution amp x void out_param void problem_param FILE outfile FILE out_param convert param to FILE fprintf outfile e x getT output t output each variable and its third derivative for int j 0 j lt x getNumVariables j fprintf outfile he he x getX j 0 x getX j 3 fprint
55. his function see 3 3 and PA 2 2 Solution path properties the DAEsolution class In the method descriptions of this section x denotes a DAEsolution object at the current point t on the solution path x t and X denotes the values of x and derivatives that x holds That is X is the irregular vector comprising Xj am the computed value at t of the kth derivative of the variable x t for j k in the index set J Gk lj 1 n 05k lt d a 2 5 where a 1 if the DAE is quasilinear a 0 if not see 112 2 The basic role of an object x of the DAEsolution class is to store the numerical values t and X that jointly describe the current point on the solution path However x stores other important attributes It records its position along the solution path using the enumerated type typedef enum SolutionState Initial InitialConsistent OnPath EndOfPath A newly created x is in state Initial and DAETS deems it an inconsistent point Its t and each component of its X is flagged uninitialized and integrate will not accept x until all these values are set When consistent values at the initial ty are found x moves into state InitialConsistent this can be seen if integrate is called with tena to to find a consistent point and then exit After the first accepted step x is in state OnPath and is deemed consistent It remains in this state until a a call to integrate fails with a htoosmall or nonconsistentpt err
56. ically increases solution time quite modestly Tf you would like MaxLength increased contact N Nedialkov who can produce a precompiled library with larger MaxLength 2 1 Controlling integration the DAEsolver class 23 setOrder method void setOrder int order throw std out_of_range order Input the order to which Taylor series are generated Constraint Pmin lt order lt Pmax or order 0 e If order 0 integrate selects the order by the heuristic 2 3 e If order is in Pin Pmax integrate uses order p order e If order does not satisfy the above constraint exception std out_of_range is thrown with a message saying that the order must be in Pin Pmax Or 0 getOrder method int getOrder const Returns the current order p getMinOrder method int getMinOrder const Returns Pmin getMaxOrder method int getMaxOrder const Returns Pmax 2 1 5 Controlling the stepsize On each integration step after the first integrate computes a trial stepsize h for the next step based on the current h and the ratio r in 2 1 such that a prediction of the estimated error is below a tolerance The stepsize is restricted by the values of the largest max and smallest Amin allowed stepsizes If the trial h satisfies A lt Ain this is a fatal error and the integration is termi nated If Al gt Amax h is replaced by Amax with the same sign as h 2 1 Controlling integr
57. iven by equation 5 3 getStructuralIndex method int getStructuralIndex const Returns The structural or Taylor index of the DAE given by equation 5 0 isIllPosed method bool isI11Posed const Returns true if fcn describes an ill posed DAE see M T else false isQuasilinear method bool isQuasilinear const Returns true if the DAE is deemed to be quasilinear see 1 2 2 else false 2 1 2 Reporting structural analysis data printDAEinfo method void printDAEinfo ostream amp s cout const s Input output stream default is cout Result Outputs into s whether the problem is quasilinear size of DAE index of DAE For the pendulum example this method prints 2 1 Controlling integration the DAEsolver class 19 DAE quasilinear printDAEtableau method void printDAEtableau ostream amp s cout const s Input output stream default is cout Result Outputs into s signature matrix of the DAE indicating the positions of a highest value transversal HVT offsets of the problem index number of degrees of freedom DOF For definitions of these items see Chapter 5 For the pendulum example this method prints TABLEAU O 1 2 lci ol 2 Ox O 1l 2 0 0 21 Ox 0 2 2 1 Controlling integration the DAEsolver class 20 printDAEpointStructure method void printDAEpointStructure ostream amp s cout const s Input output stream defaul
58. j of X Write d beneath column j for each j 2 do EN Set each c so that dj c oj holds on T Write c to the right of row i of 2 for each 1 4 Increase each d where needed by the least amount that makes dj c gt gij hold everywhere 5 until nothing changes The c and d can only increase during the loop The algorithm terminates if and only if T really is a HVT We like to show the results by a signature tableau which is X annotated with the offsets c d and the names of the functions and variables and the positions of a HVT We like to write entries as a dash for readability in larger systems typically almost all entries are oo Also these are forbidden positions a HVT of a well posed DAE cannot have an entry in a position This is illustrated below for the pendulum example which has two HVT s marked and E Y fp 070 g 2 w 0 h Loe 0 J 2 dj 2 2 0 DOF 2 It is easily seen that the above algorithm only needs one iteration to find the offsets of this signature matrix Also shown is the number F of degrees of freedom DOF given by F value of the HVT Sod gt Ci 5 3 5 1 Structural analysis 63 It is the number of independent IVs it requires which is the same as the maximum number of IVs that can be specified fixed see setX in L3 3 For n up to about 8 setting up X and finding a HVT by eye is usually easy Otherwise
59. lize x integrate x from tO up to some t1 DAEsolution y Solver copy the values stored in x and their types y updatePoint x set t and other values in y and integrate onward from tl 2 2 Solution path properties the DAEsolution class 29 2 2 3 State control The rationale of these methods is discussed at the start of 22 setFirstEntry method void setFirstEntry Result x setFirstEntry puts x into the Initial state so it can be re initialized It also resets to zero the variables that gather the statistics of setOneStepMode method void setOneStepMode bool onestep onestep Input if onestep true integrate returns after it accepts a step If onestep false integrate attempts to reach tend from current t X stored in x 2 2 4 Reporting during integration setPrintProgress method The setPrintProgress method switches a DAEsolution object x in and out of a printprogress mode where the progress of integration is displayed If x is in this mode then integrate prints at each accepted step the current t number of steps as given by getNumAccSteps in 2 2 5 stepsize and local error estimate in a fixed position on the screen without scrolling void setPrintProgress double s 0 s Input Ifs gt 0 printprogress mode is switched on In addition the integration pauses for s seconds after each output to make it easier to read what is displayed If s lt 0 printprogress
60. llow the arguments must have the same shape that is the irregular arrays X that they hold must be indexed over the same J in 2 5 If not an exception std logic_error is thrown which also prints the reason for the mis match The symbol e stands for any of and l For convenience in the notation below by X i 7 we denote the value for al in object X stored operator e method DAEpoint amp operator e const DAEpoint amp X throw std logic_error Input a DAEpoint or a DAEsolution object Returns Reference to this such that this j k e X j k for all j k J 2 4 DAETS default settings 35 operator e method DAEpoint operator e const DAEpoint amp X const DAEpoint amp Y throw std logic_error X Y Input DAEpoint or DAEsolution objects Returns A DAEpoint object containing X j k e Y i j for all i j J 2 3 5 Comparison operator method bool operator const DAEpoint amp X const DAEpoint amp Y throw std logic_error X Y Input DAEpoint or DAEsolution objects Returns true if X i j Y i j for all i j J false otherwise operator method bool operator const DAEpoint amp X const DAEpoint amp Y throw std logic_error Returns X Y is equivalent to X Y 2 3 6 Norm norm method double norm const Returns X norm returns max j p eJ Xjx where X is the irregular array stored in X
61. ly several simulation packages offer SA facilities for model checking The steps to perform SA are as follows 1 Form the n x n signature matrix X 0 of the DAE where order of the derivative to which the jth variable x oi q occurs in the ith equation fi or 5 1 oo if x does not occur in fj 2 A transversal T of X is a set of n positions in the matrix with one entry in each row and each column That is these positions can be put on the main diagonal by suitably permuting the columns or rows Its value is the sum of the oj in those positions Find a highest value transversal HVT one that makes the value as large as possible For a well posed DAE the value must be finite hence an integer gt 0 Otherwise it is oo which happens if there does not exist a transversal all of whose oj are 61 5 1 Structural analysis 62 finite The DAE is then structurally ill posed there is some error in the problem formulation 3 Find the offsets integer vectors c c1 c and d d d with all c gt 0 that satisfy d c gt dij for all 4 7 and 5 2 with equality holding on the HVT They are not unique but there exist canonical smallest offsets in the sense of a lt b if a lt b for each 7 There is a simple method to find these Algorithm Finding the offsets Input Signature matrix X and a transversal T that must be a HVT 1 The initial guess for d is the largest 0 in column
62. most of his troubles solving problems with DAETS have been C ones like forgetting the namespace line at the start of a program For specific DAE difficulties the suggestions given in 93 1 have worked quite well If you find some typical difficulties or traps in using DAETS we shall welcome tips on how to avoid them for possible inclusion in these FAQs 5 2 11 Can I have more derivatives printed We plan this in the next version One tip is that if you have a quasilinear problem so the code only prints out derivatives up to order d 1 and you would like it to do so up to order d you can do so by fooling the code into thinking it is not quasilinear For instance in the simple pendulum DAE you could add A or 10780 x A to one of the equations 5 2 12 How is the simple pendulum DAE derived Let the pendulum be made of a light rigid rod of length L with a bob which is a point of mass m at the end Let it be suspended from the origin O and take axes with x horizontal and y vertically downward Let x y be the position of the bob Let T be the force in the rod with tension being positive and compression negative Let 0 be the angle between the rod and the y axis in the direction of positive x Let G be gravity Then Newton s second law gives the equations of motion Resolving horizontally mx T sin 0 Resolving vertically my T cos mG Now sin x L and cos y L 5 2 Frequently asked questio
63. n O the theory of DAETS based on Pryce s structural analysis of DAEs and gives answers in 95 2 to some frequently asked questions Appendix A records changes in each new version of DAETS 1 1 Description DAETS advances the solution of a DAE system over a range using a variable step Taylor series TS method The system has the form fi t the x and derivatives of them 0 1 m 1 1 in terms of the unknown state variables x t j 1 n It is defined by a user supplied function that evaluates the functions fi Many DAE solvers are limited by the index of the system a non negative integer that generally indicates how difficult the system is to solve numerically see B 7 and Chapter DAETS uses a different approach that is not inherently limited by the index It interprets the user s function code symbolically and uses automatic differentiation to expand the solution in Taylor series to as many terms as required This implies 1 1 Description 2 that the function code must be built of functions that are analytic along the solution path See L3 1 for more details Unlike an ODE a DAE usually does not have a solution through each point in the solution space but only through initial values IVs that are consistent 3 1 DAETS has a structural analysis stage 95 1 that finds the underlying structure of the system and tells it just what variables and derivatives must be specified to define consistent IVs This may
64. nce There are frequent forward references to Chapter B They tell you where to go if you want to look up an item of structural analysis theory You can understand the gist of this chapter without doing so 2 1 Controlling integration the DAEsolver class The method descriptions in this section assume a DAEsolver object call it Solver has been constructed from a DAE function called fcn and that the DAE s functions and variables are f and x respectively as in equation LI Structural analysis touched on in 1 2 2 and detailed in 45 1 is done during ex ecution of the DAEsolver constructor Subsections 2 1 1 and L 1 2 describe methods that get results of this analysis or display them on standard output During integration the stepsize h changes from step to step The Taylor series TS order p is fixed over one call to integrate Either it is chosen by the code on entry to integrate as a function of the tolerance or the user program sets it before calling integrate Subsections 2 1 3 onward outline aspects of the integration algorithm alongside the methods by which the user can tune its behavior by changing strategy parameters to do with h p and the accuracy tolerance These parameters are attributes of a DAEsolver object call it Solver If one wants to solve the same problem using say different Taylor orders p one may either use one Solver and change its p between integrations or use several Solvers each having a p that
65. nd_at_load all examples EXAMPLES chemakzo chemakzoUnknown golfputt pendulumshow_sa pendulumsimple pendulumsimple_outfcn laynewatson vdp examples EXAMPLES clean rm rf out o core bak EXAMPLES This file is self explanatory We make the following comments If the linker cannot find the LAPACK and BLAS libraries you need to add to LDFLAGS the paths to them The g2c GNU FORTRAN runtime library in LDLIBS is normally present on the systems on which we have compiled DAETS Install it if it is not on your system and if you use GNU compilers If you compile e g the PathScale C 18 compiler libg2c a is not needed 4 2 3 Running the examples We recommend that you execute the programs in the examples directory After you edit the makefile in this directory type make This will create the executable files pendulumsimple 2 pendulumshow_sa 2 1 2 chemakzoUnknown 93 1 chemakzo 82 vdp E3 laynewatson 8 4 and golfputt 8 5 In subdirectory gnuplot of examples there are the files vdp gp and laynewatson gp used to generate the plots in Figures 8 3 and B2 respectively The executable vdp produces the data files vdp 1 0e 00 out vdp 1 0e 01 out and vdp 1 0e 02 out in the current directory from which you can generate the plots in Figure 3 3 by 4 2 Installation 60 gnuplot gnuplot vdp gp Similarly laynewatson produces the data file laynewatson2 out from which the plot in Figure B 2 can be
66. not be obvious to the user who must set up this data before calling the integrate method of DAETS If the user gets this wrong DAETS prints an error message listing missing and or superfluous IVs Finding a consistent initial point can be the hardest part of the solution task for highly non linear DAEs Many DAE solvers require exactly consistent values to be given This is desirable but difficult e g the value of a Lagrange multiplier or an acceleration is often not known a priori DAETS allows guesses to be given mixed with fixed values For instance one may specify The IVs are z 2 x2 0 fixed and x 5 13 0 guessed DAETS treats computation of consistent data as an optimization problem and passes this to the optimization code IPOPT 22 During this computation guessed values are allowed to be altered while fixed values are left unchanged Subject to this the consistent point that is found is the one that is closest to the point given by the user in a least squares sense The main work of DAETS is done by the integrate method 1 3 2 which ad vances the solution step by step from the found consistent point x at the initial t t up to a final t tend Each step generates and sums a Taylor series and then does a projection stage to preserve consistency within the specified accuracy tolerance Assuming one reaches tend successfully one can efficiently integrate to further tend values as required The TS order is chosen as
67. ns 68 Substituting these into the the equations of motion and defining T mL gives x x 0 y Aiy g 0 The final equation is the rigidity constraint x y L 5 2 13 I m new to DAEs why are they so different from ODEs No short explanation will be useful to everyone but here s a try Take the typical case where the independent variable t is time and the dependent state variables 21 t 2 t describe how the state of some physical system evolves in time With an ODE the rate of change of each state variable is explicitly given in terms of the state variables and possibly t If there is just one scalar state variable x so that the ODE can be written dx dt f t x one can draw a picture in the t x plane Through any point say to zo draw a line segment whose gradient is f t x This line is along the tangent to the solution of the ODE that passes through to Zo Drawing enough of these segments produces a direction field for the ODE which gives a good feel for the shape of the solutions One can join the segments up by eye to sketch solutions quite accurately This is illustrated in Figure 5 1 for the ODE dx A _t dt x It is clear from the direction field that the solutions are circles and indeed the general solution of this ODE is easily seen to be x c where cis a constant The same goes for an ODE with several state variables x t except that the line segments are now d
68. ntegrate x tend flag integrate until t tend Solver printDAEtableau if flag success printSolverExitFlag flag check the exit flag x printSolution display solution return 0 Figure 1 2 Program for simple pendulum problem DAE definition 1 1 The components are numbered 0 1 as in fcn above and setX j d v sets the dth derivative of the 7th variable to v Thus x setX 1 1 1 0 initializes the first derivative of the 1 th variable to 1 i e y 1 at t 0 here Line PI does the real work of solving the IVP using the integrate method of the DAEsolver class which moves x along the solution path to the desired end point at t 100 Usually integrate accounts for most of the CPU time of the solution process Since it is a numerical method there are accuracy tolerances involved For the default tolerance values how to change them and the choice between absolute relative and mixed accuracy control see 2 1 3 If the integration failed to reach the end point this is reported by lines B3H24 using the variable flag which is of an enumeration type described in 34 The solution at t 100 or at the last t successfully reached in case of failure is reported by the printSolution method line 25 Its output is shown in Figure ci The dth column of the jth row of the display shows the dth derivative of the jth variable with j counting from 0 as previously Thus row 0 shows the final values of x and x row
69. of the DAEs in Fortran with arrays indexed from 1 Using define statements we convert between this style and C s style of indexing from 0 Hence the Fortran code can be adapted with very little change file chemakzo cc in daets examples A program to integrate the Chemical Akzo Nobel DAE problem from the Test Set for IVP solvers include DAEsolver h using namespace daets Rossa oS a DAE DEFINITION Function for evaluating f y in formulation of Akzo Nobel DAE template lt typename T gt void fcnO T t const T x T f Since we adapt the Fortran code given in this test set it is convenient to introduce these two macros define f i f i 1 define y i x i 1 constants from problem definition double ki 18 7e0 k2 0 58e0 k3 0 09e0 k4 0 42e0 kbig 34 4e0 kla 3 3e0 ks 115 83e0 po2 0 9e0 hen 737e0 intermediate variables T ri ki pow y 1 4 sqrt y 2 r2 k2 y 3 y 4 r3 k2 kbig y 1 y 5 r4 k3 y 1 sqr y 4 r5 k4 sqr y 6 sqrt y 2 fin kla po2 hen y 2 the f y functions f 1 2e0 r1 r2 r3 r4 f 2 0 5e0 r1 r4 0 5e0 r5 fin f 3 ri r2 r3 f 4 r2 r3 2e0 r4 f 5 r2 r3 r5 6 ksty 1 y 4 y 6 undef f undef y 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 71 72 74 75 76 77 78 79 80 81 82 83
70. ontinuing the path beyond the first one found but the first occurrence of 1 on the path starting from 0 specifies a unique solution 3 5 A golf putting problem This problem was presented by S M Alessandrini i as a motivational example for teaching numerical methods for two point boundary value problems BVPs With DAETS we have to treat it as an IVP rather than a BVP but the program below can be used to explore solving the BVP by means of a manual shooting method A golf ball lies on a green at zo yo in plan view ignoring the z co ordinate for now and we wish to hit it to fall into the hole at 0 0 The green is not necessarily flat and its surface is described by a function 0 g x y z z S x y 3 6 1 2 3 5 A golf putting problem 49 The rolling dynamics of the ball is ignored it is regarded as a point sliding on a surface with coefficient of friction u It is confined to the surface so no matter how curved the green the ball can never fly into the air This is a valid assumption if and only if the motion is slow enough and the curvature small enough that the normal reaction force R of surface on ball has a positive z component In the system is modeled as an ODE via the assumption that the magnitude R of R equals the component normal to the surface of the ball s weight This gives a reasonable approximation but when one solves it using the MATLAB code BVP4C and does a rotatable
71. or putting x into the EndOfPath state or b it is returned to the Initial state by acall to setFirstEntry see Only in the Initial state can the value of t X be set by the setT and setX methods or by the updatePoint method of In all other states trying to alter the t or X values in x results in an exception std logic_error x also records for each element of X whether it has been given an initial value and for the algorithm that finds a consistent point whether this value was set fixed or 2 2 Solution path properties the DAEsolution class 26 free by the relevant setX call This is recorded through the type of an element of x declared as typedef enum VarType Uninitialized Fixed Free Hence using setX one can also mark an entry in x as uninitialized Calling integrate when X is not fully initialized causes an uninitializedpt error x also has a one stepping attribute set by setOneStepMode see When this is true x is in one step mode and integrate returns after each successful step with x On the initial step there is a special feature when returning on the first call x contains a consistent point but t is not changed when returning on the second call x contains a solution at the next selected t This simplifies writing output to a file see 43 3 for an example At present one step mode is the only way to program event location or to pro duce output at each step for plotting etc However a
72. pes of green A difficulty no doubt obvious to a golfer is that if the slope of the green near the hole is so large that gravity overcomes friction there may be no solution We give a crude program that solves the DAE as an IVP A difficulty is that n in 37 creates a singularity when the ball comes to rest If one takes no precautions over this the integration oscillates indefinitely to see why try solving the one variable ODE y sign y y 0 1 which shows essentially the same behavior with any standard ODE solver To overcome this the program solves in one step mode and terminates the integra tion when it finds y this step is so different from y at the last step that their dot product is lt 0 Another useful technique is shown output of t x y z Y 2 is sent to a file not at each step but at equal time intervals of dt so that one can estimate from a graph how the speed is changing along the path file golfputt cc in daets examples This program integrates the golf putt problem from Alessandrini 1995 000 JO AUO 3 5 A golf putting problem 50 Function output cn writes to a file Line 1 contains the parameters that define the green The rest of the file contains 1 x y 2 x y z values suitable for Gnuplot include DAEsolver h using namespace daets preneren DEFINE SURFACE OF GREEN ITS GRADIENT template lt typename T gt void dfgreen const T amp x const
73. ple and suppose its length L 10 One can fix x or y or neither but not both If one fixes x to be 8 say there are two roots for y namely y 6 and y 6 Hence an appropriate guess for y should be given and IPOPT should then find the root that is nearest it Clearly one can not necessarily impose an arbitrary fixed value if we fix x at a value gt 10 then no solutions exist If one fixes neither x nor y then IPOPT will find an approximation to the consistent point nearest to the given guess of x y For instance if guesses x y 4 3 are given the found consistent values will be close to 8 6 This happens similarly on each stage k lt 0 1 3 Basic classes and methods 15 Summary This chapter has described enough of the DAETS methods and of their underlying theory to let you write a simple program like that in Figure L2 The next chapter doc uments all the remaining methods organized by the three DAETS classes DAEsolver DAEsolution and DAEpoint that are visible to the user If you are more interested in examples of code using DAETS see Chapter B For installation matters see Chapter KM Chapter 2 All Other Methods Sections 2 1 2 2 and 2 3Jof this chapter describe the remaining methods of DAEsolver and DAEsolution and the methods of DAEpoint which is a base class of DAEsolution If a method is not found here it is in IL3 The final Section P A gives a list of default parameter values for refere
74. raint 0 h x y L and the initial velocity x y 0 1 exactly satisfies the hidden constraint 0 h 2xx 2yy which means the velocity is tangential to the circle Consistency is explained in 2 2 and in Chapter Bl When inconsistent Vs are given DAETS does not automatically report the initial consistent point it found If this is desired do a first call to integrate with tend equal to t which just sets x to a consistent point and returns Report x by printSolution and call integrate a second time to solve up to the desired tend 5 DAETS has data protection features that ensure a correct set of IVs is provided see 2 2 and make it impossible to inadvertently corrupt x in the middle of an integration see start of 2 1 2 Quick start guide 6 1 2 2 Offsets initial values quasilinearity consistency In the program of Figure L2 how do we know apart from physical knowledge that a pendulum requires an initial position and velocity to be given that x y x and y are the correct set of IVs to specify This comes from the offsets mentioned above which are crucial to the use of DAETS We explain them briefly here and more fully in El see also 6 The offsets are 2n nonnegative integers computed in the structural analysis phase d is the offset of variable x and c is the offset of equation f One way to see their role is as follows If one differentiates the ith equation c times w r t t 1 1
75. rate x x getT 1 flag outputfcn fd x nplus1 printSolverExitFlag flag x printStats Now we know lambda gt 1 at this step and was lt 1 at previous step DAETS doesn t do event location yet so hand code this here x setOneStepMode false while fabs x getX 0 0 1 0 gt tol amp amp flag success 4 x setFirstEntry tend x getT x getX 0 0 1 0 x getX 0 1 Solver integrate x tend flag 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 3 4 A continuation problem outputfcn fd x nplusi printSolverExitFlag flag x printStats x printSolution if flag success return 1 return 0 47 FUNCTION FOR OUTPUT TO FILE This function outputs into a file suitable for GNUPLOT static void outputfcn FILE fd fprintf fd 4e x getT for int i 0 i lt nplus1 i fprintf fd for int k 0 k lt 2 k fprintf fd 7e gt fprintf fd n gt The output of this const DAEsolution amp x x getX i k int nplus1 program to the terminal is EXIT STATE SUCCESS The integration has completed successfuly STATISTICS TIME 2406 406 G 18 6 9097 Sees 3 387 ORDER i oii r e a ais a gids 4 5 4 15 STE PS ias 513 TOLERANCE accepted 4 509 relative 1 0e 07 rejected 4 absolute 1 0e 07 dera 0 8 STEPSIZES smallest 3 563e 04 largest 0 667 EXI
76. rawn in a multi dimensional space Most numerical methods JTI PEE A A IRNOS LAO AO AA 22 gt gt S SNS NX Xx DVD AS 1117 10 gt SN XXX MN 5 7 J 7 1 7 7 7 gt SXN XXX W 17111717 7 gt NONO X NAVA NLD EA bb AFANS vb Y OR aa COP a VA NN MANZ Z P T Bod 5 NA N N X re SL ll XN N XN N X lt 7 7 ST VNNNNN lt gt 4 YY f E E el NANA AAA AAA 10 5 0 5 10 Figure 5 1 Direction field of dx dt t z 5 2 Frequently asked questions 69 for ODE initial value problems are essentially sophisticated ways of joining up little line segments DAEs are not like that or rather only partly like that It may help to look at this from the physical rather than the mathematical viewpoint by considering the derivation of the simple pendulum equations in FAQ 6 202 There it is clear that x and y are ODE like variables two equations define their acceleration which governs their change in velocity which governs their change in position However which is a constant factor times the force T of tension compression in the rod is physically quite different The value of T is not specified by a rate of change but rather is an instantaneous response to the current position and velocity of the bob Thus it is not surprising that the mathematics shows 2 y and entangled in that they are at any instant the solution of three simultaneous equations in which the known values of x y x y
77. re or the computation see Result This does the main work in solving the problem by advancing the integration of the DAE over a range 1 3 3 DAEsolution class A DAEsolution object x holds a point comprising the current values of the indepen dent variable t and of the x and their derivatives at t as described in 912 2 for more details see Chapter It also stores various data related to the state of the integration Constructor DAEsolution const DAEsolver amp Solver throw std logic_error Solver Input a DAEsolver object Result A DAEsolution object associated with the given Solver It is in the Initial state see start of 2 2 Before it can be used by the integrate method it must be given initial values by means of its setT and setX methods If Solver is constructed from an ill posed DAE an exception std logic_error results 1 3 Basic classes and methods 11 setT method DAEsolution amp setT double t throw std logic_error t Input scalar x setT t sets t t in a DAEsolution object Returns a reference to the updated DAEsolution object Constraint x must be in the Initial state see start of 2 2 If x is not in this state an exception std logic_error results setX method DAEsolution amp setX int index int order double value VarType type Free throw std logic_error std out_of_range index order value Input such that x setX j 1 k v sets the entry of x represen
78. roblem 52 a x b y h exp x c y d 2 v Path in velocity x y plane Figure 3 3 Golf putt plots Markers are at 0 2 sec intervals Vertical axis not to scale a Path in x y z space b Path in velocity y plane The program displays the expected structural information DAE quasilinear SIZ oe cee aha ka S 4 Indak oren ane a a ia a toe 3 TABLEAU 0 0 Ox O 2 1 2 1 1 0 2 O 1 2 1 0 31 0 1 1 2 O variable derivatives xO x1 0 1 x2 0 1 x3 0 1 The program as given uses SI units It has a planar green with an uphill slope of 1 20 in the y direction and a Gaussian circular hump of height 10cm and nominal radius 50cm centered at the point 3 1 1 m The start point o Yo 0 0 is hard 3 6 Setting an output function 53 wired in the code and the initial Zo Yo and the interval t are read from the command line The output for the graphs came from the call golfputt 2 2 1 24 2 Figure B 3 shows plots of the results which show the ball being deflected uphill as it passes near the top of the hump before it curls round and stops They strongly suggest what theory confirms that the ball s velocity is directly downhill at the moment it comes to rest The MATLAB script used to produce the plots is given below curve textscan fid n n n n n n n read rest fclose fid a b c d h v deal params 1 2 green s parameters t x y zZ
79. rresponding qe may not be set correctly For example in x setX 0 0 Free xo is set to 2 The reason is that Free is an enumerated constant with value 2 see the definition of VarType in 2 2 C converts enumerated values to int and this call sets xy to 2 printSolution method void printSolution ostream amp s cout const s Input output stream default is cout x printSolution s prints the t and z components into the output stream s t is displayed first Then for j 1 n on separate lines labeled 0 to n 1 to match C indexing the derivatives 4 are displayed k runs from 0 to dj 1 if the DAE is quasilinear and from 0 to d if not 1 3 4 SolverExitFlag type This enumerated data type provides an error indicator for the integrate method of DAEsolver 1 3 Basic classes and methods 13 Definition typedef enum SolverExitFlag success toofewdof nonconsistentpt uninitializedpt htoosmall success The integration has completed successfully toofewdof More IVs were specified fixed than the number of degrees of freedom of the DAE nonconsistentpt DAETS was unable to find an initial consistent point uninitializedpt setT is not called on the input x of integrate or setX is not called for all the required components of x htoosmall The stepsize during integration became too small to make progress printSolverExitFlag method void printSolverExitFlag SolverExitFlag state ostream amp s
80. ss reference is a hyperlink For instance you can click on the entry Obtaining DAETS above to jump to that section This also applies to page numbers in the Index and in the body of the text to chapter and section references and to equation numbers To return to where you just were use your PDF reader s Back command Chapter 1 Basic Usage DAETS is a C package for the numerical solution of initial value problems IVPs for differential algebraic equation systems DAEs This guide assumes knowledge of basic theory of DAEs as covered for instance in the texts 2 7 what DAEs are why they are in certain ways harder to solve than ordinary differential equations ODEs and the notion of index as a measure of a DAE s difficulty If you are unfamiliar with the area but still wish to read further consult first the brief explanation of DAEs in FAQ B212 and FAQ E213 Section LI of this chapter gives a short overview of the code s purpose and numer ical methods with references to fuller explanations later in this guide and in published work Section T2 introduces the DAETS user interface via a simple complete program Section L3 is a reference guide to the classes and methods used in this program The rest of the guide is organized as follows Chapter Plis a reference guide to all other user callable features Chapter B gives some more substantial applications of DAETS Chapter Mis about installation Chapter B summarizes i
81. t is cout Result Outputs into s the structure of a solution point indicating which variables and derivatives of them need to be initialized For the pendulum example this method prints POINT STRUCTURE variable derivatives x0 0 1 xi 0 1 x2 2 1 3 Accuracy control There is an accuracy tolerance tol and a type of accuracy control esttype At each step in the numerical solution an estimate st of the local error in a is computed internally for each j k in the set J as in 2 5 For the current step to be accepted the following condition must be satisfied lest 7 r lt 1 where r max J j k L 2 1 GET T X El Ta 2 1 and 7 and 7 are defined by the table esttype Tp Ta Description Abs 0 0 tol absolute accuracy control Rel tol e relative accuracy control AbsRel tol tol mixed accuracy control Here e is a modest multiple of the machine precision Error theory for T S methods is less developed for DAEs than for ODEs but there are good reasons to expect this to give a robust accuracy control If 2 1 is not satisfied the stepsize is reduced and the solution is recomputed on the current step If the user wishes to measure the error in the computed solution in terms of the number of correct decimal places then esttype should be set to Abs on entry whereas if the error requirement is in terms of the number of correct significant digits then esttype should be set to Rel For a mixed accur
82. test set double compSignificantCorrectDigits DAEsolution amp x double y 6 reference solution from the test set y O 0 1150794920661702e0 y 1 0 1203831471567715e 2 y 2 0 1611562887407974e0 y 3 0 3656156421249283e 3 y 4 0 1708010885264404e 1 y 5 0 4873531310307455e 2 Compute max norm of componentwise relative errors double error_norm 0 100 101 a a R O N H for int i 0 i lt 6 3 3 One step mode Passing data to a problem i 41 double r fabs x getX i 0 y il y i if r gt error_norm error_norm The output is return log10 error_norm r significant correct digits SOLUTION t 1 800000e 02 x x0 1 150795e 01 xi 1 203831e 03 x2 1 611563e 01 x3 3 656156e 04 x4 1 708011e 02 x5 4 873531e 03 STATISTICS TIMES 335 0 4 35 2 cho aed 02201 STEPS 4 lt cord 133 accepted ii 132 TeJected pisi pig 1 dit 0 8 STEPSIZES smallest 0 048 largest 1 752 Significant correct digits 0 gt x u 1 x 3 x x DAEso lver h file vdp cc in daets examples It solves Van Der Pol s ODE It illustrates how parameters can be passed to the DAE and how the one step mode is called include using namespace daets 2 265538e 04 1 358514e 07 1 127593e 04 1 036671e 06 1 380017e 06 ORDER 0 oie 11 TOLERANCE relative 1 0e 08 absolute 1 0e 08 8 5 3 3 One step mode Passing data to a problem We ill
83. ting i the kth derivative of the jth variable equal to v Constraint j k must be in the index set J of 2 5 Otherwise an exception std out_of_range results type Input optional default Free specifies only for the initial computation of a consistent point whether the value in question must be kept unchanged Fixed or may be altered Free see also 2 2 Constraint There are limits on how many IVs can be made Fixed This is explained in Returns A reference to the updated object Constraint x must be in the Initial state see start of 2 2 If x is not in this state an exception std logic_error results Notes Returning a reference lets one chain calls as in x setT setX For instance to specify The initial t is 1 and the IVs are x 2 fixed and z 0 x5 5 guessed the calling program could do x setT 1 0 setX 0 0 2 0 Fixed setX 1 0 0 0 setX 1 1 5 0 The following is equivalent but displays the setting of each x separately x setT 1 0 1 3 Basic classes and methods 12 x setX 0 0 2 0 Fixed x setX 1 0 0 0 x setX 1 1 5 0 A DAEsolution object knows which initial values it requires If a value is set for an unnecessary index order pair setX displays an error message to this effect If a required value is left unset integrate detects this when called and displays an error message WARNING If a value for the third parameter value is not given the co
84. ts to initialize Suppose as in this problem you know consistent initial values or reasonable guesses of them for some required components here x for j 0 to 4 counting from 0 as in the main program The program should set these values suitably and set the others here x5 and xq to x4 to an arbitrary value such as zero See setInitialValues at line 73 in the code below If you want to inspect the consistent point that results call integrate with teng equal to the initial ty This computes a consistent point and exits The result here is SOLUTION t 0 000000e 00 x x0 4 440000e 01 5 xi 1 230000e 03 x2 0 000000e 00 2 x3 7 000000e 03 3 097682e 02 372932e 02 548743e 02 916080e 06 If on inspection this appears to be correct you can continue the integration to the desired tena The results are shown in 48 2 below 3 2 Chemical Akzo Nobel problem 39 3 2 Chemical Akzo Nobel problem This IVP is a stiff non linear DAE of size 6 and index 1 It describes a chemical process in which 2 species are mixed while carbon dioxide is continuously added to produce a product It is taken from the Test Set for IVP Solvers 12 which describes more fully its formulation and origin We give a basic program to integrate this problem and check the accuracy of the solution by computing the number of significant correct digits as defined in J A useful coding technique is shown The Test Set supplies definitions
85. ubsequent call to integrate on x unless reset to zero by a call to setFirstEntry getCPUtime method double getCPUtime const Returns user CPU time in seconds used to integrate a DAEsolution object This is obtained by a call to the Fortran etime utility which accesses the high resolution clock and is typically accurate to a microsecond or better 2 2 Solution path properties the DAEsolution class 31 getNumAccSteps method int getNumAccSteps const Returns number of accepted steps getNumRejSteps method int getNumRejSteps const Returns number of rejected steps 2 2 6 Printing printSolutionState method void printSolutionState ostream amp s cout const s Input output stream default is cout Result x printSolutionState s outputs into s the current state of x along the solution path which is one of the values of SolutionState printStats method void printStats ostream amp s cout const s Input output stream default is cout Result x printStats s outputs into s the following data related to the integra tion of x Where relevant these are cumulative as in 2 2 5 user CPU time total number of steps number of accepted and number of rejected steps values of smallest and largest stepsizes current order of the Taylor series current tolerances 7 and Ta see equation CD 2 3 Arithmetic on solutions the DAEpoint class 32 2 3 Arithmetic on solutions the DAEpoint cl
86. ustrate how parameters are passed to the function for evaluating the DAE in DAETS Here we integrate the Van Der Pol s equation 3 1 with x x 3 2 at to 0 values for parameter u 1 10 and 100 and teng 20 60 and 500 respectively We also show how the one step mode is called We report the solution on each step into files and plot x versus t and 2 versus 2 We list and explain the DAETS program for integrating this problem 3 3 One step mode Passing data to a problem 42 template lt typename T gt void fcn T t const T x T f void param Convert the content at param to double double mu double param 0 20 u 1 23 0 x f 0 Diff x 0 2 mu 1 sqr x 0 Diff x 0 1 x 0 MAIN PROGRAM int main double tO 0 double mul t1 lei 1 e2 double tend 20 60 500 double hmax 0 05 0 005 0 001 Create a solver and pass the address of mu 0 as last parameter DAEsolver Solver 1 DAE_FCN fcn mu Create DAEpoint object and store initial point DAEpoint xp Solver xp setX 0 0 3 setX 0 1 2 Create DAEsolution object DAEsolution x Solver Set one step mode x setOneStepMode true Monitor the integration x setPrintProgress SolverExitFlag flag success for int i 0 i lt 3 i printf Integrating VDP with mu 1e tend 1f n mui tendl il if i gt 0 Solver passDataToDAE
87. ution path properties the DAEsolution class 28 2 2 2 Copying No assignment operator y x is provided for this class because copying the entire contents of x to y could compromise the locking security feature if either object is not in the Initial state However one often wants to use the X contents of an object to initialize another object The updatePoint method offers this without bypassing the security feature updatePoint method DAEsolution amp updatePoint const DAEpoint amp x throw std logic_error x Input a DAEpoint or a DAEsolution object Returns y updatePoint x has the effect of copying the X part of x into that of y If x is a DAEpoint object e if the type of a component of y is Uninitialized it becomes Free e otherwise the type of a component of y remains unchanged If x is a DAEsolution object the type of each element is also copied Constraints y must be in the Initial state else an exception std logic_error results x and y must have the same shape have the same index set J in 2 5 else an exception std logic_error results It is not necessary that x and y were constructed from the same Solver This feature is convenient for copying the elements of initialized objects For ex ample if an integration proves problematic near some t t and one wishes to experiment with the accuracy control or other parameter settings around this point one can do DAEsolution x Solver initia
88. ven So a variable whose offset is zero needs no IVs e Otherwise IVs for each x and its derivatives through the d th must be given Example The pendulum DAE is quasilinear and x y A have offsets 2 2 0 so IVs must be given for z x y y If it is changed to be not quasilinear then DAETS recognizes this and IVs must be given for x x y y y A One needs to know what the definition of consistent values as follows e If the DAE is quasilinear the IVs must satisfy each f 0 and its derivatives up to but not including the c th So a function whose offset is zero imposes no constraint e Otherwise the IVs must satisfy each f 0 and its derivatives up to and including the cith 1 3 Basic classes and methods 7 Example For the pendulum DAE f g h have offsets 0 0 2 so the consistency equations for x x y y are 0 h h that is 0 x y L and 0 2xx 2yy If it is changed to be not quasilinear the enlarged set of variables x x x y y y A must satisfy the extra consistency equations 0 f 0 g and 0 h 2 xx x 2 yy y As a consequence the numerical solution values held in a solution object such as x in the example program are not a flat vector as they are with an ODE DAETS uses an irregular array such as for the pendulum 2 x 0 1 1 o o 7 if not y y 7 if quasilinear or 2 2 Method setX 9
89. viewed as a point moving along the solution path DAE_FCN in line 5 is a macro explained in 1 3 2 x records the current values not only of z x y A the contents of z in fcn but also of the independent variable t Thus if as here t represents time and z is loosely thought of as space then x represents a point z t in space time Lines 7HI9 give x its IVs Method setT initializes the independent variable com ponent and setX initializes the state variable components including their deriva tives We use the name setX because the state variables are called x in the general V 0 JD 0 A OU N RP N N N N N N N N H H H H H H H H H H 4 A N H 0000 N H O 1 2 Quick start guide 4 file pendulumsimple cc in daets examples include DAEsolver h using namespace daets template lt typename T gt void fcn T t const T z T f void param f gzl z 11 z 2 are x y A const double G 9 8 L 10 0 f O Diff z 0 2 z 0 z 2 f 1 Diff z 1 2 z 1 z 2 G f 2 sqr z 0 sqr z 0 sqr z 1 sqr L int main const int n 8 size of the problem double tO 0 0 tend 100 0 DAEsolver Solver n DAE_FCN fcn create a solver analyse DAE DAEsolution x Solver create a DAE solution object x setT t0 initial value of t setX 0 0 10 0 setX 1 0 0 0 i and of x and y setX 0 1 0 0 setX 1 1 1 0 and of 2 and y SolverExitFlag flag Solver i
90. xp yp zp deal curvef curve data figure 3 plot xp yp o LineWidth 1 1 set gca FontSize 20 FontWeight bold YTick 0 0 5 2 axis equal title Path in velocity x y plane FontSize 18 FontWeight bold print depsc golfputt3 eps figure 4 clf hold on set gca FontSize 20 FontWeight bold CameraPosition 24 32 2 2 DataAspectRatio 15 15 2 colormap autumn wid 3 sqrt v fgreen 0 x y a x bey h exp x c 72 y d 72 2 v ezsurf fgreen min x c wid max x c wid min y d wid max y d wid 250 plot3 x y fgreen x y 0 001 bo LineWidth 1 1 print depsc golfputt4 eps 3 6 Setting an output function We illustrate how an output function is given to integrate so it can output solution data into a file without calling the one step mode The example from Figure L3 is augmented in the listing of FigureB 4 where line B3 sets the output function outfcn defined in lines B2H39 On each successful step integrate passes to this function the accepted solution and the void parameters Here we output t x y y A A Plots of these variables and derivatives versus t are given in Figure 6 8 3 7 Summary This chapter has aimed to show typical techniques in the application of DAETS The authors welcome suggestions for improving these examples and contributions t

Daets User Guide - Computing and Software

Contents

Download Pdf Manuals

Related Search

Related Contents