# Using periodic boundary conditions in NDSolve

I need to solve a 2D linear elliptical equation in polar coordinates using NDSolve, but I can’t seem to get Mathematica to accept the required periodic boundary conditions. My actual differential equation has complicated coefficients, but the same problem appears if I try to solve the Laplace equation, as indicated below. What am I doing wrong?

s=NDSolve[{r^2*D[V[r,phi],r,r]+r*D[V[r,phi],r]+D[V[r,phi],phi,phi]==0,V[1,phi]==0,
V[2,phi]==Sin[phi],V[r,2*Pi]==V[r,0],Derivative[0,1][V][r,2*Pi]
==Derivative[0,1][V][r,0]},V,{r,1,2},{phi,0,2*Pi}]

(* NDSolve[{(V^(0,2))[r,phi]+r (V^(1,0))[r,phi]+r^2 (V^(2,0))[r,phi]==0,V[1,phi]==0,
V[2,phi]==Sin[phi],V[r,2 π]==V[r,0],(V^(0,1))[r,2 π]
==(V^(0,1))[r,0]},V,{r,1,2},{phi,0,2 π}] *)

Evaluate[V[1.5,2]/.s]


This last line yields

ReplaceAll::reps: {NDSolve[{(V^(0,2))[r,phi]+r (V^(1,0))[r,phi]+r^2 (V^(2,0))[r,phi]==0,V[1,phi]==0,V[2,phi]==Sin[phi],V[r,2 π]==V[r,0],(V^(0,1))[r,2 π]==(V^(0,1))[r,0]},V,{r,1,2},{phi,0,2 π}]} is neither a list of replacement rules nor a valid dispatch table, and so cannot be used for replacing. >>

(* V[1.5,2]/. NDSolve[{(V^(0,2))[r,phi]+r (V^(1,0))[r,phi]+r^2 (V^(2,0))[r,phi]==0,
V[1,phi]==0,V[2,phi]==Sin[phi],V[r,2 π]==V[r,0],(V^(0,1))[r,2 π]==
(V^(0,1))[r,0]},V,{r,1,2},{phi,0,2 π}] *)

• Welcome to Mathematica.SE! I suggest the following: 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Read the faq! 3) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! – Michael E2 Jun 21 '15 at 1:37
• You can format inline code and code blocks by selecting the code and clicking the {} button above the edit window. The edit window help button ? is also useful for learning how to format your questions and answers. Also, it's better without the In/Out tags (many put comment signs (*, *) around the output). That makes it easier for those who might help you to copy the code into Mathematica. The easier it is, the more likely they will try. – Michael E2 Jun 21 '15 at 1:39
• You might try FiniteElementMethod, as described here. – bbgodfrey Jun 21 '15 at 6:11
• I'd consider it a bug that it does not print an error message. The method of lines can't solve problems where boundary values have been specified on all four boundaries. I'm not sure what the problem with the finite element method is, but halirutan has already shown a better way to set that up, unless polar coordinates are required. To use the method of lines, you can manually do a shooting method on your simple example, but I doubt it will work on your actual one, without knowing what it is. – Michael E2 Jun 22 '15 at 1:04

First thing: I have very limited experience with PDEs but what you present doesn't look like the Laplace equation. AFAIK the r and r^2 need to be in the denominator.

You try to solve this in polar coordinates. Maybe it is possible to use NDSolve if you can do your real problem in Cartesian coordinates, because it seems the new region functionality within Mathematica makes it possible.

Let's say we specify your annulus region and the boundary conditions in Cartesian coordinates, then we could simply write:

Ω = ImplicitRegion[1 <= x^2 + y^2 < 4, {x, y}];
bound = {
DirichletCondition[u[x, y] == 0, x^2 + y^2 == 1],
DirichletCondition[u[x, y] == Sin[ArcTan[x, y]], x^2 + y^2 == 4]
};
RegionPlot[Ω]


The boundary conditions are directly translated from yours. So you want u to be 0 at the circle r=1. Therefore I specified

DirichletCondition[u[x, y] == 0, x^2 + y^2 == 1]

Additionally, you want V[2,phi]==Sin[phi] which is nothing more than modulating the circle at r=2 with a sine. The angle polar phi can be written as ArcTan[x,y] and therefore, I believe your second condition is

DirichletCondition[u[x, y] == Sin[ArcTan[x, y]], x^2 + y^2 == 4]

The good thing is that this can directly be used with NDSolve.

duval = NDSolveValue[{Laplacian[u[x, y], {x, y}] == 0, bound}, u, {x, y} ∈ Ω]


and a quick plot reveals, that it seems to give indeed a solution where at r=0 the values and derivatives are periodic:

Plot3D[duval[x, y], {x, y} ∈ Ω]


Or you plot this as rectangular where one side is r and the other phi to see that the solution goes along with the answer of bbgodfrey

Plot3D[duval[r*Cos[phi], r*Sin[phi]], {r, 1, 2}, {phi, 0, 2 Pi}]


• Very nicely done. Apparently, using a region causes NDSolve to default to the finite element method. I believe that your outer boundary condition should be DirichletCondition[u[x, y] == y/2, x^2 + y^2 == 4] to reproduce the Sin[phi] variation in the Question. – bbgodfrey Jun 22 '15 at 2:32
• @bbgodfrey I have to say that I just presented the idea, I have chosen an arbitrary boundary condition, because the deq of the OP doesn't look exactly like the laplace equation in polar coordinates. I think the boundary needs to be Sin[ArcTan[x, y]] and a plot (in Cartesian) shows that it is similar to your result: i.stack.imgur.com/DVbOq.png – halirutan Jun 22 '15 at 11:00
• @bbgodfrey I have added a section about the boundary. I additionally compared your first solution with my one. Seems they are the same i.stack.imgur.com/8hetW.png – halirutan Jun 22 '15 at 11:21

