使用Numeric.Integration.TanhSinh进行N维集成

Question

我在Haskell中使用Numeric.Integration.TanhSinh进行数值积分.这定义了一个函数

parTrap :: (Double -> Double) -> Double -> Double -> [Result]

其中第一个参数是要集成的1维函数,然后是上限和下限.我有一个包装器,将此功能转换为

ttrap f xmin xmax = (ans, err)
   where
      res = absolute 1e-6 $ parTrap f xmin xmax
      ans = result res
      err = errorEstimate res

为了集成二维功能,我可以使用

 ttrap2 f y1 y2 g1 g2 = ttrap h y1 y2 -- f ylower yupper (fn for x lower) (fn for x upper) 
       where
           h y = fst $ ttrap (flip f y) (g1 y) (g2 y)

ttrap2_fixed f y1 y2 x1 x2 = ttrap2 f y1 y2 (const x1) (const x2)

ttrap2_fixed我的想法是,我现在可以在函数(Double - > Double - > Double)中进行双积分,并且边界是y1 y2 x1 x2.

使用这种模式,我可以定义更高阶的积分函数

ttrap3_fixed :: (Double -> Double -> Double -> Double) -> Double -> Double -> Double -> Double -> Double -> Double -> Double
ttrap3_fixed f z1 z2 y1 y2 x1 x2 = fst $ ttrap h z1 z2
  where
    h z  = fst $ ttrap2_fixed (f z)  x1 x2 y1 y2


ttrap4_fixed :: (Double -> Double -> Double -> Double -> Double) -> Double -> Double -> Double -> Double -> Double -> Double -> Double -> Double -> Double
ttrap4_fixed f w1 w2 z1 z2 y1 y2 x1 x2 = fst $ ttrap h w1 w2
  where
    h w  = ttrap3_fixed (f w) z1 z2 x1 x2 y1 y2


ttrap5_fixed :: (Double -> Double -> Double -> Double -> Double -> Double) -> Double -> Double -> Double -> Double -> Double -> Double -> Double -> Double -> Double -> Double -> Double
ttrap5_fixed f u1 u2 w1 w2 z1 z2 y1 y2 x1 x2 = fst $ ttrap h u1 u2
  where
    h u  = ttrap4_fixed (f u) w1 w2 z1 z2 x1 x2 y1 y2

但是,我想集成一个类型的函数

f :: [Double] -> Double

理念是函数的维度可以在程序中变化.理想情况下,我想要一个类型的函数

int_listfn :: ([Double] -> Double) -> [(Double, Double)] -> Double

我可以在一个边界元组列表上集成多维函数.作为其中的一部分,似乎我需要使用诸如延续传递样式之类的东西来构建积分器函数,但这是我遇到的问题.提前致谢.

---

f我希望整合的一个例子是

> let f [x,y,z] = x**3.0 * sin(y) + (1.0/z)

考虑边界

> let bounds = [(1.0,2.0),(2.0,4.0), (0.0,3.0)]
> int_listfn f bounds

这应该相当于计算

---

编辑:添加另一个示例函数

f1 :: Double -> Double
f1 x = 1.0 * x


fn_maker :: [Double -> Double] -> ([Double] -> Double)
fn_maker inlist = myfn
  where
    myfn xlist = product $ zipWith (\f x -> f x) inlist xlist


m = 4

f_list = fn_maker (replicate f1 m)

f_list有类型[Double] -> Double,相当于f x y z w = x * y * z * w.我认为类型[Double] -> Double是合适的,因为函数的维度m是一个参数.也许我需要改变设计.

---

列举函数到 curried 函数,@ findruk我一直在我的代码中使用它并且似乎工作,虽然我需要根据边界列表的大小跟踪调用哪个函数.

static_1 :: ([Double] -> Double) -> (Double -> Double)
static_1 f = f'
  where
    f' x = f [x]


static_2 :: ([Double] -> Double) -> (Double -> Double -> Double)
static_2 f = f'
  where
    f' x y = f [x,y]


static_3 :: ([Double] -> Double) -> (Double -> Double -> Double -> Double)
static_3 f = f'
  where
    f' x y z = f [x,y,z]


static_4 :: ([Double] -> Double) -> (Double -> Double -> Double -> Double -> Double)
static_4 f = f'
  where
    f' x y z w = f [x,y,z,w]


static_5 :: ([Double] -> Double) -> (Double -> Double -> Double -> Double -> Double -> Double)
static_5 f = f'
  where
    f' x y z w u = f [x,y,z,w,u]


