algorithm - 如何正确地将龙格误差估计规则添加到这个例子中？

Question

我有一个并行计算某个积分的算法。如果您使用多线程，此解决方案可提供非常好的时间加速。而且线程越多，计算速度越快。我测试到-N4，加速因子达到了8。即在4核上启动程序是积分计算比在1核上启动该程序快8倍。但是我想添加一个用于估计龙格误差的规则。从现在开始，为了增加积分计算的准确性，有必要增加 N。这表示我们需要打破原始段的多少部分。我该怎么做？

import Data.Time
import System.Environment
import Data.Massiv.Array as A

main = do
    begT <- getCurrentTime
    putStrLn $ show $ integrateA 100000 f 0.00005 10000
    endT <- getCurrentTime
    putStrLn $ init $ show $ diffUTCTime endT begT

f :: Double -> Double
f x = sin x * cos x*x*x

integrateA :: Int -> (Double -> Double) -> Double -> Double -> Double
integrateA n f a b =
 let step = (b - a) / fromIntegral n
  sz = size segments - 1
  segments = computeAs P $ A.map f (enumFromStepN Par a step (Sz (n + 1)))
  area y0 y1 = step * (y0 + y1) / 2
  areas = A.zipWith area (extract' 0 sz segments) (extract' 1 sz segments)
 in A.sum areas

发射示例:

Answer 1

您无需更改已提供的积分估计器中的任何内容，即可使用朗格规则为其增加精度。我认为是这样的:

-- | Returns estimated integral up to a precision, or value estimated at max
-- number of steps
rungeRule ::
     Int -- ^ Maximum number of steps as an upper bound, to prevent unbounded computation
  -> Double -- ^ ε -- precision
  -> Int -- ^ Starting value of @n@
  -> Double -- ^ Θ -- ^ Either 1/3 for trapezoidal and midpoint or 1/15 for Simpson's
  -> (Int -> Double) -- ^ Integral estimator
  -> Either Double (Int, Double)
rungeRule nMax epsilon n0 theta integralEstimator =
  go (integralEstimator n0) (2 * n0)
  where
    go prevEstimate n
      | n >= nMax = Left prevEstimate
      | theta * abs (curEstimate - prevEstimate) < epsilon =
        Right (n, curEstimate)
      | otherwise = go curEstimate (2 * n)
      where
        curEstimate = integralEstimator n

trapezoidal ::
     Double -- ^ ε -- precision
  -> (Double -> Double) -- ^ f(x) - function to integrate
  -> Double -- ^ a - from
  -> Double -- ^ b - to
  -> Either Double (Int, Double)
trapezoidal epsilon f a b =
  rungeRule 100000 epsilon 10 (1 / 3) (\n -> integrateA n f a b)

如果我们运行它，我们会得到有希望的结果:

λ> trapezoidal 0.0005 (\x -> x * x) 10 20
Right (640,2333.333740234375)
λ> trapezoidal 0.00005 (\x -> x * x) 10 20
Right (2560,2333.3333587646484)
λ> trapezoidal 0.00000005 (\x -> x * x) 10 20
Right (81920,2333.3333333581686)
λ> trapezoidal 0.000000005 (\x -> x * x) 10 20
Left 2333.3333333581686

边注:

• 您的函数 f 按照您编写它的方式表明:

  • 你期望:f x = (sin x) * (cos (x*x*x))

  • 实际情况是:f x = (sin x) * (cos x) * x * x

编辑:

上面给出的解决方案足够通用，适用于所有积分近似规则。但是在龙格规则的每次迭代中都会发生一些重复的工作，在梯形规则的情况下，每次都会重新计算一半的元素，我认为这是一种潜在的优化。接下来是 massiv 的更高级用法，因此我无法详细说明它是如何工作的，除了 segments 数组传递的事实用于访问在上一步计算的值。

trapezoidalMemoized ::
     Int
  -> Array P Ix1 Double
  -> (Double -> Double)
  -> Double
  -> Double
  -> (Double, Array P Ix1 Double)
trapezoidalMemoized n prevSegments f a b =
  let step = (b - a) / fromIntegral n
      sz = size segments - 1
      curSegments = 
        fmap f (enumFromStepN Seq (a + step) (2 * step) (Sz (n `div` 2)))
      segments =
        computeAs P $
        makeLoadArrayS (Sz (n + 1)) 0 $ \w -> do
          A.iforM_ prevSegments $ \i e -> w (i * 2) e
          A.iforM_ curSegments $ \i e -> w (i * 2 + 1) e
      area y0 y1 = step * (y0 + y1) / 2
      areas = A.zipWith area segments (extract' 1 sz segments)
   in (A.sum areas, segments)


trapezoidalRungeMemo ::
     Double -- ^ ε -- precision
  -> (Double -> Double) -- ^ f(x) - function to integrate
  -> Double -- ^ a - from
  -> Double -- ^ b - to
  -> Either Double (Int, Double)
trapezoidalRungeMemo epsilon f a b = go initEstimate initSegments 4
  where
    (initEstimate, initSegments) =
      trapezoidalMemoized 2 (A.fromList Seq [f a, f b]) f a b
    nMax = 131072 -- 2 ^ 17
    theta = 1 / 3
    go prevEstimate prevSegments n
      | n >= nMax = Left prevEstimate
      | theta * abs (curEstimate - prevEstimate) < epsilon =
        Right (n, curEstimate)
      | otherwise = go curEstimate curSegments (2 * n)
      where
        (curEstimate, curSegments) =
          trapezoidalMemoized n prevSegments f a b

使其可并行化甚至更加复杂:

-- Requires additional import: `Data.Massiv.Array.Unsafe`
trapezoidalMemoizedPar ::
     Int
  -> Array P Ix1 Double
  -> (Double -> Double)
  -> Double
  -> Double
  -> (Double, Array P Ix1 Double)
trapezoidalMemoizedPar n prevSegments f a b =
  let step = (b - a) / fromIntegral n
      sz = size segments - 1
      curSegments =
        fmap f (enumFromStepN Seq (a + step) (2 * step) (Sz (n `div` 2)))
      segments =
        computeAs P $
        unsafeMakeLoadArray Par (Sz (n + 1)) Nothing $ \scheduler _ w -> do
          splitLinearlyWith_
            scheduler
            (unSz (size prevSegments))
            (unsafeLinearIndex prevSegments) $ \i e -> w (i * 2) e
          splitLinearlyWith_
            scheduler
            (unSz (size curSegments))
            (unsafeLinearIndex curSegments) $ \i e -> w (i * 2 + 1) e
      area y0 y1 = step * (y0 + y1) / 2
      areas = A.zipWith area segments (extract' 1 sz segments)
   in (A.sum areas, segments)