计算科学 - 使用高斯雅可比求积法或其他方法计算积分 - 吾爱随笔录

我需要整合以下积分：

\begin{aligned} I = \int^{z} \frac{1 - ζ^{2}}{(1 + ζ^{2}) (ζ - ζ_{l}) (1 - ζ_{l} ζ)} \prod_{k = 2}^{n - 1} {(\frac{ζ - z_{k}}{1 - ζ z_{k}})}^{- β_{k}} d ζ \end{aligned}

$\begin{align} I = \int^z\frac{1-\zeta^2}{(1+\zeta^2)(\zeta-\zeta_l)(1-\zeta_l\zeta)}\prod_{k=2}^{n-1}\left ( \frac{\zeta-z_k}{1-\zeta z_k} \right )^{-\beta_k} d\zeta \end{align}$

沿着 argand 平面单位圆上半部分内的任何复杂路径，其中 $z_k$ , $\beta_k$ 和 $\zeta_l$ 有规定值。

我想知道是否有人可以提出一种方法来用他们自己的代码稳健地计算这个积分，或者评论我在下面所做的是否明智：

积分出现在 Kirchoff 流动问题中，是 Schwarz Christoffel 参数问题的一部分。我已经修改了来自 Toby Driscoll 的以下代码，试图计算我的积分。请注意，这假定qdat这是可用的高斯雅可比正交权重矩阵。

正如下面的评论中所解释的，这个想法是改进应用正交的间隔，使得任何奇点 $z_k$ （除了积分的终点，距离至少一半的积分间隔更远。这被称为“自适应或复合”高斯雅可比正交。由于被积函数的形式（本质上它是类型 $\int f(x) (1+x)^\alpha(1-x)^\beta dx$ ) 它适用于高斯雅可比正交。

function I = modified_hpquad(z1,z2,sing1,z,beta,qdat)
%HPQUAD (not intended for calling directly by the user)
%   Numerical quadrature for the half-plane map.

%   Copyright 1998 by Toby Driscoll.
%   $Id: hpquad.m 298 2009-09-15 14:36:37Z driscoll $

%   HPQUAD(z1,z2,sing1,z,beta,qdat)
%   z1,z2 are vectors of left and right endpoints.  sing1 is a vector of
%   integer indices which label the singularities in z1.  So if sing1(5)
%   = 3, then z1(5) = z(3).  A zero means no singularity.  z is the
%   vector of finite singularities; beta is the vector of associated
%   turning angles.  qdat is quadrature data from SCQDATA.
%
%   Make sure z and beta are column vectors.
%   
%   The integral is subdivided, if necessary, so that no singularity
%   lies closer to the left endpoint than 1/2 the length of the
%   integration (sub)interval.


%Some parameters for modified SCT:
sl = 0;


nqpts = size(qdat,1);
% Note: Here n is the total number of *finite* singularities; i.e., the
% number of terms in the product appearing in the integrand.
n = length(z);

bigz = z(:,ones(1,nqpts));
bigbeta = beta(:,ones(1,nqpts));
if isempty(sing1)
  sing1 = zeros(length(z1),1);
end

%I = zeros(size(z1));
I = 0;
%nontriv = find(z1(:)~=z2(:))';
nontriv = find(z1~=z2);

if nontriv == 1

    za = z1; 
    zb = z2;
    sng = sing1;

    if sng > n
        sng = 0; %No singularity 
    end

    % Allowable integration step, based on nearest singularity.
    dist = min(1,2*min(abs(z([1:sng-1,sng+1:n])-za))/abs(zb-za)); % z([1:sng-1,sng+1:n]) jumps over sing
    zr = za + dist*(zb-za);
    ind = rem(sng+n,n+1)+1;

    % Adjust Gauss-Jacobi nodes and weights to interval.
    nd = ((zr-za)*qdat(:,ind) + zr + za).'/2; % G-J nodes
    wt = ((zr-za)/2) * qdat(:,ind+n+1);     % G-J weights
    %terms = nd(ones(n,1),:) - bigz;
    terms = (nd(ones(n,1),:) - bigz)./(1-bigz.*nd(ones(n,1),:)); %Modified SC
    if any(terms(:)==0) 
        % Endpoints are practically coincident.
        I = 0;
    else
        % Use Gauss-Jacobi on first subinterval, if necessary.
        if sng > 0
        terms(sng,:) = terms(sng,:)./abs(terms(sng,:));
        wt = wt*(abs(zr-za)/2)^beta(sng);
        end
        %I(k) = exp(sum(log(terms).*bigbeta,1))*wt;
        I = exp(log((1-nd.^2)./((1+nd.^2).*(nd-sl).*(1-sl.*nd)))+sum(log(terms).*bigbeta,1))*wt;
        while dist < 1              
            % Do regular Gaussian quad on other subintervals.
            zl = zr;
            dist = min(1,2*min(abs(z-zl))/abs(zl-zb));
            zr = zl + dist*(zb-zl);
            nd = ((zr-zl)*qdat(:,n+1) + zr + zl).'/2;
            wt = ((zr-zl)/2) * qdat(:,2*n+2);
            %terms = nd(ones(n,1),:) - bigz;
            terms = (nd(ones(n,1),:) - bigz)./(1-bigz.*nd(ones(n,1),:));
            I = I + exp(log((1-nd.^2)./((1+nd.^2).*(nd-sl).*(1-sl.*nd)))+sum(log(terms).*bigbeta,1))*wt;
        end
    end
