# Unexpected 'switch' of partial derivatives when differentiating InterpolatingFunction

Bug introduced in 11.3 or earlier and persisting through 12.1.0 or later

In short, there is a simple initial boundary-value problem for which NDSolveValue produces an InterpolatingFunction f[x,y,t]. The issue is that Derivative[1,0,0][f] and Derivative[0,1,0][f] seem to be switched! By this I mean that for $$\partial_{x} f$$ one gets the values one expects for $$\partial_{y} f$$, and vice-versa.

It seems that the problem is related to the fact that the InterpolatingFunction is defined over a mesh for the finite element method.

I am using Mathematica 11.3.0.0 on Linux x86 (64 bit) platform (Kubuntu 18.04 on Dell Precision M4800).

# The initial boundary-value problem

The PDE giving rise to the issue is the following. Solve the heat equation $$\nabla^{2}T=\partial_{t}T$$ on the rectangle $$(x,\,y)\in [0, 1] \times [0, 2]$$ subject to the initial condition $$T(x,y,t=0)=1$$ and the boundary conditions that correspond to the state-state solution of $$T(x,y)=y/2$$. In other words, the boundary conditions are:

$$T=0$$ along the horizontal edge $$y=0$$;
$$T=1$$ along the horizontal edge $$y=2$$; and,
along the vertical edges ($$x=0$$ and $$x=1$$), we linearly connect the boundary conditions at the horizontal edges; in other words, $$T=y/2$$ along the vertical edges $$x=0$$ and $$x=1$$.

For the steady-state solution, we have $$\partial_{t}T=0$$, and so we are simply solving the Laplace's equation $$\nabla^{2}T=\partial^{2}_{x}\,T+\partial^{2}_{y}\,T=0$$ subject to the above boundary conditions. It is obvious that $$T(x,y)=y/2$$ satisfies both the boundary conditions and the Laplace's equation, and thus is the correct solution.

Let's first check that the steady-state solution works out:

Ω = Rectangle[{0, 0}, {1, 2}];
DSolveValue[{Laplacian[T[x, y], {x, y}] == 0, DirichletCondition[T[x, y] == y/2, True]},
T[x, y], {x, y} ∈ Ω]


And indeed, the output is y/2.

Now the full time-dependent solution:

tmin = 0;tmax = 10;
Ω = Rectangle[{0, 0}, {1, 2}];
f = NDSolveValue[{Laplacian[T[x, y, t], {x, y}] == D[T[x, y, t], t],
T[x, y, 0] == 1, DirichletCondition[T[x, y, t] == y/2, True]},
T, {x, y} ∈ Ω, {t, tmin, tmax}]


Here is the plot of the solution at tmax. Clearly, it corresponds well to the steady-state solution $$T(x,y)=y/2$$:

  Manipulate[
Plot3D[f[x, y, t], {x, 0, 1}, {y, 0, 2}, PlotRange -> {0, 1}],
{{t, tmax}, tmin, tmax, Appearance -> "Open" }]


So far, everything has been as expected.

# The issue

But now compute the partial derivatives in the $$x$$ and $$y$$-directions anywhere, say right in the middle of the region, $$(x,y)=(0.5,1)$$:

Derivative[1, 0, 0][f][0.5, 1, tmax]


The output should be zero, since $$T$$ has no $$x$$-dependence within the region. But in fact, the output is 0.499968 instead!

And if the compute Derivative[0, 1, 0][f][0.5, 1, tmax], we get the output 1.28073*10^-14.

This is opposite of what it should be: after all, at t=tmax, the solution is almost exactly g[x_, y_] = y/2; and Derivative[1, 0][g][0.5, 1] comes out as 0, while Derivative[0, 1][g][0.5, 1] comes out as 1/2, as it should.

One can also try this:

h[x_, y_] = f[x, y, tmax];
hx[x_, y_] = D[h[x, y], x];
hy[x_, y_] = D[h[x, y], y];


We still get that hx[0.5, 1] is 0.499968 while hy[0.5, 1] is 1.28073*10^-14.

And yet, if we try to compute these partial derivatives 'by hand',

