\documentclass{article}
%\pagestyle{empty}
\usepackage{graphicx, amssymb,amsmath,amsfonts, psfrag}
\usepackage{amsmath}
\usepackage[usenames,dvipsnames]{xcolor}
\setlength{\textheight}{9in}
\setlength{\topmargin}{0in}
\setlength{\textwidth}{6in}
\setlength{\oddsidemargin}{+0.5in}
\newcommand{\dis}{\displaystyle}
\newcommand{\bd}{\partial}
\newcommand{\x}{\times} 

\newcommand{\green}[1]{\textcolor{ForestGreen}{#1}}
\newcommand{\red}[1]{\textcolor{red}{#1}}
\newcommand{\blue}[1]{\textcolor{blue}{#1}}


\begin{document}


\begin{center}{\bf The Heat Equation for a Square Plate}\end{center}

\vspace{0.5in}

Let $u(x,y,t)$ be the temperature at $(x,y)$ at time $t$.
The two dimensional heat equation in rectangular coordinates is
\[
\alpha^2 (u_{xx} + u_{yy}) = u_t.
\]
For simplicity we will set $\alpha=1$ and work on the $1 \x 1$ square
with corners (0,0), (1,0), (1,1), (0,1). 

Suppose $u(x,y,t) = X(x)Y(y)T(t)$. Then the heat equation becomes
\[
X''YT + XY''T = XYT'.
\]
Divide through by $XYT$ to get
\[
\frac{X''Y + XY''}{XY} = \frac{T'}{T}
\]
Since the left does not change with $t$, both sides are equal to some constant; call it $-\sigma$. Now we have
\[
T'+\sigma T = 0    \:\:\: \& \:\:\: \frac{X''Y + XY''}{XY} = -\sigma.
\]
The equation in $T$ we know how to deal with. For the equation in $X$ and $Y$ we have
\[
\frac{X''}{X} = -\sigma - \frac{Y''}{Y}.
\]
Since one side depends only on $x$ and the other only on $y$, both must be a constant; call it $-\beta$.
Now we have
\[
X'' +\beta X = 0 \:\:\:\&\:\:\: Y'' + (\sigma-\beta) Y = 0.
\]
The equation in $X$ will only have nontrivial solutions if $\beta = n^2\pi^2$ for some integer $n>0$. You can show that $X(x)$ can then be any multiple of $\sin (n\pi x)$. We let 
\[
X_n(x) = \sin (n\pi x).
\]
The equation in $Y$ will only have nontrivial solutions if $\sigma-\beta = m^2\pi^2$ for some integer $m>0$. Thus $\sigma = (n^2+m^2)\pi^2.$ 
We let
\[
Y_m(y) = \sin (m\pi y),
\]
\[
T_{m,n}(t) = e^{-(m^2+n^2)\pi^2t},
\]
and
\[
u_{m,n}(x,y,t) = X_n(x)Y_m(y)T_{m,n}(t).
\]
Then $u_{m,n}$ satisfies the heat equation and the boundary condition
for any integers $n>0$, $m>0$ as would any linear combination.

Now, the theory of Fourier series can be extended to functions of two
variables. If $f(x,y)$ is odd in both variables and periodic with period $2L$ is each variable then the series
\[
\sum_{n=1}^\infty \sum_{m=1}^\infty c_{m,n} \sin \left(\frac{n\pi x}{L} \right)  \sin \left( \frac{m\pi y}{L} \right)
\]
where
\[
c_{m,n} = \frac{4}{L}\int_0^L \int_0^L \sin \left(\frac{n\pi x}{L} \right)  \sin \left(\frac{m\pi y}{L}\right)  f(x,y) \, dx\, dy
\]
converges, almost everywhere, to $f(x,y)$. [See {\em Fourier Series and Boundary Value Problems} by Brown and Churchill.]

Thus, if we are given $u(x,y,0) = f(x,y)$ we just use the above  formulas for the $c_{m,n}$. Then
\[
u(x,y,t) = \sum_{n=1}^\infty \sum_{m=1}^\infty c_{m,n}   \sin (n\pi x)  \sin (m\pi y) e^{-(m^2+n^2)\pi^2t}.
\]

\pagebreak
\begin{center}{\bf Example}\end{center}
\vspace{0.5in}

Suppose $u(x,y,0) = f(x,y) = (x-0.5)^2 + (y-0.5)^2$. Then we can compute the $c_{m,n}$'s and animate the result by the commands below.

\vspace{.2in}
{\tt
\noindent\red{> c := (m,n) -> 4*int(int(sin(n*Pi*x) * sin(m*Pi*y) * \\((x-0.5)\^{}2 + (y-0.5)\^{}2),x=0..1),y=0..1);}

\vspace{.1in}
\noindent\red{> animate(plot3d ,[sum(sum(c(m,n)*sin(n*Pi*x)*sin(m*Pi*y)*\\exp(-(m\^{}2+n\^{}2)*t),m=1..10),n=1..10),x=0..1,y=0..1],t=0..0.3,frames=100);}

}

\vspace{0.51in}

\includegraphics[height=2.5in]{HeatSquareplot3d1.eps}
%\includegraphics[height=1.8in]{HeatSquareplot3d2.eps}
\includegraphics[height=2.5in]{HeatSquareplot3d3.eps}

\vspace{.5in}

\includegraphics[height=2.5in]{HeatSquareplot3d4.eps}
\includegraphics[height=2.5in]{HeatSquareplot3d5.eps}

\vspace{.2in}

See the course website to view the animation.

\vfill
\pagebreak
\begin{center}{\bf Discussion of Heat Equation on a Metal Disk}\end{center}
\vspace{0.5in}


We consider a metal disk of radius 1 with temperature held to 0 degrees around the circumference. The initial temperature is given by $f(r,\theta)$. Let $u(r,\theta, t)$ be the temperature at $(r,\theta)$ at time $t$. Thus,
\[
u(1,\theta, t) = 0\:\:\:\&\:\:\: u(r,\theta,0) = f(r, \theta).
\]
Now, the heat equation in rectangular coordinates in 
\[
\alpha^2 (u_{xx} + u_{yy}) = u_t.
\]
Our first task is to translate this into polar coordinates. I'll show one of the steps involved and then give the final result. The main tool is the multivariable chain rule.
\[
\frac{\bd u(r,\theta,t)}{\bd x} = \frac{\bd u}{\bd r}\frac{\bd r}{\bd x} +\frac{\bd u}{\bd \theta}\frac{\bd \theta}{\bd x}
+\frac{\bd u}{\bd t}\frac{\bd t}{\bd x}.
\]
Since $t$ is independent of $x$ we know $\bd_x t =0$, so we can drop the last term. For polar coordinates we know
\[
r = \sqrt{x^2 + y^2} \:\:\:\&\:\:\: \theta = \arctan (y/x).
\]
Therefore
\[
\frac{\bd r}{\bd x} = \frac{2x}{2\sqrt{x^2+y^2}} = \frac{r\cos \theta}{r} = \cos \theta,
\]
and
\[
\frac{\bd \theta}{\bd x} = \frac{1}{1 + (y/x)^2}\left(\frac{y}{x}\right)' = \frac{1}{1 + (y/x)^2}\frac{-y}{x^2} = \frac{-y}{x^2+y^2} = \frac{-\sin \theta}{r}.
\]
Hence,
\[
u_x = u_r \cos \theta - \frac{1}{r} u_\theta \sin \theta.
\]
One repeats  this to get $u_{xx}$ and similarly for $u_{yy}$. After simplifying the result is
\[
\alpha^2 (u_{rr} + \frac{1}{r} u_r + \frac{1}{r^2} u_{\theta\theta} ) = u_t.
\]
 
\vfill
\pagebreak

Now we try to find solutions. As before we suppose $u(r,\theta,t) = R(r) \Theta(\theta) T(t)$. Plugging into the heat equation gives.
\[
\alpha^2\left(R''\Theta T + \frac{1}{r}R'\Theta T + \frac{1}{r^2}R\Theta''T\right) = R\Theta T'.
\]
Divide through by $\alpha^2 R\Theta T$ to get
\[
\frac{R''}{R}  + \frac{1}{r}\frac{R'}{R}  +   \frac{1}{r^2}\frac{\Theta''}{\Theta} =  \frac{1}{\alpha^2} \frac{T'}{T}.
\]
Since the right side depends only on $t$ and the left side is independent of $t$, both sides are constant; call it $-\sigma$.
We get
\[
r^2\frac{R''}{R}  + r\frac{R'}{R}+r^2\sigma =-  \frac{\Theta''}{\Theta}  \:\:\:\&\:\:\: T' + \alpha^2 \sigma T = 0
\]
Thus, there is a constant $-\beta$ such that 
\[
r^2\frac{R''}{R}  + r\frac{R'}{R} + r^2\sigma =-  \frac{\Theta''}{\Theta}  =\beta.
\]
Thus we get
\[
r^2R''  + rR' + (r^2\sigma-\beta)R = 0   \:\:\:\&\:\:\: \Theta'' + \beta \Theta = 0.
\]
Thus, we have there separate ODE's. 

We shall make a simplifying assumption, that we will hopefully drop later. But for now
assume $u(r,\theta,t)$ is independent of $\theta$. This means $\Theta' = 0$ and hence 
so does $\Theta''$. Thus $\beta \Theta(\theta) = 0$ So, unless the temperature is always 
zero, will require $\beta = 0$. The $R$ equation becomes
\[
r^2R''  + rR' + \sigma r^2R = 0.
\]
Let $\gamma^2 = \sigma$ and $p = \gamma r$. Also define $P(p) = R(p/\gamma) = R(r)$.
Then 
\[
r \frac{d R}{dr} = \frac{p}{\gamma} \frac{d R(p\gamma)}{dr} = \frac{p}{\gamma}\frac{dP}{dp}\frac{dp}{dr} = pP'(p).
\]
Similarly, $r^2 R(r) = p^2 P(p)$. Thus our equation for $R$ becomes 
\[
p^2P'' + p P' + p^2 P = 0.
\]
This is Bessel's Equation of Order Zero. It is solved in Section 5.7. The general solution is 
\[
P(p) = C_1 J_0(p) + C_2 Y_0(p).
\]
Switching back to $R$ and $r$ we get
\[
R(r) =  C_1 J_0(\gamma r) + C_2 Y_0(\gamma r).
\]
The functions $J_0$ and $Y_0$ are called Bessel functions of the first and second kind, respectively, both of order zero.
They are determined by certain power series.  

\vfill
\pagebreak
\begin{center}{\bf Bessel's Equations}\end{center}

Below are graphs of $J_0(x)$ and $Y_0(x)$. 

\begin{center}
\includegraphics[height=2.3in]{BesselJYgraphsplot1.eps} 
\hspace{.3in}
\includegraphics[height=2.3in]{BesselJYgraphsplot2.eps} 
\end{center}

Next we consider the boundary conditions. First, $u(1,\theta,0) = 0$ implies $R(1) = 0$. 
But this is a second ``boundary'' condition that is implicit: the temperature is bounded.
In particular $|u(0,\theta,t)| < \infty$. Thus, $|R(0)| < \infty$. But, $Y_0(\gamma r)$ is unbounded 
as $r\to 0$. Hence $C_2 = 0$.

Now, $J_0(x)$ has infinitely many zeros. Enumerate them $\gamma_1$, $\gamma_2$, $\gamma_3$, etc. There are tables that list these or 
they can be generated by a computer. Here are the first twenty as generated by Maple.

\vspace{.2in}
{\tt
\noindent\red{> evalf(BesselJZeros(0,1..20));}}

\vspace{.2in}

{\em
\noindent\blue{2.404825558, 5.520078110, 8.653727913, 11.79153444, 14.93091771,\\ 
18.07106397, 21.21163663, 24.35247153, 27.49347913, 30.63460647,\\
 33.77582021, 36.91709835, 40.05842576, 43.19979171, 46.34118837, \\
49.48260990, 52.62405184, 55.76551076, 58.90698393, 62.04846919
}}


%\pagebreak
\vspace{.2in}
Then we let
\[
R_n(r) = J_0 (\gamma_n r).
\]
Each $R_n(r)$  satisfies the boundary conditions. Let 
\[
T_n(t) = e^{-\gamma_n^2\alpha^2 t} \:\:\:\&\:\:\: u_n(r,t) = X_n(r)T_n(t).
\]

Then each $u_n$ satisfies the heat equation and the boundary conditions and so too 
would any linear combination. (We dropped $\theta$ since we are still assuming $u$ is independent of $\theta$.) 
It can be shown that this holds for 
\[
u(r,t) = \sum_{n=1}^\infty c_n J_0(\gamma_n r) e^{-\gamma_n^2\alpha^2 t}
\]
if the $c_n$'s are such that the series convergences. Thus if we can find $c_n$'s such that
\[
\sum_{n=1}^\infty c_n J_0(\gamma_n r) = u(r,0) = f(r)
\]
we are done! It looks almost like the Fourier expansion with $ J_0(\gamma_n r)$ in stead of $\sin (n\pi r)$ or $\cos(n\pi r)$.
In fact, Section 11.4 shows that this can be done and gives a formula for $c_n$.
\[
c_n = \frac{\int_0^1 r f(r) J_0(\gamma_n r)\, dr}{\int_0^1 r J_0^2 (\gamma_n r) \, dr}.
\]

\begin{center}{\bf Example}\end{center}


Suppose $f(r) = 1-r^2$. We approximate $f(r)$ with a series of Bessel functions. Then we plot a few partial sums 
to see how good the approximation is. We define $f(r)$, {\tt g(n)} is $\gamma_n$ and {\tt c(n)} is $c_m$.

\vspace{.2in}
{\tt
\noindent\red{> f := r -> 1-r\^{}2;}
 
\noindent\red{> g := n -> evalf(BesselJZeros(0,n)); \# g(n) is n-th zero of Bessel function Jo.} 
 
\noindent\red{> c := n -> (int(r*f(r)*BesselJ(0,g(n)*r) ,r=0..1))/(int(r*BesselJ(0,g(n)*r)\^{}2 ,r=0..1));}
 


\noindent\red{> N:=1:}

\noindent\red{> plot([sum(c(n)*BesselJ(0,g(n)*r),n=1..N),1-r\^{}2 ] ,r=0..1,\\thickness=2,color=[red,blue]);}
}

\begin{center}
\includegraphics[height=1.9in]{HeatDiskplot1.eps} 
\end{center}

We repeat for $N=2, 3, \& 5$. 

\begin{center}
\includegraphics[height=1.9in]{HeatDiskplot2.eps} 
\includegraphics[height=1.9in]{HeatDiskplot3.eps}
\includegraphics[height=1.9in]{HeatDiskplot4.eps}
\end{center}

Looks good! But let's plot the graph out to $r=5$; $N$ is still 5.

\begin{center}
\includegraphics[height=1.9in]{HeatDiskplot5.eps}
\end{center}

Yikes! Fortunately, we only care about $r\in [0,1]$.
Next we compute $u(r,t)$ and plot it for several values of $t$. I used $N=10$ and then did plots for $t=0.01, 0.02, 0.1, 0.2, 0.3 \& 1.0$. Only the command for $t=0.01$ is shown.

\vspace{.2in}

{\tt
\noindent\red{> N:=10:t:=0.01:}

\noindent\red{> plot(sum(c(n)*BesselJ(0,g(n)*r)*exp(-g(n)\^{}2*t),n=1..N) ,r=0..1,\\thickness=2,color=red,view=0..1);}
 }

\begin{center}
\includegraphics[height=1.9in]{HeatDiskplot6.eps}
\includegraphics[height=1.9in]{HeatDiskplot7.eps}
\includegraphics[height=1.9in]{HeatDiskplot8.eps}

\includegraphics[height=1.9in]{HeatDiskplot9.eps}
\includegraphics[height=1.9in]{HeatDiskplot10.eps}
\includegraphics[height=1.9in]{HeatDiskplot11.eps}

\includegraphics[height=1.9in]{HeatDiskplot12.eps}
\end{center}


\begin{center}{\bf Extra Credit}\end{center}

Suppose $f(r) = r(1-r^2)$. The temperature at the origin starts at 0 and in the limit as $t\to\infty$ it will tend towards 0. What is the maximum temperature of the origin and when will it occur? Approximate as best you can. 

\pagebreak
\begin{center}{\bf Example}\end{center}

Suppose $\alpha=1$ and $f(r) = r(1-r^2)\sin(3\theta)^2$. Here is a graph.

\vspace{.2in}{\tt 
\noindent\red{> plot3d(sqrt(x\^{}2+y\^{}2)*(1-x\^{}2-y\^{}2)*sin(3*arctan(y,x))\^{}2, x=-1..1,y=-1..1,\\view=0..0.5,color=sqrt(x\^{}2+y\^{}2)*(1-x\^{}2-y\^{}2)*sin(3*arctan(y,x))\^{}2,numpoints=10000);} 
}


\begin{center}
\includegraphics[height=4in]{HeatDisk2plot3d2a.eps}
\end{center}

We revisit the separated equations. 


\[r^2R''  + rR' + (r^2\sigma-\beta)R = 0   \:\:\:\:\:\: \Theta'' + \beta \Theta = 0 \]\[ T' + \sigma T = 0\]

The implied boundary condition on $\Theta$ is just $\Theta(0) = \Theta(2\pi)$.  For $\beta<0$ only the
trivial solution is periodic. For $\beta = 0$, any constant solution will do. For $\beta>0$ the general
solution is 
\[
\Theta(\theta) = C_1 \sin \sqrt{\beta}\theta + C_2 \cos  \sqrt{\beta}\theta. 
\]
Therefore, $\beta = n^2$ for $n>0$. Thus, using any linear combination of $1$, $\sin n\theta$, and $\cos n\theta$
will do. 
Let
\[
\Theta_0(\theta) = C_0 \:\:\:\&\:\:\: \Theta_n(\theta) = a_n \sin n\theta + b_n \cos n\theta,
\]
for $n=1, 2, 3, ....$. Then for each $n=0,1,2,3...$ we have that 
$\Theta_n(\theta)$ solves $\Theta'' + n^2\Theta = 0$ and $Theta_n(0) = 
Theta_n(2\pi)$.

Let $\gamma^2 = \sigma$, $\gamma>0$, and consider 
\[
r^2R''  + rR' + (\gamma^2r^2-n^2)R = 0.
\]
The change of variable $p=\gamma r$ transforms this to
\[
p^2 P'' pP' +(p^2-n^2)P = 0.
\]
This is Bessel's equation of order $n$. It is general solution
takes the form 
\[
P(p) = C_1 J_n(p) + C_2 Y_n(p).
\]
This is equivalent to 
\[
R(r) = C_1 J_n(\gamma r) + C_2 Y_n(\gamma r).
\]
Here $J_n$ is Bessel's function of the first kind of order $n$, and $Y_n$ is Bessel's function of the second kind of order $n$, as you might have guessed. They are defined have certain power series (see Section ). Again, $Y_n$ is unbounded at 0 and so $C_2=0$. Again $J_n$ has infinitely many zeros. We denote them by $\gamma_{n,m}$, $m=1, 2, 3, ...$. Their numerical values can be found in tables or via a computer. 
For each $n=0,1,2,3...$ and $m=1,2,3...$ let
\[
R_{n,m}(r) = J_n(\gamma_{n,m}r).
\]
Then $R_{n,m}(1) = 0$ and $R_{n,m}(r)$ solves $r^2R'' +r R' +(\gamma_{n,m}^2r^2 - n^2)R = 0$. 

Define
\[
T_{n,m} = e^{-\gamma_{n,m}^2 t}.
\]
Finally, define
\[
u_{n,m}(r,\theta,t) = R_{n,m}(r)T_{n,m}(t)\sin (n\pi) \:\:\:\&\:\:\: v_{n,m}(r,\theta,t) = R_{n,m}(r)T_{n,m}(t)\cos (n\pi),
\]
for integers $n\geq 0$, $m>0$. 
Each satisfies the heat equation in polar coordinates, will be 0 
when $r=1$, be bounded and periodic in $\theta$ with least period $2\pi/n$. 
So too would any linear combination. If $\{ a_{n,m} \}$ and $\{ b_{n,m}\}$ are chosen so that

\[
u(r,\theta,t) = \sum_{n=0}^\infty \sum_{m=1}^\infty a_{n,m}  R_{n,m}(r)T_{n,m}(t)\sin (n\pi) - b_{n,m} R_{n,m}(r)T_{n,m}(t)\cos (n\pi)
\]
converges, then its limit will satisfy the heat equation, be 0 when $r=1$, bounded and periodic in $\theta$ with period $2\pi$. All that if left is to find values for the  $a_{n,m}$ and  $b_{n,m}$ such that 
\[
u(r,\theta,0) = f(r,\theta) = r(1-r^2)\sin^2(3\theta).
\]
In our case we can use that 
\[
\sin^2(3\theta) = \frac{1}{2} - \frac{1}{2} \cos (6\theta).
\]
Thus 
\[
u(r,\theta,0) = r(1-r^2)\frac{1}{2} - r(1-r^2)\frac{1}{2}\cos(6\theta).
\]
We need for 
\[
\sum_{m=1}^\infty b_{0,m}  R_{0,m}(r) = \sum_{m=1}^\infty b_{0,m}J_0(\gamma_{0,m}r) = r(1-r^2)
\]
and
\[
\sum_{m=1}^\infty b_{6,m}  R_{6,m}(r) = \sum_{m=1}^\infty b_{6,m} J_6(\gamma_{6,m}r) = r(1-r^2).
\]

Each family of functions $\{ J_n(\gamma_{n,m} \}_{m=1}^\infty $ satisfies 
the orthogonality conditions needed to do series approximations. We get
\[
b_{0,m} = \frac{\int_0^1 r^2(1-r^2) J_0(\gamma_{0,m}r)  \, dr }{\int_0^1  r J_0^2 (\gamma_{0,m} r) \, dr} 
\] 
and
\[
b_{6,m} = \frac{\int_0^1 r^2(1-r^2) J_6(\gamma_{6,m}r)  \, dr }{\int_0^1  r J_6^2 (\gamma_{6,m} r) \, dr}. 
\] 

Next we plot $u(r,\theta,0)$ and compare the graph for $f(r,\theta)$ above.

{\tt
\noindent\red{> f := r -> r*(1-r\^{}2);}

\noindent\red{> g := (n,m) -> evalf(BesselJZeros(n,m));}

\noindent\red{> b0 := b0 := m -> (int(r*f(r)*BesselJ(0,g(0,m)*r)\\,r=0..1))/(int(r*BesselJ(0,g(0,m)*r)\^{}2 ,r=0..1));}

\noindent\red{> b6 := m -> (int(r*f(r)*BesselJ(6,g(6,m)*r)\\,r=0..1))/(int(r*BesselJ(6,g(6,m)*r)\^{}2 ,r=0..1));}

\noindent\red{> U := (r,theta) -> 1/2*sum( b0(m)*BesselJ(0,g(0,m)*r) ,m=1..M) - \\1/2*cos(6*theta)*sum(b6(m)*BesselJ(6,g(6,m)*r), m=1..M);}

\noindent\red{> M:=10:}

\noindent\red{> plot3d(U(sqrt(x\^{}2+y\^{}2),arctan(y,x)),x=-1..1,y=-1..1,view=0..0.5,\\color=U(sqrt(x\^{}2+y\^{}2),arctan(y,x)),numpoints=1000);}
}

\begin{center}
\includegraphics[height=4in]{HeatDisk2plot3d3.eps}
\end{center}

[I have no idea why Maple used a different color scheme.]

Now we put all the pieces together to get
\[
u(r,\theta,t) =  \frac{1}{2} \sum_{m=1}^\infty a_{0,m} J_0(\gamma_{n,m}r) e^{-\gamma_{0,m}^2 t } - \frac{1}{2}\cos(6\theta)  \sum_{m=1}^\infty a_{6,m} J_6(\gamma_{n,m}r) e^{-\gamma_{6,m}^2 t }.
\]


We plot this for a few values of $t$.

\begin{center}
\psfrag{t=0.01}{\hspace{-0.2in}$t=0.01$}
\includegraphics[width=3in]{HeatDisk2plot3d5.eps}

\psfrag{t=0.02}{\hspace{-0.2in}$t=0.02$}
\includegraphics[width=3in]{HeatDisk2plot3d6.eps}

\psfrag{t=0.05}{\hspace{-0.2in}$t=0.05$}
\includegraphics[width=3in]{HeatDisk2plot3d7.eps}
\end{center}

[Maple is re-calibrating the color scale each time. I do not know how to fix this.]

\pagebreak

Note: We took advantage of the fact the given $f(r,\theta)$ was separable. This is not always the case. 
 

General case: See {\tt http://faculty.fiu.edu/~meziani/MAP4401.html},
``NOTE12: Fourier-Bessel Series and BVP in Cylindrical Coordinates'' by Hamid Meziani, 2016. 

We have...

\[
u(r,\theta,t) = \sum_{n=0}^\infty \sum_{m=1}^\infty a_{n,m}  R_{n,m}(r)T_{n,m}(t)\sin (n\pi) - b_{n,m} R_{n,m}(r)T_{n,m}(t)\cos (n\pi) 
\]
\[
=  \sum_{n=0}^\infty \sum_{m=1}^\infty a_{n,m}  J_n(\gamma_{n,m}r)e^{-\gamma_{n,m}^2 t}\sin (n\pi) - b_{n,m} J_n(\gamma_{n,m}r)e^{-\gamma_{n,m}^2 t}\cos (n\pi). 
\]
At time $t=0$ we need
\[
u(r,\theta,0) = \sum_{n=0}^\infty \sum_{m=1}^\infty a_{n,m}  J_n(\gamma_{n,m}r)\cos (n\pi) - b_{n,m} J_n(\gamma_{n,m}r)\sin (n\pi) = f(r,\theta).
\]

For each value of $r\in [0,1]$ the function $f(r,\theta)$ is periodic with period $2\pi$ and so has a Fourier series. 
We compute the coefficients as normal, but regard them as functions of $r$. 

\[
f(r,\theta) = \frac{A_0(r)}{2} + \sum_{n=1}^\infty A_n(r) \cos (n\theta) + B_n(r) \sin (n\theta)
\]
where
\[
A_n(r) = \frac{1}{\pi}\int_{-\pi}^{\pi} f(r,\theta) \cos (n\theta) \, d\theta,
\]
and
\[
B_n(r) = \frac{1}{\pi}\int_{-\pi}^{\pi} f(r,\theta) \sin (n\theta) \, d\theta,
\]
for $n=0,1,2,3,...$ (of course $B_0 =0$ so we ignore it). 

Now each $A_n(r)$ and $B_n(r)$ has a Bessel series in terms of any $J_k$ so choose $k=n$. We have
\[
A_n(r) = \sum_{j=1}^\infty A_{nm} J_n(\gamma_{nm} r) \:\:\:\&\:\:\: B_n(r) =  \sum_{j=1}^\infty B_{nm} J_n(\gamma_{nm} r),
\]
where 
\[
A_{nm} = \frac{\int_0^1 r A_n(r) J_n(\gamma_{nm} r) \, dr}{\int_0^1 r  J_n^2(\gamma_{nm} r) \, dr} 
 \:\:\:\&\:\:\:
B_{nm} = \frac{\int_0^1 r B_n(r) J_n(\gamma_{nm} r) \, dr}{\int_0^1 r  J_n^2(\gamma_{nm} r) \, dr}. 
\]


\begin{center}{\bf Example}\end{center}


Let $f(r,\theta) = r(1-r)e^{r\sin \theta}$

\begin{center}{\bf Discontinuous Example}\end{center}

Let 
\[
f(r,\theta) = \left\{ \begin{array}{ll}
1 & 0 \leq \theta \leq \pi \& r < \sin \theta \\
0 & \mbox{ otherwise}
\end{array}\right.
\]



\end{document}
