1
$\begingroup$

I'm trying to take the partial differential of each component of an interpolation summed together. It works fine when I do it with an algebraic expression, but I can't understand why it doesn't work with the interpolation function.

Defining the magnetic field algebraically:

u = 1;

(*define magnetic field*) 

bfield[x_, y_, z_] = {y/1, (x*3^2)/1, 0};


(*calculate magnetic pressure*) 

pfield[x_, y_, z_] = {

-D[1/(2 u) (\[Sqrt]((bfield[x, y, z][[1]])^2 + (bfield[x, y, z][[
        2]])^2 + (bfield[x, y, z][[3]])^2))^2, x],

-D[1/(2 u) (\[Sqrt]((bfield[x, y, z][[1]])^2 + (bfield[x, y, z][[
        2]])^2 + (bfield[x, y, z][[3]])^2))^2, y],

-D[1/(2 u) (\[Sqrt]((bfield[x, y, z][[1]])^2 + (bfield[x, y, z][[
        2]])^2 + (bfield[x, y, z][[3]])^2))^2, z]};


VectorPlot3D[pfield[x, y, z], {x, -3, 3}, {y, -3, 3}, {z, -3, 3}]

StreamPlot[{pfield[x, y, z][[1]], pfield[x, y, z][[2]]}, {x, -3,3}, {y, -3, 3}]

But now I try the same thing with the interpolation function (I need to use the interpolation function as I will change the data eventually and thus it will have no simple algebraic function)

testdata = 

Flatten[Table[N@{x, y, z, bfield[x, y, z]}, {x, -3, 3 , 0.5}, {y, -3, 3 ,0.5}, ,{z, 0, 1, 0.5}], 2];

intf = Interpolation[testdata];

intfpressure[x_, y_, z_] := Module[{r = intf[x, y, z]}, {

(-D[1/(2 u) (\[Sqrt]((r[[1]])^2 + (r[[2]])^2 + (r[[3]])^2))^2, x]),

(-D[1/(2 u) (\[Sqrt]((r[[1]])^2 + (r[[2]])^2 + (r[[3]])^2))^2, y]),

(-D[1/(2 u) (\[Sqrt]((r[[1]])^2 + (r[[2]])^2 + (r[[3]])^2))^2, z])}];


VectorPlot3D[intfpressure[x, y, z], {x, -3, 3}, {y, -3, 3}, {z, 0, 1},

ViewPoint -> {0, 0, 200}, VectorPoints -> Fine, 

PlotLabel -> "pressure"]

StreamPlot[{intfpressure[x, y, 0][[1]], intfpressure[x, y, 0][[2]]}, {x, -3, 3}, {y, -3, 3}

And I get the correct fields from first code snippet:

enter image description here

enter image description here

Wrong from the second:

enter image description here

]

Extension of this question: Picking a dimension from an interpolating function

$\endgroup$
3
$\begingroup$

The problem with your second code block is that r[[1]], r[[2]] etc. are not the first, second etc. components of the vector field as you would expect.

If you input numerical values for x, y, z in intfpressure, e.g. intfpressure[0,0,0], then r is just a list of three numbers, so differentiating the elements of this list with respect to either of the coordinates will just give 0. On the other hand, if we call intfpressure with symbolic arguments, then r just stays unevaluated, just as if we had called intf[x,y,z]:enter image description here

This is actually what happens inside VectorPlot3D; Mathematica first tries to evaluate the expression to plot with symbolic arguments, and only later assigns numerical values to them in order to plot something. Now, taking the Parts, [[1]] or [[2]] etc., of such an expression as intf[x,y,z] is the same as taking the parts of any other expression with an arbitrary head:

dummyhead[x,y,z][[1]]
(* x *)

enter image description here

This explains your incorrect vector plot; the thing you are differentiating in intfpressure, Sqrt[r[[1]]^2 + r[[2]]^2 + r[[3]]^2]^2 is nothing but Sqrt[x^2 + y^2 + z^2]^2, which is of course not what you want.

My proposed fix is to use Derivative to produce the interpolations to the partial derivatives, since Derivative is, as the documentation says,

a functional operator which acts on functions to give derivative functions.

Here is the code. I've made a helper function intfnorm to show how one can use FunctionInterpolation to create new InterpolatingFunctions from other such objects, and to improve readability:

intfnorm = 
 FunctionInterpolation[
  intf[x, y, z].intf[x, y, z], {x, -3, 3}, {y, -3, 3}, {z, 0, 1}]

Clear[intfpressure]

intfpressure[x_, y_, 
  z_] :=
 (-1/(2 u)) {Derivative[1, 0, 0][intfnorm][x, y, z],
   Derivative[0, 1, 0][intfnorm][x, y, z],
   Derivative[0, 0, 1][intfnorm][x, y, z]}

StreamPlot[{intfpressure[x, y, 0][[1]], 
  intfpressure[x, y, 0][[2]]}, {x, -3, 3}, {y, -3, 3}]

enter image description here

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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