Mathematica Stack Exchange is a question and answer site for users of Mathematica. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

I'm attempting to calculate symbolic expressions for the (Newton's method) Jacobian of a Galerkin finite element method using Mathematica. It would seem that this is a perfect application for symbolic differentiation tools but I can't get it working (probably due to my lack of Mathematica knowledge).

Mathematical statement of the problem

Let $v(\mathbf{x})$ and $\psi_k(\mathbf{x})$ ($k=1,\ldots,N$) be simple(ish) polynomial functions of space. Let $\mathbf{m}(\mathbf{x}) = [m_1(\mathbf{x}),m_2(\mathbf{x}),m_3(\mathbf{x})] $ be a vector field.

Let $$ \mathbf{m}(\mathbf{x}) = \sum_{k=1}^N \psi_k(\mathbf{x}) \mathbf{m}_k, $$ and $$ r(\mathbf{x}) = \int_V (\mathbf{m} \times \nabla^2 \mathbf{m}) v(\mathbf{x}) dV. $$

(I have used $\nabla^2$ as the vector Laplacian $ \nabla^2 \mathbf{m} = [ \nabla^2 m_1, \nabla^2 m_2, \nabla^2 m_3]. $)

Calculate $\frac{\partial r}{\partial m_{i,k}}$ for all $i$ (directions of $\mathbf{m}$), and all $k \in {1,\ldots,N}$ (index for the $\psi$ functions).

My code

I'm sure this can be written much more elegantly but as I said before I'm not too great with Mathematica. Anyway it should be good enough to give you the idea.

psi[x_, y_, z_] := {a[x, y, z], b[x, y, z]};
(* Really we just need to declare that this is a vector somehow...*)

m[x_, y_, z_] := {m1x*psi[x, y, z][[1]] + m2x*psi[x, y, z][[2]],
  m1y*psi[x, y, z][[1]] + m2y*psi[x, y, z][[2]],
  m1z*psi[x, y, z][[1]] + m2z*psi[x, y, z][[2]]}

VectorLaplacian[f[x,y,z]_] := {
  D[f[x,y,z][[1]], {x, 2}] + D[f[x,y,z][[1]], {y, 2}] + D[f[x,y,z][[1]], {z, 2}],
  D[f[x,y,z][[2]], {x, 2}] + D[f[x,y,z][[2]], {y, 2}] + D[f[x,y,z][[2]], {z, 2}],
  D[f[x,y,z][[3]], {x, 2}] + D[f[x,y,z][[3]], {y, 2}] + D[f[x,y,z][[3]], {z, 2}]

r[x_, y_, z_] := (Cross[m[x, y, z], VectorLaplacian[m[x, y, z]]])*v[x,y,z]
 (* I've left out the integration for now to keep things simpler *)

J[x_, y_, z_] := {
  {D[r[x,y,z][[1]], m1x], D[r[x,y,z][[2]], m1x], D[r[x,y,z][[3]], m1x]},
  {D[r[x,y,z][[1]], m1y], D[r[x,y,z][[2]], m1y], D[r[x,y,z][[3]], m1y]},
  {D[r[x,y,z][[1]], m1z], D[r[x,y,z][[2]], m1z], D[r[x,y,z][[3]], m1z]},

  {D[r[x,y,z][[1]], m2x], D[r[x,y,z][[2]], m2x], D[r[x,y,z][[3]], m2x]},
  {D[r[x,y,z][[1]], m2y], D[r[x,y,z][[2]], m2y], D[r[x,y,z][[3]], m2y]},
  {D[r[x,y,z][[1]], m2z], D[r[x,y,z][[2]], m2z], D[r[x,y,z][[3]], m2z]}

J[1, 1, 1]

But all I get is:

{{0, 0, 0}, {0, 0, 0}, {0, 0, 0}, {0, 0, 0}, {0, 0, 0}, {0, 0, 0}}
share|improve this question
There are some syntax issues here. If you change all the terms like psi[[1]][x, y, z] to psi[x, y, z][[1]] then at least it makes sense (this would be a[x,y,z]) in your example. Just making this change in the m[] definition gives values for r, i.e., you can ask for r[x,y,z] and you get something fairly large. The J function has a similar issue, perhaps m[[1]] should be m[x,y,z][[1]] (presumably you mean the first element of the polynomial m[x,y,z]). – bill s Aug 15 '12 at 3:46
Thanks, I've fixed those mistakes now but I'm still not getting a real answer. – dshepherd Aug 15 '12 at 10:45
@dshepherd, note that the definiton of your VectorLaplacian is not quite right. You may want to define a function and then test it to see if it works. e.g. evaluate VectorLaplacian[m[x, y, z]] to see that it does what you want. – user21 Aug 15 '12 at 14:21
up vote 2 down vote accepted

If you want to do this for a specific N it is straightforward.

Let us load the calculus package.


For readability let us use the formatting rules

 Format[Subscript[\[Psi], a__][x, y, z]] = Subscript[\[Psi], a];
  Format[\[Nu][x, y, z]] = \[Nu];

Use cartesian coordinate for Laplacian and Cross product

  SetCoordinates[Cartesian[x, y, z]]

Let us define

 psi[k_][x_, y_, z_] = {Subscript[\[Psi], k, 1][x, y, z], Subscript[\[Psi], k, 2][x, y, z] , 
                     Subscript[\[Psi], k, 3][x, y, z] };

Now let's say we have N=2 for the sake of this example; Define

 m[x_, y_, z_] = Sum[psi[k][x, y, z]*{Subscript[m, 1, k], Subscript[m, 2, k], 
                 Subscript[m, 3, k]}, {k, 1, 2}];

Defining the integrant:

 rint = \[Nu][x, y, z] Cross[m[x, y, z], Map[Laplacian, m[x, y, z]]];

Now mm is the set of variables:

 mm = 
 Table[{Subscript[m, 1, k], Subscript[m, 2, k], Subscript[m, 3, 
 k]}, {k, 2}]

This table will be your 'Jacobian'

  jac=Map[D[rint, #] &, mm, {2}]; 
  jac // TableForm

If you want an expression for arbitrary N its also possible.


You might want to consider the case where the fields Subscript[m,i,k] depends on {x,y,z} as well.

One should then add the formatting rule

Format[Subscript[m, a__][x, y, z]] = Subscript[m, a];


m[x_, y_, z_] = Sum[psi[k][x, y, z]*{Subscript[m, 1, k][x, y, z], 
              Subscript[m, 2, k][x, y, z], Subscript[m, 3, k][x, y, z]}, {k, 1, 2}]

And deal with derivatives by integration by part (subject to proper boundary conditions...); e.g.

Subscript[\[Psi], 1, 1][x, y, z] D[D[Subscript[m, 1, 1][x, y, z], {x,2}], 
Subscript[m, 1, 1][x, y, z]] -> D[Subscript[\[Psi], 1, 1][x, y, z], {x,2}]
share|improve this answer
Of course, these days, one loads <<VectorAnalysis` instead of the <<Calculus` context... – J. M. Aug 21 '12 at 20:21
No point reminding me how old I am! :-) – chris Aug 21 '12 at 20:47
Don't worry about it, it took me a while to get used to the new way they're doing add-on packages myself... :) – J. M. Aug 21 '12 at 20:51
Thanks, I actually ended up rechecking my solutions by hand and then forgot about this question. But it will definitely be useful later on so thank you very much :) – dshepherd Sep 12 '12 at 17:43

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

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