I have not been able to obtain a solution in a general way. Nonetheless, some progress can be made. If we happen to know the value of V at phi == 0 and phi == 2 Pi, then using it of course gives the desired solution,

s = NDSolveValue[{r^2*D[V[r, phi], r, r] + r*D[V[r, phi], r] +
D[V[r, phi], phi, phi] == 0, V[1, phi] == 0,
V[2, phi] == Sin[phi], V[r, 2*Pi] == 0,
V[r, 0] == 0}, V, {r, 1, 2}, {phi, 0, 2*Pi}];
Plot3D[s[r, phi], {r, 1, 2}, {phi, 0, 2*Pi}, AxesLabel -> {r, phi, V},
AxesStyle -> Directive[Black, Bold, 12]]


A slightly more difficult problem can be solved in the same way.

s = NDSolveValue[{r^2*D[V[r, phi], r, r] + r*D[V[r, phi], r] +
D[V[r, phi], phi, phi] == 0, V[1, phi] == 0,
V[2, phi] == Sin[phi] + 1, V[r, 2*Pi] == Log[r]/Log[2],
V[r, 0] == Log[r]/Log[2]}, V, {r, 1, 2}, {phi, 0, 2*Pi}];
Plot3D[s[r, phi], {r, 1, 2}, {phi, 0, 2*Pi}, AxesLabel -> {r, phi, V},
AxesStyle -> Directive[Black, Bold, 12]]


Suppose, though, that we do not know V at the boundaries in phi. Then, a reasonable approximate answer can be obtained by guessing a boundary condition, say 0, enlarging the range of integration in phi to {phi, 0, 6*Pi} and using the central 1/3 of the resulting solution.

s2 = NDSolveValue[{r^2*D[V[r, phi], r, r] + r*D[V[r, phi], r] +
r^2 V[r, phi] + D[V[r, phi], phi, phi] == 0, V[1, phi] == 0,
V[2, phi] == Sin[phi] + 1, V[r, 6*Pi] == 0, V[r, 0] == 0},
V, {r, 1, 2}, {phi, 0, 6*Pi}];
Plot3D[s2[r, phi], {r, 1, 2}, {phi, 0, 6*Pi},
AxesLabel -> {r, phi, V}, AxesStyle -> Directive[Black, Bold, 12]]


The central 1/3 of this last plot agrees with the second plot to of order 10^-4. I look forward to a better answer. In the meantime, this may be sufficient.

If both V and its first derivative are known at one of the boundaries in r, a solution also can be obtained straightforwardly.

s = NDSolveValue[{r^2*D[V[r, phi], r, r] + r*D[V[r, phi], r] +
D[V[r, phi], phi, phi] == 0, V[1, phi] == 0, (D[V[r, phi], r] /. r -> 1) ==
.693 (Sin[phi] + 1), V[r, 2*Pi] == V[r, 0]}, V, {r, 1, 2}, {phi, 0, 2*Pi}];


which produces a result equal to that in the second plot.

• I might be wrong, but as far as I see this, the OP wants to solve the laplace equation on an annulus. There, values at the inner boundary would arise naturally since there really is no boundary in polar coordinates. I guess there is no way to specify this directly. Simply leaving them out doesn't work either: Import["http://goo.gl/NaH6rM"]["http://i.stack.imgur.com/JgpcX.png"]. Your link to the FEM is probably the best choice. – halirutan Jun 22 '15 at 0:02
• @halirutan The solution I provided is for an annulus, just not plotted in polar coordinates. I was not able to get FEM to work, although I agree that it is a more promising approach in general. – bbgodfrey Jun 22 '15 at 0:16

I'm aware to be slightly four years to late. But here I add the solution in polarcoordinates using periodic boundary conditions.

NDSolve evaluates a solution without error messages

V = NDSolveValue[{Derivative[2, 0][v][r, phi] +1/r Derivative[1, 0][v][r, phi] +1/r^2 Derivative[0, 2][v][r, phi] == 0
, DirichletCondition[v[r, phi] == 0, r == 1],
DirichletCondition[v[r, phi] == Sin[phi], r == 2]

, PeriodicBoundaryCondition[v[r, phi],
phi == \[Pi] && (1 < r < 2), TranslationTransform[{0, -2 \[Pi]}]] },
v, {r, 1, 2}, {phi, -Pi, Pi}(*,Method\[Rule]"FiniteElement"*)]


but the plot

ParametricPlot3D[{r Cos[phi], r Sin[phi], V[r, phi]}, {r, 1,2}, {phi, -Pi, Pi}]


shows a solution which differs near phi==2 Pi from the expected solution(see cartesian answer @haliutran).

Obviously the PeriodicBoundaryCondition forces a unplausible NDSolve solution in polarcoordinates!

• Perhaps, a bug? – bbgodfrey Aug 30 '19 at 4:39
• @bbgodfrey Probably yes, but I could only check it in v11.0.1. – Ulrich Neumann Aug 30 '19 at 5:53
• I checked it with 12.0. Clearly, not the correct answer. – bbgodfrey Aug 30 '19 at 12:21
• @bbgodfrey That emans no difference between v11 and v12? – Ulrich Neumann Aug 30 '19 at 12:23
• I did not check it with v11, but the problem seemed localized at -Pi and Pi, as you indicated. – bbgodfrey Aug 30 '19 at 12:25