hX[x_, y_] = (h[x + 0.01, y] - h[x, y])/0.01;
hY[x_, y_] = (h[x, y + 0.01] - h[x, y])/0.01;


we do get the expected values: hX[0.5, 1] gives -1.46022*10^-7, while hY[0.5, 1] gives 0.499968.

# The role of the mesh for the finite element method

The reason why I say that this issue seems to be related to the fact that the InterpolatingFunction is defined over the mesh for the finite element method is this:

On the one hand, so far, NDSolveValue is using the finite element method (this is automatic because I am both specifying the domain with ∈ and using DirichletCondition; see this question and this question). Indeed, if you execute

??f


you will see things like NDSolveFEMElementMesh, NDSolveFEMQuadElement, NDSolveFEMLineElement, and NDSolveFEMPointElement.

On the other hand, the issue does not arise if we specify the computational region and the boundary conditions in a way that does not trigger the use of the finite element method:

tmin = 0; tmax = 10;
v = NDSolveValue[{Laplacian[T[x, y, t], {x, y}] == D[T[x, y, t], t],
T[x, y, 0] == 1, T[x, 0, t] == 0, T[0, y, t] == y/2,
T[1, y, t] == y/2, T[x, 2, t] == 1},
T, {x, 0, 1}, {y, 0, 2}, {t, tmin, tmax}]


If you execute ??v, you'll see no sign of NDSolveFEM objects; instead, v uses DeveloperPackedArrayForm. And now the partial derivatives evaluate as expected: Derivative[1, 0, 0][v][0.5, 1, tmax] as -7.15028*10^-12, and Derivative[0, 1, 0][v][0.5, 1, tmax] as 0.499986.

However, the finite element method does not invariably lead to the 'switching' of the partial derivatives. For example, the issue does not arise if we numerically solve the steady-state problem where we both specify the domain with ∈ and use DirichletCondition:

Ω = Rectangle[{0, 0}, {1, 2}];
u = NDSolveValue[{Laplacian[T[x, y], {x, y}] == 0,
DirichletCondition[T[x, y] == y/2, True]},
T, {x, y} ∈ Ω]


Executing ??u, we again see NDSolveFEMElementMesh etc. However, Derivative[1, 0][u][0.5, 1] gives 1.80407*10^-14, while Derivative[0, 1][u][0.5, 1] gives 0.5, which are the expected results.

What is going on here? Why are the 'spatial' partial derivatives of the InterpolatingFunction f 'switched'?

More to the point, however: are there any strategies one can employ so that the issue doesn't arise? For example, it turns out (see below) that using T[t,x,y] instead of T[x,y,t] in the problem above makes the issue disappear. Is there some general principle at work here that we could employ to make sure some similar issue won't arise in other problems?

# EDIT 1

As guessed by user21 and first confirmed by andre314 (and then myself as well), the issue does not arise if instead of T[x,y,t] we use T[t,x,y]:

tmin = 0; tmax = 10;
Ω = Rectangle[{0, 0}, {1, 2}];
q = NDSolveValue[{Laplacian[T[t, x, y], {x, y}] == D[T[t, x, y], t],
T[0, x, y] == 1, DirichletCondition[T[t, x, y] == y/2, True]},
T, {x, y} ∈ Ω, {t, tmin, tmax}]


And now Derivative[0, 1, 0][q][tmax, 0.5, 1] evaluates to 1.28073*10^-14 while Derivative[0, 0, 1][q][tmax, 0.5, 1] evaluates to 0.499968, which are the expected values.

As andre314 points out, this doesn't actually answer the question of why the problem does appear if we use T[x,y,t]. But it is certainly very much worth knowing!

# EDIT 2

user21 confirms this is a bug:

Internally, the representation of interpolating functions is t first, then spatial coordinates and something seems to be going south with the t last version. Very likely during the computation of the derivatives. Now, since V12.1 will be release in the not too distant future a fix (which I still have to come up with) might not make it into 12.1. So it's best to use the t first version for the time being.

# EDIT 3

