Nonlinear finite elements/Newton method for axial bar
Newton's method applied to the axial bar problem 
Let us now return to the axial bar problem and see how the Newton-Raphson scheme can be used to solve the system of equations. The reduced system of equations is
Therefore, the residual is
The terms of the tangent stiffness matrix are given by
For our reduced stiffness matrix, we have
Now, from the stiffness matrix, we see that
The derivatives are,
Therefore, the components of are
In matrix form,
The inverse of (the reduced form shown in equation (8)) is
==== Maple code ====11 A Maple text script for doing these computations symbolically is show below.
Can the tangent stiffness matrix be assembled too? 
We have derived the tangent stiffness matrix from the assembled global stiffness matrix. If the stiffness matrix is large that does not seem to be the way to go. So the question is, can the tangent stiffness matrix be calculated on an element by element basis and assembled just like the global stiffness matrix?
Let us try to figure that out.
The element stiffness matrix has the form
The tangent stiffness matrix is given by
After going through the algebra, we get
For the two element mesh, the assembled global tangent stiffness matrix is
Applying the boundary condition , we get
This matrix is the same as that shown in equation (8). Therefore, we can actually assemble the global tangent stiffness matrix. You could also show that in more general terms. But we avoid that proof here.
Doing the actual computation 
The Newton iteration scheme is
Let us go through the process. We will have to start with some assumed values of the parameters.
Let m, m, MPa, kN. Since there are two elements, m.
For the first Newton iteration, let the initial guess be
Then, the first Newton iteration gives
The second Newton iteration gives
The third Newton iteration gives
The fourth Newton iteration gives
At this stage we stop the computation because the solution has converged. A plot of the convergence of the Newton iterations in shown in Figure 1.
Once we have computed the nodal displacements, we can compute the element stresses in the usual manner.
We write the strains as
The stresses are computed using
For the axially loaded bar with not body forces, we have a uniform state of strain (and therefore stress). For the numbers that we have used, these stresses are 10 MPa in each element.
A slightly more complicated situation 
At this stage, let us add in the body force terms and work out the solution for another test case.
Let us assume that the input parameters are , , , , and . The stress-strain curve for this material is shown in Figure 2. We use a error tolerance of as a flag to stop the Newton iterations.
Mesh with 10 elements 
If we use a mesh containing equally sized elements, we get the solutions shown in Figure 3.
Mesh with 100 elements 
If we use a mesh containing equally sized elements, we get the solutions shown in Figure 4.
Even though the displacements and strains are quite different in the linear and the nonlinear cases, the stresses are remarkably exactly the same for the linear and nonlinear cases. This is a result of the form of the modulus that we have chosen to use for this problem. Notice that the error norms show nearly quadratic convergence.
There is however one issue. Due to the shape of the stress-strain curve of the material, Newton iterations do not converge when the applied external load exceeds the amount needed for the strains to exceed 0.5. A method other than Newton-Raphson is needed to solve the problem for such situations.
The Matlab script used for the above calculations is shown below.
Matlab code