我有一个并行计算某个积分的算法。 如果您使用多线程,此解决方案可提供非常好的时间加速。而且线程越多,计算速度越快。 我测试到-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
最佳答案
您无需更改已提供的积分估计器中的任何内容,即可使用朗格规则为其增加精度。我认为是这样的:
-- | 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)
关于algorithm - 如何正确地将龙格误差估计规则添加到这个例子中?,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/56395599/