On a Newton-Type Method for Differential-Algebraic Equations

Abstract

This paper deals with the approximation of systems of differential-algebraic equations based on a certain error functional naturally associated with the system. In seeking to minimize the error, by using standard descent schemes, the procedure can never get stuck in local minima but will always and steadily decrease the error until getting to the solution sought. Starting with an initial approximation to the solution, we improve it by adding the solution of some associated linear problems, in such a way that the error is significantly decreased. Some numerical examples are presented to illustrate the main theoretical conclusions. We should mention that we have already explored, in some previous papers (Amat et al., in press, Amat and Pedregal, 2009, and Pedregal, 2010), this point of view for regular problems. However, the main hypotheses in these papers ask for some requirements that essentially rule out the application to singular problems. We are also preparing a much more ambitious perspective for the theoretical analysis of nonlinear DAEs based on this same approach.

1. Introduction

Differential-algebraic equations are becoming increasingly important in a lot of technical areas. They are currently the standard modeling concept in many applications such as circuit simulation, multibody dynamics, and chemical process engineering; see for instance [1–5] with no attempt to be exhaustive.

A basic concept in the analysis of differential-algebraic equations is the index. The notion of index is used in the theory of DAEs for measuring the distance from a DAE to its related ODE. The higher the index of a DAE is, the more difficulties are for its numerical solution. There are different index definitions, but for simple problems they are identical. On more complicated nonlinear and fully implicit systems they can be different (see [5] and the references therein.)

For simplicity, we focus our attention on problems of the form 𝑀𝑥^′(𝑡)=𝑓(𝑥(𝑡))in(0,𝑇),𝑥(0)=𝑥₀, (1.1) where 𝑀 is a given, eventually singular, matrix depending on 𝑡. More general situations can be allowed. This type of equations arises, for instance, in the functional analytic formulation of the initial value problem for the Stokes as well as for the linearized Navier-Stokes or Oseen equations [6].

For the approximation of these equations collocation-type methods are usually used. These methods are implicit, and we need to solve a nonlinear system of equations in each iteration using a Newton’s type method. Given different coefficients 𝑐_𝑖, 1≤𝑖≤𝑠, there is a (unique for ℎ sufficiently small) polynomial of collocation 𝑞(𝑡) of degree less than or equal to 𝑠 such that 𝑞𝑡₀=𝑦₀,𝑞^′𝑡₀+𝑐_𝑖ℎ=𝑓𝑡₀+𝑐_𝑖ℎ,𝑞𝑡₀+𝑐_𝑖ℎif1≤𝑖≤𝑠. (1.2) The collocation methods are defined by an approximation 𝑦(𝑡)≃𝑞(𝑡) and are equivalent to implicit RK methods of 𝑠 stages 𝑘_𝑖=𝑓𝑡₀+𝑐_𝑖ℎ,𝑦₀+ℎ^𝑠_𝑗=1𝑎_𝑖,𝑗𝑘_𝑗,𝑦₁=𝑦₀+ℎ^𝑠_𝑖=1𝑏_𝑖𝑘_𝑖, (1.3) for the coefficients 𝑎_𝑖,𝑗=^𝑐𝑖_0𝑙≠𝑗𝑢−𝑐_𝑙𝑐_𝑗−𝑐_𝑙𝑏𝑑𝑢,_𝑖=¹_0𝑙≠𝑖𝑢−𝑐_𝑙𝑐_𝑖−𝑐_𝑙𝑑𝑢. (1.4) The coefficients 𝑐_𝑖 play the role of the nodes of the quadrature formula, and the associated coefficients 𝑏_𝑖 are analogous to the weights. From (1.4), we can find implicit RK methods called Gauss of order 2𝑠, Radau IA and Radau IIA of order 2𝑠−1, and Lobatto IIIA of order 2𝑠−2. Also we can consider perturbed collocation methods like Lobato IIIC. (See [4] for more details).

