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}}
  • 1
    $\begingroup$ 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]). $\endgroup$
    – bill s
    Commented Aug 15, 2012 at 3:46
  • $\begingroup$ Thanks, I've fixed those mistakes now but I'm still not getting a real answer. $\endgroup$
    – dshepherd
    Commented Aug 15, 2012 at 10:45
  • $\begingroup$ @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. $\endgroup$
    – user21
    Commented Aug 15, 2012 at 14:21

1 Answer 1


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}]
  • $\begingroup$ Of course, these days, one loads <<VectorAnalysis` instead of the <<Calculus` context... $\endgroup$ Commented Aug 21, 2012 at 20:21
  • $\begingroup$ No point reminding me how old I am! :-) $\endgroup$
    – chris
    Commented Aug 21, 2012 at 20:47
  • $\begingroup$ Don't worry about it, it took me a while to get used to the new way they're doing add-on packages myself... :) $\endgroup$ Commented Aug 21, 2012 at 20:51
  • $\begingroup$ 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 :) $\endgroup$
    – dshepherd
    Commented Sep 12, 2012 at 17:43

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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