NDSolve error in 2-D heat equation

Question

I'm trying to solve the following PDE by Mathematica in 2-D case in the unit disk using polar cordinates,

where $\Omega$ is a bounded domain of $\mathbb{R}^n$, $\Gamma =\partial \Omega$ is the boundary of $\Omega$, $\partial_\nu$ is the normal derivative, and $\nu$ is the outer unit vector.

Thanks to pdetoode proposed by @xzczd I simplified the problem as follows:

eq = With[{u = u[t, r, z]}, D[u, t] == Laplacian[u, {r, z}, "Polar"]];
ic = u[0, r, z] == 1;
bc = With[{u = u[t, r, z]}, (eq[[1]] == eq[[2]] - D[u, r]) /. r -> 1];
tend = 1;
domain@r = {2 10^-6, 1};
domain@z = {0, 2 Pi};
points@r = points@z = 25;
difforder = 4;
(grid@# = Array[# &, points@#, domain@#]) & /@ {r, z};
ptoofunc = pdetoode[u[t, r, z], t, grid /@ {r, z}, difforder];
delbothside = #[[1 ;; -1]] &;
ode = delbothside /@ delbothside@ptoofunc@eq;
odeic = delbothside /@ delbothside@ptoofunc@ic;
odebc = MapAt[delbothside, ptoofunc@bc, {{1}, {2}}];`

When integrating the resulting ODEs with NDSolve :

sollst = NDSolveValue[{ode, odeic, odebc}, 
 Outer[u, grid@r, grid@z], {t, 0, tend}];

It generates the following error

NDSolveValue::overdet: There are fewer dependent variables
 {u[1/500000,0][t],u[1/500000,\[Pi]/12][t],u[1/500000,\[Pi]/6][t],u[
 /500000,\[Pi]/4][t],u[1/500000,\[Pi]/3][t],u[1/500000,(5 \[Pi])/12
 [t],u[1/500000,\[Pi]/2][t],u[1/500000,(7 \[Pi])/12][t],u[1/500000,(2
 \[Pi])/3][t],u[1/500000,(3 \[Pi])/4][t],u[1/500000,(5 \[Pi])/6][t],u[
 /500000,(11 \[Pi])/12][t],u[1/500000,\[Pi]][t],u[1/500000,(13 \[Pi])/12
 [t],u[1/500000,(7 \[Pi])/6][t],<<22>>,u[500023/12000000,\[Pi]
 [t],u[500023/12000000,(13 \[Pi])/12][t],u[500023/12000000,(7 \[Pi])/6
 [t],u[500023/12000000,(5 \[Pi])/4][t],u[500023/12000000,(4 \[Pi])/3
 [t],u[500023/12000000,(17 \[Pi])/12][t],u[500023/12000000,(3 \[Pi])/2
 [t],u[500023/12000000,(19 \[Pi])/12][t],u[500023/12000000,(5 \[Pi])/3
 [t],u[500023/12000000,(7 \[Pi])/4][t],u[500023/12000000,(11 \[Pi])/6
 [t],u[500023/12000000,(23 \[Pi])/12][t],u[500023/12000000,2 \[Pi]
 [t],<<575>>}, than equations, so the system is overdetermined.

I think the problem is in the part when we delete some ODEs to make place for BCs equations. I lost many days to solve this, but still unsolved.

Update:

To be more precise, the equation to solve in polar cordinates $(r,\theta)$ is

$$u_t=\frac{d^2 u}{d r^2}+\frac{1}{r} \frac{d u}{d r}+\frac{1}{r^2} \frac{d^2 u}{d \theta^2} \quad in \quad \Omega'=[0,1)\times [0,2 \pi)$$ $$u_t|_{r=1}=(\frac{d^2 u}{d \theta^2}-\frac{d u}{d r}) \Big|_{r=1}, \quad and \quad u|_{\theta=0}=u|_{\theta=2 \pi}=1, \quad (BC)$$ $$u(0,r,\theta)=1 \; in \; [0,1)\times [0,2 \pi), \quad u(0,1,\theta)=1 \; on \quad [0,2 \pi) \quad (IC).$$

Answer 1

Seems that OP still has difficulty in understanding why

1. Periodic b.c. should be added in $\theta$ direction.

2. Certain b.c. should be added at $r=0$.

If so, I suggest OP have a look at this answer first, which solves Laplace equation based on FDM from scratch i.e. without tools like pdetoode.

Anyway, the following is the fixed code.

eq = With[{u = u[t, r, z]}, r^2 D[u, t] == r^2 Laplacian[u, {r, z}, "Polar"] // Simplify];
ic = u[0, r, z] == 1;
bc = With[{u = u[t, r, z]}, (eq[[1]] == eq[[2]] - D[u, r]) /. r -> 1];
tend = 1;
domain@r = {0, 1};
domain@z = {0, 2 Pi};
points@r = points@z = 25;
difforder = 4;
(grid@# = Array[# &, points@#, domain@#]) & /@ {r, z};

Notice we can impose periodic b.c. by setting 5th argument of pdetoode:

ptoofunc = pdetoode[u[t, r, z], t, grid /@ {r, z}, difforder, {False, True}];
delbothside = #[[2 ;; -2]] &;
ode = delbothside@ptoofunc@eq;
odeic = ptoofunc@ic;
odebc = ptoofunc@bc;

At this point, we have (points@r - 2) points@z + points@z == 600 equations at hand, while points@r points@z == 625 unknown variables to be solved. Where can we find another points@z == 25 equations? Given the problem is defined in polar coordinate, we know u[0, …][t] are the same:

bcorigin = Equal @@@ Partition[u[0, #][t] & /@ grid@z, 2, 1];

This gives us points@z - 1 equations. Where can we find the final one? There're many possible ways, I'll simply utilize the idea mentioned in this page i.e. approximating $\nabla^2 u$ at $r=0$ as

$$\nabla^2 u = \frac{4 \left( u_m - u_0 \right)}{\left( \Delta r \right)^2}$$

where $u_m$ is the mean value of $u$ along $r = \Delta r$:

bcadditional = 
  With[{dr = grid@r // Differences // First}, 
   With[{um = u[dr, #][t] & /@ grid@z // Mean}, (4 (um - u[0, 0][t]))/dr^2 == 
     D[u[0, 0][t], t]]];

sollst = NDSolveValue[{ode, odeic, odebc, bcorigin, bcadditional}, 
   Outer[u, grid@r, grid@z], {t, 0, tend}];

solfunc = rebuild[sollst, grid /@ {r, z}]

Plot3D[solfunc[tend, r, z], {r, 0, 1}, {z, 0, 2 Pi}, PlotRange -> {0, 2}]

The result is trivial as expected, but the solution above should work for more general cases.

Answer 2

This problem can be solved by iterations. Here, the function u[t,r,z] is restored by the value of the function v[t,z], and the function v[t,z] by the value Derivative[0, 1, 0][u][t, 1, z]. In this example, the solution converges quickly to u[t,r,z]=1 and v[t,z]=1(I changed bc at r=r0=10^-6).

U[0][t_, r_, z_] := 1;
V[0][t_, z_] := 1; n = 5; r0 = 10^-6;
Do[U[i] = 
   NDSolveValue[{-(D[u[t, r, z], {z, 2}]/r^2) - D[u[t, r, z], r]/r - 
       D[u[t, r, z], {r, 2}] + D[u[t, r, z], t] == 0, u[0, r, z] == 1,
      u[t, r, 0] == u[t, r, 2*Pi],u[t, r0, z] == 1, 
     u[t, 1, z] == V[i - 1][t, z]}, 
    u, {t, 0, 1}, {r, r0, 1}, {z, 0, 2*Pi}]; 
  V[i] = NDSolveValue[{-D[v[t, z], {z, 2}] + D[v[t, z], t] + 
       Derivative[0, 1, 0][U[i]][t, 1, z] == 0, v[0, z] == 1, 
     v[t, 0] == v[t, 2*Pi]}, v, {t, 0, 1}, {z, 0, 2*Pi}];, {i, 1, 
   n}] 


Table[Plot3D[U[i][1, r, z], {r, r0, 1}, {z, 0, 2*Pi}, 
  PlotRange -> All], {i, 1, n}]

Table[Plot3D[V[i][t, z], {t, 0, 1}, {z, 0, 2*Pi}, 
  PlotRange -> All], {i, 1, n}]