Note for constraint equation

**Note for constraint equation **

We consoder a function

\[ f_x=\frac{1}{\varepsilon^2 +(x-x_0)^2} \]

Integrate of this function is

\[ f=-\frac{1}{\varepsilon}arctan( \frac{x_0-x}{\varepsilon}) \]

A good idear for tcest Jang's equation is : choice a very small number $\delta$

\[\begin{align} \delta f &= -\frac{\delta}{\varepsilon}arctan( \frac{x_0-x}{\varepsilon}) \\ \delta \varepsilon f &=-\delta arctran( \frac{x_0-x}{\varepsilon} ) \end{align} \]

接着要讨边界条件的设置问题：对于内侧的 $\displaystyle \frac{f_x }{ \sqrt{ 1+f_x ^2}}$ 的边界，使用的是 Newmann 边界，给出的 $\displaystyle \phi^{-2} f_x=\frac{1}{\varepsilon^2 +(x-x_0)^2} $, 所以Newmann 边界条件是这样给出来的 $$\boxed{ f_x= \frac{1}{\varepsilon^2 +(x-x_0)^2} \phi^2} $$ 对于外侧的，使用的 Dirichlet 边界条件， $$ \begin{align} \kappa(r)&=\frac{\frac{f_r}{\phi ^2} }{\sqrt{\frac{f_r^2}{\phi ^4}+1}}\\ H(r)&=-\frac{1}{r^2 \phi^6 }( r^2 \phi^4 \frac{\frac{f_r}{\phi ^2} }{\sqrt{\frac{f_r^2}{\phi ^4}+1}})_{,r}\\ &= -\frac{1}{r^2 \phi^6 }( r^2\phi^4\kappa(r))_{,r} \\ \end{align} $$ 实际上使用的是 $$ \boxed{ r^2\phi^4(r) \kappa(r)} $$ 的边界值函数。 ```matlab function T74 % This function we use Levenberger method to calculate Einstein constraint % equation which had litte change than T73.m n=600; eps=1/20; de=1/100; r0=20; r=linspace(10,400,n)'; h=r(2)-r(1); I11=[1,0;0,0]; I12=[0,1;0,0]; I21=[0,0;1,0]; I22=[0,0;0,1]; eta=0.5; phi_0=1+1./r; % Initial guess X_0=1./r; u=[phi_0;X_0]; [kappa,dkaapa,ddkappa,J,dJ]= Kappa_and_J; %================================== %=========== main part ============ %================================== 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
subplot(2,2,1)
plot(r,u(1:n).^4,'b-.')
title(['\phi ' ,'res=',num2str(nF(u)), ',  k=',num2str(k),' , d=',num2str(norm(d))])
subplot(2,2,2)
plot(r,u(n+1:end),'-.')
title('X')
%subplot(1,3,3)
%plot(r,u(1:n).^4,'b-.')
%xlim([10,50])
drawnow

end
[K,KL]=Initial_Data(u);
subplot(2,2,3)
plot(r,KL)
title('KL')
subplot(2,2,4)
plot(r,K)
title('K')
%nF **************************
function y=nF(u)
y=norm(F(u));
end
%---------------------- F--------------------------------
%* F *************************
function y=F(u)
phi=u(1:n);
X=u(n+1:2n);
phip1=circshift(phi,[-1,0]);
phim1=circshift(phi,[1,0 ]);
Xp1=circshift(X,[-1,0]);
Xm1=circshift(X,[1,0 ]);
% Boundary condition
phip1(end)=1;
phim1(1)=4;
Xp1(end)=0;
Xm1(1)=1;
F1=(36(-2phi + phim1 + phip1) + (36h.(-phim1 + phip1))./r -...
(3power(6kappa.power(phi,3).(h.phi - phim1.r + phip1.r) ...
+ r.(3dkaapa.h.power(phi,4) + 2J.(Xm1 - Xp1)),2))./...
(power(-3 + J,2).power(phi,7).power(r,2)) + (2.power(Xm1 - Xp1,2))./...
power(phi,7))./(36.power(h,2));

    F2=(-3*h*power(-3 + J,2).*power(phi,7).*r.*(Xm1 - Xp1) + 4*power(-3 + J,2).*...
           power(phi,7).*power(r,2).*(-2*X + Xm1 + Xp1)+2*(dJ.*h.*phi.*r.*(6*kappa.*...
           power(phi,3).*(h*phi - phim1.*r + phip1.*r) + r.*(3*dkaapa*h.*power(phi,4) + 2*J.*(Xm1 - Xp1))) +... 
           3*(3 - J).*(phim1 - phip1).*r.*(6*kappa.*power(phi,3).*(h*phi - phim1.*r + phip1.*r)...
           + r.*(3*dkaapa.*h.*power(phi,4) + 2*J.*(Xm1 - Xp1)))-(3 - J).*phi.*(-3*kappa.*...
           power(phi,2).*(3*power(phim1,2).*power(r,2) + 3*power(phip1,2).*power(r,2)...
           + 4*phi.*phip1.*r.*(h + r) + 2*phim1.*r.*(-2*phi.*(h - r) - 3.*phip1.*r)-2*power(phi,2).*...
            (power(h,2) + 4*power(r,2))) + r.*(4*J.*r.*(-2*X + Xm1 + Xp1) - h*(3.*power(phi,3).*...
            (-4*dkaapa.*phim1.*r + 4*dkaapa.*phip1.*r + h*phi.*(2*dkaapa + ddkappa.*r)) ...
            + 2*dJ.*r.*Xm1 - 2*dJ.*r.*Xp1)))))./(3.*power(h,2).*power(-3 + J,2).*power(phi,7).*power(r,2));   
   y=[F1;F2];

end

%--------------------------------------------------------
%Armijio***********************
function alp=Armijo(u,d)
et=1/2;
bet1=0.8;
t=0;

     while t<=200
         
         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

%---------------------- jacobi -------------------------

function y=g(u)
    y=Jac(u);
end

%*********************************************

%************** Jacobi of F(phi,X) **************
function y=Jac(u)

     phi=u(1:n);
        X=u(n+1:end);
 phip1=circshift(phi,[-1,0]);
phim1=circshift(phi,[1,0  ]);
    Xp1=circshift(X,   [-1,0]);
   Xm1=circshift(X,   [1,0  ]);

a1=(1 - h./r + (kappa.(6kappa.power(phi,3).(hphi - phim1.r + phip1.r)...
+ r.(3dkaapah.power(phi,4) + 2J.(Xm1 - Xp1))))./(power(-3 + J,2).power(phi,4).r))./power(h,2);
a2=-(72 + (36(2dkaapa.h.phi.r + kappa.(4h.phi - 3phim1.r ...
+ 3phip1.r)).(6kappa.power(phi,3).(h.phi - phim1.r + phip1.r)...
+ r.(3dkaapa.h.power(phi,4) + 2J.(Xm1 - Xp1))))./...
(power(-3 + J,2).power(phi,5).power(r,2)) - (21power(6kappa.power(phi,3).(hphi ...
- phim1.r + phip1.r) + r.(3dkaapah.power(phi,4) + 2J.(Xm1 - Xp1)),2))./...
(power(-3 + J,2).power(phi,8).power(r,2)) + (14power(Xm1...
- Xp1,2))./power(phi,8))./(36.power(h,2));
a3=(1 + h./r - (kappa.(6kappa.power(phi,3).(hphi - phim1.r + phip1.r) + r.(3dkaapah.power(phi,4)...
+ 2J.(Xm1 - Xp1))))./(power(-3 + J,2).power(phi,4).r))./power(h,2);
b1=(Xm1 - (3J.(6kappa.power(phi,3).(hphi - phim1.r + phip1.r) + r.(3dkaapa.h.power(phi,4) +...
2J.(Xm1 - Xp1))))./(power(-3 + J,2).r) - Xp1)./(9.power(h,2).power(phi,7));
b2=0;
b3=(-Xm1 + (3J.(6kappa.power(phi,3).(hphi - phim1.r + phip1.r) + r.(3dkaapa.h.power(phi,4) +...
2J.(Xm1 - Xp1))))./(power(-3 + J,2).r) + Xp1)./(9.power(h,2).power(phi,7));
c1=(2(-2kappa.power(phi,3).(-3(-3 + J).phim1.r + 3(-3 + J).phip1.r +...
phi.(-3h - 6r + dJ.h.r + J.(h + 2r))) - (-3 + J).r.(-(dkaapa.h.power(phi,4)) ...
+ 2J.(Xm1 - Xp1))))./(power(h,2).power(-3 + J,2).power(phi,7).r);
c2=(2(-2kappa.power(phi,3).(6(-3 + J).power(phim1,2).power(r,2) + 6(-3 + J).power(phip1,2).power(r,2)...
+ 3phi.phip1.r.(-3h + J.(h - 2r) + 6r + dJ.h.r)+2power(phi,2).(dJ.power(h,2).r - 3.(power(h,2) +...
4power(r,2)) + J.(power(h,2) + 4power(r,2))) - 3phim1.r.(4(-3 + J).phip1.r + phi.(-3h - 6r +...
dJh.r + J.(h + 2r))))+r.(-(hphi.(power(phi,3).(-9dkaapa.phim1.r + 9dkaapa.phip1.r +...
2h.phi.(3ddkappa.r + dkaapa.(6 + dJ.r))) + 12dJ.r.Xm1 - 12dJ.r.Xp1)) +...
2power(J,2).r.(7phim1.(Xm1 - Xp1) + 7phip1.(-Xm1 + Xp1) - 4phi.(-2X + Xm1 + Xp1)) +...
J.(3dkaapa.h.power(phi,4).phip1.r + 2power(h,2).power(phi,5).(2dkaapa + ddkappa.r) -...
3phim1.r.(dkaapah.power(phi,4) + 14Xm1 - 14Xp1) + 42phip1.r.(Xm1 - Xp1) +...
24phi.r.(-2X + Xm1 + Xp1)))))./(power(h,2).power(-3 + J,2).power(phi,8).power(r,2));
c3=(-2(2kappa.power(phi,3).(3(-3 + J).phim1.r - 3(-3 + J).phip1.r - phi.(-3h + J.(h - 2r) ...
+ 6r + dJ.h.r)) - (-3 + J).r.(-(dkaapa.h.power(phi,4)) + 2J.(Xm1 - Xp1))))./...
(power(h,2).power(-3 + J,2).power(phi,7).r);
d1=(-3hpower(-3 + J,2).power(phi,7) + 4power(-3 + J,2).power(phi,7).r + 4(dJh.J.phi ...
+ (-3 + J).(-(dJ.h) + 2J).phi - 3(-3 + J).J.(phim1 - phip1)).r)./...
(3.power(h,2).power(-3 + J,2).power(phi,7).r);
d2=(-8(-3power(phi,6) + J.(2 + power(phi,6))))./(3.power(h,2).(-3 + J).power(phi,6));
d3=(3hpower(-3 + J,2).power(phi,7) + 4power(-3 + J,2).power(phi,7).r + 4(-(dJ.h.J.phi) +...
(-3 + J).(dJh + 2J).phi + 3(-3 + J).J.(phim1 - phip1)).r)./...
(3.power(h,2).power(-3 + J,2).power(phi,7).r);
J1=sparse(diag(a2,0)+diag(a1(2:end),-1)+diag(a3(1:end-1),1));
J2=sparse(diag(b1(2:end),-1)+diag(b3(1:end-1),1));
J3=sparse(diag(c2,0)+diag(c1(2:end),-1)+diag(c3(1:end-1),1));
J4=sparse(diag(d2,0)+diag(d1(2:end),-1)+diag(d3(1:end-1),1));
Jacobi=sparse(kron(I11,J1)+kron(I12,J2)+kron(I21,J3)+kron(I22,J4));
y=Jacobi;
end

%---------------- function for kappa and J ------------------------------
function[kappa,dkappa,ddkappa,J,dJ]= Kappa_and_J
c=@(r)de./(eps^2 + (r - r0).^2);
dc=@(r)de(-2(r - r0))./power(power(eps,2) + power(r - r0,2),2);
ddc=@(r)(-2de)./power(power(eps,2) + power(r - r0,2),2) + ...
(8depower(r - r0,2))./power(power(eps,2) + power(r - r0,2),3);
ddka=@(r)(3power(c(r),3).power(dc(r),2))./power(1 + power(c(r),2),2.5) ...
- (3c(r).power(dc(r),2))./power(1 + power(c(r),2),1.5) - (power(c(r),2).ddc(r))./power(1 + power(c(r),2),1.5) +...
ddc(r)./sqrt(1 + power(c(r),2)); %<------ kappa''(x)
dka=@(r)-((power(c(r),2).dc(r))./power(1 + power(c(r),2),1.5)) + dc(r)./sqrt(1 + power(c(r),2)); %<------kappa'(x)
ka=@(r)c(r)./sqrt(1 + c(r).^2);
kappa=ka(r);
dkappa=dka(r);
ddkappa=ddka(r);
J=kappa.^2;
dJ=(-4power(de,2).(power(eps,2) + power(r - r0,2)).(r - r0))./power(power(de,2) ...
+ power(power(eps,2) + power(r - r0,2),2),2);
%---------------------------------
end

function [K,KL]=Initial_Data(u)
  phi=u(1:n);
  X=u(n+1:2*n);
  phip1=circshift(phi,[-1,0]);
  phim1=circshift(phi,[1,0 ]);
  Xp1=circshift(X,[-1,0]);
  Xm1=circshift(X,[1,0 ]);
  % Boundary condition         
   phip1(end)=1;
   phim1(1)=4;
   Xp1(end)=0;
   Xm1(1)=1;
  %-------------------------------- 
  KL=(2*kappa.*power(phi,3).*(h*phi - phim1.*r + phip1.*r)...
      + r.*(dkaapa*h.*power(phi,4) + 2*Xm1 - 2*Xp1))./(h*(-3 +J).*power(phi,6).*r);
  
  K=(6*kappa.*power(phi,3).*(h*phi - phim1.*r + phip1.*r)...
     + r.*(3*dkaapa.*h.*power(phi,4) + 2*J.*(Xm1 - Xp1)))./(h*(-3 + J).*power(phi,6).*r);
end

end

posted @ 2019-09-28 12:07 yuewen_chen 阅读(213) 评论(0) 编辑收藏举报

刷新页面返回顶部