int_listfn :: ([Double] -> Double) -> [(Double, Double)] -> Double
int_listfn f bounds  
  | f_dim == 1 = fst $ ttrap (static_1 f) (fst (bounds !! 0)) (snd (bounds !! 0))
  | f_dim == 2 = fst $ ttrap2_fixed  (static_2 f) (fst (bounds !! 0)) (snd (bounds !! 0)) (fst (bounds !! 1)) (snd (bounds !! 1))
  | f_dim == 3 = ttrap3_fixed  (static_3 f) (fst (bounds !! 0)) (snd (bounds !! 0)) (fst (bounds !! 1)) (snd (bounds !! 1)) (fst (bounds !! 2)) (snd (bounds !! 2))
  | f_dim == 4 = ttrap4_fixed  (static_4 f) (fst (bounds !! 0)) (snd (bounds !! 0)) (fst (bounds !! 1)) (snd (bounds !! 1)) (fst (bounds !! 2)) (snd (bounds !! 2)) (fst (bounds !! 3)) (snd (bounds !! 3))
  | f_dim == 5 = ttrap5_fixed  (static_5 f) (fst (bounds !! 0)) (snd (bounds !! 0)) (fst (bounds !! 1)) (snd (bounds !! 1)) (fst (bounds !! 2)) (snd (bounds !! 2)) (fst (bounds !! 3)) (snd (bounds !! 3)) (fst (bounds !! 3)) (snd (bounds !! 3))
  | otherwise = error "Unsupported integral size"
  where
    f_dim = length bounds

Answer 1

这个想法

为了方便起见，我引入了这些类型别名：

type Res   = (Double, Double)
type Bound = (Double, Double)

让我们仔细看看以下类型ttrap_fixedN：

ttrap :: (Double -> Double) -> Double -> Double -> Res
ttrap_fixed2 :: (Double -> Double -> Double) -> Double -> Double -> Double -> Double -> Res

显然，我们可以将边界配对以获得更短、更清晰的版本：

ttrap :: (Double -> Double) -> Bound -> Res
ttrap_fixed2 :: (Double -> Double -> Double) -> Bound -> Bounds -> Res

此外，我们可以Bounds一起收集N> 1：

ttrap_fixed2 :: (Double -> Double -> Double) -> (Bound, Bound) -> Res
ttrap_fixed3 :: (Double -> Double -> Double -> Double) -> (Bound, Bound, Bound) -> Res

请注意我们如何使所有ttrap_fixedN函数只接受 2 个参数。另请注意，第二个参数的类型（带有Bounds 的元组的数量）取决于第一个参数（要积分的函数的数量）。

现在应该清楚的是，通用ttrap_fixed函数将取决于给定函数的数量，我们需要类型类来实现这种多态性。除了integrate方法（这是一般方法ttrap_fixed）之外，该类还需要第二个参数的关联类型同义词：

class Integrable r where
  type Bounds r :: *
  integrate :: r -> Bounds r -> Res

解决方案草图

我们将需要以下扩展：

{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE FlexibleInstances #-}

import Numeric.Integration.TanhSinh

这是我将使用的问题中的唯一函数：

ttrap :: (Double -> Double) -> Double -> Double -> Res
ttrap f xmin xmax = (ans, err)
   where
      res = absolute 1e-6 $ parTrap f xmin xmax
      ans = result res
      err = errorEstimate res

定义Integrable类型类：

class Integrable r where
  type Bounds r :: *
  integrate :: r -> Bounds r -> Res

这是实例

instance Integrable Double where
  type Bounds Double = ()
  integrate x _ = (x, 0)

instance Integrable r => Integrable (Double -> r) where
  type Bounds (Double -> r) = (Bound, Bounds r)
  integrate f ((xmin, xmax), args) = ttrap g xmin xmax
    where
      g x = fst $ integrate (f x) args

去测试一下：

main :: IO ()
main = do
  let f :: Double -> Double -> Double -> Double
      f x y z = x ** 3.0 * sin(y) + (1.0/z)
      zbounds = (1.0, 2.0)
      ybounds = (2.0, 4.0)
      xbounds = (0.0, 3.0)
      res = integrate f (xbounds, (ybounds, (zbounds, ())))
  print res -- prints (8.9681929657648,4.732074732061164e-10)

柯里化Bounds

您可能已经注意到，我们现在没有非常方便的Bounds 嵌套元组。如果返回一个柯里化函数那就太好了integrate f，例如：

integrate (*) :: Bound -> Bound -> Res

代替

integrate (*) :: (Bound, (Bound, ())) -> Res

不幸的是，我未能找到任何简单的重构方法来Integrable实现这一点。然而，我们可以用另一种类型类黑客技术来解决这个问题。这个想法是引入多态curryBounds：

class CurryBounds bs where
  type Curried bs a :: *
  curryBounds :: (bs -> a) -> Curried bs a

这些实例很简单：

instance CurryBounds () where
  type Curried () a = a
  curryBounds f = f ()

instance CurryBounds bs => CurryBounds (b, bs) where
  type Curried (b, bs) a = b -> Curried bs a
  curryBounds f = \x -> curryBounds (\xs -> f (x, xs))

现在我们可以定义一个更好的版本integrate：

integrate' :: (Integrable r, CurryBounds (Bounds r)) => r -> Curried (Bounds r) Res
integrate' = curryBounds . integrate

例子

>>> let f x y z = x**3.0 * sin(y) + (1.0/z) :: Double
>>> integrate' f (0, 3) (2, 4) (1, 2)
(8.9681929657648,4.732074732061164e-10)