# Why can't I use Inactive[D] in NDSolve?

I have a 2D rectangular domain with -L/2 < x < L/2 and 0 < y < D0, which is divided into two parts, with a sinusoidal curve y == 1 + δ Cos[2 Pi x/L] separating the two.

I want to solve Laplace equation D[κ[x, y] D[T[x, y], x], x] + D[κ[x, y] D[T[x, y], y], y] == 0, where κ[x, y] takes different values in the two regions (equals 1 below the curve, and 0.2422 above the curve). After copying some codes in StackExchange, I have the following codes which I think give the right answer:

D0 = 3;
k = 0.2422;
L = 3;
δ = 0.2;
sol = NDSolveValue[{Inactive[
Div][-If[y <= 1 - δ Cos[2 π x/L], 1,
k].Inactive[Grad][T[x, y], {x, y}], {x, y}] == 0,
PeriodicBoundaryCondition[T[x, y], x == -L/2,
TranslationTransform[{L, 0}]],
DirichletCondition[T[x, y] == 1 + δ Cos[2 π x/L],
y == 0 && -L/2 < x < L/2],
DirichletCondition[T[x, y] == 0, y == D0 && -L/2 < x < L/2]},
T, {x, -L/2, L/2}, {y, 0, D0},
Method -> {"FiniteElement", "InterpolationOrder" -> {T -> 2},
"MeshOptions" -> {"MaxCellMeasure" -> 0.0005},
"IntegrationOrder" -> 5}];
Show[ContourPlot[sol[x, y], {x, -L/2, L/2}, {y, 0, D0},
PlotLegends -> Automatic, Contours -> 20],
Plot[1 - δ Cos[2 π x/L], {x, -L/2, L/2}]]


Now I want to solve a slightly different equation D[κ[x, y] D[T[x, y], x], x] + epsilon^2 D[κ[x, y] D[T[x, y], y], y] == 0. This is basically Laplace equation where the two spatial coordinates are scaled differently, and I want to study the effect of ignoring the second term with epsilon^2 in my analytical work.

I try to define the equation differently as Inactive[D][-If[y <= 1-δ Cos[2 π x/L], 1, k] Inactive[D][T[x, y], {x}], {x}] + (epsilon^2) Inactive[D][-If[y <= 1-δ Cos[2 π x/L], 1, k] Inactive[D][T[x, y], {y}], {y}] == 0

But it gives an error NDSolveValue::derivs: No derivatives of dependent variables were found in the equations. NDSolveValue is designed to solve differential or differential algebraic equations. Use NSolve or FindRoot to numerically solve algebraic equations.

What's the problem? In fact, I don't think I fully understand the use of Inactive. What's the best way to define the PDE (or the domain) in this case?

• Where's the epsilon in the new PDE? Have you read the "Formal Partial Differential Equations" section in the document page FEMDocumentation/tutorial/FiniteElementOverview? Mar 27, 2020 at 5:18
• @xzczd just edited. I basically don't understand why the way I write the equation (splitting it into the x and y derivatives instead of using grad and div) doesn't work. Mar 27, 2020 at 5:49
• I think this is a good question. My guess is Inactive is only for building the "formal PDE", but it's not explicitly mentioned in the document. Let's wait for user21's clarification :) . Mar 27, 2020 at 5:51

The coefficient you give for the diffusion operator is a scalar but you can give a matrix of coefficients. I think this is what you are looking for:

D0 = 3;
k = 0.2422;
L = 3;
δ = 0.2;
epsilon = 10^-1;


Note how the coeff now returns a matrix of coefficients:

coeff = If[
y <= 1 - δ Cos[2 π x/L], -{{1, 0}, {0,
epsilon^2*1}}, -{{k, 0}, {0, epsilon^2*k}}];


Well just plug that in:

sol = NDSolveValue[{Inactive[Div][
coeff.Inactive[Grad][T[x, y], {x, y}], {x, y}] == 0,
PeriodicBoundaryCondition[T[x, y], x == -L/2,
TranslationTransform[{L, 0}]],
DirichletCondition[T[x, y] == 1 + δ Cos[2 π x/L],
y == 0 && -L/2 < x < L/2],
DirichletCondition[T[x, y] == 0, y == D0 && -L/2 < x < L/2]},
T, {x, -L/2, L/2}, {y, 0, D0},
Method -> {"FiniteElement",
"MeshOptions" -> {"MaxCellMeasure" -> 0.0005}}];


As a side note, for single equations you never need to specify "InterpolationOrder"; that is only relevant for coupled PDE where you want to interpolate some of the dependent variables not at the highest level.

The new in version 12.1 HeatTransfer tutorial has a section on isotropic, orthotropic and anisotropic diffusion that explain the concepts and implementation more.

• Any comment on the Inactive[D] part? Mar 27, 2020 at 11:21
• @xzczd, basically because it is not implemented. You should be able to do everything with Inactive[Div] and Inactive[Grad], no? Mar 27, 2020 at 11:58
• Yeah, but Inactive[D] will make coding more flexible, right? :) I think the more important part is, this isn't explicitly documented. Do you mean currently only Inactive[Div] and Inactive[Grad] are supported by FiniteElement method of NDSolve? Mar 27, 2020 at 12:06
• @xzczd, yes only Inactive Div/Grad (for differential operators) are supported. And, yes, it might make sense to mention that in the documentation. Mar 27, 2020 at 12:42
• @xzczd, I tried to clarify this in the documentation and this should be in 12.1.1. Thanks. Mar 30, 2020 at 5:59