Chaire d'Analyse Appliquée                    The DNA modelling COURSE

Session 5

Introduction

To understand buckling and imperfections we will retreat to an extremely simple, finite dimensional model, which can be regarded as a finite difference approximation to the planar rod problem.
 

The Discrete Strut-problem formulation

Consider the energy $$\sum_{i=0}^n k_{i+1} (\phi_{i+1} - \phi_{i} -\beta_{i+1})^2/2 + \lambda h\, \cos\phi_{i}$$ where $\lambda$, $h$, the $\beta_{i+1}$, and the $k_{i+1}$, $i=0,\dots,n$ are given, there are $n+2$ unknown variables $\phi_{i}$, $i=0,\dots,n+1$, and we have the ``boundary conditions'' $$\phi_0=\beta_0,\quadk_{n+1}(\phi_{n+1} - \phi_{n} - \beta_{n+1})=0.$$ where $\beta_0$ is an additional parameter.

The quadratic terms are a model for the energy stored in $n$ rotational springs, with spring constant $k_{i}$ at the $i$th joint, $i=1,\cdots,n$and a minimum energy at the $i$th joint when the angles satisfy

$$\phi_{i}=\phi_{i-1} + \beta_{i}$$

The fact that $\phi_0=\beta_0$ is prescribed is analogous to the boundary condition $\phi(0)=\beta_0$ in the continuous problem, while the ``boundary condition'' involving $\phi_{n+1}$ is a discretized version of the vanishing moment condition $K_2(1)(\phi'(1)-\hatu_2(1))=0$, which says that there is no moment applied at the top of the segment with $n$th joint.

We could equally well eliminate $\phi_0$and $\phi_{n+1}$ using the boundary conditions to leave a sum with index range $i=1,\dots,n$ with $n$ unknowns (and redefining $\bar\beta_1=\beta_0+\beta_1$) and there would be only $n$ quadratic terms to be summed.

The term in the energy involving $\lambda$ and $h$ models the potential energy of a vertical load (with the convention $\lambda$ +ve acting downward) applied to the $(n+1)$st joint (or equivalently the end of the $n$th segment with each segment being of length $h$). Think of a bucket of water of weight $\lambda$ hung from the $(n+1)$st joint, then the potential energy is just the height $\sumh\,\cos\phi_i$ times the weight $\lambda$.

 In the formulation where $\phi_{n+1}$ is kept as a variable it seems a little strange that the load is not applied at the tip of the $(n+1)$st segment, but it is quite usual that in a finite difference approximation, it is necessary to introduce an additional mesh point `beyond the boundary' in order to approximate derivative boundary conditions at the boundary. In any case, in this example eliminating $\phi_{n+1}$ eliminates any confusion.
 

The Minimum Energy State

We also see that in the case of a zero external load $\lambda=0$ the minimum energy configuration is given by $$\hat\phi_i\equiv\sum_{j=0}^i\beta_k,\qquad \forall i=0,\dots,n+1$$ where we have also included indices for the boundary variables in the summation range.

Thus we can see that the $\beta_i/h$ can be identified as determining a piece-wise constant approximation of the unstressed strain $\hat u_2(s)$ in the continuous problem, or equivalently the $\beta_i$ determine a piece-wise linear approximation with corner values $\hat\phi_i$, of the continuous unstressed shape $\hat\phi(s)$determined by the property $\hat\phi'(s)=\hat u_2(s)$.
 

The Equilibrium conditions

Taking partial derivatives of the energy with respect to the variables $\phi_{i}$, $i=1,\dots,n$ yields the stationarity, or first-order necessary conditions

which is a set of $n$ nonlinear equations for the $n$ unknowns $\phi_{i}$, $i=1,\dots,n$, but the equations also contain many additional parameters, namely the single scalars $\lambda$ and $h$ and the sets of parameters $k_{i}$, $i=1,\dots,n$ and $\beta_i$, $i=0,\dots,n$.

Notice that the boundary condition relating $\phi_{n+1}$ to $\phi_n$ is exactly the additional first-order condition that is obtained by taking the partial derivative with respect to the variable $\phi_{n+1}$. This is the finite dimensional manifestation of what is called natural boundary conditions in the calculus of variations (which we will see in a little while).

The boundary condition is actually implied by the variational principle. In contrast it makes no sense to take a partial derivative with respect to the variable $\phi_0$ because the boundary condition gives $\phi_0$ a definite value in terms of the given parameter $\beta_0$. Such boundary conditions are sometimes called imposed, and are not implied by the variational principle.

These equations could have been derived directly as the balance (about the $i$th joint) of moments acting on the $i$th segment. Thus we see that there is a variational formulation or minimization principle giving the moment balance equations. Notice that following common practice we talk about a minimization principle, but the balance law comes as the first-order stationary condition of all partial derivatives vanishing (i.e. the potential is only known to be flat, but it could be a local minimum, a local maximum, or some sort of saddle point).

It is certainly of considerable interest to know which of all the equilibria (solutions of the balance laws) are actually local minima (which we expect to be stable equilibria). The type of the critical point can be classified by studying the second derivative of the potential evaluated at a particular critical point.

The Second Variation of the Strut Problem

In the case of the discrete strut potential given above, the second derivative can be regarded as the $n\times n$ symmetric matrix whose $(i,j)$th entry is the second partial derivative of the potential with respect to the variables $\phi_{i}$ and $\phi_j$ (for $i,j=1,\dots,n$).

For the particular potential above, the matrix is tri-diagonal with a typical $(3\times 3)$ block centred on the $(i,i)$entry being

$$\left( \begin{array}{ccc}(k_{i-1} + k_{i}) -\bar\lambda h\... ...+ k_{i+2}-\bar\lambda h\, \cos\bar\phi_{i+1})\end{array}\right)$$

and appropriately modified first and last blocks (after $\phi_0$ and $\phi_{n+1}$ are eliminated using the boundary conditions).

Here $\bar\phi_i$ is a particular solution of the equilibrium conditions with $\lambda=\bar\lambda$ at which the second derivative (or second variation) is being computed. The dependence of the second derivative on the the solution $(\bar\phi_1,\dots,\bar\phi_n,\bar\lambda)$ is particularly simple because the bending part of the energy is exactly quadratic.

Similarly because the terms involving the $\beta_i$ in an expansion of the energy are all linear in the $\phi_{i}$ the parameters $\beta_i$ do not appear explicitly in the second derivative (although the $\beta_i$ certainly appear explicitly in the first-order conditions, and so implicitly in the second derivative through the exact form of the equilibrium $\bar\phi_i$). Because the second derivative matrix is symmetric, it has only real eigenvalues, and we will call the solution, $(\bar\phi_1,\dots,\bar\phi_n,\bar\lambda)$non-degenerate if its second derivative matrix is nonsingular, i.e. is not an eigenvalue of the matrix.

For non-degenerate equilibria we define its index to be the number of negative eigenvalues of its second derivative matrix. In particular a Taylor series expansion of the (finite dimensional) potential shows that a non-degenerate equilibrium is a local minimum of the potential (and thus presumably stable) precisely if its index is zero.

We remark that classifying degenerate equilibria is more subtle, because the second derivative cannot be used to bound all the higher order terms. (Similarly in infinite dimensional problems the situation is not so straightforward even at non-degenerate equilibria.)

The exercise sessions will ask you to compute the second derivative at some particular solutions in the case $n=2$.

In a little while we will turn to discussing how the above discrete system, can be interpreted as a discretization of the continuous system which will be recovered in a limit $n\to\infty$, $h\to 0$. But first we will discuss methods for solving nonlinear systems of parameter-dependent algebraic equations, and understanding typical behaviours.
 

Solution of nonlinear parameter-dependent equations

We would like to understand the solution set of the discretized strut as the parameters of the problem are varied. In particular we shall consider $\lambda$ as the main parameter to be varied, with the stiffnesses $k_{i}$ and the discretized unstressed shape $\beta_i$ being for the most part regarded as given fixed parameters. In other words we want to know what are the possible equilibrium shapes as the end load is increased say from $\lambda=0$.

Generally speaking it is hopeless to find anything other than very special solutions of nonlinear equations analytically, and we must instead turn to numerical methods.

Because we know that our problem has an associated variational formulation, we could implement a numerical method exploiting that structure, say steepest descent, i.e. take an iterative method that at each step marches down the gradient of the potential. Such methods have some very desirable features. For example they can reasonably be expected to converge to a local minimum of the potential from an essentially arbitrary starting point--which is an extremely powerful statement.

On the other hand minimization by steepest descent or other methods has its draw backs, for example often there is slow convergence in the neighbourhood of a local minimum. For our purposes the bigest disadvantage of minimization methods is that they can only be expected to find local minima, and as we shall see it can be rather instructive to find other types of (presumably unstable) equilibria.

The only practical alternative to a minimization method is a root-finding technique applied directly to the system of equations that arises as the first-order necessary conditions. If for the moment we freeze all parameters in the first-order conditions, the first-order conditions assume the general form

$${\bf f}({\bf z})=0$$ of $n$ nonlinear equations in $n$ unknowns ${\bf z}$.

At this level of generality the only practical way to solve such systems numerically is some sort of iterative method starting from a sufficiently good guess for the solution. Amongst such iterative schemes, those that are some variant of Newton's method are extremely powerful.
 

Review of Newton's method

Let us therefore review Newton's method. Suppose that ${\bf z}^k$ is some good approximation to a root ${\bf z}^e$ of the governing equations. Then a Taylor series expansion gives us $${\bf f}({\bf z}^e)={\bf f}({\bf z}^k) +{\cal J}^k({\bf z}^e-{\bf z}^k) + O(\vert{\bf z}^k-{\bf z}^e\vert^2)$$ where ${\cal J}^k$ denotes the $(n\times n)$ Jacobian matrix with $(i,j)$th entry $\partial f_i/\partial z_j$ evaluated at ${\bf z}^k$.

Neglecting the higher order terms, and using the fact that by assumption ${\bf f}({\bf z}^e)=0$, we see that the next iterate ${\bf z}^{(k+1)}$ defined via

$${\bf z}^{(k+1)}= {\bf z}^k + {\bf w}^k$$ with  defined to be the solution of the inhomogeneous linear system $${\cal J}^k{\bf w}^k=-{\bf f}({\bf z}^k)$$ should be a very good estimate for . In fact Newton's method converges quadratically once the approximation is sufficiently close (which is an extremely good feature).

Newton's method had two major potential difficulties. First it may be either computationally expensive or difficult to solve the linear system for the update ${\bf w}^k$ either because the matrix ${\cal J}^k$ is very large and dense, or is ill-conditioned or actually singular. Second, to start the iteration you need a first guess ${\bf z}^0$, and often that guess has to be really quite good in order that the Newton iteration will converge.

Our particular problem has special features that can help with both of these potential difficulties. First because we have a variational principle we know that the first order conditions are of the form , but with the vector field ${\bf f}$ equal to the gradient of the scalar potential, $V$ say.

Thus $f_i({\bf z})=\partial V/\partial{\bf z}_i$. In particular the Jacobian matrix is just the symmetric matrix of second derivatives of the potential $V$, which for the planar rod problem is tri-diagonal (some such banded sparsity structure of the Jacobian is typical of problems that arise as discretizations of systems of ordinary differential equations). Linear systems involving banded matrices can be solved extremely efficiently (for example a banded Gaussian elimination), whenever they are not singular.

The second special feature that we shall exploit is the parameter dependence of the problem with an explicitly known solution at one parameter value (namely the unstressed shape when $\lambda=0$). We shall apply parameter continuation methods in our numerical computations of rod equilibria.
 

Parameter Continuation method

That is to say we shall use knowledge of previously computed solutions at nearby parameter values, to generate an accurate guess for the solution at the current parameter values, which is then corrected via a Newton-type iterative solve. The whole process can get started provided that a solution (or even just a good guess for the solution) is known at some specific set of parameter values.

Let us make this outline more concrete. Suppose therefore that we have a system of equations of the form

$${\bf f}({\bf z};\lambda)=0$$ and a known solution $({\bf z}^0,\lambda^0)$. Then provided that the Jacobian ${\cal J}({\bf z}^0,\lambda^0)$ is nonsingular the implicit function theorem guarantees that there is a neighbourhood $\vert\lambda-\lambda_0\vert<\epsilon$, such that there is a unique branch of solutions passing through $({\bf z}^0,\lambda^0)$ which can be parametrized as $({\bf z}(\lambda),\lambda)$.

How might we compute such a branch of solutions? If we seek solutions for $\lambda=\lambda^0 +\lambda^1$ with $\lambda^1$ small, we can expand the solution as

$${\bf z}(\lambda)={\bf z}^0 + \lambda^1\,{\bf z}^1 + O(\vert\lambda^1\vert^2)$$ and expand the equation $${\bf f}(\{{\bf z}^0 + \lambda^1{\bf z}^1 + O(\vert\lambda^1\vert^2)\};\{\lambda^0 +\lambda^1\}) = 0$$ as $${\bf f}({\bf z}^0;\lambda^0) + \lambda^1\,{\bf f}_{{\bf z}}(... ...{\bf f}_{\lambda}({\bf z}^0;\lambda^0)=O(\vert\lambda^1\vert^2)$$

But  as  is known to be a solution, and (as a matter of notation)  is just the Jacobian . Thus we see that by ignoring  terms we obtain a nonhomogeneous linear system that can be solved for a unique solution  provided that the Jacobian  is nonsingular:

$${\cal J}^0\,{\bf z}^1 = -{\bf f}_{\lambda}^0$$ (where the superscript on ${\cal J}$ and $-{\bf f}_{\lambda}$ indicates that the functions are to be evaluated at $({\bf z}^0,\lambda^0)$).

Accordingly we can reasonably expect that for $\lambda^1$ small a Newton iterative solve of the equations at the parameter value $\lambda=\lambda^0 +\lambda^1$ will converge if we take as initial guess

$${\bf z}^g = {\bf z}^0 + \lambda^1\,{\bf z}^1$$ where ${\bf z}^1$ is determined by a linear solve of the system with coefficients determined by the Jacobian ${\cal J}^0$ and the function ${\bf f}_{\lambda}^0$ evaluated at the previously computed solution $({\bf z}^0,\lambda^0)$.

Once the Newton iterate has converged, and if the Jacobian is nonsingular at the new solution, we may repeat the process to generate a new initial guess at another nearby value of the parameter.

Thus we can generate a number of solutions at a discrete set of values of the parameter $\lambda$ that are close to each other, and which approximate the branch of solutions $({\bf z}(\lambda),\lambda)$through the initially known solution $({\bf z}^0,\lambda^0)$.

This is an elementary (and rather naive) description of numerical continuation. As you can probably imagine there is a wealth of important detail in getting the method to work efficiently, but the above captures at least the flavour of the method that is implemented in the continuation code AUTO that we will use later on (via the graphics interface VBM).

We can see that the above method is problematic whenever we hit a solution with a singular or near singular Jacobian. The analysis of what happens at such singular points is the domain of bifurcation theory. The typical behaviour at simple singular points is that the branch of solutions persists, but there may be multiple branches of solutions, or there may be a fold point at which the parameter $\lambda$has a local extreme value along the branch of solutions.

We will illustrate these concepts by considering the case $n=1$ of the discrete strut example. The $n=1$ case can be solved in an essentially explicit fashion, but it is a little too simple to truly explain the perturbation expansion described above. You will consider the case $n=2$ in the exercise session.

This n=1 strut problem will soon be available on the web too.

last update 23.11.1998
Back to the Table of Contents