**Suggested Citation:**"8 Nonlinear Equations and Bifurcation." National Research Council. 1991.

*Research Directions in Computational Mechanics*. Washington, DC: The National Academies Press. doi: 10.17226/1909.

**8**

**NONLINEAR EQUATIONS AND BIFURCATION**

Finite-dimensional nonlinear systems of equations have been encountered in mechanics ever since the field's formalization in Newton's time. It is well known that, except for very simple cases, there are no general methods for obtaining exact solutions of such systems. This has been one of the reasons for the extensive use of linearized mechanical models throughout history. However, in recent years and with the development of electronic computers, there has been a growing trend toward extending the computational methodology to nonlinear problems.

Nonlinearities pervade nearly all parts of mechanics, and accordingly nonlinear systems of equations are found in connection with a multitude of problems. Some of the principal examples include systems arising (1) as discretizations of nonlinear equilibrium equations, (2) in connection with nonlinear optimization problems, (3) in implicit methods for the solution of dynamic processes, and (4) in the formulation of the nonlinear dynamics of constrained multibody systems.

The form and properties of these nonlinear systems are strongly affected by the problem area as well as by the different sources of nonlinearities encountered in most types of problems. The latter point can be illustrated by some typical examples in structural mechanics, including (1) geometrical nonlinearities due to nonlinear strain-displacement relations, (2) material nonlinearities due to nonlinear constitutive equations, (3) force nonlinearities due to nonlinear stress boundary conditions, and (4) kinematic constraint nonlinearities due to nonlinear displacement boundary conditions.

The choice and effectiveness of a solution technique for any nonlinear equation depend critically on factors such as (1) the form and properties of the equations, (2) the specific aims and accuracy requirements of the computation, (3) the amount of available information about the problem, and (4) the computing resources to be used.

The theoretical understanding of these and related factors has not reached, by far, a satisfactory stage. Many techniques and computer programs have been developed for

**Suggested Citation:**"8 Nonlinear Equations and Bifurcation." National Research Council. 1991.

*Research Directions in Computational Mechanics*. Washington, DC: The National Academies Press. doi: 10.17226/1909.

various types of problems, but their mathematical basis is often insufficiently understood, and there is rarely much known about the accuracy of the computed results.

Deficiencies in the numerical analysis of many problems are often directly attributable to an insufficient theoretical basis for the problems themselves. This is the case, for instance, in finite deformation plasticity where there are profound mathematical and computational difficulties in modeling phase changes, viscous effects, cracks, singularities, etc. At the other end of the spectrum, there are many novel questions relating to scientific computing, such as those concerning the effect of different types of computational resources and, especially, the development of new computer architectures on the selection, implementation, and performance of different numerical methods.

In line with the above remarks, it must be noted that the following discussion will touch only a few of the numerous open research topics in this area.

**ITERATIVE TECHNIQUES**

Mathematical theory shows that a system of n nonlinear equations in n unknowns

*f*_{i}(*x*_{1}, .... , *x*_{n}) = 0, *i* = 1, ... , *n* (8.1)

may have many solutions. In fact, there exists a theorem to the effect that for certain systems any closed subset of **R**^{n} can occur as the solution set. This also implies that any numerical algorithm must be based on specific *a priori* knowledge about the existence and type of solution that is to be computed. It also means that there cannot be any best algorithm for general-purpose applications.

In practice, a basic task relating to a system (8.1) is the development of locally convergent iterative processes for the computation of an isolated solution *x**. These are processes for which there exists an open neighborhood *U* of *x** such that for any starting point(s) in *U* the computed iterates do not leave *U* and converge to *x**.

The literature on such locally convergent methods is certainly huge, but there are still numerous open problems, especially when it comes to practical applications. Such questions as storage management, data access, vectorization, and parallel processing have a profound influence on the overall procedure, especially when the dimension *n* is large.

**Suggested Citation:**"8 Nonlinear Equations and Bifurcation." National Research Council. 1991.

*Research Directions in Computational Mechanics*. Washington, DC: The National Academies Press. doi: 10.17226/1909.

For large *n* it becomes necessary to consider the use of secondary iterative techniques for solution of the linear systems at each step. Then effective control of the combined primary and secondary iterations assumes central importance. Some challenging questions of this type arise when, as a secondary process, a multigrid technique is applied that utilizes the fact that the nonlinear system is the discretization of some differential equation.

In connection with locally convergent iterative processes, it must be noted that the attraction region for a specific solution *x** is anything but a simply shaped set. In fact, while it is usually known that this region contains some ball around *x**, the space between the attraction balls belonging to different solutions of (8.1) generally has a fractal nature. Since the size of the attraction balls is rarely accessible, this makes early detection of suspected divergence of any iterative process for (8.1) a critical issue for the overall computational cost.

