具有向量值被积函数的自适应 Gauss-Kronrod 求积

计算科学 正交 数据结构 r
2021-12-03 16:27:34

所以我正在尝试实现一个 Gauss-Kronrod 自适应正交。也就是说，我要计算

∫abf(x)dx=∑if(xi)wi

其中 f(x) 一次在多个点进行评估以提高效率。要记住的一个关键点是f也是向量值的（这里可能与向量化的双重含义和使用混淆）；在实践中，我同时整合了复数被积函数的实部和虚部，或者被积函数可以是 3D 矢量（电场）等。

我已经制定了 R 中的基本求积步骤（最终将使用 Armadillo 库移植到 C++，任何类似 Matlab 的矩阵语法都可以用于讨论），我首先将积分分成 N 部分，并评估对应于 Gauss 和 Kronrod 节点的加权和。

## nodes and weights borrowed from Octave source code

abscissa = c(-0.9914553711208126e+00, -0.9491079123427585e+00, 
             -0.8648644233597691e+00, -0.7415311855993944e+00, 
             -0.5860872354676911e+00, -0.4058451513773972e+00, 
             -0.2077849550078985e+00,  0.0000000000000000e+00, 
             0.2077849550078985e+00,  0.4058451513773972e+00, 
             0.5860872354676911e+00,  0.7415311855993944e+00, 
             0.8648644233597691e+00,  0.9491079123427585e+00, 
             0.9914553711208126e+00)

weights15 = c(0.2293532201052922e-01,  0.6309209262997855e-01, 
              0.1047900103222502e+00,  0.1406532597155259e+00, 
              0.1690047266392679e+00,  0.1903505780647854e+00, 
              0.2044329400752989e+00,  0.2094821410847278e+00, 
              0.2044329400752989e+00,  0.1903505780647854e+00, 
              0.1690047266392679e+00,  0.1406532597155259e+00, 
              0.1047900103222502e+00,  0.6309209262997855e-01, 
              0.2293532201052922e-01)

weights7  = c(0.1294849661688697e+00,  0.2797053914892767e+00, 
              0.3818300505051889e+00,  0.4179591836734694e+00, 
              0.3818300505051889e+00,  0.2797053914892767e+00, 
              0.1294849661688697e+00)

gauss_kronrod <- function(f, a, b, startN=10, dimf = 2, eps = 1e-3, ...){

  Nquad <- length(abscissa)
  intervals <- seq(a, b, length=startN)
  ## change of variable from [a, b] to [-1, 1]
  shifts <- (intervals[-length(intervals)] + intervals[-1]) / 2
  scalings <- diff(intervals) / 2

  ## scale the weights accordingly
  weightsG <- tcrossprod(weights15, scalings)
  weightsK <- tcrossprod(weights7, scalings)

  ## transformed node positions
  x.scaled <- tcrossprod(abscissa, scalings) +
    matrix(shifts, ncol=length(shifts), nrow=Nquad, byrow=TRUE)

  ## evaluate the function at all points
  fvals <- f(c(x.scaled))

  ## function evaluations are reshaped into a 3D array
  ## rows correspond to dimf values of integrand
  ## columns are the Nquad evalution points in the sub-intervals
  ## slices correspond to the (startN - 1) sub-intervals
  dim(fvals) = c(dimf, Nquad, startN-1)

  ## select which subset of fvals are used for the 7-points kronrod sum
  ikron <- rep(c(FALSE, TRUE), length.out=Nquad)


  gauss <- kronrod <- matrix(0, nrow=dimf, ncol=startN-1)
  for (ii in seq.int(startN-1)){ # integrals for each sub-interval
    gauss[,ii] <- fvals[,,ii] %*% weightsG[,ii]
    kronrod[,ii] <- fvals[,ikron,ii] %*% weightsK[,ii]
  }
  ## relative difference between the two rules 
  errors <- abs(gauss - kronrod) / abs(gauss)
  test <- which(apply(errors, 2, max) > eps)

  ## net integral, just sum the partial sums...
  rowSums(gauss)

}

# a vector-valued integrand...
f <-  function(x) rbind(exp(x/(x+2))*sin(x), cos(x))

gauss_kronrod(f, 2, 110, startN=10)

我的问题是在尝试实施自适应方案时出现的；基本上，在每次迭代中，我想比较每个子区间的gaussvs值kronrod，并从那里决定在哪里进一步细分。什么是在迭代期间存储这些中间结果的好策略，并在最终求和中跟踪哪些子区间和去哪里等？

关于这项工作的更一般的建议也将不胜感激，我通常不热衷于自己实现这样的通用操作，但它似乎比将不同的库链接在一起更容易（最初我正在考虑cuture，因为我的最终目标是 2D一体化）。

2个回答

我建议阅读 Quadpack 书籍（Quadpack, Piessens, Robert; de Doncker-Kapenga, Elise; Überhuber, Christoph W.; Kahaner, David (1983)。QUADPACK: A subroutine package for automatic integration. Springer-Verlag. ISBN 978- 3-540-12553-2）和 Pedro Gonnet 的博士论文“重新审视自适应正交（Adaptive Quadrature Revisited）（可在此处获得 pdf ）。Pedro 是 SciComp SE 的贡献者。

