\documentclass{article} %\pagestyle{empty} \usepackage{graphicx, amssymb,amsmath,amsfonts, psfrag} \usepackage{amsmath} %%% Add Color \usepackage[usenames,dvipsnames]{xcolor} \newcommand{\green}[1]{\textcolor{ForestGreen}{#1}} \newcommand{\red}[1]{\textcolor{red}{#1}} \newcommand{\blue}[1]{\textcolor{blue}{#1}} \setlength{\textheight}{9in} \setlength{\topmargin}{0.0in} \setlength{\textwidth}{6in} \setlength{\oddsidemargin}{+0.5in} \newcommand{\dis}{\displaystyle} \newcommand{\bd}{\partial} \newcommand{\x}{\times} \begin{document} \begin{center}{\bf Lecture Notes for Ch 10\\Fourier Series and Partial Differential Equations}\footnote{\copyright Michael C. Sullivan, 2017.}\end{center} \section{Introduction} Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations \tableofcontents \section{Review of Even and Odd Functions} \begin{center} \fbox{ Even: $f(-x) = f(x)$ \hspace{0.6in} Odd: $f(-x)=-f(x)$.} \end{center} \begin{itemize} \item [Even Examples:] $5x^6-3x^4+2$, $\frac{x^4+1}{x^2-7}$, $|x|$, $\cos(x)$, $\cos(x^5)$, $\sin^4(x)$. \item [Odd Examples:] $x^{1/3}$, $x^7-6x^3$, $x|x|$, $\sin(x)$, $\sin^5(x^7)$. \item [Neither:] $x^2+x$, $x+\cos(x)$, $\frac{1}{x+1}$. \begin{center} \psfrag{E}{Three Even Functions} \psfrag{O}{Three Odd Functions} \hspace*{-1in}\includegraphics[height=1.75in]{Figs/evenodd.eps} \end{center} \item [Problem:] Let $f(x)$ be an odd function and suppose it is defined at $x=0$. What is $f(0)$? Prove this! \item [Problem:] Find a function with domain the whole real line that is both even and odd. (There is only one correct answer.) \item [Fact:] If $f(x)$ is even, then $\dis \int_{-a}^a f(x)\,dx = 2\int_0^a f(x)\,dx$. \item [Fact:] If $f(x)$ is odd, then $\dis \int_{-a}^a f(x)\,dx = 0$. Draw pictures to see intuitively why these two facts hold. \item [Fact:] If $f(x)$ is even and differentiable, then $f'(x)$ is odd, and vice versa. This can be proved with the Chain Rule. Suppose $f(x)$ is even and differentiable. Then $(f(-x))' = f'(-x)(-1)$. But $f(-x)= f(x)$, so $(f(-x))'= (f(x))' = f'(x)$. Thus, $f'(-x)(-1) = f'(x)$, or $f'(-x) = -f'(x)$. You can also draw pictures of tangent lines to curves to see this intuitively. %\end{itemize} \item [Exercises:] Suppose that $f(x)$ and $g(x)$ are even and the $h(x)$ and $k(x)$ are odd. Then show that: \begin{center} \begin{tabular}{ll} (a) $f(x)g(x)$ is even. & (f) $f(g(x))$ is even. \\ (b) $f(x)h(x)$ is odd. & (g) $f(h(x))$ is even. \\ (c) $h(x)k(x)$ is even. & (h) $h(f(x))$ is even. \\ (d) $f(x)+g(x)$ is even. & (i) $h(k(x))$ is odd. \\ (e) $h(x)+k(x)$ is odd. & (j) $f(x) + h(x)$ need not be even or odd. \end{tabular} \end{center} \item [Example:] Let $p(x) = f(x)h(x)$. Then $p(-x) = f(-x)h(-x)= f(x)(-h(x)) = -f(x)h(x) = -p(x)$. Hence we have an odd function. \end{itemize} %\vfill %\vfill %\pagebreak \section{Some Damn Useful Integral Formulas} \begin{enumerate} \item $\dis \int_{-L}^L \cos \left( \frac{m\pi x}{L}\right) \cos \left(\frac{n\pi x}{L} \right) \,\, dx = \left\{ \begin{array}{rrr} 0 & \mbox{ for } & m\neq n, \\ L & \mbox{ for } & m = n. \end{array}\right. $ \vspace{.25in} \item $\dis \int_{-L}^L \cos \left( \frac{m\pi x}{L}\right) \sin \left(\frac{n\pi x}{L} \right) \,\, dx = 0,\: \mbox{for all integers $m$ and $n$} $ \vspace{.25in} \item $\dis \int_{-L}^L \sin \left( \frac{m\pi x}{L}\right) \sin \left(\frac{n\pi x}{L} \right) \,\, dx = \left\{ \begin{array}{rrr} 0 & \mbox{ for } & m\neq n, \\ L & \mbox{ for } & m = n. \end{array}\right. $ \end{enumerate} \vspace{.25in} Below and to the left are the overlaid plots of $\cos 4\pi x$ and $\cos 2\pi x$. Below and to the right is the graph of their product. Study this. Make some similar plots on your own until the formulas above make sense. \begin{center} \includegraphics[height=1.75in]{Figs/cos4xandcos2x.eps} \hspace{0.5in} \includegraphics[height=1.75in]{Figs/cos4xXcos2x.eps} \end{center} We will prove only a special case of \#1. Let $L=\pi$, $m=3$ and $n=2$. The proof uses the trig identity \[ \cos (\theta + \phi) = \cos \theta \cos \phi - \sin \theta \sin \phi \] and its corollary \[ \cos (\theta - \phi) = \cos \theta \cos \phi + \sin \theta \sin \phi \] Thus, \[ \cos \theta \cos \phi = \frac{1}{2}\left( \cos (\theta + \phi) + \cos (\theta - \phi) \right). \] Applying this to our case gives \[ \cos 3x \cos 2x = \frac{1}{2}\left( \cos 5x + \cos x\right). \] Thus, \[ \int_{-\pi}^{\pi} \cos 3x \cos 2x \, dx = \frac{1}{2}\int_{-\pi}^{\pi} \cos 5x + \cos x \, dx = 0 + 0 = 0. \] Let's consider one more special case: $L=\pi$, $m=n=2$. Then \[ \int_{-\pi}^{\pi} \cos^2 2x \, dx = \frac{1}{2}\int_{-\pi}^{\pi} \cos 4x + \cos 0 \,dx = \frac{0+2\pi}{2} = \pi. \] From these ideas you should be able to derive the three integrals formulas. If you have had linear algebra, you might notice that the integral of a product of two functions is a kind of inner product and so the cosine and sine functions used above are mutually orthogonal. %\pagebreak \section{Definition of Fourier Series} Let $f(x)$ be a piecewise continuous periodic function with period $2L$.\footnote{Technical note: It is also required that $f(x)$ is bounded and that in each period of $f(x)$ there are only finitely many extrema.} The Fourier Series of $f(x)$ is, \[ \frac{a_0}{2} + \sum_{n=1}^\infty a_n \cos \left( \frac{n\pi x}{L} \right) + \sum_{n=1}^\infty b_n \sin \left( \frac{n\pi x}{L} \right), \] where \[ a_n = \frac{1}{L} \int_{-L}^L \cos \left( \frac{n\pi x}{L} \right) f(x) \, dx, \:\:\: n = 0,1,2,3,\dots, \] and \[ b_n = \frac{1}{L} \int_{-L}^L \sin \left( \frac{n\pi x}{L} \right) f(x) \, dx, \:\:\: n = 1,2,3,\dots. \] {\bf Theorem.} The Fourier series of $f(x)$ converges to $f(x)$ if $f$ is continuous at $x$. If $f$ is discontinuous at $x$ then this must be a jump discontinuity. Let \[ f(x^+) = \lim_{c \to x^+} f(c) \hspace{.4in}\mbox{and} \hspace{.4in} f(x^-) = \lim_{c \to x^-} f(c). \] Then the Fourier series of $f$ converges to the average of these two limits, \[ \frac{f(x^+) + f(x^-)}{2}. \] If $f$ is an even function then $b_n = 0$, for $n=1, 2, 3, \dots.$ If $f$ is an odd function then $a_n = 0$, for $n=0, 1, 2, 3, \dots.$ In all cases $a_0/2$ is the average value of $f(x)$ over one period. \vspace{.2in} We will verify parts of this theorem. The full proof is covered in Math 407. \subsection{We verify the formula for $\mathbf a_0$ and show it is the average value of $\mathbf f(x)$.} Let $c_n = \cos \left(\frac{n\pi x}{L}\right)$ and $s_n = \sin \left(\frac{n\pi x}{L}\right)$. \vspace{,25in} Suppose $$f(x) = \frac{a_0}{2} + \sum_{n=1}^\infty a_n c_n + \sum_{n=1}^\infty b_n s_n.$$ We will show that \[ a_0 = \frac{1}{L} \int_{-L}^L \cos (0) f(x) \, dx. \] \begin{eqnarray*} \frac{1}{L} \int_{-L}^L \cos (0) f(x) \, dx &=& \frac{1}{L} \int_{-L}^L 1\cdot f(x) \, dx \\ && \\ && \\ &=& \frac{1}{L} \int_{-L}^L \left( \frac{a_0}{2} + \sum_{n=1}^\infty a_n c_n + \sum_{n=1}^\infty b_n s_n\right)\, dx \\ && \\ && \\ &=& \frac{1}{L}\int_{-L}^L \frac{a_0}{2}\, dx + \frac{1}{L} \sum_{n=1}^\infty a_n \int_{-L}^L c_n \, dx + \frac{1}{L} \sum_{n=1}^\infty b_n \int_{-L}^L s_n \, dx\\ && \\ && \\ &=& a_0 + \frac{1}{L} \sum_{n=1}^\infty a_n\cdot 0 + \frac{1}{L} \sum_{n=1}^\infty b_n \cdot 0 = a_0 \end{eqnarray*} since, \vspace{.2in} $\dis \int_{-L}^L \cos \left(\frac{n\pi x}{L} \right) \, dx = \frac{L}{n\pi} \sin \left( \frac{n\pi x}{L} \right) \bigg|_{-L}^{L} = 0 - 0 = 0, $ \vspace{.2in} $\dis \int_{-L}^L \sin \left(\frac{n\pi x}{L} \right) \, dx = \frac{-L}{n\pi} \cos \left( \frac{n\pi x}{L} \right) \bigg|_{-L}^{L} = \frac{-L}{n\pi} \left( \cos(n\pi) - \cos(-n\pi) \right) = 0, $ \vspace{.2in} and \vspace{.2in} $\dis \frac{1}{L} \int_{-L}^L \frac{a_0}{2}\, dx = \frac{a_0}{2L} \int_{-L}^L 1 \, dx = \frac{a_0}{2L} \cdot 2L = a_0. $ \vspace{.2in} Finally, recall that the average value of $f(x)$ over a cycle is $\dis \frac{1}{2L} \int_{-L}^L f(x) \, dx$. Thus, $a_0/2$ is the average value of $f(x)$ over a cycle. %\vfill %\pagebreak \subsection{We check the formula for $\mathbf a_7$, \\but the method is the same for all n $\mathbf \geq$ 1.} \vspace{.2in} \begin{eqnarray*} a_7 &\stackrel{?}{=} & \frac{1}{L} \int_{-L}^L \cos \left( \frac{7\pi x}{L}\right) f(x) \, dx \\ &&\\ &=& \frac{1}{L} \int_{-L}^L c_7\left( \frac{a_0}{2} + \sum_{n=1}^\infty a_n c_n + \sum_{n=1}^\infty b_n s_n\right) \, dx \\ &&\\ &=& \frac{a_0}{2L} \int_{-L}^L c_7 \, dx + \frac{1}{L} \sum_{n=1}^\infty a_n \int_{-L}^L c_n c_7 \, dx + \frac{1}{L} \sum_{n=1}^\infty b_n \int_{-L}^L s_n c_7 \, dx \\ &&\\ &=& \frac{a_0}{2L}\cdot 0 + \frac{1}{L} (0+0+0+0+0+0+a_7L+0+0+0+\cdots ) + \frac{1}{L}(0+0+0+\cdots ) \\ &&\\ &=& a_7. \end{eqnarray*} \vspace{0.5in} %\pagebreak \subsection{Next we check $\mathbf b_4$}. \vspace{.2in} \begin{eqnarray*} b_4 &\stackrel{?}{=} & \frac{1}{L} \int_{-L}^L \sin \left( \frac{4\pi x}{L}\right) f(x) \, dx \\ &&\\ &=& \frac{1}{L} \int_{-L}^L s_4\left( \frac{a_0}{2} + \sum_{n=1}^\infty a_n c_n + \sum_{n=1}^\infty b_n s_n\right) \, dx \\ &&\\ &=& \frac{a_0}{2L} \int_{-L}^L s_4 \, dx + \frac{1}{L} \sum_{n=1}^\infty a_n \int_{-L}^L c_n s_4 \, dx + \frac{1}{L} \sum_{n=1}^\infty b_n \int_{-L}^L s_n s_4 \, dx \\ &&\\ &=& \frac{a_0}{2L}\cdot 0 + \frac{1}{L} (0+0+0+\cdots ) + \frac{1}{L}(0+0+0+b_4L+0+0+0+\cdots ) \\ &&\\ &=& b_4. \end{eqnarray*} %\vfill %\pagebreak \subsection{Even and Odd Properties of Fourier Series} If $f(x)$ is an odd function then \[ a_n = \frac{1}{L} \int_{-L}^L \cos \left( \frac{n\pi x}{L} \right) f(x) \, dx = 0, \] since the product of an even function with an odd function is odd. \vspace{.2in} If $f(x)$ is an even function then \[ b_n = \frac{1}{L} \int_{-L}^L \sin \left( \frac{n\pi x}{L} \right) f(x) \, dx = 0, \] since the product of an odd function with an even function is odd. \section{An Example: The Square Wave} Let \[ f(x) = \left\{ \begin{array}{lll} -1 & \mbox{ for } & -\pi < x < 0 \\ 1 & \mbox{ for } & 0 < x < \pi, \end{array} \right. \] with $f(x+2\pi) = f(x)$ for all $x$. See the graph below. Note that the vertical lines are not part of the true graph, but are an artifact of the graphing program. I could have removed them, but they do show why this function is called a {\em Square Wave}. \vspace{.3in} \begin{center}\psfrag{Square}{Square}\psfrag{Wave}{\hspace*{0.3in}Wave} \includegraphics[width=6in]{Figs/sqplot1.eps} \end{center} {\bf Problem.} Find the Fourier series of $f(x)$. \vspace{.15in} {\bf Solution.} Since $f(x)$ is odd, $a_n=0$ for $n=0, 1, 2, \dots$. \vspace{.15in} Next, \[ b_n = \frac{1}{\pi} \int_{-\pi}^\pi \sin(nx) f(x) \\, dx = \frac{2}{\pi} \int_0^\pi \sin(nx)\cdot 1 \, dx = \] \[ \frac{-2}{n\pi} \left( \cos(n\pi) - \cos(0)\right) = \left\{\begin{array}{lll} 0 & \mbox{ for } & n \mbox{ even } \\ \frac{4}{n\pi} & \mbox{ for } & n \mbox{ odd. }\end{array} \right. \] \vspace{.25in} Thus, \[ f(x) = \sum_{k=1}^\infty \frac{4}{(2k-1)\pi} \sin \left( (2k-1)x \right) = \frac{4}{\pi} \sin x + \frac{4}{3\pi} \sin 3x + \frac{4}{5\pi} \sin 5x + \cdots. \] Recall that using $(2k-1)$ is a way to generate only odd numbers. \vspace{.3in} Next we plot a few partial sums. Let $f_N(x) = \sum_{k=1}^N \frac{4}{(2k-1)\pi} \sin \left( (2k-1)x \right)$. \begin{center} {\bf Plots of $f_N(x)$ for $N=1, 2, 3, 4, 10, 100$} \end{center} \begin{center}\psfrag{N}{$N=1$}\psfrag{=}{}\psfrag{1}{} \includegraphics[width=6in]{Figs/sqplot2.eps} \end{center} \begin{center}\psfrag{N}{$N=2$}\psfrag{=}{}\psfrag{2}{} \includegraphics[width=6in]{Figs/sqplot3.eps} \end{center} \begin{center}\psfrag{N}{$N=3$}\psfrag{=}{}\psfrag{3}{} \includegraphics[width=6in]{Figs/sqplot4.eps} \end{center} \begin{center}\psfrag{N}{$N=4$}\psfrag{=}{}\psfrag{4}{} \includegraphics[width=6in]{Figs/sqplot5.eps} \end{center}\psfrag{N}{$N=10$}\psfrag{=}{}\psfrag{10}{} \begin{center} \includegraphics[width=6in]{Figs/sqplot6.eps} \end{center} \begin{center}\psfrag{N}{$N=100$}\psfrag{=}{}\psfrag{100}{} \includegraphics[width=6in]{Figs/sqplot7.eps} \end{center} \section{Another Example! A Triangle Wave} {\bf Problem.} Find the Fourier series of the wave depicted below. \begin{center} \includegraphics[width=6in]{Figs/triplots1.eps} \end{center} {\bf Solution.} First note that since the period is 2, $L=1$. We can see that for $00$, in the table below. \begin{center}\renewcommand{\arraystretch}{2.5} \begin{tabular}{|c||c|c|c|c|c|c|c|c|} \hline $n$ & 1 & 2 & 3 & 4 & 5 & 6 & 7 & 8 \\ \hline $a_n$ & $\dis\frac{4}{\pi^2}$ & $\dis\frac{8}{4\pi^2}$ & $\dis\frac{4}{9\pi^2}$ & 0 & $\dis\frac{4}{25\pi^2}$ & $\dis\frac{8}{36\pi^2}$ & $\dis\frac{4}{49\pi^2}$ & 0 \\ \hline \end{tabular} \end{center} Therefore, \[ f(x) = \frac{1}{4} + \frac{4}{\pi^2}\cos \left(\frac{\pi x}{2}\right) + \frac{8}{4\pi^2}\cos \left(\frac{2\pi x}{2}\right) + \frac{4}{9\pi^2}\cos \left(\frac{3\pi x}{2}\right) + 0 + \frac{4}{25\pi^2}\cos \left(\frac{5\pi x}{2}\right) + \frac{8}{36\pi^2}\cos \left(\frac{6\pi x}{2}\right) + \cdots \] Plots of some partial sums are below. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \begin{center} {\bf Plots of some Partial Sums } \end{center} \begin{center} \psfrag{N=1}{$N=1$} \hspace*{-1in}\includegraphics[width=8in]{Figs/FSexampleplot1.eps} \vspace{.5in} \psfrag{N=2}{$N=2$} \hspace*{-1in}\includegraphics[width=8in]{Figs/FSexampleplot2.eps} \vspace{.5in} \psfrag{N=5}{$N=5$} \hspace*{-1in}\includegraphics[width=8in]{Figs/FSexampleplot3.eps} \vspace{.5in} \psfrag{N=20}{$N=20$} \hspace*{-1in}\includegraphics[width=8in]{Figs/FSexampleplot4.eps} \end{center} \section{An Example Done Entirely on Maple} %\vspace{.2in} {\bf Problem.} Find the Fourier series that will converge to $f(x) = 4 -x^2$ for $x \in [-2,2]$. \vspace{.1in} \begin{center} \includegraphics[height=1.5in]{Figs/ParaFplot1.eps} \end{center} \vspace{.2in} \noindent{\bf Solution.} Since $f(x)$ is even $b_n=0$. Next we compute $a_0$. \vspace{.2in} {\tt 1/2*int( 4-x\^{}2, x=-2..2);} %\vspace{.2in} \[ \frac{16}{3} \] Then we find $a_n$. \vspace{.2in} {\tt> 1/2*int( (4-x\^{}2)*cos(n*Pi*x/2), x=-2..2); } %\vspace{.2in} \[ - \frac{16(n \pi \cos(n\pi) - \sin(n\pi))}{n^3\pi^3} \] Then we plot a partial sum of the first 11 terms. Note that I simplified the expression for $a_n$. \vspace{.2in} {\tt> plot(8/3 - 16/Pi\^{}2*sum((-1)\^{}n/n\^{}2*cos(n*Pi*x/2),n=1..10),x=-10..10,thickness=2); } %\vspace{.2in} \begin{center} \includegraphics[height=1.5in]{Figs/ParaFplot2.eps} \end{center} \subsection{Extra Credit!} Looking at the graph from the last example, we notice it never quite gets to zero. We begin think, could there be a better way? The periodic extension of the parabola has very sharp corners. These make convergence ``difficult''. To get around this, use the periodic extension of a modified graph shown below. \begin{center} \includegraphics[height=2in]{Figs/paraflipplot1.eps} \end{center} The function is based on the original curve $y=4-x^2$ for $x\in[-2,2]$, but for $x\in[-4,2]$ and for $x\in[2,4]$ we have taken half of original curve flipped it twice and moved it into place. Find the equations for these pieces of the new function. Take its periodic extension, find its Fourier series and compare the convergence with the Fourier series used in the last example. \section{An Example that Requires no Computer Skills} {\bf Problem.} Find the Fourier series for $f(x) = \sin^3 x$. \begin{center} \includegraphics[height=2in]{Figs/sin3plot1.eps} \end{center} \noindent{\bf Solution.} You did this for homework already! Remember when you showed that \[ \sin^3 x = \frac{3}{4} \sin x - \frac{1}{4} \sin 3x. \] That is its Fourier series. It is odd, so all the $a_n$'s are zero. Then $b_1 = 3/4$, $b_3 = -1/4$ and all the other $b_n$'s are zero. \vspace{.2in} \noindent{\bf Conclusion.} Trigonometry is useful after all! Who would have guessed it! Call or write your trigonometry teacher and thank her or him! \section{An ODE with Square Wave Forcing Example} Let $f(x)$ be the square wave we studied earlier, shown again below. \begin{center}\psfrag{Square}{Square Wave}\psfrag{Wave}{} \includegraphics[width=6in]{Figs/sqplot1.eps} \end{center} \noindent{\bf Problem.} Solve the initial value problem $u'' + \frac{1}{4} u = f(x)$, $u(0)=u'(0)=0$. \vspace{.15in} \noindent{\bf Solution.} The solution to the corresponding homogeneous problem is $u_h = C_1 \sin (x/2) + C_2 \cos (x/2)$. To get a particular solution we recall the Fourier series we found for $f(x)$: \[ \sum_{n=1}^\infty b_n \sin (nx), \hspace{0.5in}\mbox{where }\, b_n = \left\{ \begin{array}{ll}0 & \mbox{for $n$ even,}\\ \frac{4}{n\pi} & \mbox{for $n$ odd.}\end{array}\right. \] This converges to $f(x)$ except at the jump discontinuities. The idea is, we can consider each term separately. That is for $n=1,2,3,\dots$ we have \[ u'' + \frac{u}{4} = \sin nx. \] For each $n$ we shall find a particular solution $u_n$. Let \[ u_n = A \sin nx + B \cos nx. \] Then \[ u''_n = -An^2 \sin nx - B n^2 \cos nx. \] Thus, $-An^2 + A/4 =1$ and $-Bn^2 + B/4 = 0$. Thus, $B=0$ and $A = \dis\frac{1}{\frac{1}{4} - n^2}$. Hence, \[ u_n = \frac{1}{\frac{1}{4} - n^2} \sin nx, \hspace{0.5in} \mbox{ for } n = 1, 2, 3, \dots . \] Now we let \[ u_p = \sum_{n=1}^\infty \frac{b_n}{\frac{1}{4} - n^2} \sin nx. \] It follows by linearity that $u''_p +\frac{1}{4} u_p = $ the Fourier series of $f(x)$ which equals $f(x)$, almost everywhere. We can rewrite $u_p$ as \[ u_p = \sum_{k=1}^\infty \frac{\frac{4}{(2k-1)\pi} \sin((2k-1)x)}{\frac{1}{4} - (2k-1)^2} = \frac{16}{\pi} \sum_{k=1}^\infty \frac{\sin((2k-1)x)}{(2k-1)(1-4(2k-1)^2)}. \] Now, we have that the general solution is \[ u = u_h + u_p = C_1 \sin (x/2) + C_2 \cos (x/2) + u_p. \] We still need to find $C_1$ and $C_2$ such that $u(0) = u'(0) = 0$. \[ u(0) = 0 + C_2 + \frac{16}{\pi}\sum_{k=1}^\infty 0 = C_2. \] Thus $C_2 = 0$. Now, \[ u'(x) = \frac{C_1}{2} \cos (x/2) + \frac{16}{\pi}\sum_{k=1}^\infty \frac{(2k-1)\cos((2k-1)x)}{(2k-1)(1-4(2k-1)^2)}. \] Thus, \[ u'(0) = \frac{C_1}{2} + \frac{16}{\pi}\sum_{k=1}^\infty \frac{1}{1-4(2k-1)^2}. \] It turns out $\dis \sum_{k=1}^\infty \frac{1}{1-4(2k-1)^2} = -\pi/8$. (I'll show how to do this in a bit.) Then, since $u'(0) = 0$ we have \[ C_1 = \frac{32}{\pi} \cdot \frac{\pi}{8} = 4. \] Our final solution is \[ u(x) = 4 \sin(x/2) + \frac{16}{\pi} \sum_{k=1}^\infty \frac{\sin((2k-1)x)}{(2k-1)(1-4(2k-1)^2)}. \] Next we plot some partial sums. \subsection{Graphs for ODE Example } \begin{center}\psfrag{sin(x/2)}{$4\sin x/2$} \includegraphics[height=1.5in]{Figs/ODEexampleplot1.eps} \end{center} \begin{center}\psfrag{N=1}{$N=1$} \includegraphics[height=1.5in]{Figs/ODEexampleplot2.eps} \end{center} \begin{center}\psfrag{N=2}{$N=2$} \includegraphics[height=1.5in]{Figs/ODEexampleplot3.eps} \end{center} \begin{center}\psfrag{N=3}{$N=3$} \includegraphics[height=1.5in]{Figs/ODEexampleplot4.eps} \end{center} \begin{center}\psfrag{N=10}{$N=10$} \includegraphics[height=1.5in]{Figs/ODEexampleplot5.eps} \end{center} \subsection{Back to computing that series...} Back to why $\dis \sum_{k=1}^\infty \frac{1}{1-4(2k-1)^2} = -\pi/8$. \vspace{.15in} {\bf Method one.} Use a computer. Here is the command in Maple: \[ \tt sum(1/(4*(2*k-1)^2-1), k=1..infinity); \] It returns $\pi/8$. (I switched the terms in the denominator to make it positive.) I could not find a way to do this is Maxima; it can only do numerical approximations. \vspace{.15in} {\bf Method two.} Recall that you took a course called {\em Calculus II}. Look up the Taylor series for $\arctan x$, then plug in $x=1$. This shows that \[ 1 -\frac{1}{3} + \frac{1}{5} - \frac{1}{7} + \frac{1}{9} - \cdots = \arctan (1) = \frac{\pi}{4}. \] We will apply this fact to our sum. But, first we notice that \begin{eqnarray*} \frac{1}{4(2k-1)^2 - 1} &=& \frac{1}{16k^2 -16k + 3} \\ &&\\ &=& \frac{1}{(4k-3)(4k-1)} \\ &&\\ &=& \frac{\frac{1}{2}}{4k-3} - \frac{\frac{1}{2}}{4k-1}. \end{eqnarray*} Therefore, \begin{eqnarray*} \sum_{k=1}^\infty \frac{1}{4(2k-1)^2 - 1} &=& \frac{1}{2}\sum_{k=1}^\infty \left(\frac{1}{4k-3} - \frac{1}{4k-1}\right) \\ &&\\ &=& \frac{1}{2}\left(1 -\frac{1}{3} + \frac{1}{5} - \frac{1}{7} + \frac{1}{9} - \cdots \right) \\ &&\\ &=& \frac{1}{2}\frac{\pi}{4} = \frac{\pi}{8}. \end{eqnarray*} Note: There is a subtle issue with convergence here. The radius of convergence of the Taylor series of $\arctan x$ is 1. At $x=1$ the series does converge by the alternating series test. But the proof that it converges to $\arctan 1$ requires the Weierstrass M-test which is covered in Math 352. \section{The Wave Equation for a Vibrating String} We are going to use Fourier series to solve some partial differential equations. The methods used will seem quite strange and ad hoc at first. They are in fact standard and are used to solve a number of important PDE's used in physics, chemistry and engineering. We shall start will a classic problem of modeling a vibrating string. When a string a held straight and not moving it is in equilibrium. We will start the motion but plucking it, that is we pinch a straight string in the middle and pull in up a bit a let it go. Let the string have length $L>0$ and use $x\in [0,L]$ as a as the distance to the left end point. Let $t$ be time. Let $u(x,t)$ be the height of the string away from equilibrium. If $u(x,t)$ is negative then that point of the string is below equilibrium. The PDE that is used to model this is called {\bf the wave equation}. It is \[ \frac{\bd^2 u}{\bd x^2} = k\frac{\bd^2 u}{\bd t^2}, \] where $k>0$ is a constant. This rationale is that the force at any point should be proportional to the concavity at that point. It is a simplified model; for example it does not include a term for air resistance or internal heat loss. Real strings can wear out and break, but not our ideal string. It is customary to rewrite the wave equation as \[ a^2 u_{xx} = u_{tt}, \] where $a>0$ is a constant. Next we add {\bf boundary conditions}. These simply model the fact that we shall be holding the end points fixed. Thus, for all $t$, \[ u(0,t) = u(L,t) = 0. \] Finally we add the {\bf initial conditions} or {\bf configuration}. For simplicity we shall give use $u(L/2,0) =L/2$ and assume the string is piecewise linear. Then \[ u(x,0) = f(x) = \left\{ \begin{array}{ll} x & \mbox{for } x\in [0,L/2],\\ L-x & \mbox{for } x \in (L/2,L]. \end{array}\right. \] \begin{center}\psfrag{0}{0}\psfrag{L}{$L$}\psfrag{L2}{$L/2$} \includegraphics[width=3in]{Figs/string.eps} \end{center} We encode the fact that initially the string is not moving by \[ u_t(x,0) = 0 \] for all $x \in [0,L]$. In general we could specify $u_t(x,0) = g(x)$ for some function $g(x)$. Thus, our model is as shown in the box below. \vspace{0.15in} \begin{center} \fbox{\parbox{3in}{ \[ a^2 u_{xx} = u_{tt} \] \[ u(0,t) = u(L,t) = 0 \] \[ u(x,0) = f(x) \hspace{0.4in} u_t(x,0) = 0 \] where $f(x)$ was given above.}} \end{center} \subsection{Solving the Model} Our strategy is as follows. Suppose $u(x,t)$ can be written in the form \[ u(x,t) = X(x)T(t). \] Plugging into the wave equation we get \[ a^2 X''(x) T(t) = X(x)T''(t). \] Now we separate the terms involving $x$ from the terms involving $t$ to get \[ \frac{X''(x)}{X(x)} =\frac{1}{a^2} \frac{T''(t)}{T(t)}. \] Both $x$ and $t$ are free, independent variables. But, if we change $x$ notice the right hand side will not change because it depends only on $t$. Therefore $X''/X = $ a constant. Likewise $T''/a^2T$ = the same constant. We will call this constant $-\sigma$ because we can. Thus, we now have two second order ODEs: \[ X'' + \sigma X = 0 \hspace{.25in}\&\hspace{.25in} T'' + \sigma a^2 T = 0. \] For the first we can deduce boundary conditions. \[ u(0,t) = 0 \implies X(0)T(t) = 0 \implies X(0) = 0. \] Likewise $X(L)=0$. So, let's just focus on the problem \[ X''(x) +\sigma X(x) = 0 \hspace{.25in} X(0) = X(L)=0. \] We know that the form of the solution will depend on the sign of $\sigma$. We will consider three cases, $\sigma <0$, $\sigma = 0$ and $\sigma >0$, and see what happens. It will turn out that nontrivial solutions will exist only for certain values of $\sigma$. \vspace{.2in} {\bf Case I.} Suppose $\sigma <0$. Then the general solution is \[ X(x) = C_1 e^{\sqrt{-\sigma}x} + C_2 e^{-\sqrt{-\sigma}x}. \] We impose the boundary conditions $X(0) = X(L) = 0$ to get $C_1$ and $C_2$. We get that $X(0)=0$ implies $C_1 + C_2 = 0$ and $C_1 e^{\sqrt{-\sigma}L} + C_2 e^{-\sqrt{-\sigma}L} = 0$. Thus, \[ C_1 = - C_2 \] and \[ C_1 = -e^{-2\sqrt{-\sigma}L} C_2. \] Thus we can only get a solution besides $C_1 = C_2 = 0$ if $e^{-2\sqrt{-\sigma}L} = 1$. But, this implies $\sigma = 0$. We conclude that there are no nontrivial solutions if $\sigma <0$. \vspace{.2in} {\bf Case II.} Suppose $\sigma = 0$. Now our equation is just $X''(x) = 0$. The general solution is \[ X(x) = C_1 x + C_2. \] Notice, $X(0)=0$ implies $C_2 = 0$, but then $X(L)=0$ implies $C_1L=0$ giving us that $C_1=0$. Again, there are no nontrivial solutions when $\sigma = 0$. It seems we are wasting our time! \vspace{.2in} {\bf Case III.} As a final act of desperation suppose $\sigma >0$. Now the general solution is \[ X(x) = C_1 \sin \sqrt{\sigma}x + C_2 \cos \sqrt{\sigma}x. \] Now, $X(0) = C_2$, which implies $C_2=0$. But $X(L) = C_1 \sin (\sqrt{\sigma}L)$. This implies $C_1 = 0$ unless the sine of $\sqrt{\sigma}L$ is zero. So, our only hope of finding nontrivial solutions is to suppose $\sqrt{\sigma}L$ is an integer multiple of $\pi$. That is \[ \sigma = \frac{n^2 \pi^2}{L^2} \] for some positive integer $n$. For these values \[ X(L) = C_1 \sin \left(\sqrt{\sigma}L\right) = C_1 \sin (n\pi) = 0, \] for any value of $C_1$! Thus, for any nonzero integer $n$ we have infinitely many nontrivial solutions; we will take $n$ to be positive. This is not the behavior you are used from working with initial value problems for ODEs. There we got unique solutions. But, so far we have only done the {\em boundary values} and we have not yet considered the initial shape and velocity on our string. We will record our observation of nontrivial solutions by defining \[ X_n(x) = \sin \left(\frac{n\pi x}{L}\right) \mbox{ for } n\geq 1. \] {\bf Next we study the $T(t)$ equation.} \vspace{.2in} Recall \[ T''(t) + a^2 \sigma T(t) = 0. \] Since we can only get nontrivial solutions for $X(x)$ if $ \sigma= \dfrac{n^2 \pi^2}{L^2}$ we shall assign this value to $\sigma$. Thus, we have \[ T''(t) + \frac{a^2n^2 \pi^2}{L^2}T(t) = 0. \] The general solution is \[ T(t) = C_1 \sin \left( \frac{a n \pi t}{L} \right) + C_2 \cos \left( \frac{a n \pi t}{L} \right). \] Now we shall begin looking at the initial conditions. It turns out it is easier if we study the initial velocity first. We were given \[ u_t(x,0) = 0, \] meaning that the string is not moving at the moment we let go of it. Now $u_t = \bd_t (X(x)T(t)) = X(x)T'(t)$. Thus we have \[ X(x)T'(0) = 0. \] Assuming we have a nontrivial solution for $X(x)$ this forces $T'(0) = 0$. We compute \[ T'(t) = C_1 \frac{a n \pi}{L} \cos \left( \frac{a n \pi t}{L} \right) - C_2 \frac{a n \pi}{L} \sin \left( \frac{a n \pi t}{L} \right). \] Thus, \[ T'(0) = C_1 \frac{a n \pi}{L} = 0 \implies C_1 = 0. \] There are no restrictions on $C_2$. For any value of $C_2$ \[ T(t) = C_2 \cos \left( \frac{a n \pi t}{L} \right) \] solves the differential equation in $T$ and the condition $T'(0) = 0$. We let $\dis T_n(t) = \cos \left( \frac{a n \pi t}{L} \right)$ and define $u_n(x,t) = X_n(x)T_n(t)$. \vspace{0.3in} Each $u_n$ satisfies the wave equation, the two boundary conditions, and the initial velocity condition. Furthermore, by linearity, any linear combination of the $u_n$'s will also satisfy these conditions. {\bf Check this.} \vspace{0.3in} It can be shown that if we let \[ u(x,t) = \sum_{n=1}^\infty c_n u_n(x,t) \] and the $c_n$'s are such that it converges everywhere, then $u$ satisfies the wave equation, the two boundary conditions, and the initial velocity condition. This in proven in Math 407. The final step, is that we need to choose the $c_n$'s so that $$\dis \sum_{n=1}^\infty c_n u_n(x,0) = f(x).$$ We do this next. We need for \[ f(x) = \dis \sum_{n=1}^\infty c_n u_n(x,0) = \sum_{n=1}^\infty c_n \sin \left(\frac{n\pi x}{L}\right) \cos \left( \frac{a n \pi 0}{L} \right) = \sum_{n=1}^\infty c_n \sin \left(\frac{n\pi x}{L}\right).\] But, this is just the Fourier series for the odd periodic extension of $f(x)$. \begin{center}\psfrag{L}{$L$}\psfrag{L2}{$L/2$}\psfrag{-L}{$-L$}\psfrag{-L2}{$-L/2$} \includegraphics[width=3in]{Figs/stringoddext.eps} \end{center} Thus, \begin{eqnarray*} c_n &=& \frac{1}{L} \int_{-L}^L \tilde{f_o}(x) \sin \left(\frac{n\pi x}{L}\right) \, dx \\ &=& \frac{2}{L} \int_0^L f(x) \sin \left(\frac{n\pi x}{L}\right) \, dx \\ &=& \frac{2}{L}\int_0^{L/2} x \sin \left(\frac{n\pi x}{L}\right) \, dx + \frac{2}{L}\int_{L/2}^L (-x+L) \sin\left(\frac{n\pi x}{L}\right)\, dx \\ &=& ..... \mbox{ busy work } .... \\ &=& \frac{4L}{n^2\pi^2} \sin \left( \frac{n\pi}{2} \right) \end{eqnarray*} Therefore, \[ u(x,t) = \sum_{n=1}^\infty \frac{4L}{n^2\pi^2} \sin \left( \frac{n\pi}{2} \right) \sin \left(\frac{n\pi x}{L}\right) \cos \left( \frac{a n \pi t}{L} \right). \] We can rewrite this as \[ u(x,t) = \frac{4L}{\pi^2}\sum_{k=0}^\infty \frac{(-1)^k}{(2k+1)^2} \sin \left(\frac{(2k+1)\pi x}{L}\right) \cos \left( \frac{a (2k+1) \pi t}{L} \right). \] \subsection{Graphs and an Animation} Now for the fun part. We will make an animation of $u(x,t)$. We shall take $a=1$ and $L=2$. \vspace{0.15in} {\tt \begin{tabular}{ll} > with(plots); & \# Load package of plotting commands. \\ > c:= n -> 8*sin(n*Pi/2)/(n\^{}2*Pi\^{}2): & \# Define the coefficient function.\\ > N:=100: &\# Number of terms to use.\\ > &\# Finally we execute the animation command. \end{tabular} \begin{tabular}{ll} > animate(\{sum((c(n)*sin(n*Pi*x/2)*cos(n*Pi*t/2)), n=1..N)\}, &\\ \hspace*{0.2in} x=0..2, t=0..20, frames=150, thickness=2); & \end{tabular}} \begin{center} %\includegraphics[width=5in]{../S17/a5.gif} \end{center} See the animation link on the course web site. \subsection{Summary for Wave Equation} Our model for a vibrating string is repeated below. \vspace{0.3in} \begin{center} \fbox{\parbox{3in}{ \[ a^2 u_{xx} = u_{tt} \] \[ u(0,t) = u(L,t) = 0 \] \[ u(x,0) = f(x) \hspace{0.4in} u_t(x,0) = 0 \] where $f(x)$ is given.}} \end{center} The solution is \[ u(x,t) = \sum_{n=1}^\infty c_n \sin \left(\frac{n \pi x}{L} \right) \cos \left( \frac{an\pi t}{L} \right) \] where \[ c_n = \frac{2}{L}\int_0^L \sin \left(\frac{n \pi x}{L} \right) f(x) \, dx. \] \subsection{Extra Credit!} Modify the wave equation by adding a damping term, $\dis \gamma \frac{\bd u}{\bd t}$. Use the same boundary and initial conditions as in the example we just did. Solve this model. Write this up neatly and turn it in. Make an animation, put it on web site, and send me the link. When you are at an interview for a job or internship, bring your tablet or phone, and show your animation to the interviewer. \section{The Heat Equation} We consider a metal rod of length $L$. Suppose that initially the temperature at each point $x \in [0,L]$ along the rod is given by $f(x)$. We wish to know how the temperature will evolve over time. Let \[ u(x,t) \] be the temperature at location $x$ and time $t$. Thus, $u(x,0) = f(x)$. The evolution of the temperature along to rod is modeled by the {\bf heat equation}, which is \[ \alpha^2 \frac{\bd^2 u}{\bd x^2} = \frac{\bd u}{\bd t}, \] or $\alpha^2 u_{xx} = u_t$. See Appendix in textbook. We will make an additional assumption that temperature of the ends of the metal rod are held at fixed temperatures. This gives the {\bf boundary conditions} \[ u(0,t) = T_1 \hspace{.5in} \& \hspace{.5in} u(L,t) = T_2. \] Later we will consider other types of boundary conditions. We will solve the present problem in two steps. First, we assume $T_1=T_2=0$. Then we will show how to jazz this up to the more general case. \subsection{Solving Homogeneous Case} Step 1 is to solve the following model. \begin{center} \fbox{\parbox{3in}{ \[ \alpha^2 u_{xx} = u_{t} \] \[ u(0,t) = u(L,t) = 0 \] \[ u(x,0) = f(x) \] where $f(x)$ will be given.}} \end{center} \noindent{\bf Solution.} Suppose $u(x,t) = X(x) T(t)$. Then $u_{xx}(x,t) = X''(x)T(t)$ and $u_t(x,t) = X(x)T'(t)$. Therefore, \[ \alpha^2 X''T = XT' \hspace{0.5in} \implies \hspace{0.5in} \frac{X''}{X} = \frac{1}{\alpha^2}\frac{T'}{T}. \] This implies that each side of the second equation above must be constant. We call the constant $-\sigma$. This yields two ordinary differential equations. \[ X'' + \sigma X = 0 \hspace{0.5in} \& \hspace{0.5in} T' + \sigma \alpha^2 T = 0. \] The boundary conditions $u(0,t) = u(L,t) = 0$ implies $X(0) = X(L) = 0$. Then for the equation in $X$ the analysis we did for the wave equation goes through in exactly the same way. We conclude that non trivial solution exist only when $\sqrt{\sigma} = \dfrac{n\pi}{L}$ where $n$ is an integer, and that \[ X(x) = C \sin \left( \frac{n\pi x}{L} \right) \] gives a solution for any real number $C$ and integer $n$. Let \[ X_n(x) = \sin \left( \frac{n\pi x}{L} \right), \:\:\mbox{ for }\:\: n\geq 1. \] Note: $n=0$ gives a trivial solution and negative values for $n$ are redundant. Now, consider the ordinary differential equation in $T$. Its general solution is just \[ T(t) = C e^{-\sigma\alpha^2t} = C e^{-\frac{n^2\pi^2\alpha^2}{L^2}t}. \] We let \[ T_n(t) = e^{-\frac{n^2\pi^2\alpha^2}{L^2}t} \hspace{0.5in} \& \hspace{0.5in} u_n(x,t) = X_n(x)T_n(t) \:\:\: \mbox{ for }\:\: n\geq 1. \] We have \[ u_n(x,t) = X_n(x)T_n(t) = \sin\left(\frac{n\pi x}{L}\right) e^{-\frac{n^2\pi^2\alpha^2}{L^2}t} \:\:\: \mbox{ for } n\geq 1 . \] Each $u_n$ satisfies the heat equation and the boundary conditions. So would any linear combination of $u_n$. It is shown in Math 407 that if $\{ c_n \}_{n=1}^\infty$ are such that \[ u(x,t) = \sum_{n=1}^\infty c_n u_n(x,t) \] converges, then $u(x,t)$ satisfies the heat equation (almost everywhere) and the boundary conditions. \vspace{.3in} Therefore, if we can find values for the $c_n$'s such that \[ u(x,0) = \sum_{n=1}^\infty\sin\left(\frac{n\pi x}{L}\right) = f(x) \] where $f(x)$ is the initial temperature distribution, then $u(x,t)$ will satisfy all our requirements. \subsection{An Example: Metal Rod with Ends held at Zero Degrees} Suppose $f(x) = 10$, $L=1$ and $\alpha=1$. At the instant $t=0$ the end points of the rod put into contact with degree zero heat sinks. We will model what happens next. \begin{center} \psfrag{a}{$0^o$} \psfrag{b}{$0^o$} \psfrag{10}{$10^o$} \includegraphics[height=1.0in]{Figs/hotrod.eps} \end{center} The odd periodic extension of $f(x)$ is a square wave. Let $L=1$. We compute the $c_n$'s. \[ c_n = \frac{1}{1} \int_{-1}^1 \sin(n\pi x) f(x) \, dx = 20 \int_0^1 \sin(n\pi x) \, dx = -\frac{20}{n\pi} \big[ \cos(n\pi)- \cos(0)\big] = \left\{\begin{array}{ccc}\dfrac{40}{n\pi} & \mbox{ for } & n \mbox{ odd},\\ &&\\ 0 & \mbox{ for } & n \mbox{ even}. \end{array}\right. \] Thus, \[ u(x,t) = \sum_{k=1}^\infty \frac{40}{(2k-1)\pi} \sin( (2k-1)\pi x) e^{-(2k-1)^2\pi^2 t}. \] An animation is on the web page. Below are some snapshots. \begin{center}{\bf Some Plots for the Last Example} \end{center} \begin{center} \psfrag{t=0.0}{$t=0.0$} \includegraphics[height=1.8in]{Heat2plot2.eps} \psfrag{t=0.00001}{$t=0.00001$} \includegraphics[height=1.8in]{Heat2plot3.eps} \psfrag{t=0.0001}{$t=0.0001$} \includegraphics[height=1.8in]{Heat2plot4.eps} \vspace{0.25in} \psfrag{t=0.001}{$t=0.001$} \includegraphics[height=1.8in]{Heat2plot5.eps} \psfrag{t=0.01}{$t=0.01$} \includegraphics[height=1.8in]{Heat2plot6.eps} \psfrag{t=0.1}{$t=0.1$} \includegraphics[height=1.8in]{Heat2plot7.eps} \vspace{0.25in} \psfrag{t=0.2}{$t=0.2$} \includegraphics[height=1.8in]{Heat2plot8.eps} \psfrag{t=0.3}{$t=0.3$} \includegraphics[height=1.8in]{Heat2plot9.eps} \end{center} \subsection{A Crazy Example } Suppose $L=1$, $\alpha=1$ and the initial temperature distribution is given by \[ f(x) = \left\{\begin{array}{lll} 0 & \mbox{ for } & 0\leq x < 0.2 \\ -500(x-0.2)(x-0.4) & \mbox{ for } & 0.2\leq x < 0.4 \\ 0 & \mbox{ for } & 0.4\leq x < 0.6 \\ 4 & \mbox{ for } & 0.6\leq x < 0.8 \\ 0 & \mbox{ for } & 0.8\leq x < 1.0 \end{array} \right. \] Below is a graph. \begin{center} \includegraphics[height=1.5in]{Heat2plot10.eps} \end{center} Here are the Maple commands that created it. \vspace{.15in} {\tt \noindent \red{> f := x -> piecewise(x<.2, 0, x>=.2 and x<=.4 , -500*(x-0.2)*(x-0.4), \\x>0.4 and x<0.6, 0, x>=0.6 and x<=0.8, 4, x>0.8, 0); } \noindent\red{> plot(f(x),x=0..1);}} \vspace{.15in} The coefficients are given by $c_n = 2 \dis\int_0^1 \sin(n\pi x) f(x) \, dx$. I did not bother trying to find a formula for the $c_n$, I just used the integral. Thus the solution is \[ u(x,t) = \sum_{n=1}^\infty c_n \sin(n\pi x) e^{-n^2\pi^2t}. \] I created an animation using the commands below. I took the sum out to $n=50$, but this is likely overkill. It is on the web-page. \vspace{.2in} {\tt \noindent\red{> c := n -> 2*int(sin(n*Pi*x)*f(x) ,x=0..1);} \noindent\red{> animate(sum( c(n)*sin((n)*Pi*x)*exp(-(n\^{}2*Pi\^{}2*t)) ,n=1..50), \\x=0..1,t=0..0.2,numpoints=1000,frames=100);} } \vspace{0.15in} Snapshots of the animation are below. \begin{center}{\bf Some Plots for the Last Example} \end{center} \begin{center} \psfrag{t=0.0000}{$t=0.0$} \includegraphics[height=1.8in]{Heat2plot13.eps} \psfrag{t=0.00001}{$t=0.00001$} \includegraphics[height=1.8in]{Heat2plot14.eps} \psfrag{t=0.0001}{$t=0.0001$} \includegraphics[height=1.8in]{Heat2plot15.eps} \vspace{.2in} \psfrag{t=0.001}{$t=0.001$} \includegraphics[height=1.8in]{Heat2plot16.eps} \psfrag{t=0.002}{$t=0.0002$} \includegraphics[height=1.8in]{Heat2plot17.eps} \psfrag{t=0.003}{$t=0.0003$} \includegraphics[height=1.8in]{Heat2plot18.eps} \vspace{.2in} \psfrag{t=0.005}{$t=0.005$} \includegraphics[height=1.8in]{Heat2plot19.eps} \psfrag{t=0.01}{$t=0.01$} \includegraphics[height=1.8in]{Heat2plot20.eps} \psfrag{t=0.05}{$t=0.05$} \includegraphics[height=1.8in]{Heat2plot21.eps} \end{center} \section{Heat Equation with Nonhomogeneous Boundary Conditions} Now we are working with the model \begin{center} \fbox{\parbox{3in}{ \[ \alpha^2 u_{xx} = u_{t} \] \[ u(0,t) = T_1 \hspace{0.5in} \& \hspace{0.5in} u(L,t) = T_2 \] \[ u(x,0) = f(x) \] where $f(x)$ will be given.}} \end{center} {\bf Key Idea.} As $t\to\infty$ we think that $u(x,t) \to \dfrac{T_2-T_1}{L} x + T_1$. So, we let \[ v(x) = \dfrac{T_2-T_1}{L} x + T_1 \hspace{0.5in} \mbox{and} \hspace{0.5in} w(x,t) = u(x,t) - v(x). \] Notice that $w_{xx} = u_{xx}$ and $w_t = u_t$. Therefore, $w$ satisfies the heat equation \[ \alpha^2 w_{xx} = w_t. \] Further, $w(0,t) = u(0,t) - v(0) = T_1 - T_1 = 0$ and likewise $w(L,t) = 0$. However, \[ w(x,0) = u(x,0) - v(x) = f(x) - v(x). \] Let $g(x) = f(x)- v(x)$. Then $w(x,t)$ satisfies the model below. \begin{center} \fbox{\parbox{3in}{ \[ \alpha^2 w_{xx} = w_{t} \] \[ w(0,t) = w(L,t) = 0 \] \[ w(x,0) = g(x). \] }} \end{center} Our plan is to solve this model for $w(x,t)$ and then add $v(x)$ to get $u(x,t)$. We do this in the next example. \subsection{Example: Metal Rod Ends Held at Different Temperatures} Let $L=1$, $T_1 = 10^o$ and $T_2=20^o$. We suppose the bar starts with uniform temperature $10^o$ and at time $t=0$ the ends are set to a heat bath of $T_1$ or the right and $T_2$ on the left. \begin{center} \psfrag{a}{\hspace{-0.1in}$10^o$} \psfrag{b}{$20^o$} \psfrag{10}{$10^o$} \includegraphics[height=1.0in]{Figs/hotrod.eps} \end{center} Thus, $v(x) = 10x + 10$ and $g(x)=-10$. The graph of the odd periodic extension of $g(x)$ is below along with the command that created it. Note: {\tt frac} is a Maple command that gives the fractional part of a number. \vspace{.2in} {\tt > plot(piecewise(x>=-1,-20*frac((x+1)/2)+10,x<-1,-20*frac((x+1)/2)-10), \\x=-5..5,discont=true, color=plum, thickness=2);} \begin{center} \includegraphics[height=1.8in]{Heat3plot1.eps} \end{center} Now we compute the Fourier series coefficients for the odd periodic extension of $g(x)$. \[ c_n = \frac{1}{1} \int_{-1}^1 \sin(n\pi x) g(x) \, dx = -20 \int_{0}^1 x \sin(n\pi x) \, dx = \frac{20}{n\pi}\cos(n\pi) = \frac{20}{n\pi}(-1)^n. \] Thus, \[ w(x,t) =\frac{20}{n\pi} \sum_{n=1}^\infty (-1)^n \sin(n\pi x) e^{-n^2\pi^2t}. \] Hence, \[ u(x,t) = g(x) + w(x,t) = 10 + 10x + \frac{20}{n\pi} \sum_{n=1}^\infty (-1)^n \sin(n\pi x) e^{-n^2\pi^2t}. \] Some plots are below. Animation coming soon to a website near you! \begin{center} {\bf Plots for Last Example} \end{center} \begin{center} \psfrag{t=0.0}{$t=0.0$} \includegraphics[height=1.9in]{Heat3plot3.eps} \psfrag{t=0.001}{$t=0.001$} \includegraphics[height=1.9in]{Heat3plot4.eps} \psfrag{t=0.01}{$t=0.01$} \includegraphics[height=1.9in]{Heat3plot5.eps} \vspace{0.25in} \psfrag{t=0.1}{$t=0.1$} \includegraphics[height=1.9in]{Heat3plot6.eps} \psfrag{t=0.5}{$t=0.5$} \includegraphics[height=1.9in]{Heat3plot7.eps} \end{center} \section{Heat Equation with Insulated Ends} Now we consider a metal rod where heat cannot flow out through the ends. The end are insulated. The model now becomes the following. \begin{center} \fbox{\parbox{2.5in}{ \[ \alpha^2 u_{xx} = u_{t} \] \[ u_x(0,t) = u_x(L,t) = 0 \] \[ u(x,0) = f(x), \] where $f(x)$ is given.}} \end{center} We go through the same steps are before. Suppose $u(x,t)=X(x)T(t)$. Then the heat equation becomes \[ \alpha^2 X''(x)T(t) = X(x)T'(t), \] and just as before we have \[ \frac{X''}{X} = \frac{1}{\alpha^2} \frac{T'}{T} = \mbox{ a constant } = -\sigma. \] Thus, we get two ODE's \[ X'' + \sigma X = 0 \hspace{.3in} \& \hspace{.3in} T' + \sigma\alpha^2 T = 0. \] We work with the $X$ equation first. We consider the same three cases $\sigma <0$, $\sigma = 0$ and $\sigma >0$. The boundary condition translates to \[ X'(0) = X'(L) = 0. \] {\bf Suppose $\mathbf\sigma \mathbf < \mathbf 0$.} It will turn out there are no nontrivial solutions. The general solution is \[ X(x) = C_1 e^{\sqrt{-\sigma}x} + C_2 e^{-\sqrt{-\sigma}x} \] Then \[ X'(x) = C_1 \sqrt{-\sigma}e^{\sqrt{-\sigma}x} - C_2 \sqrt{-\sigma}e^{-\sqrt{-\sigma}x} \] At the boundary points $x=0$ and $x=L$ we have \[ X'(0) = C_1\sqrt{-\sigma} - C_2\sqrt{-\sigma} = 0 \:\:\: \implies \:\:\: C_1 = C_2 \] and \[ X'(L) = C_1\sqrt{-\sigma}e^{\sqrt{-\sigma}L} - C_2\sqrt{-\sigma}e^{-\sqrt{-\sigma}L} \:\:\: \implies \:\:\: C_1 = C_2 e^{-2\sqrt{-\sigma}L} \] These only have a solution besides $C_1=C_2=0$ if $L=0$, which we do not allow. \vspace{0.25in} {\bf Suppose $\mathbf\sigma \mathbf = \mathbf 0$.} The general solution is $X(x) = C_1x + C_2$. Since $X'(x) = C_1$ it must be that $C_1=0$. Any $X(x) = C_2$ will be a solution. \vspace{.2in} {\bf Finally, suppose $\mathbf\sigma \mathbf > \mathbf 0$.} The general solution is $X(x) = C_1 \sin \sqrt{\sigma}x + C_2 \cos \sqrt{\sigma}x$. Thus, \[ X'(x) = C_1 \sqrt{\sigma}\cos \sqrt{\sigma}x - C_2 \sqrt{\sigma}\sin \sqrt{\sigma}x. \] Thus, \[ X'(0) = C_1 \sqrt{\sigma} = 0 \:\:\: \implies C_1=0. \] Then \[ X'(L) = - C_2 \sqrt{\sigma}\sin \sqrt{\sigma}L. = 0 \:\:\: \implies \:\:\: \sigma = \frac{n^2\pi^2}{L^2} \mbox{ or } C_2=0. \] Thus, if we want nontrivial solutions we require $\sigma = \frac{n^2\pi^2}{L^2}$. Then there are no restrictions on $C_2$. We let \[ X_n(x) = \cos \left(\frac{n\pi x}{L}\right) \:\:\: \mbox{ for } n = 0, 1, 2, 3, ... \] Notice $n=0$ gives us a constant. Recall the equation for $T(t)$ was $T' + \sigma\alpha^2 T = 0$. It is general solution is \[ T(t) = C e^{-\sigma\alpha^2t}. \] We let \[ T_n(t) = e^{-\frac{n^2 \pi^2 \alpha^2}{L^2}t} \:\:\: \mbox{ for } n = 0, 1, 2, 3, ... \] Let \[ u_n(x,t) = X_n(x)T_n(t)\ :\:\: \mbox{ for } n = 0, 1, 2, 3, ... \] Each $u_n(x,t)$ satisfies the heat equation and the boundary conditions on $u_x$. So too does any linear combination of them. In is shown in Math 407 that if $c_0, c_1, c_2,c_3, \dots$ are such that \[ u(x,t) = c_0 + \sum_{n=1}^\infty c_n X_n T_n \] converges, then $u(x,t)$ satisfies the heat equation and the boundary conditions on $u_x$. All that is left to do is find $c_n$'s such that $u(x,0) = f(x)$, the initial condition. But this is just the Fourier series of the {\bf even} periodic extension of $f(x)$. Thus, $c_0$ is the average value of $f(x)$ and for $n=1, 2, 3,...$ \[ c_n = \frac{1}{L}\int_{-L}^L \cos \left( \frac{n \pi x}{L}\right) f_e(x) \, dx = \frac{2}{L} \int_0^L \cos \left( \frac{n \pi x}{L}\right) f(x) \, dx. \] \subsection{Example for Heat Equation with Insulated Ends} Suppose $L=\pi$, $\alpha=1$ and $u(x,0) = f(x) = \sin(x)$. The even periodic extension of $f(x)$ is $|\sin (x)|$. \begin{center} \psfrag{Even Extension}{Even Extension} \includegraphics[height=1.5in]{Heat4plot1.eps} \end{center} \noindent But, using symmetry our integrals are only over $[0,\pi]$. Then \[ c_0 = \frac{1}{\pi} \int_0^\pi \sin x \, dx = \frac{2}{\pi} \] and for $n=1,2,3, ...$ we have \[ c_n = \frac{2}{\pi} \int_0^\pi \cos \left( \frac{n \pi x}{\pi} \right) \sin (x) = \frac{1+\cos (n\pi)}{1-n^2} =\left\{\begin{array}{lll} \frac{4/\pi}{1-n^2} & \mbox{ for } & n \mbox{ even},\\ &&\\ 0 & \mbox{ for } & n \mbox{ odd}. \end{array} \right. \] Note, the $n=1$ case should done separately to avoid division by zero, but it comes out zero. Therefore, \[ u(x,t) = \frac{2}{\pi} + \frac{4}{\pi} \sum_{k=1}^\infty \frac{1}{1-4k^2} \cos (2kx) e^{-4k^2t}. \] On the next page are plots for various values of $t$ with the sum going to $k=30$. The green lines mark $2/\pi$. An animation is on the course website. \begin{center}{\bf Plots for Last Example} \end{center} \vspace{0.25in} \begin{center} \psfrag{t=0.0}{$t=0.0$} \includegraphics[width=1.8in]{Heat4plot2.eps} \psfrag{t=0.01}{$t=0.01$} \includegraphics[width=1.8in]{Heat4plot3.eps} \psfrag{t=0.02}{$t=0.02$} \includegraphics[width=1.8in]{Heat4plot4.eps} \vspace{0.5in} \psfrag{t=0.03}{$t=0.03$} \includegraphics[width=1.8in]{Heat4plot5.eps} \psfrag{t=0.05}{$t=0.05$} \includegraphics[width=1.8in]{Heat4plot6.eps} \psfrag{t=0.1}{$t=0.1$} \includegraphics[width=1.8in]{Heat4plot7.eps} \vspace{0.5in} \psfrag{t=0.2}{$t=0.2$} \includegraphics[width=1.8in]{Heat4plot8.eps} \psfrag{t=0.3}{$t=0.3$} \includegraphics[width=1.8in]{Heat4plot9.eps} \psfrag{t=0.6}{$t=0.6$} \includegraphics[width=1.8in]{Heat4plot10.eps} \end{center} \section{Summary for Heat Equation} {\bf End points held at 0 degrees (homogeneous case).} \vspace{0.25in} \noindent Model: $\alpha^2u_{xx} = u_t$, $u(0,t)=u(L,t)=0$, $u(x,0)=f(x)$. \vspace{0.25in} \noindent Solution: \[ u(x,t) = \sum_{n=1}^\infty c_n \sin \left( \frac{n\pi x}{L}\right) e^{-\frac{n^2\pi^2\alpha^2t}{L^2}} \] where \[ c_n = \frac{2}{L} \int_0^L \sin \left( \frac{n\pi x}{L}\right) f(x)\, dx. \] \vspace{0.25in} \hrule \vspace{0.25in} {\bf End points held at fixed temperatures (Nonhomogeneous case).} \vspace{0.25in} \noindent Model: $\alpha^2u_{xx} = u_t$, $u(0,t)=T_1$, $u(L,t)=T_2$, $u(x,0)=f(x)$. \vspace{0.25in} \noindent Solution: Let $v(x) = \dfrac{T_2-T_1}{L}x + T_1$. Then \[ u(x,t) = v(x) + \sum_{n=1}^\infty c_n \sin \left( \frac{n\pi x}{L}\right) e^{-\frac{n^2\pi^2\alpha^2t}{L^2}}, \] where \[ c_n = \frac{2}{L} \int_0^L \sin \left( \frac{n\pi x}{L}\right)[f(x)-v(x)]\, dx. \] \vspace{0.25in} \hrule \vspace{0.25in} {\bf Insulated end points.} \vspace{0.25in} \noindent Model: $\alpha^2u_{xx} = u_t$, $u_x(0,t)=0$, $u_x(L,t)=0$, $u(x,0)=f(x)$. \vspace{0.25in} \noindent Solution: \[ u(x,t) = c_0 + \sum_{n=1}^\infty c_n \cos \left( \frac{n\pi x}{L}\right) e^{-\frac{n^2\pi^2\alpha^2t}{L^2}} \] where \[ c_0 = \mbox{ the average value of $f(x)$ } = \frac{1}{L} \int_0^L f(x) \, dx \] and for $n= 1, 2, 3, ...$ \[ c_n = \frac{2}{L} \int_0^L \cos \left( \frac{n\pi x}{L}\right) f(x) \, dx. \] \section{Extra Credit!} Consider a metal ring. For a given initial temperature distribution $f(\theta)$ for $\theta \in [0,2\pi]$ set up a model with the heat equation and show how to solve it. Then work through an example and display the plots for several times. If you like, create an animation. \section{ The Heat Equation for a Square Plate: OPTIONAL READING} 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}. \] \subsection{Example} 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.15in} \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. \section{Heat Equation on a Metal Disk: OPTIONAL READING} 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 location $(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 is \[ \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. \] 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 $$\dis \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 this constant $-\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. \] Now, 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(r)}{dr} = \frac{p}{\gamma} \frac{d P(p)}{dr} = \frac{p}{\gamma}\frac{dP}{dp}\frac{dp}{dr} = \frac{p}{\gamma}\frac{dP}{dp} \gamma = 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\footnote{Elementary Differential Equations \& Boundary Value Problems 10th edition, by Boyce \& DiPrima}. 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. \subsection{Bessel's Equations} Below are graphs of $J_0(x)$ and $Y_0(x)$. \begin{center} \psfrag{Jo}{$J_0$}\psfrag{Yo}{$Y_0$} \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 in increasing order by $\gamma_1$, $\gamma_2$, $\gamma_3$, etc., all positive. 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 }} \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 like the Fourier expansion with $J_0(\gamma_n r)$ in stead of $\sin (n\pi r)$ or $\cos(n\pi r)$. It is called a {\bf Fourier-Bessel series}. 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}. \] You should compare this to the formulas for the Fourier sine and cosine series. The chief difference is the integral in the denominator. For $\sin^2 (n\pi x)$ and $\cos^2(n\pi x)$ the integral was independent of $n$ and was just folded into the constants. For $J^2_0(\gamma_n r)$ this is not the case and that is why is has to be included as a normalizing term. The $r$ in the integrals is a weighting factor that comes from using polar coordinates. \subsection{Example with Axial Symmetry} Suppose $f(r) = 1-r^2$. We expand $f(r)$ as a series of Bessel functions. Then we plot a few partial sums to see how good the approximations are. We define {\tt f(r)}, {\tt g(n)} and {\tt c(n)} with the Maple commands below; {\tt g(n)} is $\gamma_n$ and {\tt c(n)} is $c_n$. \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));}} \vspace{.2in} Here is the plot for $N=1$, that is we just polar the first term in the Fourier Bessel expansion. \vspace{.2in} {\tt \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);} } \pagebreak \begin{center} $t=0.0$ \hspace{1.5in} $t=0.01$ \hspace{1.5in} $t=0.02$ \includegraphics[height=1.9in]{HeatDiskplot6.eps} \includegraphics[height=1.9in]{HeatDiskplot7.eps} \includegraphics[height=1.9in]{HeatDiskplot8.eps} \vspace{.15in} $t=0.1$ \hspace{1.5in} $t=0.2$ \hspace{1.5in} $t=0.3$ \includegraphics[height=1.9in]{HeatDiskplot9.eps} \includegraphics[height=1.9in]{HeatDiskplot10.eps} \includegraphics[height=1.9in]{HeatDiskplot11.eps} \vspace{.15in} $t=1.0$ \includegraphics[height=1.9in]{HeatDiskplot12.eps} \end{center} \subsection{Extra Credit!} 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 it as best you can. \subsection{Example without Axial Symmetry} Suppose $\alpha=1$ and $f(r,\theta) = r(1-r^2)\sin^2(3\theta)$. For the graphs here and below I used cylindrical coordinates, but Maple default plots $r$ as a function of $\theta$ and $z$. Therefore, I had to define mine own cylindrical coordinates option where $z$ is plotted as a function of $r$ and $\theta$. \vspace{.2in} {\tt \noindent\red{> addcoords(zcylindrical,[z,r,theta],[r*cos(theta),r*sin(theta),z],}} {\tt \noindent\red{[[1],[Pi],[0],[0..2,0..2*Pi,-1..1],[-4..4,-4..4,-4..4]])}} \vspace{.2in} I also created my own RGB color scheme just for fun. \vspace{.2in} {\tt \noindent\red{B:= t -> piecewise(t<0.6, 1.0, -t/0.4 + 2.5 );}} {\tt \noindent\red{R:= t -> piecewise(t<0.4, t/0.4, 1.0);}} \vspace{.2in} Here then is a graph of $z=f(r,\theta)$. \vspace{.2in} {\tt \noindent\red{> f := (r,theta) -> r*(1-r\^{}2)*sin(3*theta)\^{}2;}} {\tt\noindent\red{> plot3d(f(r,theta), r=0..1, theta=0..2*Pi,coords=zcylindrical, \\view=0..0.5,color=,numpoints=10000);} } \begin{center} \includegraphics[height=1.9in]{HeatDisk2plot3d2.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 \cos \sqrt{\beta}\theta + C_2 \sin \sqrt{\beta}\theta. \] Therefore, $\beta = n^2$ for $n>0$. Thus, $1$, $\sin n\theta$, and $\cos n\theta$, will solve the $\Theta$ equation and have an acceptable period as would any linear combination of these. Now we go back to the $R$ equation. 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$. Its 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 by certain power series (see {\em Introduction to Bessel Functions}, by Frank Bowman). Again, $Y_n$ is unbounded at 0 and so $C_2=0$. Again $J_n$ has infinitely many zeros. We denote them by $\gamma_{nm}$, $m=1, 2, 3, ...$. They are positive and numbered in increasing order. Their numerical values can be found in tables or via a computer. The command in Maple for the $m$-th of $J_n$ is {\tt\red{BesselJZeros(n,m)}}. For each $n=0,1,2,3...$ and $m=1,2,3...$ let \[ R_{nm}(r) = J_n(\gamma_{nm}r). \] Then $R_{nm}(1) = 0$ and $R_{nm}(r)$ solves $r^2R'' +r R' +(\gamma_{nm}^2r^2 - n^2)R = 0$. And of course we define \[ T_{nm} = e^{-\gamma_{nm}^2 t}. \] Now, the products $R_{0m}(r)T_{0m}(r)$, $R_{nm}(r)T_{nm}(t)\cos (n\theta)$, and $R_{nm}(r)T_{nm}(t)\sin (n\theta)$, for $n, m >0$, solve the heat equation and satisfy the boundary condition of being $0$ when $r=1$, are bounded, and have period $2\pi$ in $\theta$. So too would any linear combination or convergent infinite sum. We suppose that, \[ u(r,\theta,t) = \frac{1}{2}\sum_{m=1}^\infty a_{0m}R_{nm}(r)T_{nm}(t) + \sum_{n=1}^\infty \sum_{m=1}^\infty a_{nm} R_{n,m}(r)T_{nm}(t)\cos (n\pi) + b_{nm} R_{nm}(r)T_{nm}(t)\sin (n\pi) \] converges. All that if left is to find values for the $a_{nm}$ and $b_{nm}$ 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, only the $a_{0m}$ and $a_{6m}$ terms will be nonzero. We write \[ 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 a_{0m} R_{0m}(r) = \sum_{m=1}^\infty a_{0m}J_0(\gamma_{0m}r) = r(1-r^2) \] and \[ \sum_{m=1}^\infty a_{6m} R_{6m}(r) = \sum_{m=1}^\infty a_{6m} J_6(\gamma_{6m}r) = r(1-r^2). \] Each family of functions $\{ J_n(\gamma_{nm}r) \}_{m=1}^\infty $ satisfies the orthogonality conditions needed to do series approximations. We get \[ a_{0m} = \frac{\int_0^1 r^2(1-r^2) J_0(\gamma_{0m}r) \, dr }{\int_0^1 r J_0^2 (\gamma_{0m} r) \, dr} \] and \[ a_{6m} = \frac{\int_0^1 r^2(1-r^2) J_6(\gamma_{6m}r) \, dr }{\int_0^1 r J_6^2 (\gamma_{6m} r) \, dr}. \] Next we plot $u(r,\theta,0)$ and compare it to the graph for $f(r,\theta)$ above. \vspace{0.2in} {\tt \noindent\red{> fr := r -> r*(1-r\^{}2);} \noindent\red{> g := (n,m) -> evalf(BesselJZeros(n,m));} \noindent\red{> a0 := m -> (int(r*fr(r)*BesselJ(0,g(0,m)*r),r=0..1))\\/(int(r*BesselJ(0,g(0,m)*r)\^{}2 ,r=0..1));} \noindent\red{> a6 := m -> (int(r*fr(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( a0(m)*BesselJ(0,g(0,m)*r) ,m=1..M) - \\1/2*cos(6*theta)*sum(a6(m)*BesselJ(6,g(6,m)*r), m=1..M);} \noindent\red{> M:=10:} \noindent\red{> plot3d(U0(r,theta),r=0..1,theta=0..2*Pi,coords=zcylindrical,view=0..0.5,\\ color=[R(U0(r,theta)/0.4), 0, B(U0(r,theta)/0.4)],numpoints=10000);} } \begin{center} \includegraphics[height=1.9in]{HeatDisk2plot3d3.eps} \end{center} Now we put all the pieces together to get \[ u(r,\theta,t) = \frac{1}{2} \sum_{m=1}^\infty a_{0m} J_0(\gamma_{nm}r) e^{-\gamma_{0m}^2 t } - \frac{1}{2}\cos(6\theta) \sum_{m=1}^\infty a_{6m} J_6(\gamma_{nm}r) e^{-\gamma_{6m}^2 t }. \] We plot this for a few values of $t$. \begin{center} \psfrag{t=0.001}{\hspace{-0.2in}$t=0.001$} \includegraphics[width=1.9in]{HeatDisk2plot3d4.eps} % \psfrag{t=0.003}{\hspace{-0.2in}$t=0.003$} \includegraphics[width=1.9in]{HeatDisk2plot3d5.eps} % \psfrag{t=0.005}{\hspace{-0.2in}$t=0.005$} \includegraphics[width=1.9in]{HeatDisk2plot3d6.eps} \psfrag{t=0.01}{\hspace{-0.2in}$t=0.01$} \includegraphics[width=1.9in]{HeatDisk2plot3d7.eps} % \psfrag{t=0.02}{\hspace{-0.2in}$t=0.02$} \includegraphics[width=1.9in]{HeatDisk2plot3d8.eps} % \psfrag{t=0.05}{\hspace{-0.2in}$t=0.05$} \includegraphics[width=1.9in]{HeatDisk2plot3d9.eps} \psfrag{t=0.1}{\hspace{-0.1in}$t=0.1$} \includegraphics[width=1.9in]{HeatDisk2plot3d10.eps} \end{center} \subsection{General Case} In the last example we took advantage of the fact the given $f(r,\theta)$ was separable. This is not always the case. I found the online lecture notes ``NOTE12: Fourier-Bessel Series and BVP in Cylindrical Coordinates'' by Hamid Meziani, 2016, especially helpful at this point. \\ See {\tt http://faculty.fiu.edu/$\sim$meziani/MAP4401.html}. \[ u(r,\theta,t) = \frac{1}{2}\sum_{m=1}^\infty a_{0m} R_{0m}(r)T_{0m}(t) + \sum_{n=1}^\infty \sum_{m=1}^\infty a_{nm} R_{nm}(r)T_{nm}(t)\cos (n\theta) + b_{nm} R_{nm}(r)T_{nm}(t)\sin (n\theta) \] \[ = \sum_{n=1}^\infty \sum_{m=1}^\infty a_{nm} J_n(\gamma_{nm}r)e^{-\gamma_{nm}^2 t}\cos (n\theta) + b_{nm} J_n(\gamma_{nm}r)e^{-\gamma_{nm}^2 t}\sin (n\theta). \] At time $t=0$ we need \[ u(r,\theta,0) = \frac{1}{2}\sum_{m=1}^\infty a_{0m} J_0(\gamma_{0m}r) + \sum_{n=1}^\infty \sum_{m=1}^\infty a_{nm} J_n(\gamma_{nm}r)\cos (n\theta) + b_{nm} J_n(\gamma_{nm}r)\sin (n\theta) = 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 with respect to $\theta$. 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_{0}^{2\pi} f(r,\theta) \cos (n\theta) \, d\theta, \] and \[ B_n(r) = \frac{1}{\pi}\int_{0}^{2\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 Fourier-Bessel series in terms of any $\{J_k(\gamma_{km}\|_{m=1}^\infty$, so we 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}. \] Combining these with the integral formulas for $A_n(r)$ and $B_n(r)$ gives the following. \[ a_{nm} = \frac{1}{\pi}\frac{\int_0^1 \int_{0}^{2\pi} r f(r,\theta) J_n(\gamma_{nm} r) \cos (n\theta)\, d\theta \, dr}{\int_0^1 r J_n^2(\gamma_{nm} r) \, dr} \] and \[ b_{nm} = \frac{1}{\pi} \frac{\int_0^1 \int_{0}^{2\pi} r f(r,\theta) J_n(\gamma_{nm} r) \sin (n\theta)\, d\theta \, dr}{\int_0^1 r J_n^2(\gamma_{nm} r) \, dr}. \] \subsection{Example: Not Axially Symmetric and Not Separable} Let $f(r,\theta) = r(1-r)e^{r\sin \theta}$. Its graph is below. I did not use color because the plots took too long. \vspace{.2in} {\tt \noindent\red{>f := (r,theta) -> r*(1-r)*exp(r*sin(theta));}} {\tt \noindent\red{> plot3d(f(r,theta),r=0..1,theta=0..2*Pi,coords=zcylindrical,color=white,view=-0.001..0.5)}} \begin{center} \includegraphics[height=2.5in]{HeatDisk3plot3d1.eps} \end{center} Here are the commands I used to construct the Fourier-Bessel series. First I defines the coefficients. Here ${\tt g(n,m)}$ is the $m$-th zero of $J_n$. \vspace{.2in} {\tt \noindent\red{> g := (n,m) -> evalf(BesselJZeros(n,m)); } \noindent\red{> a := (n,m) -> evalf((1/Pi)*int(r*BesselJ(n,g(n,m)*r)*int( f(r,theta)*cos(n*theta),\\theta=0..2*Pi), r=0..1)/ int(r*BesselJ(n,g(n,m)*r)\^{}2,r=0..1)); } \noindent\red{> b := (n,m) -> evalf((1/Pi)*int(int( f(r,theta)*sin(n*theta)*r*BesselJ(n,g(n,m)*r),\\theta=0..2*Pi),r=0..1)/ int(r*BesselJ(n,g(n,m)*r)\^{}2,r=0..1));}} \vspace{.2in} Because these calculations take so long I stored the results in two data matrices, $A$ and $B$. Since $n$ starts at 0 we get $A[1,1] = a(0,1)$, and so on. The {\tt sum} does not like entries that are from a matrix, so I had to use the {\tt add} command instead. Next we plot $u(r,\theta,0)$ and compare with the plot of $f(r,\theta)$. I ended up only needing a few terms. \vspace{.2in} {\tt \noindent\red{> N:=2: M:=2:}} {\tt \noindent\red{> U0:= (r,theta) -> add(BesselJ(0,g(0,m)*r)*A[1,m] ,m=1..M)/2 + \\add(add(BesselJ(n,g(n,m)*r)*(A[n+1,m]*cos(n*theta) + B[n+1,m]*sin(n*theta)),m=1..M ), n=1..N );}} {\tt \noindent\red{> plot3d(U0(r,theta),r=0..1,theta=0..2*Pi,coords=zcylindrical,numpoints=1000,\\color=white,view=-0.001..0.5); }} \begin{center} \includegraphics[height=2.5in]{HeatDisk3plot3d2.eps} \end{center} This looks pretty close. Next I defined our solution $u(r,\theta,t)$. Then I plotted in for several values of $t$. \vspace{.2in} {\tt \noindent\red{> U := (r,theta,t) -> add(BesselJ(0,g(0,m)*r)*A[1,m]*exp(-g(0,m)\^{}2*t) ,m=1..M)/2 +\\ add(add(BesselJ(n,g(n,m)*r)*exp(-g(n,m)\^{}2*t)*(A[n+1,m]*cos(n*theta) +\\ B[n+1,m]*sin(n*theta)),m=1..M ), n=1..N );} %\noindent\red{> t:=0.01;} \noindent\red{> plot3d(U(r,theta,0.01),r=0..1,theta=0..2*Pi,coords=zcylindrical,numpoints=1000,\\color=white,view=-0.001..0.5); } } \begin{center} $t=0.01$ \hspace{2in} $t=0.02$ \includegraphics[height=2.5in]{HeatDisk3plot3d3.eps} \includegraphics[height=2.5in]{HeatDisk3plot3d4.eps} $t=0.05$ \hspace{2in} $t=0.1$ \includegraphics[height=2.5in]{HeatDisk3plot3d5.eps} \includegraphics[height=2.5in]{HeatDisk3plot3d6.eps} \end{center} Finally I did an animation. I'll put it on the web site. Here is the command that generates it. \vspace{.2in} {\tt\noindent\red{> with(plots):}} {\tt\noindent\red{> animate(plot3d ,[ U(r,theta,a),r=0..1,theta=0..2*Pi,coords=zcylindrical, \\numpoints=1000,color=white,view=-0.001..0.5 ] ,a=0.0..0.2,frames=200);}} \subsection{Example: Discontinuous Initial Condition} Suppose the initial temperature distribution is 1 inside a disk of radius 1/2 with center (0,1/2) and 0 elsewhere on the unit disk. Let \[ f(r,\theta) = \left\{ \begin{array}{ll} 1 & 0 \leq \theta \leq \pi \,\&\, r < \sin \theta \\ 0 & \mbox{ otherwise} \end{array}\right. \] See the graph below. \vspace{.2in} {\tt \noindent\red{> plot3d(f,0..1,0..2*Pi,coords=zcylindrical,numpoints=5000) \#The default coloring has no physical significance. } } \begin{center} \includegraphics[width=2.5in]{HeatDIsk4plot3d1.eps} \end{center} Next we compute the coefficients and plot $u(r,\theta, 0)$. \vspace{.2in} {\tt \noindent\red{> g := (n,m) -> evalf(BesselJZeros(n,m));} \noindent\red{> a := (n,m) -> evalf((1/Pi)*evalf(Int(evalf(Int( 1*r*BesselJ(n,g(n,m)*r)*cos(n*theta),\\ r=0..sin(theta))), theta=0..Pi))/int(r*BesselJ(n,g(n,m)*r)\^{}2,r=0..1));} \noindent\red{> b := (n,m) -> evalf((1/Pi)*evalf(Int(evalf(Int( 1*r*BesselJ(n,g(n,m)*r)*sin(n*theta),\\ r=0..sin(theta))), theta=0..Pi))/int(r*BesselJ(n,g(n,m)*r)\^{}2,r=0..1));} \noindent\red{> U0:= (r,theta) -> add(BesselJ(0,g(0,m)*r)*(a(0,m)),m=1..M)/2 + \\add(add(BesselJ(n,g(n,m)*r)*(a(n,m)*cos(n*theta) \\ + b(n,m)*sin(n*theta)),m=1..M ), n=1..N );} \noindent\red{> N:=4: M:=4:} \noindent\red{>plot3d(U0(r,theta),r=0..1,theta=0..2*Pi,coords=zcylindrical,color=white); } } \begin{center} \includegraphics[width=2.5in]{HeatDIsk4plot3d2.eps} \end{center} Unsurprisingly, convergence is much slow where the initial function is discontinuous. We plot several more partial sums: $N=M=5$, $N=M=6$, $N=8, M=6$, and $N=M=9$. These plots take up a lot of computer time. \begin{center} \includegraphics[width=2.5in]{HeatDIsk4plot3d3.eps} \includegraphics[width=2.5in]{HeatDIsk4plot3d4.eps} \includegraphics[width=2.5in]{HeatDIsk4plot3d5.eps} \includegraphics[width=2.5in]{HeatDIsk4plot3d6.eps} \end{center} Next we define $u(r,\theta,t)$ and plot it for several values of $t$. \vspace{.2in} {\tt\noindent\red{> U := (r,theta,t) -> add(BesselJ(0,g(0,m)*r)*exp(-g(0,m)\^{}2*t)*(a(0,m)),m=1..M)/2 + \\add(add(BesselJ(n,g(n,m)*r)*exp(-g(n,m)\^{}2*t)*(a(n,m)*cos(n*theta) + \\b(n,m)*sin(n*theta)),m=1..M ), n=1..N );}} {\tt\noindent\red{> N:=4: M:=4: t:=0.01;}} {\tt\noindent\red{> plot3d(U(r,theta,t),r=0..1,theta=0..2*Pi,coords=zcylindrical,color=white);}} \vspace{.2in} Also shown are the plots for $t= 0.02, 0.03$ and $0.06$. \begin{center} \includegraphics[width=2.5in]{HeatDIsk4plot3d7.eps} \includegraphics[width=2.5in]{HeatDIsk4plot3d8.eps} \includegraphics[width=2.5in]{HeatDIsk4plot3d9.eps} \includegraphics[width=2.5in]{HeatDIsk4plot3d10.eps} \end{center} \section{Discussion and Further Reading} Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations Partial Differential Equations \end{document}