本次分享的是基于Modified Newton–Raphson的非線性有限元編程小案例，將會(huì)介紹牛頓法的優(yōu)缺點(diǎn)，以及如何一步步修正，最后結(jié)合彈簧單元的例子進(jìn)行Modified Newton–Raphson數(shù)值實(shí)現(xiàn)，望大家喜歡。

牛頓法缺點(diǎn)

考慮使用N-R法求解以下非線性方程：

該方程的精確解為 ，設(shè)置收斂容差 ，初始值 。

每次迭代的切線：

牛頓法求解代碼

xdata=zeros(40,1);
ydata=zeros(40,1);
tol = 1.0e-5;
iter = 0;
u = 0.5;
uold = u;
c=0;
P = u+atan(5*u);
R = -P;
conv= R^2;
xdata(1)=u;
ydata(1)=P;
while conv > tol && iter < 20
Kt = 1+5*(cos(atan(5*u)))^2;
delu = R/Kt;
u = uold + delu;
P = u+atan(5*u);
R = -P;
conv= R^2;
uold = u;
iter = iter + 1;
xdata(2*iter)=u; ydata(2*iter)=0;
xdata(2*iter+1)=u; ydata(2*iter+1)=P;
end
%
plot(xdata,ydata);
hold on;
x=[-1:0.1:1];
y=x+atan(5*x);
plot(x,y)

以上程序在迭代次數(shù)超過預(yù)置的20次后停止迭代，每次迭代的結(jié)果皆偏離精確解，導(dǎo)致數(shù)值發(fā)散。

不收斂原因：初始取值不合理

解決方法：初值小于0.5，程序?qū)?huì)很快收斂；大于0.5，程序發(fā)散。

由此可以得出Newton–Raphson method的優(yōu)缺點(diǎn)：

1. 收斂快
2. 每次都要更新切線剛度，計(jì)算成本高
3. 初始值在接近精確值時(shí)。才能保證收斂

修正牛頓法

為克服牛頓法在求解時(shí)的不足，可以將牛頓法的求解思路稍加改進(jìn)，本節(jié)將簡要回顧一下Modi?ed Newton–Raphson method的思想。相比于Newton–Raphson method，Modified的區(qū)別：

1. 每次迭代只需 the initial  tangent stiffness matrix
2. 迭代次數(shù)多，每次迭代時(shí)間長，但計(jì)算成本相對(duì)較低

為了提高收斂性，可先迭代一些tangent stiffness matrix，然后進(jìn)行Modi?ed Newton–Raphson method。但該方法在構(gòu)建合理的tangent stiffness matrix，很難指定確定迭代次數(shù)。

非線性應(yīng)用

Example

如下圖所示，模型由兩個(gè)彈簧單元組成，共3個(gè)節(jié)點(diǎn)，每個(gè)節(jié)點(diǎn)包含1個(gè)自由度，共3個(gè)自由度，模型右端受到 的集中載荷，左端固定，彈簧剛度與單元位移量相關(guān)，其中， ， ，試求 和 的值。假設(shè)單位統(tǒng)一，不計(jì)單位。

每個(gè)彈簧元的平衡方程如下式：

整體剛度組裝，去除0左端節(jié)點(diǎn)項(xiàng)（固支條件，劃行劃列）后的整體剛度方程可表示如下：

進(jìn)一步講上式化簡為1個(gè)非線性方程組 ：

Solution

tol = 1.0e-5;
iter = 0;
u = [0.2; 0.4];
uold = u;
c=0;
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=f-P;
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);
Kt = [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 < 100
delu = Kt\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=f-P;
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;
uold = u;
iter = iter + 1;
fprintf('\n %3d %7.5f %7.5f %12.3e %7.5f',iter,u(1),u(2),conv,c);
end

由以上求解過程，可以看出，Modi?ed Newton–Raphson method的迭代次數(shù)相對(duì)于Newton–Raphson method的迭代次數(shù)較多，問題規(guī)模較為簡單，體現(xiàn)不出求解時(shí)間的差別，讀者可自行驗(yàn)證。

【注】:在本程序中，設(shè)置的初值是u = [0.2; 0.4]，讀者可設(shè)置不同的初值，探究初值對(duì)收斂次數(shù)的影響。在設(shè)置初值時(shí)，應(yīng)保持 (常識(shí)判斷)，若反向設(shè)置，將導(dǎo)致程序錯(cuò)誤。