6
$\begingroup$

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}];

Any suggestion please?

$\endgroup$
2
  • $\begingroup$ 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) $\endgroup$
    – Nasser
    Sep 29, 2019 at 3:01
  • $\begingroup$ @Nasser I took this equation and bcs from here dx.doi.org/10.1590/S0104-66322008000100004 $\endgroup$
    – zhk
    Sep 29, 2019 at 3:06

2 Answers 2

5
$\begingroup$

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

enter image description here

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

$\endgroup$
1
  • 3
    $\begingroup$ 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... $\endgroup$
    – user21
    Sep 30, 2019 at 12:47
2
$\begingroup$

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

Needs["NDSolve`FEM`"]
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}]

Figure 1

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

Needs["NDSolve`FEM`"]
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}]

Figure 2

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.