Come differenziare gli integrali con la libreria di vettori (haskell)

Question

Quando si utilizza il pacchetto vector-space per le torri derivative (vedere derivative towers), ho riscontrato la necessità di differenziare gli integrali. dalla matematica è del tutto chiaro come raggiungere questo obiettivo:

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

con una funzione

g : R -> R 

per esempio.

La derivata rispetto a x sarebbe:

f'(x) = g(x) 

Ho cercato di ottenere questo comportamento definendo prima classe "integrazione"

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

un'istanza di base è

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

con integrateQAGS da hmatrix

il problema viene con valori b che rappresentano torri di derivati:

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

NC.T da Numeric.Complex (numerico-prelude). La funzione integrateD è definito come segue (ma sbagliato):

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) 

La funzione non restituisce quello che voglio, deriva l'integrando, ma non l'integrale. Il problema è che ho bisogno di una mappa lineare che restituisca f u. Il a :> b è definito come segue:

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

Non so come definire derivative. Qualsiasi aiuto sarà apprezzato, grazie

edit:

ho dimenticato di fornire l'istanza per 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 

e posso dare un esempio di quello che voglio dire: Diciamo che ho un funzione

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

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

(AR. *) significa moltiplicazione da Algebra.Anello (numerico-preludio)

posso facilmente integrare questa funzione con la funzione di cui sopra integrateD:

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

Quando prendo uno sguardo alla derivata di f:

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

e valutare questa a 0 e pi/2 ottengo 1 e qualche valore:

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

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

Ora, quando deriva l'integrale, ottengo la derivazione della funzione f non il suo valore al limite superiore

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

ma mi aspetto:

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

Answer 1

Ho finalmente trovato una soluzione alla mia domanda. La chiave della soluzione è la funzione >-< dal pacchetto spazio vettoriale, che rappresenta la regola della catena.

Così, ho definire una funzione integrateD' come questo:

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) 

il d_one significa una variabile di derivazione e la sua derivata deve essere 1. Con questa funzione ora posso creare alcuni casi come

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

e

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

sfortunatamente Non riesco a utilizzare integrateD con valori complessi fuori dalla scatola, devo usare liftD2. Il motivo per questo sembra essere la funzione idD, non so se c'è una soluzione più elegante.

Quando guardo l'esempio nella questione, io ora ottenere la mia soluzione desiderata:

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

o utilizzando l'istanza:

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

Answer 2

'hmatrix' è legata strettamente a Doppio. Non è possibile utilizzare le sue funzioni con altri tipi di dati numerici come quelli forniti da 'vector-space' o 'ad'.

Answer 3

Se si vuole solo fare dC su una funzione che implica l'integrazione numerica, senza che il sistema di AD conosca l'integrazione per-se, dovrebbe "funzionare". Ecco un esempio. (Questa routine integrazione è piuttosto icky, da cui il nome.)

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

Gli unici distinguo sono che la routine di integrazione numerica ha bisogno di giocare bene con AD, ed accettare anche argomenti complessi.Qualcosa di ancora più naif, come

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

è costante a tratti nel limite superiore di integrazione, richiede i limiti di integrazione essere Enum e quindi non complessa, e spesso anche di valutare l'integrando esterno all'intervallo dato a causa di questo :

Prelude> last [0..9.5] 
10.0