前幾期給大家介紹了兩期有關有限元非線性求解的基本數值算法：牛頓法和修正牛頓法，本期繼續帶著大家學習新的非線性算法——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次迭代達到收斂，而且在設置初值時，不用考慮 這一物理關系，不依賴于初值的選定，具有較好的收斂性。

思考

帶上這期的推文，木木共分享了三種非線性求解算法，在學習之余，可以看出算法背后的意義，為什么會出現這么多算法？

影響非線性求解速度的兩個關鍵點：

• 迭代時所要形成矩陣的方式
• 每次求解矩陣方程的時間

由以上兩個方面可以看出，所有的非線性算法都在聚焦于如何在保證精度的前提下，加速收斂求解。明白了求解目的后，便會對其算法本身有一定的“大局觀”，才不會懼怕枯燥且抽象的公式。