5
$\begingroup$

I can't get Mathematica to solve this standard textbook PDE, which is Laplace inside a disk of some radius. One of the boundary conditions needed is that the solution is finite (bounded) in center of disk, and I do not know how specify this boundary condition. I do not know if this is why Mathematica can't solve it. Most examples in help are for NDSolve and finite elements, and I am looking to analytical solution.

This is the formal description of the PDE

Solve for $u\left( r,\theta\right) $ the PDE

\begin{align*} \frac{\partial^{2}u}{\partial r^{2}}+\frac{1}{r}\frac{\partial u}{\partial r}+\frac{1}{r^{2}}\frac{\partial^{2}u}{\partial\theta^{2}} & =0\\ 0 & \leq r\leq a\\ 0 & <\theta\leq2\pi \end{align*}

Boundary conditions

\begin{align*} u\left( a,\theta\right) & =f\left( \theta\right) \\ \left\vert u\left( 0,\theta\right) \right\vert & <\infty\\ u\left( r,0\right) & =u\left( r,2\pi\right) \\ \frac{\partial u}{\partial\theta}\left( r,0\right) & =\frac{\partial u}{\partial\theta}\left( r,2\pi\right) \end{align*}

The standard solution using seperation of variables by hand gives

\begin{align*} u\left( r,\theta\right) & =A_{0}+\sum_{n=1}^{\infty}r^{n}\left( c_{n}\cos n\theta+b_{n}\sin n\theta\right) \\ A_{0} & =\frac{1}{2\pi}\int_{0}^{2\pi}f\left( \theta\right) d\theta\\ c_{n} & =\frac{1}{\pi a^{n}}\int_{0}^{2\pi}f\left( \theta\right) \cos\left( n\theta\right) d\theta\qquad n>0\\ b_{n} & =\frac{1}{\pi a^{n}}\int_{0}^{2\pi}f\left( \theta\right) \sin\left( n\theta\right) d\theta\qquad n>0 \end{align*}

This is what I tried

ClearAll[u,theta,r,a];
pde=D[u[r,theta],{r,2}]+1/r D[u[r,theta],r]+1/r^2 D[u[r,theta],{theta,2}]==0;
bc=u[a,theta]==f[theta];
sol=DSolve[{pde,bc},u[r,theta],{r,theta},Assumptions->a<r&&a>0&& 0<theta<=2 Pi]

Mathematica graphics

Now I tried adding the bounded condition, hoping it will help. But I do not know what the syntax should be

ClearAll[u,theta,r,a];
pde=D[u[r,theta],{r,2}]+1/r D[u[r,theta],r]+1/r^2 D[u[r,theta],{theta,2}]==0;
bc={u[a,theta]==f[theta],Abs[u[0,theta]]<Infinity}
sol=DSolve[{pde,bc},u[r,theta],{r,theta},Assumptions->a<r&&a>0&& 0<theta<=2 Pi]

Mathematica graphics

Ok, so I tried now to add the periodic conditions. But this did not help

ClearAll[u,theta,r,a];
pde=D[u[r,theta],{r,2}]+1/r D[u[r,theta],r]+1/r^2 D[u[r,theta],{theta,2}]==0;

bc={u[a,theta]==f[theta],
    u[r,0]==u[r,2 Pi],
    Derivative[0,1][u][r,0]==Derivative[0,1][u][r,2 Pi]}

sol=DSolve[{pde,bc},u[r,theta],{r,theta},Assumptions->a<r&&a>0&& 0<theta<=2 Pi]

Mathematica graphics

Question is: Is Mathematica not able to solve it, because it is missing the bounded conditions part? If so, how to specify this? Any one could get Mathematica to solve this analytically?

This is a very standard PDE problem and its analytical solution is all over the net. The above solution I show is the standard text book solution.

11.3 on windows 7

$\endgroup$
  • $\begingroup$ The Dirichlet problem under consideration can be solved in Mathematica by making use of Poisson kernel (see en.wikipedia.org/wiki/Poisson_kernel ). $\endgroup$ – user64494 Mar 16 '18 at 9:55
  • 1
    $\begingroup$ @user, maybe write an answer demonstrating how to use the Poisson kernel, should you find time? $\endgroup$ – J. M. is away Mar 17 '18 at 11:26
