Bug introduced in 10 or earlier and fixed in version 10.4

I am really not sure if this is a bug or I am missing something very trivial.

QUESTION: What I am missing in order to obtain the partial derivative of the interpolating function (in this case the fem solution)?

EDIT: take a look at the comment by user21 directly underneath this question, the solution is to use t as the first variable.


I am testing the following pde with the Mathematica FEM

xmin = 0;
xmax = 10;
reg = ImplicitRegion[xmin <= x <= xmax, {x}];
(*Field equation*)
sig[x_, t_] := Ef[x, t]*eps[x, t];
eps[x_, t_] := D[u[x, t], x];
feq = D[Af[x, t]*sig[x, t], x] + nf[x, t];
E0 = 3;
Ef[x_, t_] := E0;
A0 = 7;
Af[x_, t_] := A0;
sig0 = 2;
nf[x_, t_] := 0;
F[t_] := A0*sig0*3*t;
iconds = u[x, 0] == 0;
bconds = {DirichletCondition[u[x, t] == 0, x == xmin]};
nconds = NeumannValue[-F[t], x == xmax];
(*Fem solution*)
ufem = NDSolveValue[{feq == nconds, iconds, bconds}, u, 
   Element[x, reg], {t, 0, 1}];
Plot[Table[ufem[x, ti], {ti, 0, 1, 0.2}], {x, xmin, xmax}, 
 AxesLabel -> {"x", "u(x,t)"}]


From mechanics, the solution looks fine, it is linear in space $x$ and time $t$, just as it should be in this case. Just for fun I wanted to compute the spatial derivative of the solution. But I dont get a constant spatial derivative, as you can see in the plots below, e.g., for $t=0.2$.

Plot[ufem[x, 0.2], {x, xmin, xmax}, AxesLabel -> {"x", "ufem(x,0.2)"}]
loc = Derivative[1, 0][ufem];
Plot[loc[x, 0.2], {x, xmin, xmax}, AxesLabel -> {"x", "loc(x,0.2)"}]


From the InterpolatingFunction documentation, I assumed that partial derivates are computable, see screenshot below.

QUESTION: What I am missing in order to obtain the partial derivative of the interpolating function (in this case the fem solution)?


  • 4
    $\begingroup$ It works if you use t as the first indepdent variable in your equations. This is a bug. $\endgroup$
    – user21
    Commented Jan 25, 2016 at 18:06
  • $\begingroup$ Wow, thanks, that's kind of weird to have to use t as the first variable. But that is a weirdness I can live with :D Thanks! $\endgroup$ Commented Jan 25, 2016 at 19:18
  • 1
    $\begingroup$ This is now fixed in the sources and will be available with the next release after 10.3.1. In the meanwhile put the time variable in the first position. $\endgroup$
    – user21
    Commented Jan 27, 2016 at 17:40
  • $\begingroup$ @user21 For completeness: was it a bug in FEM or in InterpolatingFunction? Or both? Introduced in 10 or earlier? $\endgroup$ Commented Jan 28, 2016 at 10:53
  • 1
    $\begingroup$ @AlexeyPopkov, it was a bug in computing the derivative of an interpolating function when the time variable was last and the interpolating function is based on an element mesh. I did not check when it started. $\endgroup$
    – user21
    Commented Jan 28, 2016 at 17:10

1 Answer 1


Partial derivatives should work on interpolation functions just fine, and I cannot see why your interpolating function should be any different. But it is. Here is a ridiculous workaround

ufem2 = Interpolation[
   Flatten[Table[{{x, t}, ufem[x, t]}, {x, 0, 10, .1}, {t, 0, 
      1, .02}], 1]];
loc2 = Derivative[1, 0][ufem2];
Plot[{ufem[x, .2], loc2[x, 0.2]}, {x, xmin, xmax}, 
 AxesLabel -> {"x", "loc(x,0.2)"}]

enter image description here

This seems like a bug to me, that you get completely the wrong derivative from the original interpolating function.

  • $\begingroup$ Yeah, it is really driving my nuts. I did not want to generate data from the fem solution, since I thought/assumed that the interpolating function from the fem conserves the elements order at their positions and the derivative reflects this information. Do you happen to know if this is true? For linear elements, piecewise constant derivative. Generating data and a sequential interpolation would kill this information, right? $\endgroup$ Commented Jan 25, 2016 at 9:00
  • $\begingroup$ What you say is true, the order of the arguments should be preserved, and reflected by the Derivative command. Look at this code above, the "new" interpolating function ufem2 should produce identical results to ufem, and take the same arguments in the same order. But it behaves properly for the derivative, while ufem does not. $\endgroup$
    – Jason B.
    Commented Jan 25, 2016 at 9:02
  • $\begingroup$ I understand the new interpolation, but assume that you have just 3 points in the fem with linear elements, i.e., the derivative in the first and second element would be constant. A sequential interpolation of these 3 points (assuming we keep this resolution) would give me a polynomial of second order and the derivate would be a line. I would only be able to get the "same" information by a sequential interpolation with higher resolution than the fem, right? $\endgroup$ Commented Jan 25, 2016 at 9:09
  • $\begingroup$ @MauricioLobos - not sure what you mean by 3 points - I sampled the same region that ufem is defined over, with $0\leq x \leq 10$ and $0 \leq t \leq 1$, taking 50 data points along each coordinate. So with 2500 data points to build an interpolating function on, you should be able to get decent derivatives. $\endgroup$
    – Jason B.
    Commented Jan 25, 2016 at 9:18
  • $\begingroup$ To be clear I know absolutely nothing about finite element methods. All I'm demonstrating here is that if you have an interpolating function that behaves well in that it gives values within its specified range, but misbehaves in that it won't give derivatives, you can build an identical interpolating function by taking all the values of the function on a rectangular grid - the resolution of this grid is up to you, your function ufem should give values for any number in its range $\endgroup$
    – Jason B.
    Commented Jan 25, 2016 at 9:19

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.