# NDSolve error in 2-D heat equation

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).$$

• Several issues I can spot at the moment: 1. What's the meaning of $\Delta_\Gamma$? Are you sure the b.c. is translated correctly? 2. Where's the b.c. at $r=0$? 3. Where's the b.c. in $z$ direction? Feb 4, 2019 at 2:42
• 1. Please explain what's Laplace-Beltrami operator. 2. You do need boundary condition for $r=0$ and in $z$ direction. Notice when analytically solve PDE in polar coordinate, though we don't use explicit b.c. for $r=0$, we always make use of implicit constraint such as "the solution is bounded in the whole domain", and you need to translate such implicit constraint to explicit boundary condition when implement FDM and its friends. Needless to say you also need b.c. for $z$, given it's the angle of the polar coordinate, the b.c. is usually periodic b.c.. Feb 4, 2019 at 9:55
• The obvious solution is u[t,r,z]=1. Therefore, we can put u[t,0,z]==1 and add periodic boundary conditions u[t,r,0]==u[t,r,2*Pi]. Feb 4, 2019 at 13:01
• @Alex Reading the follow up comments of OP, the problem seems to be circularly symmetric. If so, a possible b.c. for $r=0$ is $\frac{\partial u}{\partial r}\big|_{r=0}=0$. Deduction: Expand[(r #1 &) /@ With[{u = u[r]}, D[u, t] == Laplacian[u, {r, theta}, "Polar"]]] /. r -> 0(* 0 == Derivative[1][u][0] *). Feb 6, 2019 at 3:06
• @Alex, I think I get your idea, and it seems to be true following this answer math.stackexchange.com/a/1794380/515871 Feb 6, 2019 at 17:11

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.

• Thank you. Yes your'e right. I'm confused about adding periodic bc condition because in the equation with Cartesian coordinates this bc don't arise. And since polar coordinates yield equivalent equation it should be same bc fixed in the beginning. I wish that I explained the problem well. May be I need time to understand many steps. Mar 10, 2019 at 15:27
• @S.Cho You may understand the periodic b.c. in the following way: the property of polar coordinate should be relfected in the solving process. What property does it own? Periodicity i.e. $u(t, r, \theta)=u(t,r,\theta+2\pi),\frac{\partial u(t, r, \theta)}{\partial \theta}=\frac{\partial u(t, r, \theta+2\pi)}{\partial \theta},…$. Mar 10, 2019 at 15:38
• We have already discussed this at mathematica.stackexchange.com/questions/190360/… Mar 10, 2019 at 16:25
• @AlexTrounev Nope, OP was trying to solve a 1+1D problem in Cartesian coordinate in that question, while this one is a 2+1D problem in polar coordinate. As to rebuild, it's included in the recent update of pdetoode/pdetoae. Mar 10, 2019 at 16:35
• @Migalobe The page is still available in way back machine. (The images are broken, but still readable, I've added it to my answer. ) In short, this formula is a generalization of $$\begin{array}{r@{}l} \nabla^2 U & {}= \displaystyle\frac{u_1 + u_2 + u_3 + u_4 - 4u_0}{\left( \Delta r \right)^2} \\[4mm] &{}= \displaystyle 4 \left[ \frac{u_1 + u_2 + u_3 + u_4}{4} - u_0 \right] \bigg/ \left( \Delta r \right)^2 \end{array}$$ (The page doesn't contain a strict proof, but I think the explanation is plausible. ) Feb 7 at 6:27

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


• @S.Cho Show your solution in Cartesian coordinates using pdetoode. Mar 9, 2019 at 20:47
• @S.Cho See update, now the solution converges quickly to u[t,r,z]=1 and v[t,z]=1. Mar 10, 2019 at 0:04
• Well, with all due respect, though u[t, r0, z] == 1` leads to the correct result in this case, it's just an opportunistic choice, because genrally nobody knows the value at the origin in polar coordinate. One possible general treatment for $r=0$ can be found here: homepages.see.leeds.ac.uk/~amt6xw/Distance%20Learning/CFD5030/… (Particularly the part around formula (1.6). ) Mar 10, 2019 at 12:46
• Yeah, but generally we don't know the exact solution. Mar 10, 2019 at 13:25
• I don't think that's a good choice, the value at $r=0$ generally changes over time. The proper b.c. should be deduced from the truth "the solution is smooth at $r=0$". (Not sure if this is the only possible way though. ) Mar 10, 2019 at 14:06