写了两天的非线性程序终于可以计算结果了【附图】
最大体会是不平衡力的计算要精确！

算了一个文献上都计算的算例【悬臂梁的加载】
E=3000.;
L=5m【划分5个单元】
A=1.
I=1/12
M=60N*m【一个加载步】
算的的竖向位移为-2.6562m，水平位移是1.101m

对于UL法的编程具体实施：谈几点（以下采用两个构形：C1为最后一次计算得到的平衡构形，如初始构形；C2为计算构形）


1>计算【{K}*{X}={F}】初始得到的增量位移为总体系下的。


2>将总体系下的位移 Σ{X}变换到局部坐标系，计算局部的内力。此内力是C1构形下的单元局部节点力{F1}。此步计算单元节点力的{K}可以与第一步不同。但精确度要高于第一步的{K}。


3>用增量位移修改单元坐标。将第二步得到的{F1}旋转到新的坐标系下{F2}，从而计算得到不平衡力。


4>将不平衡力旋转到与总体{K}相同的构形中，计算下一轮的增量位移。


下图中长曲线为线性解，另为非线性解。