Levenberger-Marquardt for nonlinear elliptical system

考虑如下的方程组，测试Levenberger-Marquardt 方法：
\[ \begin{align*} \varphi_{rr}+\frac{2}{r}\varphi_{,r}+\frac{1}{8}(A)^2-\frac{1}{12}\varphi^5 &=f_1\\ A_{,r}-\frac{2}{3}\varphi^6+A &=f_2\\ f_1&=\frac{r^2}{2}-\frac{2}{r (r+1)^2}+\frac{2}{(r+1)^3}-\frac{1}{12 (r+1)^5}\\ f_2&=2 r-\frac{2}{3 (r+1)^6}+2 \end{align*} \]
使用中心差分，以下计算出他们的Jacobi 矩阵：
\[ \begin{align*} a_{i+1}&= F^1_{\varphi_{i+1}} =\frac{1}{h^2}+\frac{1}{h r(i)}\\ a_i &= F^1_{\varphi_i} =-\frac{2}{h^2}-\frac{1}{12} 5 \varphi(i)^4 \\ a_{i-1}&= F^1_{\varphi_{i-1}} = \frac{1}{h^2}-\frac{1}{h r(i)}\\ b_{i} &= F^1_{A_{i}} =\frac{A(i)}{4} \\ c_{i} &= F^2_{\varphi_{i}} = -4 \varphi(i)^5 \\ d_{i+1}&= F^2_{A_{i+1}} =\frac{1}{2h} \\ d_{i} &= F^2_{A_{i}} =1 \\ d_{i-1}&= F^2_{A_{i-1}} =-\frac{1}{2h} \\ \end{align*} \]
\[ J^1_\varphi= \left( \begin{array}{ccc} a_1 & a_2 & 0 \\ a_1 & a_2 & a_3 \\ 0 & a_{n-1} & a_n \\ \end{array} \right) \]
\[ J^1_A= \left( \begin{array}{ccc} b_1 & 0 & 0 \\ 0 & b_2 & 0 \\ 0 & 0 & b_n \\ \end{array} \right) \]
\[ J^2_\varphi= \left( \begin{array}{ccc} c_1 & 0 & 0 \\ 0 & c_2 & 0 \\ 0 & 0 & c_n \\ \end{array} \right) \]
\[ J^2_A= \left( \begin{array}{ccc} d_1 & d_2 & 0 \\ d_1 & d_2 & d_3 \\ 0 & d_{n-1} & d_n \\ \end{array} \right) \]
完整的jacobi 是以上四个构成的分块矩阵
\[ J= \left( \begin{array}{cc} J^1_\varphi & J^1_A \\ J^2_\varphi & J^2_A \end{array} \right) \]

function T70
% 本函数计算简化版的constrain equation using  global  Levenberg-Marquardt method
% yuewen_chen@qq.com
n=200;
r=linspace(1,5,n)';
h=r(2)-r(1);
r0=[r(1),r(end)];  % ghost point
phi0=1./(1+r0);
A0=2*r0;
I11=[1,0;0,0];
I12=[0,1;0,0];
I21=[0,0;1,0];
I22=[0,0;0,1];
f1=power(r,2)/2. - 1./(12.*power(1 + r,5)) + 2./power(1 + r,3) - 2./(r.*power(1 + r,2));
%f1(1)=f1(1)+2/r0(1)*(-phi0(1))/(2*h)+phi0(1)/(h^2);  %不要边界条件，问题就出在这里
%f1(end)=f1(end)+2/r0(2)*(phi0(2))/(2*h)+phi0(2)/(h^2);
f2=2 + 2*r - 2./(3.*power(1 + r,6));
%f2(1)=f2(1)+(-A0(1))/(2*h);
%f2(end)=f2(end)+A0(2)/(2*h);      % Dir Boundary condition 
eta=0.5;
phi=sin(r);  % 初始猜测任意给
A=0*r;
u=[phi;A];
%*********************************************
%*************** 迭代 ************************
%*********************************************
for k=1:500
    if norm(g(u)'*F(u))<1e-15
        fprintf('好了\n')
        break
        
    end
    
    
    lam=norm(F(u));
    L=g(u)'*g(u)+lam*speye(2*n,2*n);
    R=-g(u)'*F(u);
    d=L\R;
     
     
    if norm(F(u+d))<=eta*norm(F(u))
         
        u=u+d;
    else
        alp=Armijo(u,d);
        u=u+alp*d;
        fprintf('alpha=%d\n',alp)
    end
    plot(r,u(1:n),'b.',r,1./(1+r),'r')
    title(['res=',num2str(nF(u)), ',  k=',num2str(k),' , d=',num2str(norm(d))])
    drawnow
end
       
%*********************************************
%*****************nF **************************
  function y=nF(u)
         y=norm(F(u));
    end
%****************** F *************************
function y=F(u)
        ph=u(1:n,1);
        a=u(n+1:end,1);
        ph_p=circshift(ph,[-1,0]);
       ph_m=circshift(ph,[  1,0]);
          a_p=circshift(a,[-1,0]);
         a_m=circshift(a,[  1,0]);
         ph_p(n)=phi0(2);
         ph_m(1)=phi0(1);
         a_p(end)=A0(2);
         a_m(1)=A0(1);
         
         dphi=(ph_p-ph_m)/(2*h);
         d2phi=(ph_p-2*ph+ph_m)/(h^2);
         da=(a_p-a_m)/(2*h);
         
         F1=2./r.*dphi+d2phi+1/8*(a).^2-1/12*ph.^5-f1;
         F2=da-2/3*ph.^6+a-f2;
         y=[F1;F2];
end
%************** Jacobi of F(phi,A) **************
    function y=Jac(u)
        Phi=u(1:n);
        Aa=u(n+1:end);
             
             a_i=-2./(h.^2)-5/12*Phi.^4;
        a_im1=1./(h.^2)-1./(h.*r);
         a_ip1=1./(h.^2)+1./(h.*r);
             b_i=Aa/4;
              c_i=-4*Phi.^5;
             d_i=1*ones(size(r));
        d_im1=-1/(2*h)*ones(size(r));
         d_ip1=1/(2*h)*ones(size(r));
        J1=sparse(diag(a_i,0)+diag(a_im1(2:end),-1)+diag(a_ip1(1:end-1),1));
        J2=sparse(diag(b_i,0));
        J3=sparse(diag(c_i,0));
        J4=sparse(diag(d_i,0)+diag(d_im1(2:end),-1)+diag(d_ip1(1:end-1),1));
        Jacobi=sparse(kron(I11,J1)+kron(I12,J2)+kron(I21,J3)+kron(I22,J4));
        y=Jacobi;
    end

%***************Armijio***********************
function alp=Armijo(u,d)
        et=1/2;
        bet1=0.8;
        t=0;
        
       
         while t<=2000 
             
             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
%****************** Jaco **********************
    function y=g(u)
        y=Jac(u);
    end
%**********************************************
end

转载于:https://www.cnblogs.com/yuewen-chen/p/11525996.html