# Bug? Problem with derivative of interpolating function from FEM

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.

Context:

I am testing the following pde with the Mathematica FEM

Needs["NDSolveFEM"]
(*Domain*)
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];
(*Paremeters*)
E0 = 3;
Ef[x_, t_] := E0;
A0 = 7;
Af[x_, t_] := A0;
sig0 = 2;
(*Inhomogeneity*)
nf[x_, t_] := 0;
F[t_] := A0*sig0*3*t;
(*Conditions*)
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*)
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)?

• It works if you use t as the first indepdent variable in your equations. This is a bug. – user21 Jan 25 '16 at 18:06
• 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! – Mauricio Fernández Jan 25 '16 at 19:18
• 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. – user21 Jan 27 '16 at 17:40
• @user21 For completeness: was it a bug in FEM or in InterpolatingFunction? Or both? Introduced in 10 or earlier? – Alexey Popkov Jan 28 '16 at 10:53
• @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. – user21 Jan 28 '16 at 17:10

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)"}]


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

• 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? – Mauricio Fernández Jan 25 '16 at 9:00
• 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. – Jason B. Jan 25 '16 at 9:02
• 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? – Mauricio Fernández Jan 25 '16 at 9:09
• @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. – Jason B. Jan 25 '16 at 9:18
• 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 – Jason B. Jan 25 '16 at 9:19