弧长法求解非线性方程

说到非线性方程（组）的解法，牛顿-拉夫逊方法（Newton-Raphson method）应该是大家最为熟悉的方法。但在某些的力学问题中，牛顿法却是无能为力。面对这种困境，Riks提出了弧长法（Arc length method）来求解非线性方程组。下面在力学背景下介绍弧长法的基本思想。

载荷-位移方程

对于一个结构，常用（广义）载荷\(~F~\)与（广义）位移\(~u~\)之间的关系来描述其变形行为。考虑真实的加载过程，有两种控制形式：

载荷控制，即已知载荷位移关系为\(~F=F(u)\)，根据载荷增量\(~\Delta F\)，求位移增量\(~\Delta u\)；
位移控制，即已知载荷位移关系为\(~u=u(F)\)，根据位移增量\(~\Delta u\)，求载荷增量\(~\Delta F\)。

事实上，\(F(u)~\)与\(~u(F)~\)两者互为反函数。但对于某些问题，\(F(u)~\)或\(~u(F)~\)无法写为初等函数的形式

本文以载荷控制为例，研究\(~F(u)~\)为非线性方程时，如何根据载荷增量\(~\Delta F\)，求得位移增量\(~\Delta u~\)。

位移控制类似，这里不再赘述

非线性方程的解法

为简化方程形式，可将\(~F(u)~\)归一化为无量纲形式，记为\(~\lambda(a)\)，其中\(~\lambda~\)对应载荷\(~F\)，\(a~\)对应位移\(~u~\)。

牛顿法

牛顿法的迭代示意图如下，已知载荷增量\(~\Delta \lambda\)，求位移增量\(~\Delta a\)。

图片来自于Nick的文章Nonlinear Analysis of Structures

如下图所示，当载荷增量\(~\Delta \lambda~\)越过局部最高点时，牛顿法无法捕捉到局部最高点以下的载荷路径。这种现象被称为跳跃（Snap through），除此之外还有跳回（Snap back）现象。

事实上，牛顿法的本质决定其载荷增量只能沿一个方向，不允许反向加载。因此，牛顿法无法求解类似于加载-卸载-加载的过程。

弧长法

弧长法的迭代示意图如下，已知弧长\(~\Delta l\)，求载荷增量\(~\Delta \lambda~\)和位移增量\(~\Delta a\)。

基本思想：