A number of convergence results have been derived for these methods introducing the so-called 𝐵-convergence theory. In [7, 8] the authors extend the 𝐵-convergence theory to be valid for a class of nonautonomous weakly nonlinear stiff systems, in particular, including the linear case. As pointed out by the same authors, it is not clear if it is possible to cover, in a satisfactory way, highly nonlinear stiff problems, that is, problems where also the nonlinear terms are affected by large parameters. Moreover, any result should assume that, in each step, the associated nonlinear system is well approximated [5]. In particular, we should be able to start with a good initial guess for the iterative scheme. This might be very restrictive for many stiff problems.

On the other hand, iterative methods are the typical tool to solve nonlinear systems of equations. In these schemes we compute a sequence of approximations solving associated linear problems. In this paper, we would like to introduce a new variational approach for the treatment of DAEs where we linearize the original equations obtaining an iterative scheme. Our ideas are based on the analysis of a certain error functional of the form 𝐸(𝑥)=12^𝑇₀𝑀𝑥^′(𝑡)−𝑓(𝑥(𝑡))²𝑑𝑥 (1.5) to be minimized among the absolutely continuous paths 𝑥:(0,𝑇)→𝐑^𝑁 with 𝑥(0)=𝑥₀. Note that if 𝐸(𝑥) is finite for one such path 𝑥, then automatically 𝑀𝑥^′ is square integrable. This error functional is associated, in a natural way, with the Cauchy problem (1.1). Indeed, the existence of solutions for (1.1) is equivalent to the existence of minimizers for 𝐸 with vanishing minimum value. This is elementary.

In this initial contribution, we want to concentrate on the approximation issue through this perspective. We will place ourselves under the appropriate hypotheses so that there are indeed solutions for (1.1), that is, there are minimizers for the error with vanishing minimum value. In addition, we would like to guarantee that the main ingredients for the iterative approximating scheme to work are valid. More explicitly, our approach for the numerical approximation of such problems relies on three main analytical hypotheses that we take for granted here. (1) The Cauchy problem (1.1) admits a unique solution for every feasible initial condition 𝑥₀ (the definition of feasible path should depend on the index of the equation). (2) The linearization around any feasible, absolutely continuous, path 𝑥(𝑡) with 𝑥(0)=𝑥₀, 𝑀𝑦^′(𝑡)−∇𝑓(𝑥(𝑡))𝑦(𝑡)=𝑓(𝑥(𝑡))−𝑀𝑥^′(𝑡)in(0,𝑇),𝑦(0)=0, (1.6) always has a unique solution, and moreover, for some constant 𝐿>0 depending on 𝑀, 𝑓, 𝑥 and its derivatives, 𝑦²_{𝐿∞()0,𝑇}≤𝑇𝐿𝑓(𝑥(𝑡))−𝑀𝑥^′(𝑡)_{𝐿2(0,𝑇)}. (1.7)(3)The only solution of the problem 𝑑𝑀𝑑𝑡^⊤𝑧(𝑡)+∇𝑓(𝑥(𝑡))^⊤𝑧(𝑡)=0in(0,𝑇),𝑀^⊤𝑧(𝑇)=0, (1.8) is 𝑧≡0, for every feasible, absolutely continuous, path 𝑥(𝑡) with 𝑥(0)=𝑥₀.Here the superscript ⊤ indicates transpose.

These requirements depend on the index of the equation and on some regularity on the pair (𝑀,∇𝑓(𝑥(𝑡))). They should be more restrictive for equations with high index. More details can be found, for example, in [9, Theorem 3.9], where the authors consider DAEs transferable into standard canonical form. More precise information is outside of the scope of this paper. In any case, the equations verifying our hypotheses are, in general, a subclass of all analytically solvable systems.