For completeness, here is an analytic solution of the problem. Since the equilibrium solution is known ($$T_{\text{eq.}}(x,y)=y/2$$), following the usual procedure, we set $$u(x,y,t)=T(x,y,t)-T_{\text{eq.}}(x,y)$$ and notice that $$u$$ is the solution the heat equation subject to homogeneous Dirichlet boundary conditions (i.e. $$u=0$$ on the boundary of the rectangle $$(x,\,y)\in [0, 1] \times [0, 2]$$), with the initial condition $$u(x,y,t=0)=1-y/2$$. This problem can be solved by expanding $$u$$ in terms of the eigenfunctions of the 2D Laplacian subject to the specified boundary conditions; once $$u$$ is found, the full solution is given by $$T(x,y,t)=T_{\text{eq.}}(x,y)+u(x,y,t)$$.

We make the ansatz that $$u$$ is a sum of 'modes' $$u_{\lambda}$$, where each mode satisfies the heat equation and the boundary conditions, but not necessarily the initial condition (to satisfy the initial condition we will need to sum the modes with appropriate prefactors). To proceed, we also make the ansatz that each mode can be factorized: $$u_{\lambda}(x,y,t)=X(x)Y(y)T(t)$$; we find that $$X''=a X$$, $$Y''=b Y$$, and $$T'=(a+b) T$$, where $$a$$ and $$b$$ are constants. These are eigenvalue problems (indeed, Sturm-Liouville problems); taking account of the boundary conditions, the orthonormal eigenfunctions of the $$X$$ equation are $$X_{m}(x)=\sqrt{2}\sin(\pi m x)$$ for $$m=1,\,2,\,\ldots$$, and those of the $$Y$$ equation are $$Y_{n}(y)=\sin(\frac{\pi}{2} n y)$$, $$n=1,\,2,\,\ldots$$. The corresponding eigenfunctions of the $$T$$ equation are $$T_{mn}(t)=c\exp[-\pi^{2}(m^{2}+n^{2}/4)\,t]$$, where $$c$$ is a constant. We will pick $$c=1$$ so that $$T(0)=1$$. Thus $$u(x,y,t)=\sum_{m,n=1}^{\infty}A_{mn}\,X_{m}(x)Y_{n}(y)T_{mn}(t)$$. The coefficients $$A_{mn}$$ are obtained by demanding that $$u(x,y,0)=1-y/2$$: $$A_{mn}=\left(\int_{0}^{1}\,X_{m}(x)\,dx\right)\left(\int_{0}^{2}\,(1-y/2)Y_{n}(y)\,dy\right)\,.$$ The final result is $$T(x,y,t)=\frac{y}{2}+\frac{8}{\pi^{2}}\sum_{\underset{\scriptstyle\text{m odd}}{m=1}}^{\infty}\sum_{n=1}^{\infty}\frac{1}{mn}\sin(m\pi x)\sin\left(\frac{n}{2}\pi y\right) e^{-\pi^{2}\left(m^{2}+\frac{n^{2}}{4}\right)\,t}\,.$$

• Plot[Evaluate[f[x, 1, tmax]], {x, 0, 1}] shows that it is not zero along x. Very small change exist of order 10^(-6). That is why the derivative is not zero. May be you need more accuracy in the numerical solution. btw, Did you try DSolve on this? – Nasser Mar 8 at 8:31
• @Nasser 1. It's not just that the partial derivative does not evaluate to strictly zero; it is the fact that it evaluates to almost exactly 0.5. That is not a matter of insufficient precision. At the same time, 0.5 is what the partial derivative in the y-direction should be evaluating to, but that partial derivative instead evaluates to almost exactly zero. As I said, they are 'switched'. – linguisticturn Mar 8 at 8:47
• @Nasser 2. I did try DSolve; it is returned unevaluated, even when I change the initial condition to the steady-state value T[x,y,0]==y/2. But I am not sure how the analytical solution, if it were available, would help as far as answering my question, which, after all, is about the numerical solution. We already know the steady-state solution exactly, and it is indeed reproduced by the numerical solution. – linguisticturn Mar 8 at 8:47
• Try to put t as a first independent variable: t, x, y in stead of x, y, t and see if the problem persists? – user21 Mar 8 at 12:20
• Thank you for reporting this. – user21 Mar 9 at 6:25