Some general problems concerning any such iterative process for (8.1) are the design of effective and reliable tests for termination of the process either when, as mentioned, divergence is suspected or when convergence within given error tolerances can be declared.

Various approaches have been proposed for control of round-off errors. In particular, interval arithmetic has been shown to be capable of providing excellent guarantees for computed solutions of nonlinear equations. But at present these methods can be expected, at best, to become applicable for the large problems arising in computational mechanics only when the needed operations are available in hardware.

**PARAMETERIZED PROBLEMS**

In general, nonlinear problems in computational mechanics involve a number of parameters that influence the behavior of the physical system and that have to be taken into account in the computation. For example, in structural mechanics these parameters may characterize intrinsic quantities such as material properties, geometric dimensions, or extrinsic influences such as load points, directions, and intensities. Thus, the system of nonlinear equations (8.1) has. the more general form

*f*_{i}(*x*_{1}, ... , *x*_{n}, λ_{1}, ... , λ_{m}) = 0, *i* = 1, ... , *n* (8.2)

involving a vector of parameters λ = (λ_{1}, ... , λ_{m}) ∈ **R**^{m}.

**Suggested Citation:**"8 Nonlinear Equations and Bifurcation." National Research Council. 1991.

*Research Directions in Computational Mechanics*. Washington, DC: The National Academies Press. doi: 10.17226/1909.

In general, it is rarely desired to compute only some specific state vector *x* = (*x*_{1}, ... , *x*_{n}) which, together with a specific choice of the parameter vector λ, solves (8.2). Instead, interest centers on assessing the changes of the state vectors when the parameters are varied and, in particular, determining parameter values where a change in the physical solution behavior is expected, such as a loss of stability for certain equilibrium problems.

Under natural conditions, the set of all solutions

M = {(*x*, λ), *f*_{i} (*x*, λ) = 0, *i* = 1, ... , *n*} (8.3)

forms a submanifold of dimension *m* in the product space **R**^{n} × **R**^{m} of the state and parameter variables. Then the question is to determine the shape and properties of this solution manifold.

The classical technique for analyzing the solution manifold (8.2) of (8.3) is the continuation method. It consists in the *a priori* selection of some path on the manifold M and in computing a sequence of points along that path starting from some known point. In practice, selection of the continuation path usually consists of specification of the parameters as some function λi = λi(μ), *i* = 1, ... *m* of a single variable (μ). The solution step from a current point then consists of determining an appropriate predicted point along the path and application of a local iterative procedure to a suitable augmented nonlinear system *n* + 1 equations in as many unknowns. Various forms of continuation processes have been proposed that differ in choice of the predicted point, form of the augmented equations, and type of local iterative process.

Thus, all the errors for local iterative processes arise also in connection with continuation methods. In addition, it becomes important how to determine accurately the required ''singular points'' on the path where the solution behavior is expected to change. In this connection it is also critical to decide when the process leaves a particular component of the solution path. These questions of continuation method accuracy are as yet largely unresolved.

For a numerical analysis of a given solution manifold, it appears to be advantageous to consider continuation methods in a broader sense as a collection of numerical procedures for a variety of tasks, such as (1) following numerically any curve on the manifold and providing estimates of the accuracy of the computed points; (2) on any such curve determining the location of target points where a given variable has a specified

**Suggested Citation:**"8 Nonlinear Equations and Bifurcation." National Research Council. 1991.

*Research Directions in Computational Mechanics*. Washington, DC: The National Academies Press. doi: 10.17226/1909.

value and, if desired, at such a target point switching to the trace of a different *a priori* specified path; and (3) on any solution path determining the location and properties of the singular points where the solution behavior changes and, depending on the type of the detected point, taking appropriate action.

For specific applications, additional tasks may arise. For instance, if (8.3) represents the equilibrium equations of a system of differential equations, there may be interest in locating Hopf bifurcation points on a solution path where periodic solutions of the dynamic system branch off from the static equilibrium.

In recent years much research has been devoted to these various tasks, but numerous questions remain. Library-quality software is available only for some tasks. For specific applications, such as fluid dynamics problems or structural problems involving plasticity, creep, and viscoelastic effects, the situation is even more open.

The continuation processes, and more generally the above tasks, mean that a picture of an *m*-dimensional manifold in *n* + *m* space is to be developed solely from computed information along *a priori* specified paths. But, in many problems in mechanics, it is essential to assess the sensitivity of a solution (*x*, λ) of (8.2) under potential changes of several parameters. In other words, the interest is to determine the curvature behavior of the manifold M at the particular point. For this purpose, in recent years, some methods have been proposed that allow for computation of triangulations of the manifold in the neighborhood of a given point and for subsequent determination of some measure of the curvature. These multiparameter methods were in their infancy, but are no longer so. In particular, they have not yet been implemented in production codes nor applied to large practical problems.

