(************** Content-type: application/mathematica **************
                     CreatedBy='Mathematica 4.2'

                    Mathematica-Compatible Notebook

This notebook can be used with any Mathematica-compatible
application, such as Mathematica, MathReader or Publicon. The data
for the notebook starts with the line containing stars above.

To get the notebook into a Mathematica-compatible application, do
one of the following:

* Save the data starting with the line of stars above into a file
  with a name ending in .nb, then open the file inside the
  application;

* Copy the data starting with the line of stars above to the
  clipboard, then use the Paste menu command inside the application.

Data for notebooks contains only printable 7-bit ASCII and can be
sent directly in email or through ftp in text mode.  Newlines can be
CR, LF or CRLF (Unix, Macintosh or MS-DOS style).

NOTE: If you modify the data for this notebook not in a Mathematica-
compatible application, you must delete the line below containing
the word CacheID, otherwise Mathematica-compatible applications may
try to use invalid cache data.

For more information on notebooks and Mathematica-compatible 
applications, contact Wolfram Research:
  web: http://www.wolfram.com
  email: info@wolfram.com
  phone: +1-217-398-0700 (U.S.)

Notebook reader applications are available free of charge from 
Wolfram Research.
*******************************************************************)

(*CacheID: 232*)


(*NotebookFileLineBreakTest
NotebookFileLineBreakTest*)
(*NotebookOptionsPosition[     20797,        457]*)
(*NotebookOutlinePosition[     21875,        491]*)
(*  CellTagsIndexPosition[     21831,        487]*)
(*WindowFrame->Normal*)