end

上面的代码是代码的修改版本，它完全符合我对积分的要求

\begin{aligned} I = \int^{z} \prod_{k = 1}^{n - 1} {(ζ - z_{k})}^{- β_{k}} d ζ \end{aligned}

$\begin{align} I = \int^z\prod_{k=1}^{n-1}\left ( \zeta-z_k \right )^{-\beta_k} d\zeta \end{align}$

我已将该代码复制下来以供参考：

function I = hpquad(z1,z2,varargin)
%HPQUAD (not intended for calling directly by the user)
%   Numerical quadrature for the half-plane map.

%   Copyright 1998 by Toby Driscoll.
%   $Id: hpquad.m 298 2009-09-15 14:36:37Z driscoll $

%   HPQUAD(z1,z2,sing1,z,beta,qdat)
%   z1,z2 are vectors of left and right endpoints.  sing1 is a vector of
%   integer indices which label the singularities in z1.  So if sing1(5)
%   = 3, then z1(5) = z(3).  A zero means no singularity.  z is the
%   vector of finite singularities; beta is the vector of associated
%   turning angles.  qdat is quadrature data from SCQDATA.
%
%   Make sure z and beta are column vectors.
%   
%   HPQUAD integrates from a possible singularity at the left end to a
%   regular point at the right.  If both endpoints are singularities,
%   you must break the integral into two pieces and make two calls, or
%   call HPQUAD(z1,z2,sing1,sing2,z,beta,qdat) and accept an automatic
%   choice. 
%   
%   The integral is subdivided, if necessary, so that no singularity
%   lies closer to the left endpoint than 1/2 the length of the
%   integration (sub)interval.

if nargin==7
  % Break into two pieces with recursive call.
  [sing1,sing2,z,beta,qdat] = deal(varargin{:});
  mid = (z1+z2)/2;
  mid = mid + 1i*abs(mid);
  I1 = hpquad(z1,mid,sing1,z,beta,qdat);
  I2 = hpquad(z2,mid,sing2,z,beta,qdat);
  I = I1-I2;
  return
else
  [sing1,z,beta,qdat] = deal(varargin{:});
end




nqpts = size(qdat,1);
% Note: Here n is the total number of *finite* singularities; i.e., the
% number of terms in the product appearing in the integrand.
n = length(z);
bigz = z(:,ones(1,nqpts));
bigbeta = beta(:,ones(1,nqpts));
if isempty(sing1)
  sing1 = zeros(length(z1),1);
end

I = zeros(size(z1));
nontriv = find(z1(:)~=z2(:))';



for k = nontriv
  za = z1(k);
  zb = z2(k);
  sng = sing1(k);

  % Allowable integration step, based on nearest singularity.
  dist = min(1,2*min(abs(z([1:sng-1,sng+1:n])-za))/abs(zb-za));
%%  if isempty(dist), dist=1; end
  zr = za + dist*(zb-za);
  ind = rem(sng+n,n+1)+1;
  % Adjust Gauss-Jacobi nodes and weights to interval.
  nd = ((zr-za)*qdat(:,ind) + zr + za).'/2; % G-J nodes
  wt = ((zr-za)/2) * qdat(:,ind+n+1);   % G-J weights
  terms = nd(ones(n,1),:) - bigz;
  if any(terms(:)==0) 
    % Endpoints are practically coincident.
    I(k) = 0;
  else
    % Use Gauss-Jacobi on first subinterval, if necessary.
    if sng > 0
      terms(sng,:) = terms(sng,:)./abs(terms(sng,:));
      wt = wt*(abs(zr-za)/2)^beta(sng);
    end
    I(k) = exp(sum(log(terms).*bigbeta,1))*wt;
    while dist < 1              
      % Do regular Gaussian quad on other subintervals.
      zl = zr;
      dist = min(1,2*min(abs(z-zl))/abs(zl-zb));
      zr = zl + dist*(zb-zl);
      nd = ((zr-zl)*qdat(:,n+1) + zr + zl).'/2;
      wt = ((zr-zl)/2) * qdat(:,2*n+2);
      terms = nd(ones(n,1),:) - bigz;
      I(k) = I(k) + exp(sum(log(terms).*bigbeta,1)) * wt;
    end
  end
end