In addition to the basic facts just stated on existence and uniqueness of solutions for our problems, the analysis of the approximation scheme, based on a minimization of the error functional 𝐸, requires one main basic assumption on the nonlinearity 𝑓:𝐑^𝑁→𝐑^𝑁. It must be smooth, so that ∇𝑓:𝐑^𝑁→𝐑^𝑁×𝑁 is continuous and globally Lipschitz with constant 𝐾>0 (|∇𝑓|≤𝐾). Moreover, the main result of this paper demands a further special property on the map 𝑓: for every positive 𝐶>0 and small 𝜖>0, there is 𝐷_𝐶,𝜖>0 so that 𝑓𝑥+𝑦−𝑓(𝑥)−∇𝑓(𝑥)𝑦≤𝐷_𝐶,𝜖𝑦²,|𝑥|≤𝐶,𝑦≤𝜖. (1.9) This regularity is somehow not surprising as our approach here is based on regularity and optimality. On the other hand, that regularity holds for most of the important problems in applications. It certainly does in all numerical tests performed in this work. Our goal here is placed on the fact that this optimization strategy may be utilized to set up approximation schemes based on the minimization of the error functional. Indeed, we provide a solid basis for this approximation procedure. One very important and appealing property of our approach states that typical minimization schemes like (steepest) descent methods will work fine as they can never get stuck in local minima and converge steadily to the solution of the problem, no matter what the initialization is.

We should mention that we have already explored, in some previous papers, this point of view. Since the initial contribution [10], we have also treated the reverse mechanism of using first discretization and then optimality [11]. The perspective of going through optimality and then discretization has already been indicated and studied in [12], though only for the steepest descent method, and without going through any further analytical foundation for the numerical procedure. However, the main hypotheses in these papers ask for some requirements that essentially rule out the application to singular problems. We will however address shortly [13] a complete treatment of DAEs with no a priori assumptions on existence and uniqueness. Rather, we will be interested in showing existence and uniqueness from scratch by examining the fundamental properties of the error functional 𝐸.

On the other hand, variational methods have been used also before in the context of ODEs. See [14, 15], where numerical integration algorithms for finite-dimensional mechanical systems that are based on discrete variational principles are proposed and studied. This is one approach to deriving and studying symplectic integrators. The starting point is Hamilton’s principle and its direct discretization. In those references, some fundamental numerical methods are presented from that variational viewpoint where the model plays a prominent role.

The rest of the paper is divided in three sections. In Section 2 we introduce our variational approach and develop a convergence analysis. Section 3 introduces the numerical procedure. We analyze some numerical results in Section 4. Finally, we present the main conclusions in Section 5.

2. A Main Descent Procedure

We start with a fundamental fact for our approach.

Proposition 2.1. Let 𝑥 be a critical point for the error 𝐸. Then 𝑥 is the solution of the Cauchy problem (1.1).

Proof. The proof is elementary. Based on the smoothness and bounds assumed on the mapping 𝑓, we conclude that if 𝑥≡𝑥 is a critical point for the error 𝐸, then 𝑥 ought to be a solution of the problem −𝑑𝑀𝑑𝑡^⊤𝑀𝑥^′(𝑡)−𝑓(𝑥(𝑡))−∇𝑓(𝑥(𝑡))^⊤𝑀𝑥^′(𝑡)(𝑥(𝑡(),𝑥(0)−𝑓))=0in0,𝑇=𝑥₀,𝑀^⊤𝑀𝑥^′(𝑇)−𝑓(𝑥(𝑇))=0. (2.1) The vector-valued map 𝑦(𝑡)=𝑀𝑥^′(𝑡)−𝑓(𝑥(𝑡)) is then a solution of the linear, nondegenerate problem 𝑀^⊤𝑦^′(𝑡)+∇𝑓(𝑥(𝑡))^⊤𝑦(𝑡)=0in(0,𝑇),𝑀^⊤𝑦(𝑇)=0. (2.2) The only solution of this problem, by our initial conditions on uniqueness of linearizations, is 𝑦≡0, and so 𝑥 is the solution of our Cauchy problem.

On the other hand, suppose we start with an initial crude approximation 𝑥⁽⁰⁾ to the solution of our basic problem (1.1). We could take 𝑥⁽⁰⁾=𝑥₀ for all 𝑡 or 𝑥⁽⁰⁾(𝑡)=𝑥₀+𝑡𝑓(𝑥₀). We would like to improve this approximation in such a way that the error is significantly decreased. We have already pointed out that descent methods can never get stuck on anything but the solution of the problem, under global lipschitzianity hypotheses.

