# 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) – Nasser Sep 29 '19 at 3:01
• @Nasser I took this equation and bcs from here dx.doi.org/10.1590/S0104-66322008000100004 – zhk Sep 29 '19 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... – user21 Sep 30 '19 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. – xzczd Sep 29 '19 at 13:15
• @xzczd you're right, corrected, matches your answer. – Alex Trounev Sep 29 '19 at 15:10
• @AlexTrounev How to construct a prolate or oblate region inside a square and then mesh around it? – zhk Oct 1 '19 at 4:46
• @zhk Are you interested in solving this equation on such a mesh or in general? – Alex Trounev Oct 1 '19 at 7:50
• @AlexTrounev yes I want to solve this equation on such a mesh – zhk Oct 1 '19 at 8:11