Notebook[{
Cell["\<\
Cell 1.  Modules to form the internal force vector of the  \
geometrically nonlinear, C0 plane beam.\
\>", "Text",
  CellFrame->True,
  CellMargins->{{18, 43}, {Inherited, Inherited}},
  Evaluatable->False,
  CellLabelMargins->{{9, Inherited}, {Inherited, Inherited}},
  AspectRatioFixed->True,
  Background->RGBColor[0, 1, 1]],

Cell[TextData[{
  StyleBox["FormIntForceC0TwoNodePlaneBeam[XY1_,XY2_,u1_,u2_,S0_,z0_]:=Module[\
\n {X1,Y1,X2,Y2,X21,Y21,L0,uX1,uY1,theta1,uX2,uY2,theta2,x21,y21,\n  \
uX0,uY0,thetam,ctheta,stheta,Lcpsi,Lspsi,N0,V0,M0,\n  \
EA0,GA0,EI0,cphi,sphi,cm,sm,Nm,Vm,Mm,Bm,kappa,p},\n  {X1,Y1}=XY1; \
{X2,Y2}=XY2; X21=X2-X1; Y21=Y2-Y1;\n  {uX1,uY1,theta1}=u1; \
{uX2,uY2,theta2}=u2;\n   x21=X21+uX2-uX1; y21=Y21+uY2-uY1; \n   \
L0=PowerExpand[Sqrt[X21^2+Y21^2]]; thetam=(theta1+theta2)/2; \n   \
ctheta=Cos[thetam]; stheta=Sin[thetam];\n   \
Lcpsi=Simplify[(X21*x21+Y21*y21)/L0]; \n   \
Lspsi=Simplify[(X21*y21-Y21*x21)/L0];\n   em= \
(ctheta*Lcpsi+stheta*Lspsi)/L0-1;\n   gm=-(stheta*Lcpsi-ctheta*Lspsi)/L0;\n   \
kappa=(theta2-theta1)/L0;\n   {N0,V0,M0}=z0; {EA0,GA0,EI0}=S0;\n   \
Nm=Simplify[N0+EA0*em];  Vm=Simplify[V0+GA0*gm];  \n   \
Mm=Simplify[M0+EI0*kappa];\n   cphi=X21/L0; sphi=Y21/L0; \n   \
cm=ctheta*cphi-stheta*sphi; sm=stheta*cphi+ctheta*sphi;\n   Bm= \
(1/L0)*{{-cm,-sm, L0*gm/2,     cm,sm, L0*gm/2    }, \n               { \
sm,-cm,-L0*(1+em)/2,-sm,cm,-L0*(1+em)/2},\n               {  0, 0,  -1,       \
   0,  0,          1 }};\n   p=L0*Transpose[Bm].{{Nm},{Vm},{Mm}};\n   \
Return[Simplify[p]] ];\n      ",
    CellFrame->True,
    CellMargins->{{18, 43}, {Inherited, Inherited}},
    InitializationCell->True,
    CellLabelMargins->{{9, Inherited}, {Inherited, Inherited}},
    AspectRatioFixed->True,
    Background->RGBColor[1, 0.738933, 0.130022]],
  StyleBox[" \n  \
p=FormIntForceC0TwoNodePlaneBeam[{0,0},{L/Sqrt[2],L/Sqrt[2]},\n         \
{0,0,Pi/2},{-2*L/Sqrt[2],0,Pi/2},\n         {EA,GA0,EI},{0,V0,0}];\n  \
Print[\"p=\",p];",
    CellFrame->True,
    CellMargins->{{18, 43}, {Inherited, Inherited}},
    InitializationCell->True,
    CellLabelMargins->{{9, Inherited}, {Inherited, Inherited}},
    AspectRatioFixed->True,
    FontColor->RGBColor[1, 0, 0],
    Background->RGBColor[1, 0.738933, 0.130022]]
}], "Input",
  CellFrame->True,
  CellMargins->{{18, 43}, {Inherited, Inherited}},
  InitializationCell->True,
  CellLabelMargins->{{9, Inherited}, {Inherited, Inherited}},
  AspectRatioFixed->True,
  Background->RGBColor[1, 0.738933, 0.130022]],

Cell[TextData[StyleBox[
"Cell 2.  Module to form the material component of the tangent stiffness \
matrix of the \ngeometrically nonlinear, C0 plane beam.",
  CellFrame->True,
  CellMargins->{{18, 43}, {Inherited, Inherited}},
  CellLabelMargins->{{9, Inherited}, {Inherited, Inherited}},
  AspectRatioFixed->True,
  FontFamily->"Times",
  FontWeight->"Plain",
  Background->RGBColor[0, 1, 1]]], "Input",
  CellFrame->True,
  CellMargins->{{18, 43}, {Inherited, Inherited}},
  CellLabelMargins->{{9, Inherited}, {Inherited, Inherited}},
  AspectRatioFixed->True,
  Background->RGBColor[0, 1, 1]],

Cell[TextData[{
  StyleBox[
  "FormMatStiffC0TwoNodePlaneBeam[XY1_,XY2_,u1_,u2_,S0_,z0_]:=Module[\n  \
{X1,Y1,X2,Y2,X21,Y21,L0,uX1,uY1,theta1,uX2,uY2,theta2,x21,y21,\n   \
uX0,uY0,thetam,ctheta,stheta,Lcpsi,Lspsi,N0,V0,M0,\n   \
EA0,GA0,EI0,cphi,sphi,cm,sm,a1,KM},\n  {X1,Y1}=XY1; {X2,Y2}=XY2; X21=X2-X1; \
Y21=Y2-Y1;\n  {uX1,uY1,theta1}=u1; {uX2,uY2,theta2}=u2;\n   x21=X21+uX2-uX1; \
y21=Y21+uY2-uY1; \n   L0=PowerExpand[Sqrt[X21^2+Y21^2]]; \
thetam=(theta1+theta2)/2; \n   ctheta=Cos[thetam]; stheta=Sin[thetam];\n   \
Lcpsi=Simplify[(X21*x21+Y21*y21)/L0]; \n   \
Lspsi=Simplify[(X21*y21-Y21*x21)/L0];\n   em= \
(ctheta*Lcpsi+stheta*Lspsi)/L0-1;\n   gm=-(stheta*Lcpsi-ctheta*Lspsi)/L0;\n   \
cphi=X21/L0; sphi=Y21/L0; \n   cm=ctheta*cphi-stheta*sphi; \
sm=stheta*cphi+ctheta*sphi;\n   {N0,V0,M0}=z0; {EA0,GA0,EI0}=S0; a1=1+em;\n   \
KM = (EA0/L0)*{\n   { cm^2,cm*sm,-cm*gm*L0/2,-cm^2,-cm*sm,-cm*gm*L0/2},\n   { \
cm*sm,sm^2,-gm*L0*sm/2,-cm*sm,-sm^2,-gm*L0*sm/2},\n   \
{-cm*gm*L0/2,-gm*L0*sm/2,gm^2*L0^2/4,\n     \
cm*gm*L0/2,gm*L0*sm/2,gm^2*L0^2/4},\n   \
{-cm^2,-cm*sm,cm*gm*L0/2,cm^2,cm*sm,cm*gm*L0/2},\n   \
{-cm*sm,-sm^2,gm*L0*sm/2,cm*sm,sm^2,gm*L0*sm/2},\n   \
{-cm*gm*L0/2,-gm*L0*sm/2,gm^2*L0^2/4,\n     \
cm*gm*L0/2,gm*L0*sm/2,gm^2*L0^2/4}}+\n   (EI0/L0)*{{0,0,0,0,0,0},  \
{0,0,0,0,0,0},\n             {0,0,1,0,0,-1}, {0,0,0,0,0,0},\n             \
{0,0,0,0,0,0},  {0,0,-1,0,0,1}}+\n   (GA0/L0)*{\n   \
{sm^2,-cm*sm,-a1*L0*sm/2,-sm^2,cm*sm,-a1*L0*sm/2},\n   \
{-cm*sm,cm^2,cm*a1*L0/2,cm*sm,-cm^2,cm*a1*L0/2},\n   {-a1*L0*sm/2,cm*a1*L0/2, \
a1^2*L0^2/4,\n     a1*L0*sm/2,-cm*a1*L0/2, a1^2 L0^2/4}, \n   {-sm^2,cm*sm, \
a1*L0*sm/2,sm^2,-cm*sm,a1*L0*sm/2},\n   { \
cm*sm,-cm^2,-cm*a1*L0/2,-cm*sm,cm^2,-cm*a1*L0/2},\n   \
{-a1*L0*sm/2,cm*a1*L0/2,a1^2*L0^2/4,\n     a1*L0*sm/2,-cm*a1*L0/2, \
a1^2*L0^2/4}}; \n   Return[KM] ];\n   \n\
FormGeoStiffC0TwoNodePlaneBeam[XY1_,XY2_,u1_,u2_,S0_,z0_]:=Module[\n  \
{X1,Y1,X2,Y2,X21,Y21,L0,uX1,uY1,theta1,uX2,uY2,theta2,x21,y21,\n   \
uX0,uY0,thetam,ctheta,stheta,Lcpsi,Lspsi,N0,V0,M0,\n   \
EA0,GA0,EI0,cphi,sphi,cm,sm,Nm,Vm,KG},\n   {X1,Y1}=XY1; {X2,Y2}=XY2; \
X21=X2-X1; Y21=Y2-Y1;\n   {uX1,uY1,theta1}=u1; {uX2,uY2,theta2}=u2;\n   \
x21=X21+uX2-uX1; y21=Y21+uY2-uY1; \n   L0=PowerExpand[Sqrt[X21^2+Y21^2]]; \
thetam=(theta1+theta2)/2; \n   ctheta=Cos[thetam]; stheta=Sin[thetam];\n   \
Lcpsi=Simplify[(X21*x21+Y21*y21)/L0]; \n   \
Lspsi=Simplify[(X21*y21-Y21*x21)/L0];\n   em= \
(ctheta*Lcpsi+stheta*Lspsi)/L0-1;\n   gm=-(stheta*Lcpsi-ctheta*Lspsi)/L0;\n   \
kappa=(theta2-theta1)/L0;\n   {N0,V0,M0}=z0; {EA0,GA0,EI0}=S0;\n   \
Nm=Simplify[N0+EA0*em];  Vm=Simplify[V0+GA0*gm];  \n   cphi=X21/L0; \
sphi=Y21/L0; \n   cm=ctheta*cphi-stheta*sphi; sm=stheta*cphi+ctheta*sphi;   \n\
   KG = Nm/2*{{ 0, 0, sm, 0, 0, sm}, {0, 0, -cm, 0, 0,-cm},\n              \
{sm, -cm, -(L0/2)*(1+em), -sm, cm,-(L0/2)*(1+em)},\n              {0, 0, -sm, \
0, 0, -sm}, {0, 0, cm, 0, 0, cm},\n              {sm, -cm, -(L0/2)*(1+em), \
-sm, cm,-(L0/2)*(1+em)}}+\n        Vm/2*{{0, 0, cm, 0, 0, cm}, {0, 0, sm, 0, \
0, sm},\n              {cm, sm, -(L0/2)*gm, -cm, -sm,-(L0/2)*gm},\n           \
   {0, 0, -cm, 0, 0, -cm}, {0, 0, -sm, 0, 0, -sm},\n              {cm, sm, \
-(L0/2)*gm, -cm, -sm,-(L0/2)*gm}};\n   Return[KG] ];\n   \n\
FormTanStiffC0TwoNodePlaneBeam[XY1_,XY2_,u1_,u2_,S0_,z0_]:=\n  Module[{},\n  \
KM=FormMatStiffC0TwoNodePlaneBeam[XY1,XY2,u1,u2,S0,z0];\n  \
KG=FormGeoStiffC0TwoNodePlaneBeam[XY1,XY2,u1,u2,S0,z0];\n  Return[KM+KG]];\n",
    
    CellFrame->True,
    CellMargins->{{18, 43}, {Inherited, Inherited}},
    InitializationCell->True,
    CellLabelMargins->{{9, Inherited}, {Inherited, Inherited}},
    AspectRatioFixed->True,
    Background->RGBColor[1, 0.738933, 0.130022]],
  StyleBox[
  "\nClearAll[EA,GA,EI];\n\
KM=FormMatStiffC0TwoNodePlaneBeam[{0,0},{10,0},{0,0,0},{0,0,0},\n    \
{EA,GA,EI},{0,0,0}];\n Print[KM];\n(* Print[Chop[Eigensystem[N[KM]]]];*)\n  \n\
KG=FormGeoStiffC0TwoNodePlaneBeam[{0,0},{4,3},{0,0,0},{0,0,0},\n    \
{1,1,1},{10,30,20}];\n Print[N[KG]];\n Print[Chop[Eigenvalues[N[KG]]]];",
    CellFrame->True,
    CellMargins->{{18, 43}, {Inherited, Inherited}},
    InitializationCell->True,
    CellLabelMargins->{{9, Inherited}, {Inherited, Inherited}},
    AspectRatioFixed->True,
    FontColor->RGBColor[1, 0, 0],
    Background->RGBColor[1, 0.738933, 0.130022]],
  StyleBox["  ",
    CellFrame->True,
    CellMargins->{{18, 43}, {Inherited, Inherited}},
    InitializationCell->True,
    CellLabelMargins->{{9, Inherited}, {Inherited, Inherited}},
    AspectRatioFixed->True,
    Background->RGBColor[1, 0.738933, 0.130022]],
  StyleBox["   ",
    CellFrame->True,
    CellMargins->{{18, 43}, {Inherited, Inherited}},
    InitializationCell->True,
    CellLabelMargins->{{9, Inherited}, {Inherited, Inherited}},
    AspectRatioFixed->True,
    FontColor->RGBColor[1, 0, 0],
    Background->RGBColor[1, 0.738933, 0.130022]]
}], "Input",
  CellFrame->True,
  CellMargins->{{18, 43}, {Inherited, Inherited}},
  InitializationCell->True,
  CellLabelMargins->{{9, Inherited}, {Inherited, Inherited}},
  AspectRatioFixed->True,
  Background->RGBColor[1, 0.738933, 0.130022]],

Cell[TextData[
"Cell 3.  This is a Finite Difference check (by Central FDs)  that K = \
partial p /partial u"], "Text",
  CellFrame->True,
  CellMargins->{{18, 43}, {Inherited, Inherited}},
  Evaluatable->False,
  CellLabelMargins->{{9, Inherited}, {Inherited, Inherited}},
  AspectRatioFixed->True,
  Background->RGBColor[0, 1, 1]],

Cell[TextData[
"ClearAll[EA,GA,EI,N0,V0,M0];\n XY1={-3,7}; XY2={-8,11}; u1={-.6,-.9,0}; \
u2={1.5,2.1,0};\n{X1,Y1}=XY1; {X2,Y2}=XY2; X21=X2-X1; \
Y21=Y2-Y1;L0=Sqrt[X21^2+Y21^2];\n u1={0.1,0.25,-0.3}; u2={0.45,0.58,-0.6};  \n\
 EA=100; GA=240; EI=300;  (* GA=12*EI/L0^2; RBF *)\n S={EA,GA,EI}; \
z0={10,-20,30}; eps=0.001;\n u0=Flatten[{u1,u2}]; KFD=Table[0,{6},{6}];\n \
(*p0=FormIntForceC0TwoNodePlaneBeam[XY1,XY2,u1,u2,S,z0];*)\n For [i=1, i<=6, \
i++,\n   u=u0; u[[i]]=u[[i]]+eps; u1=Take[u,3]; u2=Take[u,-3];\n   pplus= \
FormIntForceC0TwoNodePlaneBeam[XY1,XY2,u1,u2,S,z0];\n   u=u0; \
u[[i]]=u[[i]]-eps; u1=Take[u,3]; u2=Take[u,-3];\n   \
pminus=FormIntForceC0TwoNodePlaneBeam[XY1,XY2,u1,u2,S,z0];\n   \
KFD[[i]]=Flatten[Simplify[(pplus-pminus)/(2*eps)]] ];\n KFD=N[KFD];\n \
Print[\"K by finite differences:\"];Print[KFD//MatrixForm];\n u=u0; \
u1=Take[u,3]; u2=Take[u,-3];\n \
K=FormTanStiffC0TwoNodePlaneBeam[XY1,XY2,u1,u2,S,z0]; K=N[K];\n Print[\"K by \
stiffness modules:\"]; Print[K//MatrixForm];\n Print[\"This diff should be \
zero or small:\"];\n Print[Chop[K-KFD,0.0001]//MatrixForm];\n \
Print[Chop[Eigenvalues[K]]];\n Print[Chop[Eigenvalues[KFD]]];"], "Input",
  CellFrame->True,
  CellMargins->{{18, 43}, {Inherited, Inherited}},
  CellLabelMargins->{{9, Inherited}, {Inherited, Inherited}},
  AspectRatioFixed->True,
  Background->RGBColor[0, 1, 0]],

Cell[TextData[
"Cell 4.  These are two utility modules to merge element forces and stiffness \
into master p and master K."], "Text",
  CellFrame->True,
  CellMargins->{{18, 43}, {Inherited, Inherited}},
  Evaluatable->False,
  CellLabelMargins->{{9, Inherited}, {Inherited, Inherited}},
  AspectRatioFixed->True,
  Background->RGBColor[0, 1, 1]],

Cell[TextData[
"MergeElemIntoMasterIntForce[pe_,eftab_,pm_]:=   \n  \
Module[{i,ii,nf=Length[eftab],p}, p=pm; \n       For[i=1, i<=nf, i++, \
ii=eftab[[i]];\n           If [ii>0,p[[ii,1]]+=pe[[i,1]]] \n          ]; \
Return[p]\n  ];\n \nMergeElemIntoMasterStiff[Ke_,eftab_,Km_]:=   \n  \
Module[{i,j,ii,jj,nf=Length[eftab],K}, K=Km; \n       For[i=1, i<=nf, i++, \
ii=eftab[[i]]; If[ii==0,Continue[]];\n         For[j=i, j<=nf, j++, \
jj=eftab[[j]];\n           If [ii>0 && jj>0, \n           \
K[[jj,ii]]=K[[ii,jj]]+=Ke[[i,j]]] \n            ]\n          ]; Return[K]\n   \
     ];"], "Input",
  CellFrame->True,
  CellMargins->{{18, 43}, {Inherited, Inherited}},
  InitializationCell->True,
  CellLabelMargins->{{9, Inherited}, {Inherited, Inherited}},
  AspectRatioFixed->True,
  Background->RGBColor[1, 0.703273, 0.110018]],

Cell[TextData[StyleBox[
"Cell 5. The following  modules assemble the internal force vector and the \
master stiffness matrix of\na compressed cantilever beam-column  of total \
length L, clamped on the left, and discretized\nwith Ne elements.  The \
structure is subjected to an applied axial end force P, positive\n o the \
right (so that P<0 means the beam-column is compressed)",
  CellFrame->True,
  CellMargins->{{18, 43}, {Inherited, Inherited}},
  CellLabelMargins->{{9, Inherited}, {Inherited, Inherited}},
  AspectRatioFixed->True,
  FontFamily->"Times",
  FontWeight->"Plain",
  Background->RGBColor[0, 1, 1]]], "Input",
  CellFrame->True,
  CellMargins->{{18, 43}, {Inherited, Inherited}},
  CellLabelMargins->{{9, Inherited}, {Inherited, Inherited}},
  AspectRatioFixed->True,
  Background->RGBColor[0, 1, 1]],

Cell[TextData[
"AssembleMasterIntForceOfCantBeam[{L_,P_,Ne_},{Em_,A0_,I0_},u_]:=\n  \
Module[{e,numnod,numdof,u1,u2,u3,Le,X1,X2,S0,z0,eftab,\n          pe,p},\n   \
numnod=Ne+1; numdof=3*numnod; Le=L/Ne;\n   p=Table[0,{numdof},{1}];\n   \
z0={P,0,0}; S0={Em*A0,12*Em*I0/Le^2,Em*I0}; (* RBF *)\n   For [e=1, e<=Ne, \
e++,\n     X1=(e-1)*Le; X2=e*Le; \n     If [e==1,u1={0,0,0}; \
u2=Take[u,{3*e-2,3*e}];\n         eftab={0,0,0,3*e-2,3*e-1,3*e}];\n     If \
[e>1, u1=Take[u,{3*e-5,3*e-3}]; u2=Take[u,{3*e-2,3*e}];\n         \
eftab={3*e-5,3*e-4,3*e-3,3*e-2,3*e-1,3*e}];\n     \
pe=FormIntForceC0TwoNodePlaneBeam[{X1,0},{X2,0},u1,u2,S0,z0];\n     \
(*Print[\"pe=\",pe];*)\n     p=MergeElemIntoMasterIntForce[pe,eftab,p];\n   \
];\n   Return[Simplify[p]]\n ];\n \n\
AssembleMasterStiffOfCantBeam[{L_,P_,Ne_},{Em_,A0_,I0_},u_]:=\n  \
Module[{e,numnod,numdof,u1,u2,u3,Le,X1,X2,S0,z0,eftab,\n          pe,p},\n   \
numnod=Ne+1; numdof=3*numnod-3; Le=L/Ne;\n   K=Table[0,{numdof},{numdof}];\n  \
 z0={P,0,0}; S0={Em*A0,12*Em*I0/Le^2,Em*I0}; (* RBF *)\n   For [e=1, e<=Ne, \
e++,\n     X1=(e-1)*Le; X2=e*Le; \n     If [e==1,u1={0,0,0}; \
u2=Take[u,{3*e-2,3*e}];\n         eftab={0,0,0,3*e-2,3*e-1,3*e}];\n     If \
[e>1, u1=Take[u,{3*e-5,3*e-3}]; u2=Take[u,{3*e-2,3*e}];\n         \
eftab={3*e-5,3*e-4,3*e-3,3*e-2,3*e-1,3*e}];\n     \
Ke=FormTanStiffC0TwoNodePlaneBeam[{X1,0},{X2,0},u1,u2,S0,z0];\n     \
(*Print[\"Ke=\",Ke];*)\n     K=MergeElemIntoMasterStiff[Ke,eftab,K];\n   ];\n \
  Return[Simplify[K]]\n ];\n "], "Input",
  CellFrame->True,
  CellMargins->{{18, 43}, {Inherited, Inherited}},
  InitializationCell->True,
  CellLabelMargins->{{9, Inherited}, {Inherited, Inherited}},
  AspectRatioFixed->True,
  Background->RGBColor[1, 0.703273, 0.110018]],

Cell[TextData[StyleBox[
" Cell 6. The following script calculates the buckling load factor of the \
compressed beam-column for an\n increasing number of elements: 1,2, 4,... \
etc, using several of the  foregoing modules in cells 1,2,4,5.   \n The \
calculation is brute force, looking for the roots of the exact characteristic \
polynomial of K\n  with NSolve.  Method breaks  down at about 16 or 32 \
elements, depending on computational power, \n  because the coefficients of \
det[K] grow very rapidly.   \n \n This is the topic of Ex 9.3.",
  CellFrame->True,
  CellMargins->{{18, 43}, {Inherited, Inherited}},
  CellLabelMargins->{{9, Inherited}, {Inherited, Inherited}},
  AspectRatioFixed->True,
  FontFamily->"Times",
  FontWeight->"Plain",
  Background->RGBColor[0, 1, 1]]], "Input",
  CellFrame->True,
  CellMargins->{{18, 43}, {Inherited, Inherited}},
  CellLabelMargins->{{9, Inherited}, {Inherited, Inherited}},
  AspectRatioFixed->True,
  Background->RGBColor[0, 1, 1]],

Cell[TextData[
"ClearAll[L,P,Em,A0,I0]; Em=A0=I0=L=1;  Ne=1;\nFor [k=1, k<=6, k++, \n   \
Print[\"Number of elements=\",Ne];  \n   LPNe={L,P,Ne}; EAI={Em,A0,I0}; \
u=Table[0,{3*Ne+3}];  \n   K=AssembleMasterStiffOfCantBeam[LPNe,EAI,u]; \n   \
(*Print[\"K=\",K//MatrixForm];*)\n   detK=Det[K]; \n   \
(*Print[\"det(K)=\",detK//InputForm];*)\n   roots=NSolve[detK==0,P];\n   \
Print[\"roots of stability det=\",roots];  Ne=2*Ne;\n ];\nPrint[\"exact \
buckling load coeff is \",-N[Pi^2/4]];"], "Input",
  CellFrame->True,
  CellMargins->{{18, 43}, {Inherited, Inherited}},
  CellLabelMargins->{{9, Inherited}, {Inherited, Inherited}},
  AspectRatioFixed->True,
  Background->RGBColor[1, 0.703273, 0.110018]],

Cell[TextData[
"Cell 7. This is another, quicker and more stable method, to get the critical \
buckling load factor. \nNeeds modules in cells 1,2,4,5.   On my Mac it worrks \
up to 128 elements. \n\nIt is the subject of Exercise 9.4."], "Text",
  CellFrame->True,
  CellMargins->{{18, 43}, {Inherited, Inherited}},
  Evaluatable->False,
  CellLabelMargins->{{9, Inherited}, {Inherited, Inherited}},
  AspectRatioFixed->True,
  Background->RGBColor[0, 1, 1]],

Cell[TextData[{
  StyleBox[
  "ClearAll[L,P,Em,A0,I0]; Em=A0=I0=L=1;  Ne=1;\nFor [k=1, k<=8, k++, \n   \
Print[\"Number of elements=\",Ne];  \n   LPNe={L,P,Ne}; EAI={Em,A0,I0}; \
u=Table[0,{3*Ne+3}];  \n   K=AssembleMasterStiffOfCantBeam[LPNe,EAI,u]; \n   \
KM=Coefficient[K,P,0];  KG=Coefficient[K,P,1]; \n   \
SG=LinearSolve[N[KM],N[KG]];\n ",
    CellFrame->True,
    CellMargins->{{18, 43}, {Inherited, Inherited}},
    CellLabelMargins->{{9, Inherited}, {Inherited, Inherited}},
    AspectRatioFixed->True,
    Background->RGBColor[1, 0.738933, 0.130022]],
  StyleBox[
  "  (*Print[\"KM=\",KM//MatrixForm];  Print[\"KG=\",KG//MatrixForm];*)\n   \
(*Print[\"SG=\",SG//MatrixForm];*) ",
    CellFrame->True,
    CellMargins->{{18, 43}, {Inherited, Inherited}},
    CellLabelMargins->{{9, Inherited}, {Inherited, Inherited}},
    AspectRatioFixed->True,
    FontColor->RGBColor[1, 0, 0],
    Background->RGBColor[1, 0.738933, 0.130022]],
  StyleBox[
  "\n   emax=-Max[Eigenvalues[SG]]; \n   Print[\"FEM lambda \
cr=\",(1/emax)//InputForm];\n   Ne=2*Ne;\n ];\nPrint[\"exact buckling lambda \
coeff is \",-N[Pi^2/4]//InputForm];",
    CellFrame->True,
    CellMargins->{{18, 43}, {Inherited, Inherited}},
    CellLabelMargins->{{9, Inherited}, {Inherited, Inherited}},
    AspectRatioFixed->True,
    Background->RGBColor[1, 0.738933, 0.130022]]
}], "Input",
  CellFrame->True,
  CellMargins->{{18, 43}, {Inherited, Inherited}},
  CellLabelMargins->{{9, Inherited}, {Inherited, Inherited}},
  AspectRatioFixed->True,
  Background->RGBColor[1, 0.738933, 0.130022]],

Cell[TextData[
"Extrapolation applying the Shanks transformation yields -Pi^2/4 to 9 \
digits."], "Text",
  CellFrame->True,
  CellMargins->{{18, 43}, {Inherited, Inherited}},
  Evaluatable->False,
  CellLabelMargins->{{9, Inherited}, {Inherited, Inherited}},
  AspectRatioFixed->True,
  Background->RGBColor[0, 1, 1]],

Cell["\<\
Shanks[Z_]:= Module[{m,n,shankedZ},
  If [Length[Dimensions[Z]]<=1, n=Length[Z];  
      If [n<=2, Return[{}]]; shankedZ=Table[0,{n-2}];
      Do [
      shankedZ[[j]]=( Z[[j+2]]*  Z[[j]]  -Z[[j+1]]^2)/
                     (Z[[j+2]]-2*Z[[j+1]]+Z[[j]]),
      {j,1,n-2}]; Return[shankedZ];
     ];
  {m,n}=Dimensions[Z]; 
  If [n<=2, Return[{{}}]];
  shankedZ=Table[0,{m},{n-2}];
  Do [Do [
      shankedZ[[i,j]]=( Z[[i,j+2]]*  Z[[i,j]]  -Z[[i,j+1]]^2)/
                       (Z[[i,j+2]]-2*Z[[i,j+1]]+Z[[i,j]]),
   {j,1,n-2}], {i,1,m}];
  Return[shankedZ]];
  
Z={-3.165151389911680,-2.603908889389283,-2.499688345543354,
   -2.475364052205859,-2.469385123708907,-2.467896687873336,
   -2.467524972077594,-2.467432067643745};
Print[\"Shanked sequence=\",Shanks[Z]//InputForm,
     \" vs. exact: -2.467401100272339\"];
Print[\"Re-Shanked sequence=\",Shanks[Shanks[Z]]//InputForm,
     \" vs. exact: -2.467401100272339\"];
Print[\"Re-Shanked sequence=\",Shanks[Shanks[Shanks[Z]]]//InputForm,
     \" vs. exact: -2.467401100272339\"];\
\>", "Input",
  CellFrame->True,
  CellMargins->{{18, 43}, {Inherited, Inherited}},
  CellLabelMargins->{{9, Inherited}, {Inherited, Inherited}},
  AspectRatioFixed->True,
  Background->RGBColor[0.768933, 0.939986, 0.206806]]
},
FrontEndVersion->"4.2 for Macintosh",
ScreenRectangle->{{0, 1920}, {0, 1180}},
AutoGeneratedPackage->None,
WindowToolbars->{},
CellGrouping->Manual,
WindowSize->{950, 1100},
WindowMargins->{{12, Automatic}, {Automatic, 13}},
PrivateNotebookOptions->{"ColorPalette"->{RGBColor, -1}},
ShowCellLabel->True,
ShowCellTags->False,
RenderingOptions->{"ObjectDithering"->True,
"RasterDithering"->False},
MacintoshSystemPageSetup->"\<\
00<0001804P000000]P2:?oQon82n@960dL5:0?l0080001804P000000]P2:001
0000I00000400`<300000BL?00400@00000000000000060801T1T00000000000
00000000000000000000000000000000\>"
]

(*******************************************************************
Cached data follows.  If you edit this Notebook file directly, not
using Mathematica, you must remove the line containing CacheID at
the top of  the file.  The cache data will then be recreated when
you save this file from within Mathematica.
*******************************************************************)

(*CellTagsOutline
CellTagsIndex->{}
*)

(*CellTagsIndex
CellTagsIndex->{}
*)

(*NotebookFileOutline
Notebook[{
Cell[1754, 51, 337, 9, 46, "Text",
  Evaluatable->False],
Cell[2094, 62, 2173, 43, 463, "Input",
  InitializationCell->True],
Cell[4270, 107, 595, 14, 56, "Input"],
Cell[4868, 123, 5150, 95, 1273, "Input",
  InitializationCell->True],
Cell[10021, 220, 332, 8, 47, "Text",
  Evaluatable->False],
Cell[10356, 230, 1368, 22, 598, "Input"],
Cell[11727, 254, 346, 8, 64, "Text",
  Evaluatable->False],
Cell[12076, 264, 827, 15, 268, "Input",
  InitializationCell->True],
Cell[12906, 281, 821, 17, 118, "Input"],
Cell[13730, 300, 1750, 28, 598, "Input",
  InitializationCell->True],
Cell[15483, 330, 984, 19, 126, "Input"],
Cell[16470, 351, 704, 12, 208, "Input"],
Cell[17177, 365, 457, 9, 94, "Text",
  Evaluatable->False],
Cell[17637, 376, 1563, 35, 238, "Input"],
Cell[19203, 413, 318, 8, 46, "Text",
  Evaluatable->False],
Cell[19524, 423, 1269, 32, 418, "Input"]
}
]
*)



(*******************************************************************
End of Mathematica Notebook file.
*******************************************************************)