множественные интегралы и матрицы (сложные элементы) из числовой прелюдии - PullRequest
Новая
       21

множественные интегралы и матрицы (сложные элементы) из числовой прелюдии

1 голос
/ 06 марта 2012

Мне нужно вычислить интегралы вида 

/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, я получаю ту же ошибку.

Что я могу сделать, чтобы решить эту проблему?

Спасибо.

1 Ответ

1 голос
/ 06 марта 2012

Ты мог сделать 

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] ]

но я согласен, что это довольно некрасиво, хотя лень Хаскелла должна гарантировать, что не все элементы матрицы будут вычисляться на каждом шаге интеграции.

...