Jak rozróżnić całki z biblioteką przestrzeni wektorowej (haskell)

Question

Podczas korzystania z pakietu vector-space dla wież pochodnych (patrz derivative towers), natknąłem się na potrzebę różnicowania całek. Z matematyką jest całkiem jasne, jak to osiągnąć:

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

z funkcją

g : R -> R 

na przykład.

Pochodna względem x byłoby:

f'(x) = g(x) 

Próbowałem dostać ten problem, najpierw zdefiniować klasę "integracja"

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

podstawowym instancji jest

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

z integrateQAGS z hmatrix

problem jest związany z wartościami B, które stanowią wieże pochodnych:

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

NC.T wynosi od Numeric.Complex (numerycznych-wstęp). Funkcja integrateD jest zdefiniowana w następujący sposób (ale źle):

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) 

Funkcja nie zwraca to, co chcę, to wywodzi się całki, ale nie całki. Problem polega na tym, że potrzebuję liniowej mapy, która zwraca f u. a :> b jest zdefiniowany następująco:

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

nie wiem jak zdefiniować derivative. Każda pomoc będzie mile widziane, dzięki

edit:

zapomniałem podać instancję 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 

i mogę dać przykład, co mam na myśli: Powiedzmy mam funkcja

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

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

(AR. *) oznacza mnożenie z algebry.Pierścień (numeryczna-wstępem)

można łatwo zintegrować z tej funkcji w wyżej funkcji integrateD:

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

Gdy spojrzeć na pochodnej f:

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

i oceny tej pod numerami: 0 i pi/2 Otrzymuję 1 i niektóre wartości:

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

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

Teraz, kiedy pochodzący całki, otrzymuję wyprowadzenie funkcji f nie jej wartość w górnej granicy

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

Ale spodziewam:

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

Answer 1

W końcu znalazłem rozwiązanie mojego pytania. Kluczem do rozwiązania jest funkcja >-< z pakietu vector-space, która oznacza regułę łańcucha.

Więc zdefiniować funkcję integrateD' tak:

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 oznacza zmienną pochodzenie i jego pochodna, musi być 1. Za pomocą tej funkcji można teraz tworzyć pewne instancje jak

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

i

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

Niestety Nie mogę użyć integrateD na skomplikowanych wartościach po wyjęciu z pudełka, muszę użyć liftD2. Powodem tego wydaje się być funkcja idD, nie wiem, czy istnieje bardziej eleganckie rozwiązanie.

Kiedy patrzę na przykład w pytaniu, ja teraz dostać mój pożądany rozwiązanie:

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

lub używając instancję:

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

Answer 2

„hmatrix” jest związana bardzo ściśle Podwójnie. Nie możesz korzystać z jego funkcji z innymi numerycznymi typami danych, na przykład z "przestrzeni wektorowej" lub "reklamy".

Answer 3

Jeśli chcesz po prostu zrobić reklamę w funkcja polegająca na integracji numerycznej, bez znajomości systemu AD wiedząc o integracji per se, powinna "po prostu działać". Oto przykład. (Procedura ta integracja jest dość icky, stąd nazwa.)

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

Jedyne zastrzeżenia to, że rutynowe numeryczny integracja musi grać dobrze z AD, a także zaakceptować złożone argumenty.Coś jeszcze naiwne, jak

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

jest stała odcinkowo w górnej granicy całkowania wymaga ograniczenia integracji być wyliczeniowego oraz tym samym nie złożone i często ocenia całki poza podanym zakresie ze względu na to :

Prelude> last [0..9.5] 
10.0