非线性有限元编程 | Incremental Secant Method

前几期给大家介绍了两期有关有限元非线性求解的基本数值算法：牛顿法和修正牛顿法，本期继续带着大家学习新的非线性算法——Incremental Secant Method，也就是常说的割线法。

该算法也是Kim教授书中第二章最后的非线性算法介绍了，至于弧长法和别的非线性算法可能后期会加以补充，后续将会介绍力控制加载与位移加载的区别与联系，再后来将会带来各种非线性有限元案例，敬请期待~

割线法

什么是割线法呢？

顾名思义，就是在迭代求解过程中，使用割线矩阵进行更新。

在上图中，我们可以看到，第一次迭代时，使用的是切线矩阵（与牛顿法相同），第2次迭代以后用的是割线矩阵进行更新。这样做的优点：计算过程中只需计算一次切线矩阵（求导），后续的过程不涉及公式的求导，迭代速度自然加快许多。接下来将割线法分别应用到单变量问题、多变量问题。

单变量求解

考虑使用Incremental Secant Method求解以下非线性方程：

该方程的精确解为 ，设置收敛容差 ，初始值 。

第一次迭代的切线：

每次迭代时的割线：

求解代码

xdata=zeros(40,1);
ydata=zeros(40,1);
tol = 1.0e-5; iter = 0; c = 0;
u = 2.0;
uold = u;
P = u+atan(5*u); Pold = P;
R = -P;
conv= R^2;
fprintf('\n iter u conv c');
fprintf('\n %3d %7.5f %12.3e %7.5f',iter,u,conv,c);
Ks = 1+5*(cos(atan(5*u)))^2;
while conv > tol && iter < 20
    delu = R/Ks;
    u = uold + delu;
    P = u+atan(5*u);
    R = -P;
    conv= R^2;
    c = abs(u)/abs(uold)^2;
    Ks = (P - Pold)/(u - uold);
    uold = u;
    Pold = P;
    iter = iter + 1;
    xdata(2*iter)=u; ydata(2*iter)=0;
    xdata(2*iter+1)=u; ydata(2*iter+1)=P;
    fprintf('\n %3d %7.5f %12.3e %7.5f',iter,u,conv,c);
end

plot(xdata,ydata);
hold on;
x=[-2:0.1:2];
y=x+atan(5*x);
plot(x,y)

在该程序中，设置的初值为2.0，相对较大于精确值0，但是程序在经过6次迭代后达到收敛，由此可知，Incremental Secant Method对初值的依赖性远小于之前介绍的两种方法。

多变量求解

多变量问题在迭代时，与单变量有所不同，具体在于切线矩阵和割线矩阵的形式不同。Broyden提出在使用割线法进行有限元求解时，仅第一迭代使用Jacobian matrix（切线矩阵），在后续的迭代中使用的secant stiffness matrix（割线矩阵），形式如下：

至于为什么是上式的形式，我也说不明白，记住就行了，详情可翻阅Kim教授的教材或Broyden的文献（A class of methods for solving nonlinear simultaneous equations）。接下来还是以弹簧问题，讲述割线法如何应用于有限元问题求解。

有限元问题

如下图所示，模型由两个弹簧单元组成，共3个节点，每个节点包含1个自由度，共3个自由度，模型右端受到 的集中载荷，左端固定，弹簧刚度与单元位移量相关，其中， ， ，试求 和 的值。假设单位统一，不计单位。

每个弹簧元的平衡方程如下式：

整体刚度组装，去除0左端节点项（固支条件，划行划列）后的整体刚度方程可表示如下：

进一步讲上式化简为1个非线性方程组 ：

求解代码

tol = 1.0e-5; iter = 0; c = 0;
u = [0.1; 0.1]; uold = u;
f = [0; 100];
P = [300*u(1)^2+400*u(1)*u(2)-200*u(2)^2+150*u(1)-100*u(2)
200*u(1)^2-400*u(1)*u(2)+200*u(2)^2-100*u(1)+100*u(2)];
R=P-f; Rold = R;
conv= (R(1)^2+R(2)^2)/(1+f(1)^2+f(2)^2);
fprintf('\n iter u1 u2 conv c');
fprintf('\n %3d %7.5f %7.5f %12.3e %7.5f',iter,u(1),u(2),conv,c);
Ks = [600*u(1)+400*u(2)+150 -400*u(2)+400*u(1)-100
400*u(1)-400*u(2)-100 400*u(2)-400*u(1)+100]; % 雅可比矩阵
while conv > tol && iter < 20
    delu = -Ks\R;
    u = uold + delu;
    P = [300*u(1)^2+400*u(1)*u(2)-200*u(2)^2+150*u(1)-100*u(2);
    200*u(1)^2-400*u(1)*u(2)+200*u(2)^2-100*u(1)+100*u(2)];
    R=P-f;  
    conv= (R(1)^2+R(2)^2)/(1+f(1)^2+f(2)^2);
    c = abs(0.9-u(2))/abs(0.9-uold(2))^2;
    delR = R - Rold;
    Ks = Ks + (delR-Ks*delu)*delu'/norm(delu)^2; % 割线矩阵
    uold = u; Rold = R;
    iter = iter + 1;
    fprintf('\n %3d %7.5f %7.5f %12.3e %7.5f',iter,u(1),u(2),conv,c);
end

该程序经过5次迭代达到收敛，而且在设置初值时，不用考虑 这一物理关系，不依赖于初值的选定，具有较好的收敛性。

思考

带上这期的推文，木木共分享了三种非线性求解算法，在学习之余，可以看出算法背后的意义，为什么会出现这么多算法？

影响非线性求解速度的两个关键点：

• 迭代时所要形成矩阵的方式
• 每次求解矩阵方程的时间

由以上两个方面可以看出，所有的非线性算法都在聚焦于如何在保证精度的前提下，加速收敛求解。明白了求解目的后，便会对其算法本身有一定的“大局观”，才不会惧怕枯燥且抽象的公式。