量子三体问题: 高精度数值计算、数值积分

目的

三体问题计算中需要高精度积分计算，高精度特征值求解，因此汇总一些现有的高精度数值计算工具和软件。并做些实际计算问题

网络资源

• 任意精度计算工具库 .最早接触的就是David H. Bailey的高精度计算库，2000年左右就很稳定，此牛对浮点计算研究功力很深，当年GPU仅支持float浮点类型，依据此人方法，用float-float方式模拟了GPU上double计算。
• DE积分方法，也简称IMT积分方法，是日本学者1974发现的，和Tanh-Sinh积分有异曲同工之妙，对有限区间上平滑函数(在端点处可以有可积奇性)精度很高，matlab程序提供$[0,1]$区间4096点DE积分的位置和权重

计算实例

• 计算积分: $$Q={\int }_0^1\exp (-\frac{1}{t}-\frac{1}{1-t})dt=0.007029858406609656239241270530354$$

     %matlab函数vpaintegral计算任意精度积分结果
     syms t x y;
     y=exp(-1/t-1/(1-t);
     %32 digits有效数字
     Q=vpaintegral(y,t,0,1,'RelTol',1e-32,'AbsTol',0);
  

  使用matlab的vpaintegral函数，可以指定计算精度有效位数，vpaintegral不受digits函数影响
• 计算DE积分公式的积分点和权重(Double Exponetial transformation). 记 $\phi (\tau )=\frac{1}{Q}{\int }_0^{\tau }\exp (-\frac{1}{t}-\frac{1}{1-t})dt$, 有积分变换公式: $${\int }_0^{1}f(t)dt={\int }_0^{1}f(\phi (\tau )){\phi }^{\prime }(\tau )d\tau ={\int }_0^{1}g(\tau )d\tau ,g(\tau )=f(\phi (\tau )){\phi }^{\prime }(\tau )$$
  $g(x)$在积分端点$0,1$处无限可微，且任意阶导数为0，采用Trapzoidal积分公式可获高精度积分结果(依据Euler-Maclaurin公式,因 $g(x)$在积分端点处的超级平滑性质，Trapzoidal积分公式中的误差项都会消失)。现计算4096个积分点Trapzoidal积分中的位置和权重:

     syms t x y;
     y=exp(-1/t-1/(1-t));
     %32 digits有效数字
     Q=vpaintegral(y,t,0,1,'RelTol',1e-32,'AbsTol',0);  %使用RelTol参数控制有效计算位数
     
     bias=[]; % 积分位置
     ws=[]; %权重
     digits(32); %设置有效位数，注意vpaintegral不受digits函数影响
     
     for k=1:4095
     		t=vpa(k/4096);
     		w=exp(-1/t-1/(1-t))/Q;
     		ws=[ws;w];
     		y=exp(-1/x-1/(1-x))/Q;
     		b=vpaintegral(y,0,t,'RelTol',1e-32, 'AbsTol',0);  %使用RelTol参数控制有效计算位数
     		bias=[bias;b];
     end 
  

  依据4096个点的积分位置和权重，可依次得到2048, 1024, 512等积分点位置和对应权重！
• 计算函数：$$U_m(\xi )={\int }_0^{\pi }\frac{\cos (m\beta )}{\sqrt{1+\xi \cos \beta }}d\beta ,m=0,1,...$$ 此函数出现在三体问题库伦势的矩阵元计算中。其中$\xi$ 位于$[0,1]$区间，$\xi$取值对应DE积分公式的积分点位置。给定 { ξ , m } \{\xi, m\} {ξ,m}, 利用vpaintegral计算$U_m(\xi )$, 根据 m = { 4 k , 4 k + 1 , 4 k + 2 , 4 k + 3 } m=\{4k,4k+1,4k+2,4k+3\} m={4k,4k+1,4k+2,4k+3}, 可以把积分区间缩小到 { 0 , π / 2 } \{0,\pi/2\} {0,π/2},进而加速计算。

---

积分精度测试

计算测试积分: ${\int }_0^1\frac{cos(10x)}{\sqrt{x}}dx$

  !-----------------------------------------------------
  !matlab,vpaintegral结果
  !   y=cos(10*x)/sqrt(x)
  !   vpaintegral(y,0,1,'RelTol',1e-48,'AbsTol,0)
  !   0.346366232384436488606080435492338630630888057352
  !-----------------------------------------------------

Fortran测试程序

!----------------------------------------
!
!4096个DE(IMT)积分点坐标和对应权重:
!   坐标文件: imt4096_bias.txt
!   权重文件: imt4096_weight.txt
!---------------------------------------

program main

  integer,parameter::wp=kind(1.0_16)
  
  integer,parameter:: NIMT=4096

  !4096-Point DE Quadrature
  real(wp),dimension(1:NIMT-1)::IMT4096_Bias,IMT4096_Weight
  !2048-Point DE Quadrature
  real(wp),dimension(1:NIMT/2-1)::IMT2048_Bias,IMT2048_Weight
  !1024-Point DE Quadrature
  real(wp),dimension(1:NIMT/4-1)::IMT1024_Bias,IMT1024_Weight
  !512-Point DE Quadrature
  real(wp),dimension(1:NIMT/8-1)::IMT512_Bias,IMT512_Weight

  !real(wp),dimension(511)::y,t
  real(wp), dimension(1023)::y,t

  call read_imt_file

 
  !t=IMT512_Bias;
  !y=cos(10*t)/sqrt(t)
  !write(*,*) sum(y*IMT512_Weight)/512  !0.346366232384436488606081234949409695 

  !积分变换后这个更准确
  !y=20*sqrt(t)*sin(10*t)
  !write(*,*) 2*cos(10.0_wp)+sum(y*IMT512_Weight)/512 !0.346366232384436488606080435492338269
  
  t=IMT1024_Bias;
  y=cos(10*t)/sqrt(t)
  write(*,*) sum(y*IMT1024_Weight)/1024 !0.346366232384436488606080435492339328
 
 stop

contains
  subroutine read_imt_file()
    implicit none 
    
    !读取bias
    open(1,file='imt4096_bias.txt',status='old')
    read(1,*) IMT4096_Bias(1:NIMT/2)
    close(1)

    !读取权重
    open(1,file='imt4096_weight.txt',status='old')
    read(1,*) IMT4096_Weight(1:NIMT/2)
    close(1)

    !利用对称性
    forall(i=1:NIMT/2-1)
       IMT4096_Bias(NIMT/2+i)=1.0_wp-IMT4096_Bias(NIMT/2-i)
       IMT4096_Weight(NIMT/2+i)=IMT4096_Weight(NIMT/2-i)
    end forall
   
    IMT2048_Bias=IMT4096_Bias(2::2)
    IMT1024_Bias=IMT4096_Bias(4::4)
    IMT512_Bias=IMT4096_Bias(8::8)
    
    IMT2048_Weight=IMT4096_Weight(2::2)
    IMT1024_Weight=IMT4096_Weight(4::4)
    IMT512_Weight=IMT4096_Weight(8::8)

  end subroutine read_imt_file
 
end program main

---

库伦矩阵元计算: 辅助函数计算

库伦势能矩阵元计算中，需要计算中间辅助函数，其中变量$\xi$取DE积分点坐标，每个下标$m$，计算得到一组数据
$$U_m(\xi )={\int }_0^{\pi }\cos (m\beta )\sqrt{\frac{1-\xi }{1+\xi \cos \beta }}d\beta$$

4096个DE积分点计算matlab程序。其中bias数组是DE积分点坐标。积分中的因子$\sqrt{1-\xi }$是为了减弱边缘奇性引入的(DE积分点边缘处的坐标非常小，或非常接近1)

syms t xi tt x y yy;

U=[];

tt=vpa(pi);
Deta=tt/1000000; !10^-6

for m=0:10
    m
    tU=[];
    for k=1:2048
        k
        xi=bias(k);
        y=cos(m*x)*sqrt(1-xi)/sqrt(1+xi*cos(x));
        t=vpaintegral(y,0,tt,'RelTol',1e-36, 'AbsTol',0, 'MaxFunctionCalls',Inf);
        tU=[tU;t];
    end

    for k=2049:3900
        k
        xi=bias(4096-k);
        y=cos(m*x)*sqrt(xi)/sqrt(2*cos(x/2)^2-xi*cos(x));
        t=vpaintegral(y,0,tt,'RelTol',1e-36, 'AbsTol',0, 'MaxFunctionCalls',Inf);
        tU=[tU;t];
    end

    for k=3901:4090
        k
        xi=bias(4096-k);
        y=cos(m*x)*sqrt(xi)/sqrt(2*cos(x/2)^2-xi*cos(x));
        t=vpaintegral(y,0,tt-Deta,'RelTol',1e-36, 'AbsTol',0, 'MaxFunctionCalls',Inf);
  
        %[PI-Deta,PI]之间,将cos(x)展开
        yy=cos(m*(tt-x))*sqrt(xi)/sqrt(xi+(1-xi)*x^2/2*(1-x^2/12+x^4/360));
        t=t+vpaintegral(yy,0,Deta,'RelTol',1e-36, 'AbsTol',0, 'MaxFunctionCalls',Inf);
        tU=[tU;t];
    end

    for k=4091:4095
        t=0;
        tU=[tU;t];
    end

    U=[U tU];
end