1

我需要计算形式的积分

/t     /x
| P(x)*| P(y) dydx
/t0    /t0

其中 P 是一个函数R -> C^(nxn)，通常是一个矩阵，我想在 Haskell 中执行此操作。我已经为标量函数实现了这一点：

import Numeric.GSL.Integration
import Data.Complex
import Data.List

prec :: Double
prec = 1E-9

integrate :: (Double -> Double) -> Double -> Double -> Double
integrate f a b = fst $ integrateQNG prec f a b 

integrateC :: (Double -> Complex Double) -> Double -> Double -> Complex Double
integrateC cf a b = (integrate (\x -> realPart (cf x)) a b :+ integrate (\x -> imagPart (cf x)) a b)

multipleIntegration :: Int -> (Double -> Complex Double) -> Double -> (Double -> Complex Double)
multipleIntegration n f a =  foldl' (\ acc g' -> (\ x -> integrateC (g'*acc) a x)) (\_ -> 1:+0) (replicate n f)

到目前为止，这很有效，尽管对于 n>5 来说速度很慢。

现在我需要将此计算扩展到矩阵，我用数字前奏尝试了这个，因为我可以将函数作为矩阵的元素。我能够积分 的矩阵Double -> Complex Double，但我在积分内乘以矩阵的实际目标失败了，首先是我的代码：

import MathObj.Matrix as Mat
import Algebra.Ring as AR
import Control.Applicative
import qualified Prelude as P
import Prelude hiding ((*))
import Number.Complex as NC
import Numeric.GSL.Integration
import Data.List

type Complex a = NC.T a

prec :: Double
prec = 1E-9

testMat :: Mat.T (Double -> Complex Double)
testMat = Mat.fromRows 2 2 [[\x-> 0.5 +: 2*x,\y-> cos y +: sin y],[\x-> 0.1*x +:x,\_-> 1 +: 1]]

integrateC :: (Double -> Complex Double) -> Double -> Double -> Complex Double
integrateC cf a b = (integrate (\x -> real (cf x)) a b +: integrate (\x -> imag (cf x)) a b)

integrate :: (Double -> Double) -> Double -> Double -> Double
integrate f a b = fst $ integrateQNG prec f a b 

integrateCMat' :: Mat.T (Double -> Complex Double) -> Double -> Mat.T (Double -> Complex Double)
integrateCMat' cmf a =  ((\f -> integrateC f a ) <$> cmf)

multipleIntegrationMat :: Int -> Mat.T (Double -> Complex Double) -> Double -> Mat.T (Double -> Complex Double)
multipleIntegrationMat n mf a =  integrateCMat' ( testMat * (integrateCMat' testMat a)) a

这里，multipleIntegrationMat只是一个测试功能，我没有使用折叠，所以 n 是多余的。错误信息是：

matmul.hs:27:59:
No instance for (C (Double -> Complex Double))
  arising from a use of `*'
Possible fix:
  add an instance declaration for (C (Double -> Complex Double))
In the first argument of `integrateCMat'', namely
  `(testMat * (integrateCMat' testMat a))'
In the expression:
  integrateCMat' (testMat * (integrateCMat' testMat a)) a
In an equation for `multipleIntegrationMat':
    multipleIntegrationMat n mf a
      = integrateCMat' (testMat * (integrateCMat' testMat a)) a
Failed, modules loaded: none.

我知道没有函数乘法的实例。对于这种情况，最好的方法是什么？另一方面，在标量示例中，乘法有效，尽管复数数据类型取自Data.Complex. 当我尝试使用标量示例时，Number.Complex我得到了同样的错误。

我能做些什么来解决这个问题？

谢谢你。

4

1 回答 1

1

你可以做

integrateCMat' :: (Double -> Mat.T (Complex Double))
                -> Double -> Double
                 -> Mat.T (Complex Double)
integrateCMat' cmf a b = Mat.fromRows n m integratedRows 
     where (n,m) = Mat.dimension(cmf undefined)
           integratedCell idx = integrateC (cellFunction idx) a b
           cellFunction idx = (\(Mat.Cons arr) -> arr ! idx) . cmf
           integratedRows = [ [ integratedCell(i,j) | i<-[1..n] ] | j<-[1..m] ]

但我同意它相当难看，尽管 Haskell 的懒惰应该确保在每个集成步骤中不会计算所有矩阵条目。

于 2012-03-06T15:46:44.630 回答