# 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[] == eq[] - 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? – xzczd Feb 4 at 2:42
• $\Delta_\Gamma$ is the Laplace-Beltrami operator on the boundary $\Gamma$ (as Riemannian manifold) which is in this case the unit cercle (corresponding to $r=1$ in polar cordinates). So I think we don't need any other BCs to translate the Bcs in the equation. – S. Maths Feb 4 at 8: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.. – xzczd Feb 4 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]. – Alex Trounev Feb 4 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[u] *). – xzczd Feb 6 at 3:06

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[] == eq[] - 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. – S. Maths Mar 10 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},…$. – xzczd Mar 10 at 15:38
• We have already discussed this at mathematica.stackexchange.com/questions/190360/… – Alex Trounev Mar 10 at 16:25
• @xzczd I can not find the definition of the function rebuild[] – Alex Trounev Mar 10 at 16:28
• @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. – xzczd Mar 10 at 16:35

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[t_, r_, z_] := 1;
V[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. – Alex Trounev Mar 9 at 20:47
• @S.Cho See update, now the solution converges quickly to u[t,r,z]=1 and v[t,z]=1. – Alex Trounev Mar 10 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). ) – xzczd Mar 10 at 12:46
• Yeah, but generally we don't know the exact solution. – xzczd Mar 10 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. ) – xzczd Mar 10 at 14:06