Как правильно добавить правило оценки ошибок Рунге в этот пример?

Question

У меня есть алгоритм для параллельного вычисления определенного интеграла.Это решение дает очень хорошее ускорение времени, если вы используете несколько потоков.И чем больше потоков, тем быстрее расчет.Я тестировал до -N4 и коэффициент ускорения достигал 8. Т.е. запуск программы на 4 ядрах - это расчет интеграла в 8 раз быстрее, чем запуск этой программы на 1 ядре.Но я хочу добавить правило для оценки ошибки Рунге.Поскольку теперь, чтобы повысить точность вычисления интеграла, необходимо увеличить 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

Edit :

Представленное выше решение достаточно общее, чтобы работать для всех интегральных правил аппроксимации.Но есть некоторая дублирующая работа, происходящая на каждой итерации правила Рунге, в случае трапецеидального правила половина элементов каждый раз пересчитывается каждый раз, что я рассматривал как потенциальную оптимизацию.Далее будет более продвинутое использование 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)