It is straightforward to find the Gâteaux derivative of 𝐸 at a given feasible 𝑥 in the direction 𝑦 with 𝑦(0)=0. Namely 𝐸^′(𝑥)𝑦=^𝑇₀𝑀𝑥^′(𝑡)−𝑓(𝑥(𝑡))⋅𝑀𝑦^′(𝑡)−∇𝑓(𝑥(𝑡))𝑦(𝑡)𝑑𝑡. (2.3) This expression suggests a main possibility to select 𝑦 from. Choose 𝑦 such that 𝑀𝑦^′(𝑡)−∇𝑓(𝑥(𝑡))𝑦(𝑡)=𝑓(𝑥(𝑡))−𝑀𝑥^′(𝑡)in(0,𝑇),𝑦(0)=0. (2.4) In this way, it is clear that 𝐸^′(𝑥)𝑦=−2𝐸(𝑥), and so the (local) decrease of the error is of the size 𝐸(𝑥). Finding 𝑦 requires solving the above linear problem which is assumed to have a unique solution by our main hypotheses in the introduction. In some sense, this is like a Newton method with global convergence.

Suppose 𝑥⁽⁰⁾ is a feasible path in the interval (0,𝑇) so that 𝑥⁽⁰⁾(0)=𝑥₀, 𝑀(𝑥⁽⁰⁾)^′ is square integrable, |𝑥⁽⁰⁾(𝑡)|≤𝐶 for a fixed constant 𝐶, all 𝑡∈(0,𝑇), and the quantity 𝐸𝑥⁽⁰⁾=12^𝑇₀𝑀𝑥^(0)′(𝑡)−𝑓𝑥⁽⁰⁾(𝑡)²𝑑𝑡 (2.5) measures how far such 𝑥⁽⁰⁾ is from being a solution of our problem.

Theorem 2.2. For 𝑇 sufficiently small, the iterative procedure 𝑥^(𝑗)=𝑥^(𝑗−1)+𝑦^(𝑗), starting from the above feasible 𝑥⁽⁰⁾ and defining 𝑦^(𝑗) as the solution of the linear problem 𝑀𝑦^(𝑗)′(𝑡)−∇𝑓𝑥^(𝑗−1)(𝑡)𝑦^(𝑗)(𝑡)=𝑓𝑥^(𝑗−1)(𝑡)−𝑀𝑥^{(𝑗−1)′}(𝑡)𝑖𝑛(0,𝑇),𝑦^(𝑗)(0)=0, (2.6) converges strongly in 𝐿^∞(0,𝑇) to the solution of (1.1).

