2012-07-11 18 views
5

Khi sử dụng gói vector-space cho tháp phái sinh (xem derivative towers) Tôi thấy cần phải phân biệt các tích phân. Từ toán học nó là khá rõ ràng làm thế nào để đạt được điều này:Cách phân biệt tích phân với thư viện không gian véc tơ (haskell)

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

với một chức năng

g : R -> R 

ví dụ.

Các phái sinh đối với x sẽ là:

f'(x) = g(x) 

tôi đã cố gắng để có được hành vi này bằng cách đầu tiên xác định một lớp "Hội nhập"

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

một ví dụ cơ bản là

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

với integrateQAGS từ hmatrix

vấn đề đi kèm với giá trị b đại diện cho tháp của các dẫn xuất:

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

NC.T là từ Numeric.Complex (số-khúc dạo đầu). Chức năng integrateD được định nghĩa như sau (nhưng sai):

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) 

Chức năng không trả lại những gì tôi muốn, nó xuất phát integrand, nhưng không phải là không thể thiếu. Vấn đề là, tôi cần một bản đồ tuyến tính trả về f u. a :> b được định nghĩa như sau:

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

Tôi không biết cách xác định derivative. Bất kỳ trợ giúp sẽ được đánh giá cao, nhờ

chỉnh sửa:

Tôi quên để cung cấp các ví dụ cho 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 

và tôi có thể đưa ra một ví dụ về những gì tôi có nghĩa là: Hãy nói rằng tôi có một chức năng

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

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

(AR. *) có nghĩa là nhân từ Đại số.Vòng (số-khúc dạo đầu)

tôi có thể dễ dàng tích hợp chức năng này với các chức năng trên integrateD:

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

Khi tôi hãy nhìn vào đạo hàm của f:

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

và đánh giá này tại 0pi/2 Tôi nhận được 1 và một số giá trị:

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

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

Bây giờ, khi phát sinh sự không thể thiếu, tôi nhận được nguồn gốc của hàm f không giá trị của nó ở phía trên bên ràng buộc

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

Nhưng tôi mong đợi:

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

Trả lời

0

Cuối cùng tôi đã tìm được giải pháp cho câu hỏi của mình. Chìa khóa của giải pháp là hàm >-< từ gói không gian vectơ, nó là viết tắt của quy tắc chuỗi.

Vì vậy, tôi định nghĩa một hàm integrateD' như thế này:

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) 

các d_one nghĩa một biến nguồn gốc và đạo hàm của nó phải 1. Với chức năng này ngay bây giờ tôi có thể tạo ra một số trường hợp như

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

thật không may Tôi không thể sử dụng integrateD về các giá trị phức tạp ra khỏi hộp, tôi phải sử dụng liftD2. Lý do cho điều này có vẻ là chức năng idD, tôi không biết nếu có một giải pháp thanh lịch hơn.

Khi tôi nhìn vào ví dụ trong câu hỏi, tôi bây giờ có được giải pháp mong muốn của tôi:

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

hoặc bằng cách sử dụng các ví dụ:

*Main> derivAtBasis (integrate f 0 (pi AF./ 2))() 
D 23.140692632779267 ... 
0

'hmatrix' được gắn rất chặt chẽ với Gấp đôi. Bạn không thể sử dụng các hàm của nó với các kiểu dữ liệu số khác như các kiểu dữ liệu được cung cấp bởi 'vectơ-không gian' hoặc 'quảng cáo'.

+0

Đúng vậy. Nhưng nó hoạt động, khi tôi sử dụng hàm 'powVal' trên một giá trị' Double:> Double'. Nhưng tất nhiên, tôi mất thông tin về đạo hàm. Tôi phải cung cấp thông tin này một cách rõ ràng, và đó là nơi tôi bị mắc kẹt: ( – TheMADMAN

1

Nếu bạn chỉ muốn làm AD trên chức năng liên quan đến tích hợp số, không có hệ thống AD biết về tích hợp mỗi-se, nó sẽ "chỉ hoạt động". Đây là một ví dụ. (Thói quen tích hợp này là khá icky, do đó tên.)

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

Các hãy cẩn thận chỉ được rằng thói quen hội nhập số cần phải chơi tốt với AD, và cũng có thể chấp nhận lập luận phức tạp.Một cái gì đó thậm chí ngây thơ hơn, như

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

là piecewise liên tục trong giới hạn trên của hội nhập, đòi hỏi phải có các giới hạn của hội nhập là Enum và do đó không phức tạp, và cũng thường đánh giá tích phân bên ngoài phạm vi nhất định vì điều này :

Prelude> last [0..9.5] 
10.0 
Các vấn đề liên quan