# Multiply list of pure functions by list of scalars?

I am trying to write a function that returns a list of pure functions. The last step in the function, is to multiply the list of pure functions by a list of scalars (i.e. matrix multiplication). When I try to evaluate the resulting functions, it does not work.

Code

calcShapeFunctions[nnpe_] := Module[
{m, xiCoord, etaCoord, p, c},
m = Array[
Function[{xi,
eta}, {1, xi, eta, xi*eta, xi^2, eta^2, xi^2*eta, xi*eta^2,
xi^2*eta^2}[[#]]] &, nnpe];
xiCoord = {-1, 1, 1, -1, 0, 1, 0, -1, 0};
etaCoord = {-1, -1, 1, 1, -1, 0, 1, 0, 0};
p = m[[1 ;; nnpe]];
c = Array[Through[p[xiCoord[[#]], etaCoord[[#]]]] &, nnpe];
Return[p.Inverse[c]];
];


Incorrect Output

Through[calcShapeFunctions[3][x, y]]


returns:

{(1/2 Function[{xi, eta}, {1, xi, eta, xi eta, xi^2, eta^2, xi^2 eta, xi eta^2, xi^2 eta^2}[[1]]] - 1/2 Function[{xi, eta}, {1, xi, eta, xi eta, xi^2, eta^2, xi^2 eta, xi eta^2, xi^2 eta^2}[[2]]])[x, y], (1/2 Function[{xi, eta}, {1, xi, eta, xi eta, xi^2, eta^2, xi^2 eta, xi eta^2, xi^2 eta^2}[[2]]] - 1/2 Function[{xi, eta}, {1, xi, eta, xi eta, xi^2, eta^2, xi^2 eta, xi eta^2, xi^2 eta^2}[[3]]])[x, y], (1/2 Function[{xi, eta}, {1, xi, eta, xi eta, xi^2, eta^2, xi^2 eta, xi eta^2, xi^2 eta^2}[[1]]] + 1/2 Function[{xi, eta}, {1, xi, eta, xi eta, xi^2, eta^2, xi^2 eta, xi eta^2, xi^2 eta^2}[[3]]])[x, y]}

Desired Output

Through[calcShapeFunctions[3][x, y]]


should return:

{1/2 - x/2, x/2 - y/2, 1/2 + y/2}

Further Thoughts

Assuming we can get this to work, how would I go about obtaining the derivatives of the pure function with respect to one of its independent variables (xi and eta)?

As with the last question on which I believe this one is based, it may be simpler to use a single function rather than a list of functions.

To understand my code you will need to know that:

I chose to use Slot (# and #2) because it makes this manipulation much simpler, though it would be possible with named parameters with more work.

Other notes:

• I made output of the calc function to be a single function as well, simplifying both the code and the application of it
• I eliminated p as it appeared to be a copy of m
• With is use to evaluate certain parts of the output Function

The code:

calcFn2[nnpe_] :=
Module[{m, xiCoord, etaCoord, p, c, xi, eta},
m = {1, #, #2, # #2, #^2, #2^2, #^2 #2, # #2^2, #^2 #2^2} &;
m = m[[{1}, ;; nnpe]];
xiCoord = {-1, 1, 1, -1, 0, 1, 0, -1, 0};
etaCoord = {-1, -1, 1, 1, -1, 0, 1, 0, 0};
c = Array[m[xiCoord[[#]], etaCoord[[#]]] &, nnpe];
With[{ic = Inverse[c], mev = m}, mev[##].ic &]
]

calcFn2[3][x, y]

{1/2 - x/2, x/2 - y/2, 1/2 + y/2}

• This is a perfect solution, exactly what I was looking for. My only question is could you explain what's going on when you call m = m[[{1}, ;; nnpe]];? I'm confused on how Part is working there. Apr 5 '13 at 1:03
• More perspicacious than mine. But for Derivative to work on xi and eta (#1 and #2), won't you need something like Evaluate[mev[#1, #2].ic] & in the body of With? Apr 5 '13 at 1:09
• @MichaelE2 D[calcFn2[3][x, y], y] does what I need it to do. Apr 5 '13 at 1:13
• @Matthew OK, if that's what you want. D[.., y] calculates the derivative of the composition with respect to y, not the derivative of calcFn2[3] with respect to eta per se though, which would be Derivative[0,1][calcFn2[3]]. For instance, Derivative[0,1][calcFn2[3]][1,2] or Derivative[0,1][calcFn2[3]][x-y,2y] could not be done by simply substituting D[.., y]. Of course, in each case, it's possible to make D work one way or another. Apr 5 '13 at 1:24
• @Matthew I added an explanation to my answer - too long for a comment. Apr 9 '13 at 3:04

First, to evaluate the functions embedded in calcShapeFunctions[3], you could do the following:

calcShapeFunctions[3] /. f_Function :> f[x, y]

{1/2 - x/2, x/2 - y/2, 1/2 + y/2}


But to find the derivatives, you need to alter how calcShapeFunctions function is defined so that the bodies of your Functions are evaluated and the part of the list of formulas is extracted (line 4):

calcShapeFunctions[nnpe_] :=
Module[{m, xiCoord, etaCoord, (*p,*) c},
m = Array[Function[{xi, eta},
Evaluate@{1, xi, eta, xi*eta, xi^2, eta^2, xi^2*eta, xi*eta^2, xi^2*eta^2}[[#]]] &,
nnpe];
xiCoord = {-1, 1, 1, -1, 0, 1, 0, -1, 0};
etaCoord = {-1, -1, 1, 1, -1, 0, 1, 0, 0};
(*p = m[[1 ;; nnpe]];*) (* same as p = m *)
c = Array[Through[m[xiCoord[[#]], etaCoord[[#]]]] &, nnpe];
Return[m.Inverse[c]];];


Then the derivatives can be done as follows:

calcShapeFunctions[3] /. f_Function :> Derivative[1, 0][f][x, y]

{-(1/2), 1/2, 0}

calcShapeFunctions[3] /. f_Function :> Derivative[0, 1][f][x, y]

{0, -(1/2), 1/2}


(Note: In calcShapeFunctions, the variable p was the same as m, so I replaced it.)

Response to comment

In response to the OP's request as to how I would modify Mr. Wizard's answer, I offer the following. As I suggested in a comment, I would Evaluate the function body inside With:

calcFn2[nnpe_] :=
Module[{m, xiCoord, etaCoord, p, c, xi, eta},
m = {1, #, #2, # #2, #^2, #2^2, #^2 #2, # #2^2, #^2 #2^2} &;
m = m[[{1}, ;; nnpe]];
xiCoord = {-1, 1, 1, -1, 0, 1, 0, -1, 0};
etaCoord = {-1, -1, 1, 1, -1, 0, 1, 0, 0};
c = Array[m[xiCoord[[#]], etaCoord[[#]]] &, nnpe];
With[{ic = Inverse[c], mev = m}, Evaluate[mev[#1, #2].ic] &]]

calcFn2[4]
(* {1/4 - #1/4 - #2/4 + (#1 #2)/4, 1/4 + #1/4 - #2/4 - (#1 #2)/4,
1/4 + #1/4 + #2/4 + (#1 #2)/4, 1/4 - #1/4 + #2/4 - (#1 #2)/4} & *)


One can take the Derivative of the Function:

Derivative[0, 1][calcFn2[4]]
(* {-(1/4) + #1/4, -(1/4) - #1/4, 1/4 + #1/4, 1/4 - #1/4} & *)


And plug in expressions:

Derivative[0, 1][calcFn2[4]][x, y]
(* {-(1/4) + x/4, -(1/4) - x/4, 1/4 + x/4, 1/4 - x/4} *)


Or one can apply D to the expression calcFn2[3][x, y]:

D[calcFn2[4][x, y], y]
(* {-(1/4) + x/4, -(1/4) - x/4, 1/4 + x/4, 1/4 - x/4} *)


You could define your own operator that maps a scalar and a function to a function i.e. $\otimes: \mathbb{R}\times\frak{F} \to \frak{F}$

CircleTimes[c_, f_] := (c f[##] &)


and use it like c Escc*Esc f to get:

f[x_]:=x
10 ⊗ f
(* 10 f[##] & *)


And you can Thread over lists:

constants = RandomReal[1, 3];
functions = {Sin[#1]Cos[#2]&, Cos, Exp};
(* {0.845413 (Sin[#1] Cos[#2] &)[##1] &, 0.835235 Cos[##1] &, 0.103793 Exp[##1] &} *)


You can easily get derivatives with D by naming the arguments:

D[result[[1]][x, y], x]
(* 0.845413 Cos[x] Cos[y] *)

D[result[[2]][x], x]
(* -0.835235 Sin[x] *)


Derivative is a little more picky, it doesn't handle ## (SlotSequence) like one would expect:

Derivative[1, 0][(Sin[#1] Cos[#2] &)[##1] &][x, y]
(* 0 *)

Derivative[1, 0][(Sin[#1] Cos[#2] &)[#1, #2] &][x, y]
(* Cos[x] Cos[y] *)


I'd go with D that symbolically evaluates the nested functions to avoid further surprises.