5

在将向量空间包用于导数塔时（请参阅导数塔），我遇到了对积分进行微分的需要。从数学上很清楚如何实现这一点：

f(x) = int g(y) dy from 0 to x

有一个功能

g : R -> R

例如。

关于 x 的导数是：

f'(x) = g(x)

我试图通过首先定义一个“集成”类来获得这种行为

class Integration a b where
--standard integration function
integrate :: (a -> b) -> a -> a -> b

一个基本的例子是

instance  Integration Double Double where
  integrate f a b = fst $ integrateQAGS prec 1000 f a b

integrateQAGS来自hmatrix _

问题在于代表衍生品塔的值 b：

instance Integration Double (Double :> (NC.T Double)) where
  integrate = integrateD

NC.T来自 Numeric.Complex (numeric-prelude)。该函数integrateD定义如下（但错误）：

integrateD ::(Integration a b, HasTrie (Basis a), HasBasis a, AdditiveGroup b) =>  (a -> a :> b) -> a -> a -> (a :> b)
integrateD f l u = D (integrate (powVal . f) l u) (derivative $ f u)

该函数不返回我想要的，它导出被积函数，但不是积分。问题是，我需要一个返回的线性映射f u。a :> b定义如下：

data a :> b = D { powVal :: b, derivative :: a :-* (a :> b) }

不知道怎么定义derivative。任何帮助将不胜感激，谢谢

编辑：

我忘了提供实例Integration Double (NC.T Double)：

instance  Integration Double (NC.T Double) where
  integrate f a b = bc $ (\g -> integrate g a b) <$> [NC.real . f, NC.imag . f]
      where bc (x:y:[]) = x NC.+: y

我可以举一个例子来说明我的意思：假设我有一个函数

f(x) = exp(2*x)*sin(x)

>let f = \x -> (Prelude.exp ((pureD 2.0) AR.* (idD x))) * (sin (idD x)) :: Double :> Double

(AR.*) 表示从 Algebra.Ring (numeric-prelude) 乘法

我可以轻松地将此功能与上述功能集成integrateD：

>integrateD f 0 1 :: Double :> Double
D 1.888605715258933 ...

当我看一下 f 的导数时：

f'(x) = 2*exp(2*x)*sin(x)+exp(2*x)*cos(x)

并对此进行评估0，pi/2我得到1了一些价值：

> derivAtBasis (f 0.0) ()
D 1.0 ...

> derivAtBasis (f (pi AF./ 2)) ()
D 46.281385265558534 ...

现在，在推导积分时，我得到的是函数的推导，而f不是它在上限的值

> derivAtBasis (integrate f 0 (pi AF./ 2)) ()
D 46.281385265558534 ...

但我期望：

> f (pi AF./ 2)
D 23.140692632779267 ...
4

3 回答 3

1

如果您只想对涉及数字积分的功能进行 AD，而 AD 系统本身不知道积分，它应该“正常工作”。这是一个例子。（这个集成例程非常棘手，因此得名。）

import Numeric.AD
import Data.Complex

intIcky :: (Integral a, Fractional b) => a -> (b -> b) -> b -> b -> b
intIcky n f a b = c/n' * sum [f (a+fromIntegral i*c/(n'-1)) | i<-[0..n-1]]
  where n' = fromIntegral n
        c = b-a

sinIcky t = intIcky 1000 cos 0 t
cosIcky t = diff sinIcky t

test1 = map sinIcky [0,pi/2..2*pi::Float]
 -- [0.0,0.9997853,-4.4734867e-7,-0.9966421,6.282018e-3]
test2 = map sin [0,pi/2..2*pi::Float]
 -- [0.0,1.0,-8.742278e-8,-1.0,-3.019916e-7]
test3 = map cosIcky [0,pi/2..2*pi::Float]
 -- [1.0,-2.8568506e-4,-0.998999,2.857402e-3,0.999997]
test4 = map cos [0,pi/2..2*pi::Float]
 -- [1.0,-4.371139e-8,-1.0,1.1924881e-8,1.0]
test5 = diffs sinIcky (2*pi::Float)
 -- [6.282019e-3,0.99999696,-3.143549e-3,-1.0004976,3.1454563e-3,1.0014982,-3.1479746e-3,...]
test6 = diffs sinIcky (2*pi::Complex Float)
 -- [6.282019e-3 :+ 0.0,0.99999696 :+ 0.0,(-3.143549e-3) :+ 0.0,(-1.0004976) :+ 0.0,...]

唯一需要注意的是，数值积分例程需要与 AD 很好地配合使用，并且还需要接受复杂的参数。更天真的东西，比如

intIcky' dx f x0 x1 = dx * sum [f x|x<-[x0,x0+dx..x1]]

在积分的上限中是分段常数，要求积分的限制是 Enum 并且因此是非复数的，并且由于这个原因，还经常评估给定范围之外的被积函数：

Prelude> last [0..9.5]
10.0
于 2013-04-10T11:53:20.350 回答
0

'hmatrix' 与 Double 密切相关。您不能将其函数与其他数字数据类型一起使用，例如“向量空间”或“广告”提供的那些。

于 2012-07-11T15:32:19.133 回答
0

我终于找到了我的问题的解决方案。解决方案的关键是>-<向量空间包中的函数，它代表链式法则。

所以，我定义了一个integrateD'这样的函数：

integrateD' :: (Integration a b, HasTrie (Basis a), HasBasis a, AdditiveGroup b , b ~ Scalar b, VectorSpace b) => (a -> a :> b) -> a -> a -> (a:>b) -> (a :> b)
integrateD' f l u d_one =  ((\_ -> integrate (powVal . f) l  (u)) >-< (\_ ->  f u)) (d_one)

这d_one意味着派生变量及其导数必须为1。使用此函数，我现在可以创建一些实例，例如

instance Integration Double (Double :> Double) where
integrate f l u = integrateD' f l u (idD 1)

和

instance Integration ( Double) (Double :> (NC.T Double)) where
integrate f l u = liftD2 (NC.+:) (integrateD' (\x -> NC.real <$>> f x) l u (idD 1.0 :: Double :> Double)) (integrateD' (\x -> NC.imag <$>> f x) l u (idD 1.0 :: Double :> Double))

不幸的是，我不能使用integrateD开箱即用的复杂值，我必须使用liftD2. 原因似乎是idD功能，不知道有没有更优雅的解决方案。

当我查看问题中的示例时，我现在得到了我想要的解决方案：

*Main> derivAtBasis (integrateD' f 0 (pi AF./ 2) (idD 1.0 :: Double :> Double )) ()
D 23.140692632779267 ...

或通过使用实例：

*Main> derivAtBasis (integrate f 0 (pi AF./ 2)) ()
D 23.140692632779267 ...
于 2012-08-09T13:25:04.790 回答