通常做的是优先队列数据结构：

  1. 您在整个初始区间上应用 Gauss-Kronrod 并估计误差
  2. 您将间隔对称地分成两部分，然后在间隔上再次应用 Gauss-Kronrod，估计错误并将一对（或元组）“间隔，错误”添加到优先级队列
  3. 如果总误差小于要求的容差（或者如果达到最大细分数），则停止
  4. 如果不是，则从队列中弹出顶部元素（最高错误），将其一分为二并应用 Gauss-Kronrod，将这两个新元素添加到优先级队列并返回 3。

好吧，我想到了一个策略，如下所示。首先，我创建了一个函数来计算区间列表中的部分总和，并检查哪些足以保留。其次，对于那些不好的间隔，辅助函数将它们一分为二。最后，在每一步中，好的贡献都会被添加到最终结果中。这是R中的一个实现，

gk_segments <- function(f, segments=list(c(2, 50), c(50, 100)),  dimf = 2, 
                        eps = 1e-3, ...){

  Nquad <- length(abscissa)
  Nsegments <- length(segments)
  shifts <- sapply(segments, sum) / 2
  scalings <- sapply(segments, diff) / 2
  weightsG <- tcrossprod(weights15, scalings)
  weightsK <- tcrossprod(weights7, scalings)

  x.scaled <- tcrossprod(abscissa, scalings) +
    matrix(shifts, ncol=length(shifts), nrow=Nquad, byrow=TRUE)

  fvals <- f(c(x.scaled), ...)

  ## function evaluations are reshaped into a 3D array
  ## rows correspond to dimf values of integrand
  ## columns are the Nquad evalution points in the sub-intervals
  ## slices correspond to the Nsegments sub-intervals
  dim(fvals) = c(dimf, Nquad, Nsegments)

  ikron <- rep(c(FALSE, TRUE), length.out=Nquad)

  gauss <- kronrod <- matrix(0, nrow=dimf, ncol=Nsegments)
  for (ii in seq.int(Nsegments)){ # integrals for each sub-interval
    gauss[,ii] <- fvals[,,ii] %*% weightsG[,ii]
    kronrod[,ii] <- fvals[,ikron,ii] %*% weightsK[,ii]
  }

  errors <- abs(gauss - kronrod) / abs(gauss)
  test <- which(apply(errors, 2, max) > eps)

  if(!length(test))
    return( list(converged = rowSums(gauss),
                 todo = list()))

  list(converged = rowSums(gauss[ ,-test]),
       todo = segments[test])
}

# f <-  function(x) rbind(exp(x/(x+2))*sin(x), cos(x))
# gk_segments(f, segments=split_interval(c(2, 110), 10))

# split a given interval [x0, x1] into a list of N subintervals
split_interval <- function(x, N=2){

  breaks <- seq(x[1], x[2], length=N+1)
  lapply(seq.int(N), function(ii) c(breaks[ii], breaks[ii+1]))
}

gk_adaptive <- function(fun, a, b, startN=10, dimf = 2, 
                         eps = 1e-3, maxiter = 100, ...){

  ## first step
  segments <- split_interval(c(a, b), startN)
  tmp <- gk_segments(f=fun, segments=segments, dimf=dimf, eps=eps, ...)

  result <- tmp$converged # partial sums that have already converged
      segments <- unlist(lapply(tmp$todo, split_interval, N=2), recursive=FALSE)
  iter <- 1

  ## we stop if there's no subinterval left to split
  ## or we've reached the maximum number of iterations
  finished <- (!length(tmp$todo)) || (iter > maxiter)

  while(!finished){
    message("iteration #", iter)
    ## work on the new subsegments
    tmp <- gk_segments(f=fun, segments=segments, dimf=dimf, eps=eps, ...)
    result <- result + tmp$converged ## add the good chunks
        segments <- unlist(lapply(tmp$todo, split_interval, N=2), recursive=FALSE)

    iter <- iter + 1
    stopping <- iter > maxiter
    if(stopping)
      message("Maximum number of iteration reached, result is likely wrong")
    finished <- (!length(tmp$todo)) || stopping
  }

    return(result)
}

f <-  function(x) rbind(sin(x), cos(x))
gk_adaptive(f, 2, 110)
gk_adaptive(f, 2, 110, startN=10, eps=1e-10)
gk_adaptive(f, 2, 110, startN=3, eps=1e-10, maxiter=2)

我认为将其推广到 2 维或更多维并不容易，因为选择正确的正方形 / 立方体以迭代地关注似乎更难编码。作为一种效率较低的替代方案，我可以在一维中链接多个积分，一次保持除一个变量之外的所有变量不变。

其它你可能感兴趣的问题
上一篇Matlab中的准蒙特卡罗 下一篇使用流函数围绕非对称障碍物的势流