**BIFURCATION WITH OR WITHOUT SYMMETRY**

As noted earlier, for a deeper understanding of the behavior of a mechanical system, it is necessary to analyze the shape of the full solution manifold and, in particular, to determine the singular points where the solution behavior is expected to change. Here singularity may be characterized by the condition that the system's Jacobian matrix is not subjective with respect to the stated variable. Of course, in this context the choice of the parameters entering into the definitions of the equations is critically important.

**Suggested Citation:**"8 Nonlinear Equations and Bifurcation." National Research Council. 1991.

*Research Directions in Computational Mechanics*. Washington, DC: The National Academies Press. doi: 10.17226/1909.

Typical computational tasks arising in connection with bifurcation phenomena include accurate determination of the singular points along paths. Usually, this is accomplished via a characterization of the expected singular point as a solution of an augmented system. In their present form the available iterative methods for such augmented equations are local in nature and often require for their convergence the availability of exceptionally close initial guesses. Moreover, currently used augmentations generally assume knowledge of the dimension of the null space. It would be of practical importance to design improved algorithms where both these constraints are relaxed as much as possible.

**Calculation of the Bifurcated Branches**

Singular points along a path are the only ones where several solution branches may intersect and, hence, where multiple options arise for the path continuation. Different branches typically correspond to different physical behavior of the underlying mechanical system, although some branches may have little or no physical reality. However, it is generally necessary to calculate parts of a solution branch before it becomes possible to decide its physical relevance. While the cases of simple turning points and simple bifurcation points have been treated exhaustively in the theoretical and numerical literature, there is as yet little material available regarding calculation of all the branches emanating from multiple bifurcation points.

**Calculation of Paths of Singular Points**

Singularity theory has provided some rigorous basis for the intuitive idea that a hierarchy exists among singularities. This in turn has given some substance to the statement that a singular point is more or less singular than another one. Singularities of high rank in the hierarchy, especially organizing centers, are of practical interest because they admit a large number of qualitatively nonequivalent perturbations, and hence their analysis may reveal unsuspected phenomena. While both the dimension *n* of the state variable and the number *m* of the parameters are important factors for a given type of singularity to be generically present in an equation system, a singularity of high rank can hardly be detected by following a few randomly chosen paths. In fact, high-rank singularities are rather scarce and, hence, lie only on selected

**Suggested Citation:**"8 Nonlinear Equations and Bifurcation." National Research Council. 1991.

*Research Directions in Computational Mechanics*. Washington, DC: The National Academies Press. doi: 10.17226/1909.

and *a priori* unknown solution paths. In other words, detection of such singularities requires special strategies, such as following paths of singular points along which higher singularities may be found. In practice, paths of singular points may be obtained by solving the original system augmented with additional relations that ensure singularity. This technique has been implemented in some of the simpler cases, but a satisfactory general procedure is not available at this time.

It should be pointed out that one ingredient entering the solution of all three problems mentioned above is the availability of cost-effective and accurate techniques for the numerical evaluation of higher-order mapping derivatives. Some recent methods for calculating the curvature of the solution manifold have also provided interesting possibilities for the numerical calculation of second-order derivatives, but further contributions to this area are certainly needed.

Bifurcations are often closely connected with a loss or "breaking" of symmetry. Several problems, both theoretical and numerical, are worthy of examination. Certainly, it is important to be able to predict as accurately as possible what symmetries will be preserved along some bifurcated branches when a bifurcation point is encountered along a path of "full" symmetry. For this, bifurcation theory provides a certain number of criteria but by no means a comprehensive answer. On the positive side, it may be noted that the case of the dihedral or orthogonal groups and their subgroups governing symmetry in many problems of mechanics is now much better understood than the general case.

Among the many questions related to numerical aspects is, for instance, the preservation of symmetries of the original "continuous" (e.g. partial differential equation) version of the problem. Clearly, symmetry can be utilized in the computation only if the discretized equations have inherited at least part of the original problem's symmetry. This often requires specific constraints in the mesh design, etc., that can be met, in practice, only with some difficulties. Another important question concerns numerical identification of the invariant subspaces and, in finite element methods, of suitable bases of such subspaces. Here it may be noted that the required symmetry properties of the basis functions may significantly enlarge the size of their supports, which in turn may have a negative effect on the sparsity of the matrices involved in the calculations. It has not yet been examined how much this inconvenience balances the gain in the reduction of the dimension of the state space or whether the symmetries induce properties in

**Suggested Citation:**"8 Nonlinear Equations and Bifurcation." National Research Council. 1991.

