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)
$$
```matlab
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
