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,

enter image description here

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? — xzczd, Feb 04 '19 at 02:42
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 04 '19 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] *). — xzczd, Feb 06 '19 at 03:06
@Alex, I think I get your idea, and it seems to be true following this answer https://math.stackexchange.com/a/1794380/515871 — S. Maths, Feb 06 '19 at 17:11
Hello @xzczd and Alex. It have been a long time and my question still unaswered. I updated my post with the precise quation based on Alex coments. Thanks for any help. — S. Maths, Mar 09 '19 at 12:39

xzczd · Accepted Answer · 2022-02-07T06:14:26.937

2

Seems that OP still has difficulty in understanding why

Periodic b.c. should be added in $\theta$ direction.
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}]

Mathematica graphics

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

edited Feb 07 '22 at 06:14

answered Mar 10 '19 at 14:27

xzczd

65,995
9
163
468

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 '19 at 15:27
1

@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 '19 at 15:38
1

We have already discussed this at https://mathematica.stackexchange.com/questions/190360/boundary-conditions-involving-time-derivative/190438#190438 – Alex Trounev Mar 10 '19 at 16:25
@xzczd I can not find the definition of the function rebuild[] – Alex Trounev Mar 10 '19 at 16:28
1

@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 '19 at 16:35
@xzczd Can you give a link to the update? – Alex Trounev Mar 10 '19 at 16:44
@Alex It's just there :) . – xzczd Mar 10 '19 at 16:47
@xzczd Thank you! When executing the code, the system displays a message NDSolveValue::chknic: Structural analysis indicates that 601 initial conditions are needed to fix the state of the system. Currently only 575 initial conditions are specified. NDSolve may return one of a family of solutions. – Alex Trounev Mar 10 '19 at 16:56
@Alex Oops, there should be no delbothside@ in the definition of odeic, fixed. Thx for pointing out. (BTW this warning doesn't exist in v9, I believe it's added attempting to resolve the problem mentioned here and here, but sadly the problem is still pending. ) – xzczd Mar 10 '19 at 17:09
@xzczd thank you, now the code works. – Alex Trounev Mar 10 '19 at 17:23
@xzczd what is the idea behind $\nabla^2 u = \frac{4 \left( u_m - u_0 \right)}{\left( \Delta r \right)^2}$? (the page you mentioned is not working). – Migalobe Feb 06 '22 at 21:10
1

@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. ) – xzczd Feb 07 '22 at 06:27
@xzczd many thanks – Migalobe Feb 07 '22 at 10:45

Alex Trounev · Answer 2 · 2019-03-10T00:02:10.893

1

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

edited Mar 10 '19 at 00:02

answered Mar 09 '19 at 18:21

Alex Trounev

44,369
3
48
106

Nice idea but still not accurate enough. I think the use of pdetoode by xzczd should give better results. Is it necessary to add periodic BCs ? Because the equation in Cartesian coordinates don't contain this extra assumption and it is well-posed. Why we need to add this BCs ? – S. Maths Mar 09 '19 at 20:43
1

@S.Cho Show your solution in Cartesian coordinates using pdetoode. – Alex Trounev Mar 09 '19 at 20:47
It is difficult to solve the pde in Cartesian coordinates, because the Laplace-Beltrami and normal derivative are not easy to calculate. They are defined in the framework of differential geometry, in which case we use parametrisation and the simple case is to work with polar coordinates. – S. Maths Mar 09 '19 at 20:55
@S.Cho If the equation is solved in polar coordinates, then the solution on the z coordinate is, by definition, periodic. – Alex Trounev Mar 09 '19 at 21:32
1

@S.Cho See update, now the solution converges quickly to u[t,r,z]=1 and v[t,z]=1. – Alex Trounev Mar 10 '19 at 00:04
1

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: http://homepages.see.leeds.ac.uk/~amt6xw/Distance%20Learning/CFD5030/node10.html (Particularly the part around formula (1.6). ) – xzczd Mar 10 '19 at 12:46
@xzczd In this case, we know that the initial data u[0,r,z]=1 is the exact solution u[t,r,z]=1. Therefore we can put u[t,r0,z]=1. – Alex Trounev Mar 10 '19 at 13:08
1

Yeah, but generally we don't know the exact solution. – xzczd Mar 10 '19 at 13:25
@xzczd We always know the initial data, so we can use u[t,r0,z]==u[0,r0,z] if no other offers. – Alex Trounev Mar 10 '19 at 13:37
1

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 '19 at 14:06

NDSolve error in 2-D heat equation

2 Answers2

Linked