4
$\begingroup$

Once again, since DSolve doesn't work well at the moment, I'd like to post an answer based on finite Fourier transform:

(* Definition of finiteFourierTransform etc. are not included in this post, 
   please find them in the link above. *)
pde = D[u[r, theta], {r, 2}] + 1/r D[u[r, theta], r] + 
    1/r^2 D[u[r, theta], {theta, 2}] == 0;

bc = {u[a, theta] == f[theta], u[r, 0] == u[r, 2 Pi], 
   Derivative[0, 1][u][r, 0] == Derivative[0, 1][u][r, 2 Pi]};

Format@finiteFourierTransform[f_, __] := ℱ[f]
Format@theta := θ

finiteFourierTransform[{pde, bc[[1]]}, {theta, 0, 2 Pi}, n]

Mathematica graphics

We substitute the periodic b.c.s in, and solve the resulting system:

% /. Rule @@@ bc[[2 ;; 3]]

tset = % /. HoldPattern@finiteFourierTransform[f_ /; ! FreeQ[f, u], __] :> f
tsol = DSolve[tset, u[r, theta], r][[1, 1, -1]]

Mathematica graphics

One constant C[1] remains, but we still have the bounded condition / finiteness condition at hand. With a bit of transforming:

tsolcollect = Collect[Simplify[tsol // TrigToExp, n ∈ Integers], r^_]

Mathematica graphics

It becomes clear that the coefficient of r^-n should be 0, because $\lim_{r\to 0^+} r^{-n}=\infty$ for $n>0$:

tsolfinal = tsolcollect /. Flatten@Solve[tsolcollect[[1]] == 0, C[1]] // Simplify

The last step is to transform back:

sol = inverseFiniteFourierTransform[tsolfinal, n, {theta, 0, 2 Pi}]
sol // transformToIntegrate

Mathematica graphics

Though it's troublesome to make Mathematica simplify sol further, it's not hard to notice this answer is equivalent to the one in your question.

sol can be used for further calculation of course. For example, when $f=\sin(3\theta)$:

ReleaseHold[sol /. f -> (Sin[3 #] &) /. C -> 7 /. a -> 2] // ComplexExpand // Simplify
(* 1/8 r^3 Sin[3 theta] *)
RevolutionPlot3D[%, {r, 0, 2}, {theta, 0, 2 Pi}, PlotRange -> All]

Mathematica graphics

$\endgroup$
  • $\begingroup$ Sorry, but C -> 7 (as well as several other steps) is built on the sand. $\endgroup$ – user64494 Mar 16 '18 at 21:02
  • 2
    $\begingroup$ @user64494 Er… what do you mean by "build on the sand"? Taking the first several terms of a series to approximate its value is a common practice. (f -> (Sin[3 #] & happens to be a special case that only finiteFourierTransform[Sin[3 theta], {theta, 0, 2 Pi}, 3] is non-zero so the final solution is not a series though. ) $\endgroup$ – xzczd Mar 17 '18 at 3:44
2
$\begingroup$

This is very strange. I suspect Mathematica desktop is being updated/patched from the cloud when one is connected to the internet, because I just tried again the first attempt I showed in my question above and now DSolve return the solution to the Laplace PDE !

I do have internet setting in preference to allow Wolfram access to do updates. But I thought this is for documentations and things like that only. I do not understand how this could have happened. But Mathematica 11.3 can actually solve this PDE.

Mathematica graphics

code again (same as in my question)

ClearAll[u,theta,r,a];
pde=D[u[r,theta],{r,2}]+1/r D[u[r,theta],r]+1/r^2 D[u[r,theta],{theta,2}]==0;
bc=u[a,theta]==f[theta];
sol=DSolve[{pde,bc},u[r,theta],{r,theta},Assumptions->a<r&&a>0&&0<theta<=2 Pi]
sol=sol/.K[1]->n

The third attempt I had in my question, still does not work, but the first works.

I assume someone tried that also and it did not work for them, else they would have screamed at me saying that it works. Very strange. But happy that 11.3 can solve this.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.