Suppose I have an expression similar to this: y=23.23*(h[x]^2)*(D[h[x],x])*(D[h[x],{x,3}])/(13.2+Cos[0.245x]) The actual expression is a lot more complicated, but it still consists of the function h[x], its derivatives, and some common functions of x (like Sin[x], Cos[x], etc).

Now I have obtained the values of h[x] and the corresponding x through other programs, and I want to substitute x, h[x] and its derivatives (h[x] is reasonable smooth, so I numerically differentiate h[x] to obtain its derivatives) into y. If this is in Matlab or Python, then I basically have x and h[x] as arrays of length N, and I can just apply the function y(x,h,h2,...) to the array, and Matlab/Python(numpy) will vectorize everything. At the end, I will get y as an array with the same length N as my input arrays.

How do I achieve the same thing in Mathematica? I tried something like y/.{h[x]->u1, h'[x]->u2, h''[x]->u3, h'''[x]->u4, x->mesh], where u1,u2,u3,u4,mesh are lists of length N. However, I end up with a NxN matrix instead of a vector (or list) of length N. What's the correct way to do the 'element-wise' substitution with numerical data?

  • 2
    $\begingroup$ Your approach works for me when I generate sample data arrays. Give us a few lines of your sample data so we can all work on the same problem. I used 'mesh = Subdivide[1.0, 10]; u1 = Sin[mesh]; u2 = Cos[mesh]; u3 = -u1; u4 = -u2;` to generate my sample data. Are your arrays really simple one-dimensional lists of numeric data of the same length? $\endgroup$ – LouisB Jan 22 at 6:21

I do not exactly know what you mean. Maybe one or more of the symbols that you use in the replacement have already definitions attached to them?

If you prefer to do it by vectorization (which indeed should generally be prefered over ReplaceAll for performance reasons), you can do it like this:

Define a function f; Block takes care of shielding the symbols from potentially existing global definitions:

Block[{u1, u2, u3, u4, mesh, h, x},
  f[u1_, u2_, u3_, u4_, mesh_] = 
   y /. {h[x] -> u1, h'[x] -> u2, h''[x] -> u3, h'''[x] -> u4, x -> mesh}

Also notice the usage of Set and not of SetDelayed because this enforces that the replacement is performed only once. Basically

f[u1_, u2_, u3_, u4_, mesh_] = ...

is equivalent to

f[u1_, u2_, u3_, u4_, mesh_] := Evaluate[...]


f = {u1, u2, u3, u4, mesh}  \[Function] Evaluate[...]

(at least for most purposes).

Now we can generate data in vectors:

{U1, U2, U3, U4, M} = data = RandomReal[{-1, 1}, {5, 20}];

And f will tread over these vectors when they are used as arguments:

f[U1, U2, U3, U4, M] // Dimensions


However, the reason for vectorization is that the constitutents of f have the attribute Listable, not f itself. (So strictly speaking, f does not thread over lists.)

If you data lies as rows in a matrix data, you have to Apply f instead:

f @@ data // Dimensions


Again, this relies on the fact that all functions in the definition of f are listable. To be on the safe side, just do

MapThread[f,data] // Dimensions


| improve this answer | |

Not the answer you're looking for? Browse other questions tagged or ask your own question.