Internally, the FEM-generated interpolating function always stores t first. This means that an input where we have an independent variable ordering of x, y, t is internally reordered to t,x,y. There is a chance that something goes wrong during the derivative computation. To try this, change your ordering from x,y,t to t,x,y and see if this fixes the issue. Should that indeed work, then the issue you reported is a bug.

Update:

OK, the proper workaround for this is to put t first as in t,x,y. However I'd like to show another way to check this:

Compute the solution:

tmin = 0; tmax = 10;
Ω = Rectangle[{0, 0}, {1, 2}];
f = NDSolveValue[{Laplacian[T[x, y, t], {x, y}] == D[T[x, y, t], t],
T[x, y, 0] == 1, DirichletCondition[T[x, y, t] == y/2, True]},
T, {x, y} ∈ Ω, {t, tmin, tmax}];


Compute and evaluate the derivative:

fx = Derivative[1, 0, 0][f];
fx[0.5, 1, tmax]
0.4999678532261518


Inspect the interpolating data structure, as explained here:

fx[[2]]
{5, 12417, 1, {63, 1394, 0}, {4, 3, 3}, {0, 0, 1}, 0, 0, 0,
Indeterminate &, {}, {}, False}


Note the 6th position:

fx[[2, 6]]
{0, 0, 1}


This gives the derivative order that the interpolating function represents. When we know that the internal representation is t,x,y then {0,0,1} means fy and not fx. Let's change that:

fx[[2, 6]] = {0, 1, 0}
{0, 1, 0}


And....

fx[0.5, 1, tmax]
1.4206200904431552*^-14


But again: the proper workaround is to use the t,x,y ordering and not the x,y,t one. I found the issue in the code and fixed it. But I am hesitant to include this for 12.1 since the release is probably close and even though all tests pass fiddling with this so short of the release makes me nervous. Should there be a 12.1.1 (and that is not clear at this point in time) I'll consider adding it there. Else it will have to wait for 12.2. Since this has an easy fix I think it's not the end of the world to not have it in 12.1. However, I apologize for my mistake. Sorry.

Update: Should there be a 12.1.1 then the fix for this issue will be in it.

Note an answer; just too long for a comment

@user21 's suggestion (taking t,x,y instead of x,y,t as variables) works.
Here is the code :

Ω = Rectangle[{0, 0}, {1, 2}];
tmin = 0; tmax = 10;
Ω = Rectangle[{0, 0}, {1, 2}];
f = NDSolveValue[{
Laplacian[T[t, x, y], {x, y}] == D[T[t, x, y], t]
, T[0, x, y] == 1
, DirichletCondition[T[t, x, y] == y/2, True]
}
, T
, {x, y} ∈ Ω
, {t, tmin, tmax}
(*, Method\[Rule]{"MethodOfLines","TemporalVariable"\[Rule] t,
"SpatialDiscretization"\[Rule]{"FiniteElement"}}*)]

Plot3D[f[tmax, x, y], {x, 0, 1}, {y, 0, 2}, PlotRange -> {0, 1},
AxesLabel -> {"X", "Y"}]

Echo[Derivative[0, 1, 0][f][tmax, 0.5, 1],
"derivative along X axis : "];


• I think it's worth mentionning my progress on the subject because the OP may find it difficult to try things. The trial and error approach may be very time consuming for those who are not acustomed to the Mathematica + FEM system. – andre314 Mar 8 at 14:10
• .. and the OP has allready done a great job in the wording of the question. – andre314 Mar 8 at 14:24
• Thanks Andre. I'll file it tomorrow then. – user21 Mar 8 at 15:05
• @user21 By 'filing it tomorrow', do you mean reporting this as a bug to Wolfram? (This is a bug, correct?) – linguisticturn Mar 8 at 18:16
• @linguisticturn, yes, I meant to say, file a bug report. Internally, the representation of interpolating functions is t first, then spatial coordinates and something seems to be going south with the t last version. Very likely during the computation of the derivatives. Now, since V12.1 will be release in the not too distant future a fix (which I still have to come up with) might not make it into 12.1. So it's best to use the t first version for the time being. Sorry for the inconvenience. – user21 Mar 9 at 6:20