Proof. Choose 𝜖>0 and 0<𝛼<1 so that 𝜖1−√𝛼𝑓(𝑧)(𝑧)𝑦𝑦≤𝐶,𝑧+𝑦−𝑓−∇𝑓≤𝐷²,𝑦≤𝜖,|𝑧|≤2𝐶, (2.7) for some constant 𝐷>0 (see the main hypotheses in Section 1). We then solve for 𝑦⁽⁰⁾ as the solution of the nonautonomous linear problem 𝑀𝑦^′(𝑡)−∇𝑓(𝑥(𝑡))𝑦(𝑡)=𝑓(𝑥(𝑡))−𝑀𝑥^′(𝑡)in(0,𝑇),𝑦(0)=0, (2.8) and pretend to update 𝑥⁽⁰⁾ to 𝑥⁽⁰⁾+𝑦⁽⁰⁾ in such a way that the error for 𝑥⁽⁰⁾+𝑦⁽⁰⁾ be less than the error for the current iteration 𝑥⁽⁰⁾. Note that 𝐸𝑥⁽⁰⁾+𝑦⁽⁰⁾=12^𝑇₀𝑓𝑥⁽⁰⁾(𝑡)+𝑦⁽⁰⁾(𝑡)−𝑓𝑥⁽⁰⁾(𝑡)−∇𝑓𝑥⁽⁰⁾(𝑡)𝑦⁽⁰⁾(𝑡)²𝑑𝑡, (2.9) where we have used the differential equation satisfied by 𝑦⁽⁰⁾ and the definition of 𝐸(𝑥). By our assumption on 𝑓 above, 𝑓𝑥⁽⁰⁾(𝑡)+𝑦⁽⁰⁾(𝑡)−𝑓𝑥⁽⁰⁾(𝑡)−∇𝑓𝑥⁽⁰⁾(𝑡)𝑦⁽⁰⁾(𝑡)≤𝐷𝑦⁽⁰⁾(𝑡)²,𝑡∈(0,𝑇), (2.10) provided that |𝑦⁽⁰⁾(𝑡)|≤𝜖. Since we know that 𝑦⁽⁰⁾ is the solution of a certain linear problem, by the upper bound assumed in Section 1 on the size of these solutions, 𝑦⁽⁰⁾(𝑡)²≤𝑇𝐿𝐸𝑥⁽⁰⁾∀𝑡∈[0,𝑇],𝐿∈ℝ⁺. (2.11) Assume that we select 𝑇>0 so small that 𝐸_0,𝑇𝑥⁽⁰⁾≡𝐸𝑥⁽⁰⁾≤𝜖²,𝑇𝐿 (2.12) and then |𝑦⁽⁰⁾(𝑡)|≤𝜖 for all 𝑡∈[0,𝑇]. By (2.9), (2.10), and (2.11), we can write 𝐸𝑥⁽⁰⁾+𝑦⁽⁰⁾≤𝐷²2^𝑇₀𝑦⁽⁰⁾(𝑡)⁴𝐷𝑑𝑡≤²2𝐿²𝑇³𝐸𝑥⁽⁰⁾². (2.13) If, in addition, we demand, by making 𝑇 smaller if necessary, 𝐸𝑥⁽⁰⁾≤2𝛼𝐷²𝑇³𝐿², (2.14) then 𝐸(𝑥⁽⁰⁾+𝑦⁽⁰⁾)≤𝛼𝐸(𝑥⁽⁰⁾). Moreover, for all 𝑡∈(0,𝑇), 𝑥⁽⁰⁾(𝑡)+𝑦⁽⁰⁾(𝑡)≤𝐶+𝜖≤2𝐶. (2.15) All these calculations form the basis of a typical induction argument, verifying ^𝑗−1_𝑖=0𝑥^(𝑖)(𝑡)≤𝐶+𝜖^𝑗−2_𝑖=0√𝛼^𝑖),𝑥(≤2𝐶_(𝑗−1)(𝑡)≤𝜖√𝛼^𝑗−2[],𝐸∀𝑡∈0,𝑇_{𝑗−1𝑖=0}𝑥^(𝑖)≤𝛼^𝑗−1𝐸𝑥⁽⁰⁾≤𝐸𝑥⁽⁰⁾. (2.16)
It is therefore clear that the sum ^∞_𝑖=0𝑥^(𝑖)(𝑡) (2.17) converges strongly in 𝐿^∞(0,𝑇) to the solution of our initial Cauchy problem in a small interval (0,𝑇).

Since the various ingredients of the problem do not depend on 𝑇, we can proceed to have a global approximation in a big interval by successively performing this analysis in intervals of appropriate small size. For instance, we can always divide a global interval (0,𝑇) into a certain number 𝑛 of subintervals of small length ℎ (𝑇=𝑛ℎ) with 𝐸_0,𝑇𝑥⁽⁰⁾𝐷²𝐿²≤12𝛼ℎ³, (2.18) according to (2.14).

3. Numerical Procedure

Since our optimization approach is really constructive, iterative numerical procedures are easily implementable.(1)Start with an initial approximation 𝑥⁽⁰⁾(𝑡) compatible with the initial conditions (e.g. 𝑥⁽⁰⁾(𝑡)=𝑥₀+𝑡𝑓(𝑥₀)).(2)Assume we know the approximation 𝑥^(𝑗)(𝑡) in [0,𝑇].(3)Compute its derivative 𝑀(𝑥^(𝑗))^′(𝑡).(4)Compute the auxiliar function 𝑦^(𝑗)(𝑡) as the numerical solution of the problem 𝑀𝑦^′(𝑡)−∇𝑓𝑥^(𝑗)(𝑡)𝑦(𝑡)=𝑓𝑥^(𝑗)(𝑡)−𝑀𝑥^(𝑗)′(𝑡)in(0,𝑇),𝑦(0)=0, (3.1) by making use of a numerical scheme for DAEs with dense output (like collocation methods). (5)Change 𝑥^(𝑗) to 𝑥^(𝑗+1) by using the update formula 𝑥^(𝑗+1)(𝑡)=𝑥^(𝑗)(𝑡)+𝑦^(𝑗)(𝑡). (3.2)(6)Iterate (3), (4), and (5), until numerical convergence. In practice, we use the stopping criterium max𝑦^(𝑗)_∞,2𝐸𝑥^(𝑗)≤𝑇𝑂𝐿. (3.3)

In particular, one can implement, in a very easy way, this numerical procedure using a problem-solving environment like MATLAB [16].

4. Some Experiments

In this section, we approximate some problems well known in the literature for a different index [5, 17, 18]. High-order accuracy and stability are major areas of interest in this type of simulations. We do not perform an analysis of the convergence conditions imposed in the above section. We are only interested to test numerically our approach.

In our approach we only need to approximate, with at least order one, the associated linear system for 𝑦^(𝑗), in order to obtain the convergence of our scheme (see Theorem 2.2). The stability can be ensured by the fact that we approximate a linear problem using specific implicit methods [4]. This is not the case with a general nonlinear problem [19], where we need to approximate well (with a Newton-type iterative method) the nonlinear system related to the implicitness of the scheme (see the above section). This approximation should be a difficult task due to the local (nonglobal) convergence of any iterative scheme for nonlinear problems.

In this section, we consider the convergent Lobatto IIIC method [18] valid for indexes 1–3, in order to approximate the associated linear problem for 𝑦^(𝑗) in each iteration. This method can be considered as a perturbation collocation method. The final error depends only on the stopping criterium. In the following examples, we stop the algorithm when max𝑦^(𝑗)_∞,2𝐸𝑥^(𝑗)≤10⁻⁶ (4.1) and plot the solution and the approximation given by our approach.(i)Index 1 [5] 𝑦^′(𝑡)(𝑡),𝑦(𝑡)=𝑧₂+𝑧(𝑡)²𝑦(0)(0)=√=1,=𝑧22. (4.2) The solution of this problem is (sin(𝑥+𝜋/4),cos(𝑥+𝜋/4)).(ii)Index 2 [17] 𝑦^′₁(𝑡)=⁵_𝑖=1𝑓_𝑖𝑦₁(𝑡),𝑦₂(𝑡)(𝑡),𝑦,𝑧_′2(𝑡)=⁵_𝑖=1𝑔_𝑖𝑦₁(𝑡),𝑦₂(𝑡)(𝑡),𝑦,𝑧₁(𝑡)²𝑦₂(𝑡)𝑦=1,₁(0)=𝑦₂(0)𝑧(0)=1,=1, (4.3) where 𝑓₁𝑦₁(𝑡),𝑦₂(𝑡),𝑧(𝑡)=𝑦₂(𝑡)−2𝑦₁(𝑡)²𝑦₂(𝑡)+𝑦₁(𝑡)𝑦₂(𝑡)²𝑧(𝑡)²+2𝑦₁(𝑡)𝑦₂(𝑡)²−2𝑒^−2𝑡𝑦₁(𝑡)𝑦₂(𝑡),𝑓₂𝑦₂(𝑡),𝑧(𝑡)=−𝑦₂(𝑡)²𝑧(𝑡)+2𝑦₂(𝑡)²𝑧(𝑡)²,𝑔₁𝑦₁(𝑡),𝑦₂(𝑡)=−𝑦₁(𝑡)²+𝑦₁(𝑡)²𝑦₂(𝑡)²,𝑔₂𝑦₁(𝑡),𝑦₂(𝑡),𝑧(𝑡)=−𝑦₁(𝑡)+𝑒^−𝑡𝑧(𝑡)−3𝑦₂(𝑡)²𝑧(𝑡)+𝑧(𝑡). (4.4) The solution of this problem is (𝑒^𝑡,𝑒^−2𝑡,𝑒^2𝑡).(iii) Index 3 [18] 𝑦^′₁(𝑡)=2𝑦₁(𝑡)𝑦₂(𝑡)𝑧₁(𝑡)𝑧₂(𝑡),𝑦^′₂(𝑡)=−𝑦₁(𝑡)𝑦₂(𝑡)𝑧₂(𝑡)²,𝑧^′₁(𝑡)=𝑦₁(𝑡)𝑦₂(𝑡)+𝑧₁(𝑡)𝑧₂(𝑡𝑢(𝑡),𝑧_′2)(𝑡)=−𝑦₁(𝑡)𝑦₂(𝑡)²𝑧₂(𝑡)²𝑢(𝑡),𝑦₁(𝑡)𝑦₂(𝑡)²𝑦=1,₁(0)=𝑦₂(0)𝑧=1,₁(0)=𝑧₂(0)𝑢(0)=1,=1. (4.5) The solution of this problem is (𝑒^2𝑡,𝑒^−𝑡,𝑒^2𝑡,𝑒^−𝑡,𝑒^𝑡).

In Figures 1, 2, 3, 4, 5, 6, 7, 8, 9, and 10, we compare the solution of the corresponding three problems with the approximations given by our approach. The results are very satisfactory in all cases, obtaining always the convergence to the true solution. In a first look, the exact and computed solutions are indistinguishable, since after convergence the error is smaller than the tolerance (=10⁻⁶) used in the stopping criterium. A more systematic and careful analysis of the numerical possibilities of the method will be pursued in the future.

Figure 1 

Index 1, 𝑇=2𝜋, the 𝑦-coordinate, “o”-original, “+”-approximation.

Figure 2 

Index 1, 𝑇=2𝜋, the 𝑧-coordinate, “o”-original, “+”-approximation.

Figure 3 

Index 2, 𝑇=2, the 𝑦1-coordinate, “o”-original, “+”-approximation.

Figure 4 

Index 2, 𝑇=2, the 𝑦2-coordinate, “o”-original, “+”-approximation.

Figure 5 

Index 2, 𝑇=2, the 𝑧-coordinate, “o”-original, “+”-approximation.

Figure 6 

Index 3, 𝑇=2, the 𝑦1-coordinate, “o”-original, “+”-approximation.

Figure 7 

Index 3, 𝑇=2, the 𝑦2-coordinate, “o”-original, “+”-approximation.

Figure 8 

Index 3, 𝑇=2, the 𝑧1-coordinate, “o”-original, “+”-approximation.

Figure 9 

Index 3, 𝑇=2, the 𝑧2-coordinate, “o”-original, “+”-approximation.

Figure 10 

Index 3, 𝑇=2, the 𝑢-coordinate, “o”-original, “+”-approximation.

5. Conclusions

A new variational approach to the analysis and numerical implementation of regular ODEs has been recently introduced in [10, 20]. Because of its flexibility and simplicity, it can easily be extended to treat other types of ODEs like differential-algebraic equations (DAEs). This has been precisely the main motivation for this paper: to explore how well those ideas can be adapted to this framework. In particular, extending to this context some of the analytical results and performing various numerical tests that confirm that indeed the variational perspective is worth pursuing. One remarkable feature is that this point of view only requires to count on good numerical schemes for linear problems, and this is the reason why it fits so well in other scenarios. Because of the many good qualities of this viewpoint, it can be considered and implemented in essentially all fields where differential equations are relevant. There is, then, a long way to go.

Acknowledgments

S. Amat and M. J. Légaz supported by MINECO-FEDER MTM2010-17508 (Spain) and by 08662/PI/08 (Murcia). P. Pedregal supported by MINECO-FEDER MTM2010-19739 (Spain).