Derivatives of a multivariant function and the corresponding Jacobian with a vector-space package

Question

I have a problem with the vector-space package. I got a very useful answer from @mnish in a recent post, but there I was dealing only with a function that depends on only one variable. What happens when I have, for example, a function that displays from polar coordinates to Cartesians

f:(0,oo) x [0,2pi] -> R² (r,phi) -> (r*cos(phi),r*sin(phi)) 

which depends on two variables.

I tried this using a rather naive approach:

 polar :: Double -> Double -> ((Double,Double) :~> (Double,Double)) polar r phi = \(r,phi) -> (((idD) r)*cos( idD phi),((idD) r)*sin( idD phi)) 

I get the following error:

 Couldn't match expected type `(Double, Double) :> (Double, Double)' with actual type `(t0, t1)' In the expression: (((idD) r) * cos (idD phi), ((idD) r) * sin (idD phi)) In the expression: \ (r, phi) -> (((idD) r) * cos (idD phi), ((idD) r) * sin (idD phi)) In an equation for `polar': polar r phi = \ (r, phi) -> (((idD) r) * cos (idD phi), ((idD) r) * sin (idD phi)) 

For one component

 polarx :: Double -> Double -> ((Double,Double) :~> Double) polarx r phi = \(r,phi) -> ((idD) r)*cos( idD phi) 

I get

 Couldn't match expected type `Double' with actual type `(Double, Double)' Expected type: (Double, Double) :> Double Actual type: (Double, Double) :> (Double, Double) In the return type of a call of `idD' In the first argument of `(*)', namely `((idD) r)' 

Apparently, there is a certain type of disorder, but I cannot understand what is wrong.

Another question arises when I want to calculate the Jacobian of such a mapping. As the name implies, it is somehow connected with line cards, which, of course, are covered by the package, in fact it is based on these cards. But then again, my knowledge of Haskell is not enough to get a solution on my own.

Answer 1

I finally found a solution to my problem, it was not so difficult, but still it took me a while to figure it out. In case anyone else is interested, I present the details.

First, here is my code for the polar case:

 polarCoordD :: ((Double,Double) :~> (Double,Double)) polarCoordD = \(r,phi) -> pairD (polarx (r,phi), polary (r,phi)) where polarx :: (Double,Double) :~> Double polarx = \(r,phi) -> (fst . unpairD $ (idD) (r,phi))*cos( snd . unpairD $ idD (r, phi)) polary :: (Double,Double) :~> Double polary = \(r,phi) -> (fst . unpairD $ (idD) (r,phi))*sin( snd . unpairD $ idD (r, phi)) 

The key was to make the "differentiation variable" (idD) aware of the tuple (r, phi) , which contains two variables that I want to differentiate. Then I have to unzip the tuple via unpairD and select the first and second part of the resulting pair (in polarx and polary ). Both are packaged again in pairs. Maybe there is a more elegant way to do this, but the way I finally understood it.

From here it is easy to pass to cylindrical coordinates or, in fact, to any other curved system of orthogonal coordinates. For cylindrical coordinates we get:

 cylCoordD :: (Vec3 Double :~> Vec3 Double) cylCoordD = \(r,phi,z) -> tripleD (cylx (r,phi,z), cyly (r,phi,z),cylz (0,0,z)) where cylx :: (Double,Double,Double) :~> Double cylx = \(r,phi,z) -> (fst' . untripleD $ (idD) (r,phi,z))*cos( snd' . untripleD $ idD (r, phi,z)) cyly :: (Double,Double,Double) :~> Double cyly = \(r,phi,z) -> (fst' . untripleD $ (idD) (r,phi,z))*sin( snd' . untripleD $ idD (r, phi,z)) cylz :: (Double,Double,Double) :~> Double cylz = \(_,_,z) -> third . untripleD $ idD (0,0,z) fst' :: (a,b,c) -> a fst' (x,_,_) = x snd' :: (a,b,c) -> b snd' (_,y,_) = y third :: (a,b,c) -> c third (_,_,z) = z 

where Vec3 Double belongs to type Vec3 a = (a, a, a) . Now we can even build the transformation matrix:

 let transmat = \(r,phi,z) -> powVal $ liftD3 (,,) (normalized $ derivAtBasis (cylCoordD (r,phi,z)) (Left ())) (normalized $ derivAtBasis (cylCoordD (r,phi,z)) (Right (Left ()))) (normalized $ derivAtBasis (cylCoordD (r,phi,z)) (Right (Right ()))) *Main> transmat (2, rad 0, 0) ((1.0,0.0,0.0),(0.0,1.0,0.0),(0.0,0.0,1.0)) *Main> transmat (2, rad 90, 0) ((6.123233995736766e-17,1.0,0.0),(-1.0,6.123233995736766e-17,0.0),(0.0,0.0,1.0)) 

rad is a convenience feature

 rad :: Double -> Double rad = (pi*) . (*recip 180) 

Now it would be interesting to convert this “matrix” into the matrix type Numeric Prelude and / or hmatrix , but I'm not sure if that would be even useful. But still, this would be a good example of using a vector-space package.

I still need to figure out the use and especially the use of line cards.

Answer 2

Just saw this next question. I'm not sure what you want:

• Jacobi matrix
• product of Jacobi vector
• product of Jacobian-transposed vectors

In such a low-dimensional system, I will take the first. (The rest come in handy mainly when the system is high enough that you do not want to store or calculate the Jacobian in your own way, but instead treat it like a generalized sparse matrix.) In any case:

 Prelude> :m + Numeric.AD Prelude Numeric.AD> let f [r,phi] = map (r*) [cos phi, sin phi] Prelude Numeric.AD> jacobian f [2,3::Float] [[-0.9899925,-0.28224],[0.14112,-1.979985]]