以初始点\(~(a_i,\lambda_i)~\)为圆心、弧长\(~\Delta l~\)为半径作圆；
过初始点\(~(a_i,\lambda_i)~\)做切线；切线与圆相交，记交点的横坐标为\(~a'\)，并计算得到\(~\lambda '=\lambda(a')\)；
若收敛，则结束；若不收敛，将\(~(a',\lambda ')~\)作为新的初始点，返回步骤2。

圆方程是在上述过程中是不变的直线方程与圆方程联立后有两个解，需要进行判断

将初始点记为\(~(a_0,\lambda_0)\)，第\(~i~\)次迭代的点记为\(~(a_i,\lambda_i)\)，下面给出圆与直线的联立方程组 \[ \left\{ \begin{array}{ccl} \Delta l^2 &=& (\lambda-\lambda_0)^2+(a-a_0)^2 \\[3pt] \lambda &=& k_i (a-a_i)+\lambda_i \end{array} \right. \]

式中，\(k_i=\left.\displaystyle\frac{\mathrm{d} \lambda}{\mathrm{d} a}\right|_{a=a_i}\)。将上式整理可得 \[ (1+k^2_i)a^2+2\big[k_i(\lambda_i-\lambda_0-k_ia_i)-a_0\big] a + \big[(\lambda_i-\lambda_0-k_ia_i)^2+a^2_0-\Delta l^2\big]=0 \]

简单判断方法：取沿加载方向的解

算例分析

模型

如图所示，由两根原长为\(~L_0\)、刚度系数为\(~k~\)的杆组成一自由度的桁架结构，初始夹角为\(~\theta_0\)，载荷\(~P~\)作用在铰点上，载荷\(~P~\)的垂直位移记为\(~u\)，材料处于线弹性状态，且不考虑结构失稳。

不难得到载荷\(~P~\)与位移\(~u~\)的关系，即 \[ \frac{P}{2kL_0}=\Bigg( \frac{1}{\sqrt{1-2\frac{u}{L_0}\sin\theta_0+\big(\frac{u}{L_0}\big)^2}} -1 \Bigg) \Big(\sin\theta_0 -\frac{u}{L_0} \Big) \]

Q：上式反映了载荷\(~P~\)与位移\(~u~\)的非线性关系，可材料为线弹性材料？

A：此为几何非线性（非线性分三种：材料、几何、边界条件），材力教材例题13-3与该模型相似

令\(~\lambda=P\big/ 2kL_0\)、\(a=u\big/ L_0\)，归一化可得 \[ \lambda(a)=\Big( \frac{1}{\sqrt{1-2a\sin\theta_0+a^2}} -1 \Big) \Big(\sin\theta_0 -a \Big) \]

该桁架结构存在三个载荷\(~P=0~\)的平衡位置，分别对应\(~u=0,~L_0\sin\theta_0,~2L_0\sin\theta_0\)（\(~a=0,~\sin\theta_0,~2\sin\theta_0\)），即下图所示

其中\(~a=0,~2\sin\theta_0~\)对应稳定平衡位置，\(a=\sin\theta_0~\)对应不稳定平衡位置

结果

为了对比牛顿法和弧长法的求解效果，编写了MATLAB计算程序，令初始夹角\(~\theta_0=60^\circ\)，结果如下所示

显然对于该问题，牛顿法无法准确跟踪加载历程。

程序

Nick也编写了Python计算程序：牛顿法、弧长法

载荷\(~\lambda~\)与位移\(~a~\)的关系：

% out = [y,dy/dx]
function [y,dy]=fun(x)

    st=sin(pi/3);
    temp=1-2*st*x+x^2;
    y=(1/sqrt(temp)-1)*(st-x);
    dy=(st-x)^2/temp^(3/2)-(1/sqrt(temp)-1);
    
end

主函数：

%version - R2018b
%%
clc;clear;close all;

pt=12;%字体大小
tol=1e-10;%误差限
num=35;%点的个数
time=3;%绘图时间
dy=0.8/num;%载荷增量
x_newton=zeros(num,1);y_newton=x_newton;
dl=0.1;%弧长
x_arc=zeros(num,1);y_arc=x_arc;


x=0:0.01:3;
[y,~]=arrayfun(@fun,x);
h=figure;figset([20,7,0.3,0.3]);
subplot(1,2,1);plot(x,y,'k','linewidth',1.5);hold on;grid on;
set(gca,'ticklabelinterpreter','latex','fontsize',pt);
xlabel('$a$','interpreter','latex','fontsize',pt);
ylabel('$\lambda$','interpreter','latex','fontsize',pt,'Rotation',0);
title('Newton''s Method','interpreter','latex','fontsize',pt);

subplot(1,2,2);plot(x,y,'k','linewidth',1.5);hold on;grid on;
set(gca,'ticklabelinterpreter','latex','fontsize',pt);
xlabel('$a$','interpreter','latex','fontsize',pt);
ylabel('$\lambda$','interpreter','latex','fontsize',pt,'Rotation',0);
title('Arc Length Method','interpreter','latex','fontsize',pt);

for i=2:num
    [x_newton(i),y_newton(i)]=newtons_method(x_newton(i-1),dy,tol);
    subplot(1,2,1);
    plot(x_newton(i),y_newton(i),'ob','linewidth',1.5,'markersize',7);
    
    [x_arc(i),y_arc(i)]=arc_method(x_arc(i-1),dl,tol);
    subplot(1,2,2);
    plot(x_arc(i),y_arc(i),'ob','linewidth',1.5,'markersize',7);
    pause(time/num);
end

%%
%设置图片尺寸与位置
function figset(parameter1,parameter2)

    %电脑屏幕尺寸
    set(0,'units','centimeters')
    cm_ss=get(0,'screensize');
    W=cm_ss(3);%电脑屏幕长度，单位cm
    H=cm_ss(4);%电脑屏幕宽度，单位cm

    %设置figure在screen中的比例位置与大小
    temp1=[parameter1(3),parameter1(4),parameter1(1)/W,parameter1(2)/H];
    set(gcf,'units','normalized','position',temp1);
    if nargin==2
        %设置axis在figure中的比例位置与大小
        temp2=[parameter2(3),parameter2(4),parameter2(1),parameter2(2)];
        set(gca,'position',temp2);
    end

end

弧长法：

% 弧长法
% input: 初始点x0, 弧长dl, 误差限tol
%output: 结束点[x1,y1]
function [x1,y1]=arc_method(x0,dl,tol)

    %tol=1e-10;%误差限
    [y0,dy0]=fun(x0);%初始点
    xi=x0;yi=y0;dyi=dy0;
    %dl=0.01;%弧长

    Nmax=1000;num=1;
    while 1
        temp_t=yi-dyi*xi-y0;
        temp_a=1+dyi^2;
        temp_b=2*dyi*temp_t-2*x0;
        temp_c=temp_t^2+x0^2-dl^2;
        delta=temp_b^2-4*temp_a*temp_c;
        if delta<0
            disp('delta<0');
            break;
        end
        xj1=(-temp_b+sqrt(delta))/2/temp_a;
        xj2=(-temp_b-sqrt(delta))/2/temp_a;
        xj=max(xj1,xj2);%取增大项
        [yj,dyj]=fun(xj);
        num=num+1;
        if (abs(yj-yi)<tol) || (num==Nmax)
            %disp(['err=',num2str(abs(yj-yi))]);
            %disp(['radius=',num2str(sqrt((xj-x0)^2+(yj-y0)^2))]);
            break;
        end
        xi=xj;yi=yj;dyi=dyj;
    end
    x1=xi;y1=yi;
    
end

牛顿法：

% 牛顿法
% input: 初始点x0, 载荷增量dy, 误差限tol
%output: 结束点[x1,y1]
function [x1,y1]=newtons_method(x0,dy,tol)
    
    xi=x0;
    [yi,dyi]=fun(xi);%初始点
    %dy=0.01;%载荷增量
    y1=yi+dy;%结束点

    Nmax=1000;num=1;
    while 1
        xj=xi+(y1-yi)/dyi;
        [yj,dyj]=fun(xj);
        xi=xj;yi=yj;dyi=dyj;
        num=num+1;
        if (abs(yi-y1)<tol) || (num==Nmax)
            break;
        end
    end
    x1=xi;y1=yi;

end