# NDSolve::ntdvdae: Cannot solve to find an explicit formula for the derivatives

I am facing these error messages,

NDSolve::ntdvdae: Cannot solve to find an explicit formula for the derivatives. NDSolve will try solving the system as differential-algebraic equations.

plus

Power::infy: Infinite expression 1/0. encountered.

The equation in question is,

PDE = D[M[t, x, y], t] == 1/(x^2 + y^2)*D[(x^2 + 1)*D[M[t, x, y], x], x] +
1/(x^2 + y^2)*D[(-y^2 + 1)*D[M[t, x, y], y], y]

nv1 = NeumannValue[0, y == 1];
nv2 = NeumannValue[0, y == 0];
nv3 = NeumannValue[0, x == 0];
nv4 = NeumannValue[-Sqrt[(x^2 + y^2)/(x^2 + 1)]*Bi*M[t, x, y], x == xf];

Bi = 0.5; xf = 5;

sol = NDSolve[{PDE == nv1 + nv2 + nv3 + nv4, M[0, x, y] == 1},
M, {x, 0, xf}, {y, 0, 1}, {t, 0, 10}];


• what does this pde represent physically? I am asking because having all Neumann BC. and no Dirichlet boundary conditions looks very strange. Most of the time this leads to no unique solution (depending on the physics of the pde) Sep 29, 2019 at 3:01
• @Nasser I took this equation and bcs from here dx.doi.org/10.1590/S0104-66322008000100004
– zhk
Sep 29, 2019 at 3:06

First of all, PDE == nv1 + nv2 + nv3 + nv4 is obviously wrong, because there already exists a == in your PDE. This is easy to fix of course.

What's confusing is the Power::infy warning. I'm not sure why this pops up, maybe NDSolve fails to notice FiniteElement should be chosen in this case, while it should be able to. Anyway, specifying the method explicitly fixes the problem:

solution = NDSolve[{Subtract @@ PDE == nv4, M[0, x, y] == 1},
M, {x, 0, xf}, {y, 0, 1}, {t, 0, 10},
Method -> {"MethodOfLines", "SpatialDiscretization" -> "FiniteElement"}];

Table[ContourPlot[M[t, x, y] /. solution, {x, 0, xf}, {y, 0, 1}, Contours -> 20,
ColorFunction -> "TemperatureMap", PlotLegends -> Automatic,
PlotLabel -> Row[{"t = ", t}]], {t, 1, 10, 3}]


nv1, nv2, nv3 are omitted because zero Neumann value is the default setting of FiniteElement.

• Yes, this seems like a short coming in the detection. Hm, I'll need to dig into the code to see what happens - or better what does not happen in this case... Sep 30, 2019 at 12:47

Your equation is not in the correct form for FEM. Homogeneous Neumann conditions are applied automatically.

Needs["NDSolveFEM"]
Bi = 0.5; xf = 5; reg = Rectangle[{0, 0}, {xf, 1}];
mesh = ToElementMesh[reg, MaxCellMeasure -> 0.0001];
PDE = D[M[t, x, y], t] - 1/(x^2 + y^2)*D[(x^2 + 1)*D[M[t, x, y], x], x] -
1/(x^2 + y^2)*D[(-y^2 + 1)*D[M[t, x, y], y], y]

nv4 = NeumannValue[-Sqrt[(x^2 + y^2)/(x^2 + 1)]*Bi*M[t, x, y],
x == xf];

sol = NDSolve[{PDE == nv4, M[0, x, y] == 1},
M, {x, y} \[Element] mesh, {t, 0, 10}]

Table[ContourPlot[M[t, x, y] /. sol, {x, y} \[Element] mesh,
Contours -> 20, ColorFunction -> "TemperatureMap",
PlotLegends -> Automatic, PlotLabel -> Row[{"t = ", t}]], {t, 1, 10,
3}]


Consider a solution in a region with a hole (homogeneous Neumann conditions are applied automatically at the edge of the hole)

Needs["NDSolveFEM"]
Bi = 0.5; xf = 5; reg = Rectangle[{0, 0}, {xf, 1}]; reg1 =
Rectangle[{1, 1/3}, {4, 2/3}];
mesh = ToElementMesh[RegionDifference[reg, reg1],
MaxCellMeasure -> 0.0001];
PDE = D[M[t, x, y], t] -
1/(x^2 + y^2)*D[(x^2 + 1)*D[M[t, x, y], x], x] -
1/(x^2 + y^2)*D[(-y^2 + 1)*D[M[t, x, y], y], y];

nv4 = NeumannValue[-Sqrt[(x^2 + y^2)/(x^2 + 1)]*Bi*M[t, x, y],
x == xf];

sol = NDSolve[{PDE == nv4, M[0, x, y] == 1},
M, {x, y} \[Element] mesh, {t, 0, 10}]

Table[ContourPlot[M[t, x, y] /. sol, {x, y} \[Element] mesh,
Contours -> 20, ColorFunction -> "TemperatureMap",
PlotLegends -> Automatic, PlotLabel -> Row[{"t = ", t}]], {t, 1, 10,
3}]


• This PDE isn't equivalent to the original one, you forgot to change the nv4 accordingly. Sep 29, 2019 at 13:15
• @xzczd you're right, corrected, matches your answer. Sep 29, 2019 at 15:10
• @AlexTrounev How to construct a prolate or oblate region inside a square and then mesh around it?
– zhk
Oct 1, 2019 at 4:46
• @zhk Are you interested in solving this equation on such a mesh or in general? Oct 1, 2019 at 7:50
• @AlexTrounev yes I want to solve this equation on such a mesh
– zhk
Oct 1, 2019 at 8:11