kill(all);

assume(EA>0,EJ>0);

assume(l>0,dx>0);

y,v A 
    |
    |
 \\||              F
 \\|O==============O ---->   ---> x,u
 \\|              /_\ 
               ---------
                ///////
x,y global coordinates
u,v displacement according to the global reference system

two motion modes are considered.
a) an axial elongation, in the absence of tranverse displacements
b) a tranverse motion, sinusoidal with respect to x

motion a

define( u(x)    , (a + da)*x);

motion b

 define(v(x)    , (0+db) * sin(1*%pi*x/l) );

derivatives of the two motion modes

define( up(x)    ,diff(u(x),x));

define(vp(x)   , diff(v(x),x)          );

define(vpp(x)  , diff(vp(x),x)         );

auxiliary functions for the nonlinear calculations:
- length of an A->B segment

l(XA,YA,XB,YB) := sqrt((XB-XA)^2 + (YB-YA)^2);

- cosinus of the angle comprised between an A->B segment, and the x axis

c(XA,YA,XB,YB):=(XB-XA)/l(XA,YA,XB,YB);

- nonlinear definition for the beam axial strain

tmp: fullratsimp( 
    sqrt( 
          (     v(x+dx)  -v(x))^2 
        + (x+dx+u(x+dx)-x-u(x))^2 

    )/dx -1
);

define(eps_nl(x) ,  limit(tmp,dx,0,plus) );

define(eps_nl(x) ,  expand(taylor(eps_nl(x), [a,0,1], [da,0,2] , [db,0,2] )) );

I cut&paste only the terms up to the second order in the "small" "a,da,db" variables

define(eps_nl(x) ,(%pi^2*db^2*cos((%pi*x)/l)^2)/(2*l^2)+da+a);

define(eps_nl(x) ,taylor((sqrt(1+vp(x)^2)-1),db,0,2)+da+a);

nonlinear energetic term

Ub_nl:  integrate( 1/2 * EJ * vpp(x)^2    ,x,0,l);
Un_nl:  integrate( 1/2 * EA * eps_nl(x)^2 ,x,0,l);

we calculate the hessian matrix in the neighborhood of the 
preloaded equilibrium condition, i.e. 
-  da = db = 0;
-  a  = F/EA;
-  b  = 0; 
such hessian matrix becomes the tangential stiffness matrix, i.e.
the matrix that describes the linearized relation between 
generalized load and displacement configurations. 

hessian( Ub_nl + Un_nl,[da,db]);
K : ev(%,[da=0,db=0]);

K : ev(K,a=F/EA);

if the following equation holds, the tangential stiffness matrix is singular;
such condition paves the way for the bending solution branches to emerge.

linsolve(K[2,2]=0,F);

the relation above defines the critical load for which solutions appear
other than the linear elastic one that was in continuity with the
unloaded/undeformed condition.
Also, at the critical load the tangent stiffness matrix loses its positive
definiteness; for loads beyond the critical value, the linear elastic solution
defines an unstable equilibrium condition.

we separate the stiffness matrix into the two following contributions

material stiffness matrix

K0:ev(K,F=0);

preload related, "geometric" stiffness matrix

K1:diff(K,F)*F;

this second contribution disappears if the strain derivation from the generalized
displacements is linearized.

The preload effect is the amplified by a lambda factor, and such condition is 
assumed equivalent to the effect of an amplified preload. 

Kt : K0 + lambda * K1;

both the Kt terms are left-multiplied by 1/lambda*invert(K0), 
then the variable substition

xi = -1/lambda is performed

to obtain the standard eigenproblem

( A - xi*I).x = 0

A : invert(K0).K1;

I : invert(K0).K0;

we now solve the eigenproblem

vals : eigenvalues
mult : eigenvalue multiplicities
vecs : eigenvectors

[[vals,mult],vecs]:eigenvectors(A);

we focus on the nonzero eigenvalue

vals[1];

the related eigenvector denotes the bending solution branch that suddenly appears,

[da,db] = vecs[1][1];