Levenberg-Marquardt method for nonlinear elliptical equation

使用 Levenberg-Marquardt 方法测试 $$u_{xx}+u^2=-sin(x)+sin(x)^2$$, 初值选取 $u(x)=cos(x)sin(x)$ 参考文献：《非线性方程组数值方法》，袁亚湘； 1. 给出初值 $x_0 \in \mathbb{R}^n; k=1,\eta \in(0,1)$ 2. 若 $\|J_k^{T} F_k\|=0$, 停; 求解 $\displaystyle (J_k^TJ_k+\lambda_k I)d=-J_k^T F_k,\lambda_k=\|F_k\|$, 得到搜索方向 $d_k$； 3. 若$d_k$ 满足 $\displaystyle \|F(x_k+d_k) \| \leq \eta \|F_k\|$, 则 $x_{k+1}=x_k+d_k$; 否则： $x_{k+1}=x_k+\alpha_k d_k$, 这里的 $\alpha_k $ 由Armijio 线搜索得到； 注：这里的Armijio 线搜索如下： $\alpha_k=\xi^t,\xi \in(0,1)$, $t\in \mathbb{Z}^+$ 是满足如下不等式的最小非负整数： $$ {\| F(x_k+\xi^t d_k) \|}^2 \leqslant {\|F(x_k) \|}^2+\beta_1 \xi^t F^T_k J_k d_k $$ 4. $k=k+1$, 转2. ```matlab function T69 % 全局二次收敛的 Levenberg-Marquardt % yuewen_chen@qq.com

n=300;
r=linspace(-pi5,pi5,n)';
h=r(2)-r(1);
rb=[r(1)-h;r;r(end)+h];
f=-sin(r)+sin(r).^2;
x0=sin(rb);
u=cos(r).*sin(r); % initial guess

% ************ D2********************
cn1=-1ones(n-1,1);
cp1=1ones(n-1,1);
c=-2ones(n,1);
D2=sparse((diag(c,0)+diag(cp1,1)+diag(-cn1,-1))/(h^2));
%*****************************************
M=D2;
eta=0.9;

for k=1:200
lam=norm(F(u));
L=g(u)'g(u)+lamspeye(n,n);
R=-g(u)'F(u); %之前转置写掉了；
d=L\R;
if norm(F(u+d))<=norm(etaF(u))
u=u+d;
else
alp=Armijo(u,d);
u=u+alp*d;
end
plot(r,u,'b.',r,sin(r),'r')
title(['res=',num2str(nF(u)), ', k=',num2str(k),' , \alpha=',num2str(norm(d))])
drawnow
end

function y=F(u)
y=M*u+u.^2-f;
y(1)=y(1)+x0(1)/h^2;
y(end)=y(end)+x0(end)/h^2;
end

function y=g(u)
J=M+diag(2*u,0);
y=J;
end

function y=nF(u)
     y=norm(F(u));
end
function alp=Armijo(u,d)
    et=1/5;
    bet1=0.5;
    t=0;
     while t<=500
         
         if (nF(u+et^t*d)^2<=nF(u)^2+bet1*et^t*F(u)'*g(u)*d)
              alp=et^t;
              break
         else
             t=t+1;
         end
         alp=et^t;
     end
end

end