____________________________________________________________________________________
Taken from Mathematica DSolve help pages
Solve the heat equation in polar coordinates for \(u(r,t)\) \[ \frac { \partial u}{\partial t}= \frac { \partial ^2 u}{\partial r^2} + \frac {1}{r} \frac { \partial u}{\partial r} \] For \(0<r<1\) and \(t>0\). The boundary conditions are \begin {align*} u(1,t) &= 0 \end {align*}
Initial condition is \(u(r,0)=1-r\).
Mathematica ✓
ClearAll[r, t, u]; pde = D[u[r, t], t] == D[u[r, t], {r, 2}] + (1*D[u[r, t], r])/r; ic = u[r, 0] == 1 - r; bc = u[1, t] == 0; sol = AbsoluteTiming[TimeConstrained[DSolve[{pde, ic, bc}, u[r, t], {r, t}], 60*10]]; sol = sol /. K[1] -> n;
\[ \left \{\left \{u(r,t)\to \sum _{n=1}^{\infty }\frac {2 e^{-t \text {BesselJZero}(0,n)^2} \text {BesselJ}(0,r \text {BesselJZero}(0,n)) \left (\frac {\text {BesselJ}(1,\text {BesselJZero}(0,n))}{\text {BesselJZero}(0,n)}-\frac {1}{3} \text {HypergeometricPFQ}\left (\left \{\frac {3}{2}\right \},\left \{1,\frac {5}{2}\right \},-\frac {1}{4} \text {BesselJZero}(0,n)^2\right )\right )}{\text {BesselJ}(0,\text {BesselJZero}(0,n))^2+\text {BesselJ}(1,\text {BesselJZero}(0,n))^2}\right \}\right \} \]
Maple ✓
r:='r'; u:='u'; t:='t'; pde:=diff(u(r,t),t)= diff(u(r,t),r$2)+ 1/r*diff(u(r,t),r); ic:=u(r,0)=1-r; bc := u(1,t) =0; cpu_time := timelimit(60*10,CodeTools[Usage](assign('sol',pdsolve([pde,ic,bc],u(r,t),HINT =boundedseries(r=0))),output='realtime'));
\[ u \left ( r,t \right ) ={\it invlaplace} \left ( \BesselJ \left ( 0,\sqrt {-s}r \right ) {\it \_F2} \left ( s \right ) ,s,t \right ) -{\it invlaplace} \left ( {\frac {\BesselY \left ( 0,\sqrt {-s}r \right ) \BesselJ \left ( 0,1/2\,\StruveH \left ( 0,\sqrt {-s} \right ) \pi +\LommelS 1 \left ( 2,0,\sqrt {-s} \right ) \right ) {\it \_F2} \left ( s \right ) }{\BesselY \left ( 0,1/2\,\StruveH \left ( 0,\sqrt {-s} \right ) \pi +\LommelS 1 \left ( 2,0,\sqrt {-s} \right ) \right ) }},s,t \right ) -1/2\,\pi \,{\it invlaplace} \left ( {\frac {\BesselY \left ( 0,\sqrt {-s}r \right ) \StruveH \left ( 0,1/2\,\StruveH \left ( 0,\sqrt {-s} \right ) \pi +\LommelS 1 \left ( 2,0,\sqrt {-s} \right ) \right ) \LommelS 1 \left ( 2,0,1/2\,\StruveH \left ( 0,\sqrt {-s} \right ) \pi +\LommelS 1 \left ( 2,0,\sqrt {-s} \right ) \right ) }{\BesselY \left ( 0,1/2\,\StruveH \left ( 0,\sqrt {-s} \right ) \pi +\LommelS 1 \left ( 2,0,\sqrt {-s} \right ) \right ) {s}^{2}}},s,t \right ) -{\it invlaplace} \left ( {\frac {\BesselY \left ( 0,\sqrt {-s}r \right ) \left ( \LommelS 1 \left ( 2,0,1/2\,\StruveH \left ( 0,\sqrt {-s} \right ) \pi +\LommelS 1 \left ( 2,0,\sqrt {-s} \right ) \right ) \right ) ^{2}}{\BesselY \left ( 0,1/2\,\StruveH \left ( 0,\sqrt {-s} \right ) \pi +\LommelS 1 \left ( 2,0,\sqrt {-s} \right ) \right ) {s}^{2}}},s,t \right ) -{\it invlaplace} \left ( {\frac {\BesselY \left ( 0,\sqrt {-s}r \right ) }{\BesselY \left ( 0,1/2\,\StruveH \left ( 0,\sqrt {-s} \right ) \pi +\LommelS 1 \left ( 2,0,\sqrt {-s} \right ) \right ) s}},s,t \right ) +{\it invlaplace} \left ( {\frac {\LommelS 1 \left ( 2,0,\sqrt {-s}r \right ) }{ \left ( -s \right ) ^{3/2}}},s,t \right ) +1 \] But has unresolved inverse Laplace transforms
____________________________________________________________________________________
Added Nov 24, 2018.
Problem 8.3.5 from Richard Haberman applied partial differential equations book, 5th edition
Solve for \(u(r,t)\) \[ \frac { \partial u}{\partial t}= k \nabla ^2 u + f(r,t) \] Inside the circle (\(r<a\)) with \(u=0\) at \(r=a\) and initially \(u=0\).
One of the problems here, is how to tell CAS the implicit condition when solving this which is that \(u(0,t)<\infty \).
Mathematica ✗
ClearAll[r, t, u, k, a]; pde = D[u[r, t], t] == k*(D[u[r, t], {r, 2}] + D[u[r, t], r]/r) + f[r, t]; ic = u[r, 0] == 0; bc = u[a, t] == 0; sol = AbsoluteTiming[TimeConstrained[DSolve[{pde, ic, bc}, u[r, t], {r, t}, Assumptions -> r < a], 60*10]];
\[ \text {Failed} \]
Maple ✗
r:='r'; u:='u'; t:='t';a:='a';k:='k';f:='f'; pde:=diff(u(r,t),t)= k*(diff(u(r,t),r$2)+ 1/r*diff(u(r,t),r)) + f(r,t); ic:=u(r,0)=0; bc := u(a,t) =0; #do not use HINT=boundedseries below, Maple will not solve it then cpu_time := timelimit(60*10,CodeTools[Usage](assign('sol',pdsolve([pde,ic,bc],u(r,t)) assuming r<a),output='realtime'));
\[ \text { sol=() } \]
Hand solution
Since this problem has homogeneous B.C. but has time dependent source (i.e. non-homogenous in the PDE itself), then we will use the method of eigenfunction expansion. In this method, we first find the eigenfunctions \(\phi _{n}\left ( x\right ) \) of the associated homogenous PDE without the source being present. Then use these \(\phi _{n}\left ( x\right ) \) to expand the source \(f\left ( x,t\right ) \) as generalized Fourier series. We now switch to the associated homogenous PDE in order to find the eigenfunctions. \(u\equiv u\left ( r,t\right ) \). There is no \(\theta \). Hence\begin {align} \frac {\partial u\left ( r,t\right ) }{\partial t} & =k\left ( \frac {\partial ^{2}u}{\partial r^{2}}+\frac {1}{r}\frac {\partial u}{\partial r}\right ) \tag {1}\\ u\left ( a,t\right ) & =0\nonumber \\ \left \vert u\left ( 0,t\right ) \right \vert & <\infty \nonumber \\ u\left ( r,0\right ) & =0\nonumber \end {align}
We need to solve the above in order to find the eigenfunctions \(\phi _{n}\left ( r\right ) \). Let \(u=R\left ( r\right ) T\left ( t\right ) \). Substituting this back into (1) gives\[ T^{\prime }R=k\left ( R^{\prime \prime }T+\frac {1}{r}R^{\prime }T\right ) \] Dividing by \(RT\)\[ \frac {1}{k}\frac {T^{\prime }}{T}=\frac {R^{\prime \prime }}{R}+\frac {1}{r}\frac {R^{\prime }}{R}=-\lambda \] Where \(\lambda \) is the separation constant. The above gives\[ T^{\prime }+k\lambda T=0 \] And\[ rR^{\prime \prime }+R^{\prime }+\lambda rR=0 \] This is a singular Sturm-Liouville ODE. Standard form is\[ \left ( rR^{\prime }\right ) ^{\prime }=-\lambda rR \] Hence\begin {align*} p & =r\\ q & =0\\ \sigma & =r \end {align*}
The ODE\(\ rR^{\prime \prime }+R^{\prime }+\lambda rR=0\) is Bessel ODE whose solution is\[ R\left ( r\right ) =C_{1}\operatorname {BesselJ}\left ( 0,\sqrt {\lambda }r\right ) +C_{2}\operatorname {BesselY}\left ( 0,\sqrt {\lambda }r\right ) \] Since \(\operatorname {BesselY}\left ( 0,\sqrt {\lambda }r\right ) \) blows up at \(r=0\), then \(C_{2}=0\) and the solution becomes \[ R\left ( r\right ) =C_{1}\operatorname {BesselJ}\left ( 0,\sqrt {\lambda }r\right ) \] At \(r=a\) the above becomes \(0=C_{1}\operatorname {BesselJ}\left ( 0,\sqrt {\lambda }a\right ) \). Non trivial solution requires that \(\sqrt {\lambda }a\) are the zeros of \(\operatorname {BesselJ}\left ( 0,x\right ) \). Let the zeros be called \(\Lambda _{n},n=1,2,3,\cdots \). Therefore \(\sqrt {\lambda _{n}}a=\Lambda _{n}\) or\[ \lambda _{n}=\left ( \frac {\Lambda _{n}}{a}\right ) ^{2}\qquad n=1,2,3,\cdots \] The corresponding eigenfunctions are \(R_{n}\left ( r\right ) =\operatorname {BesselJ}\left ( 0,\frac {\Lambda _{n}}{a}r\right ) \). Now that the eigenfunctions for the homogeneous PDE are found, eigenfunction expansion is used to find the general solution. Let \begin {equation} u\left ( r,t\right ) =\sum _{n=1}^{\infty }a_{n}\left ( t\right ) \operatorname {BesselJ}\left ( 0,\frac {\Lambda _{n}}{a}r\right ) \tag {2} \end {equation} Where \(a_{n}\left ( t\right ) \) is function of time since it includes the time solution in it. Substituting the above back into the original nonhomogeneous PDE \begin {align} u_{t} & =k\nabla ^{2}u+f\left ( r,t\right ) \tag {3}\\ & =k\left ( \frac {\partial ^{2}u}{\partial r^{2}}+\frac {1}{r}\frac {\partial u}{\partial r}\right ) +f\left ( r,t\right ) \nonumber \end {align}
Where \(\nabla ^{2}u=-\lambda ru\). Substituting (2) into (3), and using \(f\left ( r,t\right ) =\sum _{n=1}^{\infty }b_{n}\left ( t\right ) \operatorname {BesselJ}\left ( 0,\frac {\Lambda _{n}}{a}r\right ) \) gives
\[ \sum _{n=1}^{\infty }a_{n}^{\prime }\left ( t\right ) \operatorname {BesselJ}\left ( 0,\frac {\Lambda _{n}}{a}r\right ) =ka_{n}\left ( t\right ) \left ( \sum _{n=1}^{\infty }\operatorname {BesselJ}^{\prime \prime }\left ( 0,\frac {\Lambda _{n}}{a}r\right ) +\frac {1}{r}\operatorname {BesselJ}^{^{\prime }}\left ( 0,\frac {\Lambda _{n}}{a}r\right ) \right ) +\sum _{n=1}^{\infty }b_{n}\left ( t\right ) \operatorname {BesselJ}\left ( 0,\frac {\Lambda _{n}}{a}r\right ) \]
But \(\sum _{n=1}^{\infty }\operatorname {BesselJ}^{\prime \prime }\left ( 0,\frac {\Lambda _{n}}{a}r\right ) +\frac {1}{r}\operatorname {BesselJ}^{^{\prime }}\left ( 0,\frac {\Lambda _{n}}{a}r\right ) =-\lambda _{n}\operatorname {BesselJ}\left ( 0,\frac {\Lambda _{n}}{a}r\right ) \). The above becomes
\begin {align*} \sum _{n=1}^{\infty }a_{n}^{\prime }\left ( t\right ) \operatorname {BesselJ}\left ( 0,\frac {\Lambda _{n}}{a}r\right ) & =-ka_{n}\left ( t\right ) \sum _{n=1}^{\infty }\left ( \frac {\Lambda _{n}}{a}\right ) ^{2}\operatorname {BesselJ}\left ( 0,\frac {\Lambda _{n}}{a}r\right ) +\sum _{n=1}^{\infty }b_{n}\left ( t\right ) \operatorname {BesselJ}\left ( 0,\frac {\Lambda _{n}}{a}r\right ) \\ \sum _{n=1}^{\infty }\left ( a_{n}^{\prime }\left ( t\right ) +k\left ( \frac {\Lambda _{n}}{a}\right ) ^{2}a_{n}\left ( t\right ) \right ) \operatorname {BesselJ}\left ( 0,\frac {\Lambda _{n}}{a}r\right ) & =\sum _{n=1}^{\infty }b_{n}\left ( t\right ) \operatorname {BesselJ}\left ( 0,\frac {\Lambda _{n}}{a}r\right ) \end {align*}
The above simplifies to\[ a_{n}^{\prime }\left ( t\right ) +k\left ( \frac {\Lambda _{n}}{a}\right ) ^{2}a_{n}\left ( t\right ) =b_{n}\left ( t\right ) \] The solution is\[ a_{n}\left ( t\right ) =e^{-k\left ( \frac {\Lambda _{n}}{a}\right ) ^{2}t}\int _{0}^{t}b_{n}\left ( \tau \right ) e^{k\left ( \frac {\Lambda _{n}}{a}\right ) ^{2}\tau }d\tau +a_{n}\left ( 0\right ) e^{-k\left ( \frac {\Lambda _{n}}{a}\right ) ^{2}t}\] Hence the solution (2) becomes\[ u\left ( r,t\right ) =\sum _{n=1}^{\infty }\left ( e^{-k\left ( \frac {\Lambda _{n}}{a}\right ) ^{2}t}\left ( \int _{0}^{t}b_{n}\left ( \tau \right ) e^{k\left ( \frac {\Lambda _{n}}{a}\right ) ^{2}\tau }d\tau \right ) +a_{n}\left ( 0\right ) e^{-k\left ( \frac {\Lambda _{n}}{a}\right ) ^{2}t}\right ) \operatorname {BesselJ}\left ( 0,\frac {\Lambda _{n}}{a}r\right ) \] To find \(a_{n}\left ( 0\right ) \), putting \(t=0\) in the above gives\[ 0=\sum _{n=1}^{\infty }a_{n}\left ( 0\right ) \operatorname {BesselJ}\left ( 0,\frac {\Lambda _{n}}{a}r\right ) \] Hence \(a_{n}\left ( 0\right ) =0\). Therefore \(a_{n}\left ( t\right ) \) becomes.\[ a_{n}\left ( t\right ) =e^{-k\left ( \frac {\Lambda _{n}}{a}\right ) ^{2}t}\int _{0}^{t}b_{n}\left ( \tau \right ) e^{k\left ( \frac {\Lambda _{n}}{a}\right ) ^{2}\tau }d\tau \] Hence the solution from (2) now becomes\[ u\left ( r,t\right ) =\sum _{n=1}^{\infty }\left ( e^{-k\left ( \frac {\Lambda _{n}}{a}\right ) ^{2}t}\int _{0}^{t}b_{n}\left ( \tau \right ) e^{k\left ( \frac {\Lambda _{n}}{a}\right ) ^{2}\tau }d\tau \right ) \operatorname {BesselJ}\left ( 0,\frac {\Lambda _{n}}{a}r\right ) \] And finally, to find \(b_{n}\left ( t\right ) \), which is the generalized Fourier coefficient of the expansion of the source in (3) above, orthogonality is used as follows\begin {align*} \int _{0}^{a}f\left ( r,t\right ) \operatorname {BesselJ}\left ( 0,\frac {\Lambda _{n}}{a}r\right ) rdr & =b_{n}\left ( t\right ) \int _{0}^{a}\operatorname {BesselJ}^{2}\left ( 0,\frac {\Lambda _{n}}{a}r\right ) rdr\\ b_{n}\left ( t\right ) & =\frac {\int _{0}^{a}f\left ( r,t\right ) \operatorname {BesselJ}\left ( 0,\frac {\Lambda _{n}}{a}r\right ) rdr}{\int _{0}^{a}\operatorname {BesselJ}^{2}\left ( 0,\frac {\Lambda _{n}}{a}r\right ) rdr} \end {align*}
Summary of solution \begin {align*} u\left ( r,t\right ) & =\sum _{n=1}^{\infty }a_{n}\left ( t\right ) \operatorname {BesselJ}\left ( 0,\frac {\Lambda _{n}}{a}r\right ) \\ & =\sum _{n=1}^{\infty }\left ( \int _{0}^{t}b_{n}\left ( \tau \right ) e^{k\left ( \frac {\Lambda _{n}}{a}\right ) ^{2}\tau }d\tau \right ) e^{-k\left ( \frac {\Lambda _{n}}{a}\right ) ^{2}t}\operatorname {BesselJ}\left ( 0,\frac {\Lambda _{n}}{a}r\right ) \end {align*}
Where\[ b_{n}\left ( t\right ) =\frac {\int _{0}^{a}f\left ( r,t\right ) \operatorname {BesselJ}\left ( 0,\frac {\Lambda _{n}}{a}r\right ) rdr}{\int _{0}^{a}\operatorname {BesselJ}^{2}\left ( 0,\frac {\Lambda _{n}}{a}r\right ) rdr}\]
____________________________________________________________________________________
Added December 20, 2018.
Example 6.9.1 from Partial differential equations and boundary value problems with Maple/George A. Articolo, 2nd ed :
We seek the temperature distribution \(u(r,\theta ,t)\) in a thin circular plate over the two-dimensional domain \(D=\{(r,\theta ) | 0<r<1, 0< \theta < \frac {\pi }{2}\}\).
The lateral surfaces of the plate are insulated. The edges \(r = 1\) and \(\theta = 0\) are at a fixed temperature of \(0\), and the edge \(\theta = \frac {\pi }{2}\) is insulated. The initial temperature distribution \(u(r, \theta , 0) = f(r, \theta )\) is \(u(r,\theta ,0)=(r-r^3)\sin (\theta )\).
The thermal diffusivity is \(k = \frac {1}{50}\). Solve for \(u(r,\theta ,t)\) the heat PDE \[ u_t= k \left ( u_{rr} + \frac {1}{r} u_r +\frac {1}{r^2} u_{\theta \theta } \right ) \] With boundary conditions \begin {align*} |u(0,\theta ,t)| &< \infty \\ u(1,\theta ,t) &= 0\\ u(r,0,t) &=0\\ \frac {\partial u}{\partial \theta }(1,\frac {\pi }{2},t) &= 0 \end {align*}
Mathematica ✗
ClearAll[r, t, theta, u, k]; k = 1/50; pde = D[u[r, theta, t], t] == k*(D[u[r, theta, t], {r, 2}] + D[u[r, theta, t], r]/r + (1*D[u[r, theta, t], {theta, 2}])/r^2); bcOnR = u[1, theta, t] == 0; bcOnTheta = {u[r, 0, t] == 0, Derivative[0, 1, 0][u][r, Pi/2, t] == 0}; ic = u[r, theta, 0] == (r - (1*r^3)/3)*Sin[theta]; sol = AbsoluteTiming[TimeConstrained[DSolve[{pde, bcOnR, bcOnTheta, ic}, u[r, theta, t], {r, theta, t},Assumptions -> t > 0], 60*10]];
\[ \text {Failed} \]
Maple ✗
unassign('r,u,t,theta,k'); k:=1/50; pde := diff(u(r, theta, t), t) = k*( diff(u(r, theta, t), r$2) + 1/r* diff(u(r, theta, t), r) + 1/r^2 * diff(u(r ,theta, t), theta$2) ); bc_on_r:= u(1,theta,t)=0; bc_on_theta:= u(r,0,t)=0, eval(diff(u(r,theta,t),theta),theta=Pi/2)=0; ic := u(r,theta,0)=(r-1/3*r^3)*sin(theta); cpu_time := timelimit(60*10,CodeTools[Usage](assign('sol',pdsolve([pde, bc_on_r, bc_on_theta, ic], u(r, theta, t),HINT = boundedseries(r = [0])) assuming r<a,t>0),output='realtime'));
\[ \text { sol=() } \]
____________________________________________________________________________________
Added December 20, 2018.
Example 6.9.2 from Partial differential equations and boundary value problems with Maple/George A. Articolo, 2nd ed :
We seek the temperature distribution in a thin circular plate over the two-dimensional domain \( D = \{(r,\theta ) | 0 < r < 1, 0 < \theta < \pi \}\). The lateral surfaces of the plate are insulated. The sides \(\theta = 0\) and \(\theta = \pi \) are at a fixed temperature of \(0\), and the edge \(r = 1\) is insulated. The initial temperature distribution is \(u(r, \theta , 0) = \left ( r - \frac {r^3}{3} \right )\sin \theta \).
The thermal diffusivity is \(k = \frac {1}{25}\). Solve for \(u(r,\theta ,t)\) the heat PDE \[ u_t= k \left ( u_{rr} + \frac {1}{r} u_r +\frac {1}{r^2} u_{\theta \theta } \right ) \] With boundary conditions \begin {align*} |u(0,\theta ,t) &< \infty \\ u(1,\theta ,t) &= 0\\ u(r,0,t) &=0\\ u(r,\pi ,t) & 0 \end {align*}
Mathematica ✗
ClearAll[r, t, theta, u, k]; k = 1/25; pde = D[u[r, theta, t], t] == k*(D[u[r, theta, t], {r, 2}] + D[u[r, theta, t], r]/r + (1*D[u[r, theta, t], {theta, 2}])/r^2); bcOnR = Derivative[1, 0, 0][u][1, theta, t] == 0; bcOnTheta = {u[r, 0, t] == 0, u[r, Pi, t] == 0}; ic = u[r, theta, 0] == (r - (1*r^3)/3)*Sin[theta]; sol = AbsoluteTiming[TimeConstrained[DSolve[{pde, bcOnR, bcOnTheta, ic}, u[r, theta, t], {r, theta, t}], 60*10]];
\[ \text {Failed} \]
Maple ✓
r:='r'; u:='u'; t:='t';theta:='theta';k:='k'; k:=1/25; pde := diff(u(r, theta, t), t) = k*( diff(u(r, theta, t), r$2) + 1/r* diff(u(r, theta, t), r) + 1/r^2 * diff(u(r, theta, t), theta$2) ); bc_on_r:= eval(diff(u(r,theta,t),r),r=1)=0; bc_on_theta:= u(r,0,t)=0, u(r,Pi,t)=0; ic := u(r,theta,0)=(r-1/3*r^3)*sin(theta); cpu_time := timelimit(60*10,CodeTools[Usage](assign('sol',pdsolve([pde, bc_on_r, bc_on_theta, ic], u(r, theta, t),HINT = boundedseries(r = [0]))),output='realtime'));
\[ u \left ( r,\theta ,t \right ) =\mbox {{\tt casesplit/ans}} \left ( \sum _{n=0}^{\infty }-4/3\,{\frac {\BesselJ \left ( 1,\lambda _{{n}}r \right ) \sin \left ( \theta \right ) {{\rm e}^{-1/25\,{\lambda _{{n}}}^{2}t}} \left ( \BesselJ \left ( 0,\lambda _{{n}} \right ) {\lambda _{{n}}}^{3}-\BesselJ \left ( 1,\lambda _{{n}} \right ) {\lambda _{{n}}}^{2}+4\,\lambda _{{n}}\BesselJ \left ( 0,\lambda _{{n}} \right ) -8\,\BesselJ \left ( 1,\lambda _{{n}} \right ) \right ) }{{\lambda _{{n}}}^{3} \left ( \left ( \BesselJ \left ( 0,\lambda _{{n}} \right ) \right ) ^{2}\lambda _{{n}}+ \left ( \BesselJ \left ( 1,\lambda _{{n}} \right ) \right ) ^{2}\lambda _{{n}}-2\,\BesselJ \left ( 0,\lambda _{{n}} \right ) \BesselJ \left ( 1,\lambda _{{n}} \right ) \right ) }}, \left \{ {\it And} \left ( -\BesselJ \left ( 1,\lambda _{{n}} \right ) +\BesselJ \left ( 2,\lambda _{{n}} \right ) \lambda _{{n}}=0,0<\lambda _{{n}} \right ) \right \} \right ) \]
____________________________________________________________________________________
Added Feb 24, 2019.
Problem 8.2.5 from from Richard Haberman applied partial differential equations book, 5th edition.
Solve the initial value problem for a two-dimensional heat equation inside a circle (of radius \(a\)) \(u_t = k \nabla ^2 u\) with time-independent boundary conditions: \begin {align*} u(a,\theta ,t) &= g(\theta ) \end {align*}
And initial conditions \(u(r,\theta ,0) =f(r,\theta )\). There is an implied periodic boundary conditions on \(\theta \) \begin {align*} u(r,-\pi ,t) &= u(r,\pi ,t)\\ \frac {\partial u}{\partial \theta }(r,-\pi ,t) &= \frac {\partial u}{\partial \theta }(r,\pi ,t) \end {align*}
Mathematica ✗
ClearAll[r, t, theta, u, a, k, g, f]; pde = D[u[r, theta, t], t] == k*(D[u[r, theta, t], {r, 2}] + D[u[r, theta, t], r]/r + D[u[r, theta, t], {theta, 2}]/r^2); bcOnR = u[a, theta, t] == g[theta]; bcOnTheta = {u[r, -Pi, t] == u[r, Pi, t], Derivative[0, 1, 0][u][r, -Pi, t] == Derivative[0, 1, 0][u][r, Pi, t]}; ic = u[r, theta, 0] == f[r, theta]; sol = AbsoluteTiming[TimeConstrained[DSolve[{pde, bcOnR, bcOnTheta, ic}, u[r, theta, t], {r, theta, t}, Assumptions -> {a > 0, a < r, k > 0}], 60*10]];
\[ \text {Failed} \]
Maple ✗
unassign('r,u,t,theta,g,f'); pde := diff(u(r,theta,t),t)=k*(diff(u(r,theta,t),r$2) + 1/r*diff(u(r,theta,t),r)+1/r^2*diff(u(r,theta,t),theta$2)); bcOnR:= u(a,theta,t)=g(theta); bcOnTheta:= u(r,-Pi,t)=u(r,Pi,t),eval(diff(u(r,theta,t),theta),theta=-Pi)=eval(diff(u(r,theta,t),theta),theta=Pi); ic := u(r,theta,0)=f(r,theta); cpu_time := timelimit(60*10,CodeTools[Usage](assign('sol',pdsolve([pde, bcOnR,bcOnTheta, ic], u(r, theta, t), HINT = boundedseries(r = [0])) assuming a>0,a<r,k>0),output='realtime'));
\[ \text { sol=() } \]
Hand solution
Solve
\begin {align*} \frac {\partial u\left ( r,\theta ,t\right ) }{\partial t} & =k\left ( \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}}\right ) \\ \left \vert u\left ( 0,\theta ,t\right ) \right \vert & <\infty \\ u\left ( a,\theta ,t\right ) & =g\left ( \theta \right ) \\ u\left ( r,-\pi ,t\right ) & =u\left ( r,\pi ,t\right ) \\ \frac {\partial u}{\partial \theta }\left ( r,-\pi ,t\right ) & =\frac {\partial u}{\partial \theta }\left ( r,\pi ,t\right ) \end {align*}
With initial conditions \(u\left ( r,\theta ,0\right ) =f\left ( r,\theta \right ) \).
Since the boundary conditions are not homogenous, and since there are no time dependent sources, then in this case we look for \(u_{E}\left ( r,\theta \right ) \) which is solution at steady state which needs to satisfy the nonhomogeneous B.C., where \(u\left ( r,\theta ,t\right ) =\overset {\text {transient}}{\overbrace {v\left ( r,\theta ,t\right ) }}+\overset {\text {steady state}}{\overbrace {u_{E}\left ( r,\theta \right ) }}\) and \(v\left ( r,\theta ,t\right ) \) solves the PDE but with homogenous B.C. Therefore, we need to find equilibrium (steady state) solution for Laplace PDE on disk, that only needs to satisfy the nonhomogeneous B.C.\begin {align*} \nabla ^{2}u_{E} & =0\\ \frac {\partial ^{2}u_{E}}{\partial r^{2}}+\frac {1}{r}\frac {\partial u_{E}}{\partial r}+\frac {1}{r^{2}}\frac {\partial ^{2}u_{E}}{\partial \theta ^{2}} & =0 \end {align*}
With boundary condition\begin {align*} \left \vert u_{E}\left ( 0,\theta \right ) \right \vert & <\theta \\ u_{E}\left ( a,\theta \right ) & =g\left ( \theta \right ) \\ u_{E}\left ( r,-\pi \right ) & =u_{E}\left ( r,\pi \right ) \\ \frac {\partial u_{E}}{\partial \theta }\left ( r,-\pi \right ) & =\frac {\partial u_{E}}{\partial \theta }\left ( r,\pi \right ) \end {align*}
Let \[ u_{E}\left ( r,\theta \right ) =R\left ( r\right ) \Theta \left ( \theta \right ) \] Where \(R\left ( r\right ) \) is the solution in radial dimension and \(\Theta \left ( \theta \right ) \) is solution in angular dimension. Substituting \(u_{E}\left ( r,\theta \right ) \) in the PDE gives\[ R^{\prime \prime }\Theta +\frac {1}{r}R^{\prime }\Theta +\frac {1}{r^{2}}\Theta ^{\prime \prime }R=0 \] Dividing by \(R\left ( r\right ) \Phi \left ( \theta \right ) \)\begin {align*} \frac {R^{\prime \prime }}{R}+\frac {1}{r}\frac {R^{\prime }}{R}+\frac {1}{r^{2}}\frac {\Theta ^{\prime \prime }}{\Theta } & =0\\ r^{2}\frac {R^{\prime \prime }}{R}+r\frac {R^{\prime }}{R} & =-\frac {\Theta ^{\prime \prime }}{\Theta } \end {align*}
Hence each side is equal to constant, say \(\lambda \) and we obtain\begin {align*} r^{2}\frac {R^{\prime \prime }}{R}+r\frac {R^{\prime }}{R} & =\lambda \\ -\frac {\Theta ^{\prime \prime }}{\Theta } & =\lambda \end {align*}
Or\begin {align} r^{2}R^{\prime \prime }+rR^{\prime }-\lambda R & =0\tag {1}\\ \Theta ^{\prime \prime }+\lambda \Theta & =0 \tag {2} \end {align}
We start with \(\Phi \) ODE. The boundary conditions on (3) are \begin {align*} \Theta \left ( -\pi \right ) & =\Theta \left ( \pi \right ) \\ \frac {\partial \Theta }{\partial \theta }\left ( -\pi \right ) & =\frac {\partial \Theta }{\partial \theta }\left ( \pi \right ) \end {align*}
case \(\lambda =0\) The solution is \(\Phi =c_{1}\theta +c_{2}\). Hence we obtain, from first initial conditions\begin {align*} -\pi c_{1}+c_{2} & =\pi c_{1}+c_{2}\\ c_{1} & =0 \end {align*}
Second boundary conditions just says that \(c_{2}=c_{2}\), so any constant will do. Hence \(\lambda =0\) is an eigenvalue with constant being eigenfunction.
case \(\lambda >0\) The solution is \[ \Theta \left ( \theta \right ) =c_{1}\cos \sqrt {\lambda }\theta +c_{2}\sin \sqrt {\lambda }\theta \] The first boundary conditions gives\begin {align} c_{1}\cos \left ( -\sqrt {\lambda }\pi \right ) +c_{2}\sin \left ( -\sqrt {\lambda }\pi \right ) & =c_{1}\cos \left ( \sqrt {\lambda }\pi \right ) +c_{2}\sin \left ( \sqrt {\lambda }\pi \right ) \nonumber \\ c_{1}\cos \left ( \sqrt {\lambda }\pi \right ) -c_{2}\sin \left ( \sqrt {\lambda }\pi \right ) & =c_{1}\cos \left ( \sqrt {\lambda }\pi \right ) +c_{2}\sin \left ( \sqrt {\lambda }\pi \right ) \nonumber \\ 2c_{2}\sin \left ( \sqrt {\lambda }\pi \right ) & =0\tag {3} \end {align}
From second boundary conditions we obtain\[ \Theta ^{\prime }\left ( \theta \right ) =-\sqrt {\lambda }c_{1}\sin \sqrt {\lambda }\theta +c_{2}\sqrt {\lambda }\cos \sqrt {\lambda }\theta \] Therefore\begin {align} -\sqrt {\lambda }c_{1}\sin \left ( -\sqrt {\lambda }\pi \right ) +c_{2}\sqrt {\lambda }\cos \left ( -\sqrt {\lambda }\pi \right ) & =-\sqrt {\lambda }c_{1}\sin \left ( \sqrt {\lambda }\pi \right ) +c_{2}\sqrt {\lambda }\cos \left ( \sqrt {\lambda }\pi \right ) \nonumber \\ \sqrt {\lambda }c_{1}\sin \left ( \sqrt {\lambda }\pi \right ) +c_{2}\sqrt {\lambda }\cos \left ( \sqrt {\lambda }\pi \right ) & =-\sqrt {\lambda }c_{1}\sin \left ( \sqrt {\lambda }\pi \right ) +c_{2}\sqrt {\lambda }\cos \left ( \sqrt {\lambda }\pi \right ) \nonumber \\ \sqrt {\lambda }c_{1}\sin \left ( \sqrt {\lambda }\pi \right ) & =-\sqrt {\lambda }c_{1}\sin \left ( \sqrt {\lambda }\pi \right ) \nonumber \\ 2c_{1}\sin \left ( \sqrt {\lambda }\pi \right ) & =0\tag {4} \end {align}
Both (3) and (4) are satisfied if\begin {align*} \sqrt {\lambda }\pi & =n\pi \qquad n=1,2,3,\cdots \\ \lambda & =n^{2}\qquad n=1,2,3,\cdots \end {align*}
Therefore\begin {equation} \Theta _{n}\left ( \theta \right ) =\overset {\lambda =0}{\overbrace {A_{0}}}+\sum _{n=1}^{\infty }A_{n}\cos \left ( n\theta \right ) +B_{n}\sin \left ( n\theta \right ) \tag {5} \end {equation} Now we go back to the \(R\) ODE (1) given by \(r^{2}R^{\prime \prime }+rR^{\prime }-\lambda _{n}R=0\) and solve it. This is Euler ODE whose solution is found by substituting \(R\left ( r\right ) =r^{\alpha }\). The solution comes out to be\begin {equation} R_{n}\left ( r\right ) =c_{0}+\sum _{n=1}^{\infty }c_{n}r^{n}\tag {6} \end {equation} Combining (5,6) we now find \(u_{E}\) as\begin {align} u_{E_{n}}\left ( r,\theta \right ) & =R_{n}\left ( r\right ) \Theta _{n}\left ( \theta \right ) \nonumber \\ u_{E}\left ( r,\theta \right ) & =A_{0}+\sum _{n=1}^{\infty }r^{n}\left ( A_{n}\cos \left ( n\theta \right ) +B_{n}\sin \left ( n\theta \right ) \right ) \tag {7} \end {align}
Where \(c_{0}\) was combined with \(A_{0}\). Now the above equilibrium solution needs to satisfy the non-homogenous B.C. \(u_{E}\left ( a,\theta \right ) =g\left ( \theta \right ) \). Using orthogonality on (7) to find \(A_{n},B_{n}\) gives\[ g\left ( \theta \right ) =A_{0}+\sum _{n=1}^{\infty }a^{n}\left ( A_{n}\cos \left ( n\theta \right ) +B_{n}\sin \left ( n\theta \right ) \right ) \] For \(n=0\)\begin {align*} \int _{0}^{2\pi }g\left ( \theta \right ) d\theta & =A_{0}\int _{0}^{2\pi }d\theta \\ A_{0} & =\frac {1}{2\pi }\int _{0}^{2\pi }g\left ( \theta \right ) d\theta \end {align*}
For \(n>0\), applying orthogonality using cosine to find \(A_{n}\) gives\begin {align*} \int _{0}^{2\pi }g\left ( \theta \right ) \cos \left ( n\theta \right ) d\theta & =A_{n}\int _{0}^{2\pi }\cos ^{2}\left ( n\theta \right ) a^{n}d\theta \\ A_{n} & =\frac {1}{\pi }\int _{0}^{2\pi }g\left ( \theta \right ) \cos \left ( n\theta \right ) d\theta \end {align*}
Similarly, applying orthogonality using sin to find \(B_{n}\) gives\[ B_{n}=\frac {1}{\pi }\int _{0}^{2\pi }g\left ( \theta \right ) \sin \left ( n\theta \right ) d\theta \] Therefore, we have found \(u_{E}\left ( r,\theta \right ) \) completely now. It is given by (7)\begin {align} u_{E}\left ( r,\theta \right ) & =A_{0}+\sum _{n=1}^{\infty }r^{n}\left ( A_{n}\cos \left ( n\theta \right ) +B_{n}\sin \left ( n\theta \right ) \right ) \tag {7A}\\ A_{0} & =\frac {1}{2\pi }\int _{0}^{2\pi }g\left ( \theta \right ) d\theta \nonumber \\ A_{n} & =\frac {1}{\pi }\int _{0}^{2\pi }g\left ( \theta \right ) \cos \left ( n\theta \right ) d\theta \nonumber \\ B_{n} & =\frac {1}{\pi }\int _{0}^{2\pi }g\left ( \theta \right ) \sin \left ( n\theta \right ) d\theta \nonumber \end {align}
Now, since \(u\left ( r,\theta ,t\right ) =v\left ( r,\theta ,t\right ) +u_{E}\left ( r,\theta \right ) \), then we need to solve now for \(v\left ( r,\theta ,t\right ) \) with homogeneous boundary conditions \begin {align} v_{t}\left ( r,\theta ,t\right ) & =k\left ( \frac {\partial ^{2}v}{\partial r^{2}}+\frac {1}{r}\frac {\partial v}{\partial r}+\frac {1}{r^{2}}\frac {\partial ^{2}v}{\partial \theta ^{2}}\right ) \tag {8}\\ \left \vert v\left ( 0,\theta ,t\right ) \right \vert & <\theta \nonumber \\ v\left ( a,\theta ,t\right ) & =0\nonumber \\ v\left ( r,-\pi ,t\right ) & =v\left ( r,\pi ,t\right ) \nonumber \\ \frac {\partial v}{\partial \theta }\left ( r,-\pi ,t\right ) & =\frac {\partial v}{\partial \theta }\left ( r,\pi ,t\right ) \nonumber \end {align}
Let \(v\left ( r,\theta ,t\right ) =R\left ( r\right ) \Theta \left ( \theta \right ) T\left ( t\right ) \). Substituting into (8) gives\[ T^{\prime }R\Theta =k\left ( R^{\prime \prime }T\Theta +\frac {1}{r}R^{\prime }T\Theta +\frac {1}{r^{2}}\Theta ^{\prime \prime }RT\right ) \] Dividing by \(R\left ( r\right ) \Theta \left ( \theta \right ) T\left ( t\right ) \neq 0\) gives\[ \frac {1}{k}\frac {T^{\prime }}{T}=\frac {R^{\prime \prime }}{R}+\frac {1}{r}\frac {R^{\prime }}{R}+\frac {1}{r^{2}}\frac {\Theta ^{\prime \prime }}{\Theta }\] Let first separation constant be \(-\lambda \), hence the above becomes\begin {align*} \frac {1}{k}\frac {T^{\prime }}{T} & =-\lambda \\ \frac {R^{\prime \prime }}{R}+\frac {1}{r}\frac {R^{\prime }}{R}+\frac {1}{r^{2}}\frac {\Theta ^{\prime \prime }}{\Theta } & =-\lambda \end {align*}
Or\begin {align*} T^{\prime }+\lambda kT & =0\\ r^{2}\frac {R^{\prime \prime }}{R}+r\frac {R^{\prime }}{R}+r^{2}\lambda & =-\frac {\Theta ^{\prime \prime }}{\Theta } \end {align*}
We now separate the second equation above using \(\mu \) giving\begin {align*} r^{2}\frac {R^{\prime \prime }}{R}+r\frac {R^{\prime }}{R}+r^{2}\lambda & =\mu \\ -\frac {\Theta ^{\prime \prime }}{\Theta } & =\mu \end {align*}
Or\begin {align} R^{\prime \prime }+\frac {1}{r}R^{\prime }+R\left ( \lambda -\frac {\mu }{r^{2}}\right ) & =0\tag {9}\\ \Theta ^{\prime \prime }+\mu \Theta & =0\tag {10} \end {align}
Equation (9) is Sturm-Liouville ODE with boundary conditions \(R\left ( a\right ) =0\) and bounded at \(r=0\) and (10) has periodic boundary conditions as was solved above. The solution to (10) is given in (5) above, no change for this part.\begin {equation} \Theta _{n}\left ( \theta \right ) =\overset {\lambda =0}{\overbrace {\alpha _{0}}}+\sum _{n=1}^{\infty }\alpha _{n}\cos \left ( n\theta \right ) +\beta _{n}\sin \left ( n\theta \right ) \tag {11} \end {equation} Therefore (9) becomes \(R^{\prime \prime }+\frac {1}{r}R^{\prime }+R\left ( \lambda -\frac {n^{2}}{r^{2}}\right ) =0\) with \(n=0,1,2,\cdots \). We found the solution to this Sturm-Liouville before, it is given by\begin {equation} R_{nm}\left ( r\right ) =J_{n}\left ( \sqrt {\lambda _{nm}}r\right ) \qquad n=0,1,2,\cdots ,m=1,2,3,\cdots \tag {12} \end {equation} Where \(\sqrt {\lambda _{nm}}=\frac {a}{z_{nm}}\) where \(a\) is the radius of the disk and \(z_{nm}\) is the \(m^{th}\) zero of the Bessel function of order \(n\). This is found numerically. We now just need to find the time solution from \(T^{\prime }+\lambda _{nm}kT=0\). For This has solution \begin {equation} T_{nm}\left ( t\right ) =e^{-k\lambda _{nm}t}\tag {13} \end {equation} Now we combine (11,12,13) to find solution for \(v\left ( r,\theta ,t\right ) \), and combining constants gives\begin {align} v_{nm}\left ( r,\theta ,t\right ) & =\Theta _{n}\left ( \theta \right ) R_{nm}\left ( r\right ) T_{nm}\left ( t\right ) \nonumber \\ v\left ( r,\theta ,t\right ) & =\alpha _{0,1}J_{0}\left ( \sqrt {\lambda _{0,1}}r\right ) +\sum _{n=1}^{\infty }\sum _{m=1}^{\infty }J_{n}\left ( \sqrt {\lambda _{nm}}r\right ) e^{-k\lambda _{nm}t}\left ( \alpha _{nm}\cos \left ( n\theta \right ) +\beta _{nm}\sin \left ( n\theta \right ) \right ) \nonumber \\ & =\sum _{n=0}^{\infty }\sum _{m=1}^{\infty }J_{n}\left ( \sqrt {\lambda _{nm}}r\right ) \left ( \alpha _{nm}\cos \left ( n\theta \right ) +\beta _{nm}\sin \left ( n\theta \right ) \right ) \tag {14} \end {align}
We now need to find \(\alpha _{0},\alpha _{nn},\beta _{nm}\), which are found from initial conditions on \(v\left ( r,\theta ,0\right ) \) which is given by\begin {align*} v\left ( r,\theta ,0\right ) & =u\left ( r,\theta ,0\right ) -u_{E}\left ( r,\theta \right ) \\ & =f\left ( r,\theta \right ) -u_{E}\left ( r,\theta \right ) \end {align*}
Hence from (14), at \(t=0\)\begin {equation} f\left ( r,\theta \right ) -u_{E}\left ( r,\theta \right ) =\sum _{n=0}^{\infty }\sum _{m=1}^{\infty }J_{n}\left ( \sqrt {\lambda _{nm}}r\right ) \left ( \alpha _{nm}\cos \left ( n\theta \right ) +\beta _{nm}\sin \left ( n\theta \right ) \right ) \tag {15} \end {equation}
For each \(n\), inside the \(m\) sum, \(\cos \left ( n\theta \right ) \) and \(\sin \left ( n\theta \right ) \) will be constant. So we need to apply orthogonality twice in order to remove both sums. Multiplying (15) by \(\cos \left ( n^{\prime }\theta \right ) \) and integrating gives\begin {align*} \int _{-\pi }^{\pi }\left ( f\left ( r,\theta \right ) -u_{E}\left ( r,\theta \right ) \right ) \cos \left ( n^{\prime }\theta \right ) d\theta & =\int _{-\pi }^{\pi }\sum _{n=0}^{\infty }\left ( \sum _{m=1}^{\infty }\alpha _{nn}J_{n}\left ( \sqrt {\lambda _{nm}}r\right ) \right ) \cos \left ( n\theta \right ) \cos \left ( n^{\prime }\theta \right ) d\theta \\ & +\int _{-\pi }^{\pi }\sum _{n=1}^{\infty }\left ( \sum _{m=1}^{\infty }\beta _{nm}J_{n}\left ( \sqrt {\lambda _{nm}}r\right ) \right ) \sin \left ( n\theta \right ) \cos \left ( n^{\prime }\theta \right ) \end {align*}
The second sum in the RHS above goes to zero due to \(\int _{-\pi }^{\pi }\sin \left ( n\theta \right ) \cos \left ( n^{\prime }\theta \right ) d\theta \) and we end up with\[ \int _{-\pi }^{\pi }\left ( f\left ( r,\theta \right ) -u_{E}\left ( r,\theta \right ) \right ) \cos \left ( n\theta \right ) d\theta =\alpha _{nn}\int _{-\pi }^{\pi }\cos ^{2}\left ( n\theta \right ) \sum _{m=1}^{\infty }J_{n}\left ( \sqrt {\lambda _{nm}}r\right ) d\theta \] We now apply orthogonality again, but on Bessel functions and remembering to add the weight \(r\), the above becomes \begin {align*} \int _{0}^{a}\int _{-\pi }^{\pi }\left ( f\left ( r,\theta \right ) -u_{E}\left ( r,\theta \right ) \right ) \cos \left ( n\theta \right ) J_{n}\left ( \sqrt {\lambda _{nm^{\prime }}}r\right ) rd\theta dr & =\alpha _{nn}\int _{0}^{a}\int _{-\pi }^{\pi }\cos ^{2}\left ( n\theta \right ) \sum _{m=1}^{\infty }J_{n}\left ( \sqrt {\lambda _{nm}}r\right ) J_{n}\left ( \sqrt {\lambda _{nm^{\prime }}}r\right ) rd\theta dr\\ & =\alpha _{nn}\int _{0}^{a}\int _{-\pi }^{\pi }\cos ^{2}\left ( n\theta \right ) J_{n}^{2}\left ( \sqrt {\lambda _{nm^{\prime }}}r\right ) rd\theta dr \end {align*}
Therefore\[ \alpha _{nn}=\frac {\int _{0}^{a}\int _{-\pi }^{\pi }\left ( f\left ( r,\theta \right ) -u_{E}\left ( r,\theta \right ) \right ) \cos \left ( n\theta \right ) J_{n}\left ( \sqrt {\lambda _{nm}}r\right ) rd\theta dr}{\int _{0}^{a}\int _{-\pi }^{\pi }\cos ^{2}\left ( n\theta \right ) J_{n}^{2}\left ( \sqrt {\lambda _{nm}}r\right ) rd\theta dr}\qquad n=0,1,2,\cdots ,m=1,2,3,\cdots \] We will repeat the same thing to find \(\beta _{nm}\). The only difference now is to use \(\sin n\theta \). repeating these steps gives\[ \beta _{nm}=\frac {\int _{0}^{a}\int _{-\pi }^{\pi }\left ( f\left ( r,\theta \right ) -u_{E}\left ( r,\theta \right ) \right ) \sin \left ( n\theta \right ) J_{n}\left ( \sqrt {\lambda _{nm}}r\right ) rd\theta dr}{\int _{0}^{a}\int _{-\pi }^{\pi }\sin ^{2}\left ( n\theta \right ) J_{n}^{2}\left ( \sqrt {\lambda _{nm}}r\right ) rd\theta dr}\qquad n=0,1,2,\cdots ,m=1,2,3,\cdots \] This complete the solution.
Summary of solution\begin {align*} u\left ( r,\theta ,t\right ) & =v\left ( r,\theta ,t\right ) +u_{E}\left ( r,\theta \right ) \\ & =u_{E}\left ( r,\theta \right ) +\sum _{n=0}^{\infty }\sum _{m=1}^{\infty }J_{n}\left ( \sqrt {\lambda _{nm}}r\right ) e^{-k\lambda _{nm}t}\left ( \alpha _{nn}\cos \left ( n\theta \right ) +\beta _{nm}\sin \left ( n\theta \right ) \right ) \end {align*}
Where\begin {align*} u_{E}\left ( r,\theta \right ) & =A_{0}+\sum _{n=1}^{\infty }r^{n}\left ( A_{n}\cos \left ( n\theta \right ) +B_{n}\sin \left ( n\theta \right ) \right ) \\ A_{0} & =\frac {1}{2\pi }\int _{0}^{2\pi }g\left ( \theta \right ) d\theta \\ A_{n} & =\frac {1}{\pi }\int _{0}^{2\pi }g\left ( \theta \right ) \cos \left ( n\theta \right ) d\theta \\ B_{n} & =\frac {1}{\pi }\int _{0}^{2\pi }g\left ( \theta \right ) \sin \left ( n\theta \right ) d\theta \end {align*}
And\[ \alpha _{nn}=\frac {\int _{0}^{a}\int _{-\pi }^{\pi }\left ( f\left ( r,\theta \right ) -u_{E}\left ( r,\theta \right ) \right ) \cos \left ( n\theta \right ) J_{n}\left ( \sqrt {\lambda _{nm}}r\right ) rd\theta dr}{\int _{0}^{a}\int _{-\pi }^{\pi }\cos ^{2}\left ( n\theta \right ) J_{n}^{2}\left ( \sqrt {\lambda _{nm}}r\right ) rd\theta dr}\qquad n=0,1,2,\cdots ,m=1,2,3,\cdots \] And\[ \beta _{nm}=\frac {\int _{0}^{a}\int _{-\pi }^{\pi }\left ( f\left ( r,\theta \right ) -u_{E}\left ( r,\theta \right ) \right ) \sin \left ( n\theta \right ) J_{n}\left ( \sqrt {\lambda _{nm}}r\right ) rd\theta dr}{\int _{0}^{a}\int _{-\pi }^{\pi }\sin ^{2}\left ( n\theta \right ) J_{n}^{2}\left ( \sqrt {\lambda _{nm}}r\right ) rd\theta dr}\qquad n=0,1,2,\cdots ,m=1,2,3,\cdots \] Where \(\sqrt {\lambda _{nm}}=\frac {a}{z_{nm}}\) where \(a\) is the radius of the disk and \(z_{nm}\) is the \(m^{th}\) zero of the Bessel function of order \(n\).