Multiply list of pure functions by list of scalars?

Question

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

Answer 1

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

Answer 2

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:

• You can manipulate a Function with Part just as you would any other expression
• Part always wraps element sequences with the original head of the expression

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}

Answer 3

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};
result = Thread[CircleTimes[constants, functions]]
(* {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.