*Research Directions in Computational Mechanics*. Washington, DC: The National Academies Press. doi: 10.17226/1909.

the matrices that may be exploited to make up for the deterioration of their sparsity.

Finally, it has to be emphasized that all computations can work only with singular points on the solution manifold of some finite dimensional nonlinear equations (8.2). Since most of these equations in computational mechanics represent discretizations of certain differential equations, the question arises whether any computed critical phenomena, such as a particular bifurcation point and its branching behavior, correspond to similar phenomena of the original problem and, if so, what errors are encountered. These questions have as yet only partial answers and represent considerable research challenges. In this connection it must be noted that for nonlinear problems the solution manifold of the discretized equations often shows different features from that of the original equations. In particular, there may be connected components of the approximate manifold that do not occur in the exact solution manifold. Such components have been called spurious and are being observed more and more often in various applications, but their study has only recently begun.

This entire range of questions about the discretization errors of multiparameter nonlinear problems is only in a very embryonic stage. The same is true for the design of efficient and reliable *a posteriori* error estimates, although results for certain specific classes of problems and discretizations have shown that estimators for the discretization errors at computed points on the solution manifold are computationally feasible. There is certainly a need for concentrated research in this area.

**DIFFERENTIAL ALGEBRAIC EQUATIONS**

Differential algebraic equations (DAEs) may be characterized as problems of the form

*F*(*t*,*x*,*x*') = 0, (8.4)

where F : **R** × **R**^{n} × **R**^{n} ⇒ **R**^{n} and the derivative *D*_{p}*F* with respect to the third variable has constant but not full rank *r* < *n*. Thus, in contrast to the case *r* = *n*, equation (8.4) cannot be transformed into an explicit ordinary differential equation (ODE) in **R**^{n} in the vicinity of any given (*t*_{0}, *x*_{0}, *p*_{0}) with *F*(*t*_{0}, *x*_{0}, *p*_{0}) = 0.

An important class of mechanical examples involving DAEs arises in the study of the dynamics of constrained multibody systems. Inherent here is the need to account for

**Suggested Citation:**"8 Nonlinear Equations and Bifurcation." National Research Council. 1991.

*Research Directions in Computational Mechanics*. Washington, DC: The National Academies Press. doi: 10.17226/1909.

algebraic kinematic constraints due to mechanical couplings that restrict relative and absolute motion of the components in a multibody system. Thus, if *x* Є **R**_{n} represents a vector of generalized coordinates, the resulting DAE consists of the algebraic constraint equations

Φ(*x*,*t*) = 0 (8.5)

together with the constrained equations of motion

M(*x*,*t*)*x"* + Φx(*x*,*t*)T*z* = *Q*(*x*,*x'*,*t*) (8.6)

where *M* denotes the mass matrix, *z* is the vector of Lagrange multipliers, and *Q* is the vector of generalized applied forces.

The formulation (8.5) – (8.6) has been known for centuries, but it was only recognized during the past decade that DAEs have different properties than ODEs and that their numerical solution involves particular difficulties. While there has been considerable progress in the development of specific computational methods for DAEs in general, at present there exist no broadly applicable and robust integrators suitable for automated application of multibody dynamic simulation methods in a computer-aided engineering environment. Recent developments of methods of generalized coordinate partitioning methods for local parametrization of the kinematic constraints appear to offer particular promise for the effective solution of (8.5) – (8.6), but much still remains to be done.

Other recent advances have shown that a satisfactory existence and uniqueness theory for initial value problems associated with DAEs can be worked out, which extends the standard ODE theory. By showing how DAEs locally reduce to ODEs on implicitly defined manifolds, these results suggest new techniques of solution and identify new challenges. For instance, it appears that the design of higher-order schemes for the treatment of ODEs on manifolds may have to include curvature considerations.

**SINGULAR DIFFERENTIAL EQUATIONS**

Singular differential equations differ from DAEs in that noninvertibility of the derivative *D*_{p}*F* in (8.4) occurs only at the points of a set with an empty interior. In other words, in a singular ODE, singularity may be termed an "accident." Such equations occur in a variety of physical problems, including, for instance, plasticity theory. They have been studied much less than DAEs despite their first appearance in the scalar case

**Suggested Citation:**"8 Nonlinear Equations and Bifurcation." National Research Council. 1991.

*Research Directions in Computational Mechanics*. Washington, DC: The National Academies Press. doi: 10.17226/1909.

*n* = 1 well over a hundred years ago. In fact, a rigorous analysis of singular ODEs in the general case has begun only recently but has already shown that singularities may be responsible for some quite peculiar dynamic behavior. Much theoretical and numerical work remains to be done in this little explored area.