Numerical computation of normal form coefficients of bifurcations of odes in MATLAB

Normal form coefficients of codim-1 and codim-2 bifurcations of 
equilibria of ODEs are important since their sign and size determine the bifurcation 
scenario near the bifurcation points. Multilinear forms with derivatives 
up to the fifth order are needed in these coefficients. So far, in the Matlab 
bifurcation software MatCont for ODEs, these derivatives are computed either 
by finite differences or by symbolic differentiation. However, both approaches 
have disadvantages. Finite differences do not usually have the required accuracy 
and for symbolic differentiation the Matlab Symbolic Toolbox is needed. 
Automatic differentiation is an alternative since this technique is as accurate as 
symbolic derivatives and no extra software is needed. In this paper, we discuss 
the pros and cons of these three kinds of differentiation in a specific context 
by the use of several examples.


1.
Introduction. The computation of the normal form coefficients of equilibrium bifurcations of dynamical systems requires multilinear forms of order up to five in order to compute directional derivatives. Derivatives of first or second order can be approximated by finite difference methods reasonably well if sophisticated heuristics are used; for higher order derivatives this is a hopeless task. Calculation of the multilinear forms by use of symbolic derivatives (SD) is not always possible in Matlab since the Symbolic Toolbox is not always available. A third option for the computation of the directional derivatives is automatic differentiation (AD). Every function f expressible by a computer program is built from a finite set of elementary functions (this includes the basic arithmetic operations). The basic principle of automatic differentiation, see Griewank [10] and Griewank and Walther [11], is to use the known formulae for differentiating elementary functions together with the chain rule, to build up the needed derivatives of an arbitrary function.
We assume f is a vector function y = f (x) with n inputs, or independent variables, x = (x 1 , . . . , x n ) and m outputs y = (y 1 , . . . , y m ). The code for f may contain branches and loops. However, each evaluation of f at given inputs x can be written as a code list, which is a finite sequence of assignments of the simple form v i = e i (previously defined v j 's, or constants), i = 1, 2, . . . , l+m, where each e i is one of the elementary functions. The v i are called variables. In (1) it is convenient to use v 1−n , . . . , v 0 as aliases for the inputs x 1 , . . . , x n and v l+1 , . . . , v l+m as aliases for the outputs y 1 , . . . , y m , following the notation of [10]. The remaining variables v 1 , . . . , v l are called intermediate. The forward mode of automatic differentiation is the simplest, and is appropriate to our application. Each v i is represented by an object v i of a data type that holds not just the value but some needed set of derivatives. For our purposes: • There is (at any point) just one variable being treated as independent: we call it t. Thus each variable v i is regarded as a function of t. • The data structure holds some Taylor coefficients, namely the coefficients of the truncated Taylor series of v i , up to some order p, expanded about some point t = a. When we change to a new independent variable s = t − a, we get We call the data type adtayl. It is helpful to think of an adtayl object as representing an infinite power series which is only known up to the order p term. Details about the implementation are provided in [13]; here we just give a few basic ideas. When evaluating a code list, each elementary operation on real arguments is replaced by the corresponding operation on adtayl arguments. Consider first the four basic arithmetic operations. Suppose that a holds (a 0 , a 1 . . .), and so on for other named variables. Let a be defined to order p, and b to order q. Then c = a + b and d = a × b are defined to order r = min(p, q) by and similarly for a − b, and for a ÷ b provided b 0 = 0. For applying the standard functions exp, cos, . . . to power series, there are various formulas in the literature. We have aimed to choose ones that can be made reasonably fast in Matlab, especially when the argument is a vector of power series, not just a single one. We give a few examples.
First, multiplication of power series is a convolution, which can be realised by the very fast built-in filter function of Matlab.
Second, the method for exp(a) exploits the fact that if c(t) = exp(a(t)) then c (t) = a (t)c(t), which is converted to the integral form c(t) = c 0 + t 0 a (s)c(s)ds where c 0 = e a0 . In terms of the coefficients this reduces to a triangular linear system, which is fast in Matlab.
Third, for sin(a) and cos(a), it is convenient to compute both simultaneously as the real and imaginary parts of exp(i a), and to record both results, along with the argument a, in persistent storage. Each time a new a is given, it is checked against the recorded one. If they are equal, the result can be returned at once. Since cos and sin at the same argument often occur together in applications, this reduces the cost for these functions by nearly half.
Experimentally, we found in our Matlab implementation that the run-time for computing a Taylor series up to order p is approximately proportional to p 2.7 . Since this involves the computation of all coefficients up to order p, this is compatible with the known fact (see e.g. [11], Section 13.1) that the complexity of computing the d-th coefficient is proportial to d 2 .
Regarding previous work to provide AD facilities in Matlab, Rich and Hill [14] worked out a limited facility that enabled AD of simple expressions defined by a character string. The first significant work was that of Coleman and Verma  [3,4], who produced an operator-overloading AD package ADMAT with facilities for forward and reverse mode for first and second derivatives, and run-time Jacobian sparsity detection. These authors interfaced ADMAT with ADMIT [5], a package for efficient sparse Jacobian calculation, and recently the ADiMat hybrid source-transformation/operator-overloading AD tool [15] has been developed. Comparisons [2] have shown its forward mode to be more efficient than that of ADMAT. Bischof et al. [1] proposed a novel software tool which transforms a given Matlab program into another Matlab program to change the semantics of the program in a fashion based on automatic differentation. The crucial ingredient of the tool is a combination of source-to-source transformation and operator overloading. For first derivatives, currently the most efficient and comprehensive tool is probably Forth's MAD package [8]. However, none of these tools handles Taylor series, although the subject is treated in [11], Chapter 13. MatCont is a Matlab numerical bifurcation toolbox for ODEs [6,7]. So far, AD is not provided in MatCont. But it was already introduced in CL_MatContM, see [13], the pendant to MatCont for iterated maps. However, in CL_MatContM the computations of the multilinear forms A(q 1 ), B(q 1 , q 2 ), C(q 1 , q 2 , q 3 ), D(q 1 , q 2 , q 3 , q 4 ), E(q 1 , q 2 , q 3 , q 4 , q 5 ), where A is the Jacobian, B the tensor of second derivatives etc., were only implemented in the case of real vectors q i . We have extended this to the case of complex vectors. In [13] is shown that SD and AD have the same accuracy. Moreover, in the case of iterated maps, AD is faster for high values of the iteration number J. When using AD, the elapsed time grows linearly with the iteration number J, but for SD it grows more like J d for d-th derivatives.
In Section 2 we give an overview of the maximal order of derivatives needed in the computation of the normal form coefficients of equilibrium bifurcations of codimension 1 and 2 and we motivate our choice of examples. In Section 3 we investigate three examples to demonstrate the advantages and disadvantages of the three kinds of derivatives. We give a comparison of time complexity and accuracy using finite differences (FD), symbolic derivatives and automatic differentiation. In Section 4 we draw some conclusions.
2. Preliminary remarks. In Table 1 we list all codimension 1 and 2 bifurcations of equilibria of ODEs and the corresponding maximum order of derivative appearing in the normal form coefficient. The most complicated normal form coefficients are those of the Generalized Hopf points and the Hopf-Hopf points, since these computations require the evaluation of several multilinear forms of order 4 and 5.
Therefore, we will consider in particular the AB-neuron model and the Morris Lecar example, which contain several of these bifurcation points.
We note that the Jacobian is needed during the continuation of a curve. AD could also be used for these computations. However, a very accurate Jacobian is not needed since errors in the Jacobian can be remedied by the Newton corrections. Therefore, in this case it is preferable to avoid AD since in our particular implementation AD is more time-consuming, as we will see later on, and there is a huge amount of Jacobian computations during the continuation. This would slow down the speed of the continuation unnecessarily. So the options of SD and FD are in this case sufficient.
Symmetry is not taken into account for the computation of the multilinear forms with either SD or AD. In the case of SD, we note that the speed of Matlab depends strongly on the way in which the program deals with vectors. Matlab is designed to perform vector and matrix operations efficiently. By making use of the symmetry this advantage would be lost. When using AD, the matrix of derivatives itself is never computed but only the product with certain vectors. Therefore, it is not quite clear how symmetry could be dealt with (see also [11], Section 13.3).
3. Test Cases. In this section we compare SD, AD and FD in a number of test cases. In the case of FD a multilinear form of order n is computed by a symmetric finite difference formula with O(h 2 ) accuracy that uses function evaluations in n + 1 points of the form x+ihq where x is the base point, q is the unit directional vector, i is an integer, h = ( m ) 1 n+2 is the discretization step and m is the machine precision.
3.1. AB-neuron model. A good example to compare the three kinds of differentiation is the AB-neuron model. The equations and fixed parameter values of this model can be found in [9]. It is a highly complicated nonlinear 6-dimensional system where gA and gKCa are the parameters of interest to us. This model contains several GH-points and two HH-points. The investigated bifurcation points are shown in Figure 1. There is also a BT-point outside the range, with gKCa = −0.1240 and gA = 263.3610.
The first column of Table 2 shows the bifurcation, the others give the elapsed time in four cases. The second column is the case where the bifurcation point and the normal form coefficient are computed symbolically, for the third the bifurcation point is computed with SD and the normal form coefficient is computed with AD, in the fourth and fifth column the bifurcation point is computed with FD, with the normal form coefficient computed with AD in the fourth column and with FD in the last column.
Remarkable is the fact that in all cases the elapsed time needed for AD is much higher than for SD. The elapsed times for SD don't include the generation of the symbolic derivatives by Matlab, but are the precise computation times for the specific normal form coefficients. Especially the computation of the HH-points requires much more calculation time, the elapsed time needed for AD is 600 to 900 times that for SD. In the AD case the computation takes about 4 minutes, but be aware that the normal form coefficient of a HH-point consists of 5 coefficients. For the GH-points, SD is about 300 times faster than AD. The elapsed time for a H-point is with SD 400 to 500 times faster than with AD, though the normal form coefficient of a H-point doesn't require complex calculations. In the other cases the symbolic derivatives also give faster results, but the ratios are not that impressive.  In fact, the elapsed times when using AD are in all cases still acceptable (less or around one second), except in the GH-and HH-cases.
The above results are of course very dependent on the processor speed and computing power. We guess that the high time complexity in the AD case is due to the overhead caused by the implementation of Taylor series as an object type in Matlab rather than to the number of floating point computations.  Table 4. Relative distance between the bifurcation points computed by SD and FD, ABneuron model.
Another remarkable fact is that the elapsed time in SD cases seems to be erratic while in the AD cases the computation times are more consistent. Table 3 gives an indication of the accuracy of the computed normal form coefficients. We take the normal form coefficient in the second column where both bifurcation point and normal form coefficient are computed symbolically as reference value. The last three columns give the relative error compared to this case. We recall that a BT point consists of 2 normal form coefficients, a ZH point of three and a HH point of 5 coefficients.
The normal form coefficients computed by SD and AD with symbolically computed bifurcation point are equal to within machine precision. Indeed, in the calculation of the normal form coefficients, systems are solved with quite high condition number (of the order 1e + 5). This explains why the results of the second column case are not identical to the results of the third column case.
When the bifurcation point is computed with FD and the normal form coefficient with AD, the results are still acceptable, taking the relative distance between the symbolically and the with FD computed bifurcation points into account, see Table  4. We note that the bifurcation point value computed by SD is not necessarily more accurate than the one obtained by FD. However, when the normal form coefficient is computed with FD (last column in Table 3), very inaccurate results are obtained. In the computation of the second GH-point not even one digit is guaranteed.
Since the computation of the fourth and fifth order derivatives is found to cause a main slow-down in computation time, it is interesting to look at the results of Table 5. Here we have considered the case where the first to third order derivatives are computed with AD, and the fourth and fifth order derivatives by the use of first and second order finite differences on the third order tensor. The bifurcation point is computed with FD. The elapsed time given in the second column in the case of the two HH-points is the computation time for the 5 normal form coefficients together. The elapsed times should be compared with the results from the last but one column of Table 2. We obtain a time improvement of 2.5 times, both for the GH-and the HH-points. However, in comparison with a symbolically computed normal form coefficient, this is still much slower. The last column of Table 5 lists the relative error in comparison with the case where AD is used for all the derivatives. The results for the GH-points are acceptable, but for the HH-points big relative errors are obtained. So the gain of calculation time doesn't compensate for the loss of accuracy. Therefore, this is not an option to be implemented in MatCont.   [12], a quite complex system with only two state variables: The investigated bifurcations are presented in Figure 2.   Table 6 lists the elapsed time in the four interesting cases. Again, the elapsed time needed for the computation of the normal form coefficients with SD is much shorter than with AD. For the GH-points, SD is 200 to 1700 times faster than AD. Although only two multilinear forms of order 2 are required for the BT-points, SD is 4 to 70 times faster than AD. We again note the more consistent elapsed times for SD. From Table 7 we can derive analogous conclusions concerning the accuracy as for the AB-neuron model. Despite the high condition number of the linear systems to be solved (again order 1e + 5), the results of the normal form coefficients with symbolically computed bifurcation point and normal form coefficient computed either with SD or AD are very similar. Moreover, the multilinear forms are identical when using symbolic or automatic derivatives. The obtained relative errors of the last but one column are acceptable, taking the relative distances between the bifurcation points into account, which are of the same order as in the AB-neuron model. We obtain inaccurate results when using finite differences.
3.3. The Brusselator. Our last example is the 1-dimensional Brusselator: We discretize this system using an equidistant mesh of n intermediate points in the interval [0, 1], with the length of the mesh intervals equal to h = 1 n+1 and with function values X 1 , Y 1 , . . . , X n , Y n . At the beginning point and end point we take X equal to A and Y equal to B/A. The discretization is then given bẏ The point (X i , Y i ) = (A, B/A), i = 1, . . . n, is an equilibrium independently of the value of L. The parametervalues are D X = 10, D Y = 1, A = 2, B = 6 and L = 1. We start the continuation from this equilibrium, with L as free system parameter.
The original system (2) has a Hopf point for L = 10.4194808 = π DX +DY B−1−A 2 . We detect the Hopf point of the discretizations in (3) for several values of n. But for n = 25, MatCont is not able to compute the normal form coefficient by use of symbolic derivatives, an "Out of memory error" is given. With AD, there are no problems. So AD then forms the only option for the computation of the normal form coefficient of the Hopf point (FD is inaccurate). For n = 20, MatCont is able to compute the symbolic derivatives in the Symbolic Toolbox up to order three, but not up to order five, although the fourth and fifth derivatives are identically zero. So in this case restrictions of the software force us to use automatic differentiation. As an alternative one might consider to introduce sparse symbolic derivatives in MatCont (but we note that Matlab does not support sparse multidimensional arrays).

4.
Conclusion. An accurate computation of the normal form coefficients of bifurcations of ODEs is required since these coefficients are used in the detection of higher codimension bifurcations, in switching to branches of lower codimension objects etc. As supported by the examples, SD and AD are equally accurate, while FD give often inaccurate results. Therefore, when the Symbolic Toolbox is not available or, when this toolbox is available but MatCont is not capable to store the symbolic derivatives, automatic differentiation is preferable to finite differences. However, the computation time of the normal form coeficients of GH-and HH-points is often prohibitive. Therefore, we propose that in these two cases the user is given the option to choose whether the normal form coefficient is computed or not. The normal form coefficients of the other equilibrium bifurcations should always be computed and AD should be used if SD is not available.