%% 353labs.tex -- a LaTeX source file Spring 1996. %% "linalg.sty" is a special style file. Get it with this source! %% Usually location: /u/ma/gustafson/tex/texinput/linalg.sty %% \documentstyle[11pt,linalg,multicol,topics]{report} %\solutionsfalse% DON'T PRINT SOLUTIONS \solutionstrue% UNCOMMENT to print solutions. Don't change previous! \def\SCHEDULE{% \begin{itemize} % % WARNING: Look at a printout before trying to choose the problems. % The numbers below are historic, not the ones of a printout. % Usually, only collect for 9 weeks, I-IX only. % Read below for making changes to this file!!!!!! \small\sf \item[] Take-home I [1,2,3,4]: Tuesday of week 2 \item[] Take-home II [5,6,7,8]: Tuesday of week 3 \item[] Take-home III [9,10,11,12] Tuesday of week 4 \item[] Take-home IV [14,15,16,17]: Tuesday of week 5 \item[] Take-home V [18,19,20,21,22]: Tuesday of week 6 \item[] Take-home VI [23,24,25,26,27] Tuesday of week 7 \item[] Take-home VII [28,29,30,31,32]: Tuesday of week 8 \item[] Take-home VIII [33,34,35,36]: Tuesday of week 9 \end{itemize} } %% %% INSTRUCTIONS: %% To modify a problem, change both the problem and its solution, %% which are located in the order PROBLEM+SOLUTION, with no other %% problems intervening. %% %% To move a problem, move both the problem and solution to the new %% location. All numbering is taken care of by the macros. If you %% want to keep the current numbering scheme, then a new problem %% should be moved into the place previously occupied. %% %% To demote a problem to OPTIONAL, with numbering like 1a, 2c, etc, %% change the "\prob{...}{...}" to "\probx{..}{..}". That's right, %% add "x". %% %% To promote an optional problem to a REQUIRED problem, remove the %% "x", i.e., "\probx{..}{..}" becomes "\prob{...}{...}". %% %% NOTES are added using the "\NOTES{...}" macro. They will not be %% printed unless "\notestrue" appears at the top of the source. %% The definition made at the top of this source allows the NOTES to %% contain verbatim material. %% % % BOOLEAN set in "linalg.sty" for switching solutions in/out % Normally, this is FALSE to print just the problems. % \ifsolutions\long\gdef\NOTES{}\else\long\gdef\NOTES#1{}\fi \pagestyle{plain} \begin{document} \thispagestyle{empty} \setcounter{page}{1} \ifsolutions\else\begin{multicols}{2}\fi \begin{center} {\Large\bf Math 353}\\ {\Large\bf Applied Partial Differential Equations}\\ %{\Large\bf \ifsolutions Problem Notes on \fi Take-Home Exams I to VIII} {\Large\bf \ifsolutions Problem Notes on \fi Take-Home Exams} \end{center} % Collection of the take-home exam will be made in weekly increments. % WARNING: Look at a printout before trying to choose the problems. % The numbers below are historic, not the ones of a printout. % Usually, only collect for 8 weeks, I-VIII only. % \ifsolutions\else \SCHEDULE \fi \bigskip {\large\bf Periodic Functions}. Let the function $f_0(x)$ be defined on $[-\pi,\pi]$ by the formula $$ f_0(x)=\left\{ \begin{array}{rl} \displaystyle 0,&\text{~~for~~} \displaystyle -\pi \le x < 0,\\ \displaystyle \sin(x),&\text{~~for~~} \displaystyle 0 \le x \le \pi. \end{array}\right. $$ \begin{exercises} \prob{1.}{(Periodic Extension)} Define $f(x)$ to be the $2\pi$--periodic extension of $f_0(x)$ to the whole real line. Give a formula for $f(x)$ in terms of $f_0(x)$ valid for $[-3\pi,5\pi]$. The formula should depend only upon the symbols $f_0$ and $2\pi$ and not upon the actual piecewise definition of the function $f_0$. \SOL{1.}{ The formula is a piecewise defined function which equals $f_0(x)$ on $[-\pi,\pi]$ and on the intervals $[-3\pi,-\pi]$, $[\pi,3\pi]$, $[3\pi,5\pi]$ takes the function values $f_0(x+2\pi)$, $f_0(x-2\pi)$, $f_0(x-4\pi)$, respectively. All of this is verified graphically. An expected and essential part of the solution is to chase down a reference to {\em translations} $f_0(x+c)$ in a calculus textbook. According to the theory of rigid motions, each {\em panel} (a xerox copy of the base function's graph) has a corresponding equation $y=f_0(x+c)$ for some $c$. The technical part is to decide on the values of $c$, then re-assemble the equations, one for each panel, into a piecewise defined function. } \end{exercises} \NOTES{ {\large\bf Periodic Extension}. The idea of a {\bf periodic extension} is fundamental to the discussion of Fourier series because the answer produced by a Fourier series computation is itself a periodic extension. A function $f_0$ that generates a periodic extension $f$ is called {\bf a base function}. A given $T$-periodic function $f$ has an infinite number of base functions, but given one, it will generate the graph $f$. Conversely, given only a $T$-periodic function $f$, we can always take a base function $f_0$ to be the graph of $f$ restricted to the interval $0 \le t \le T$. {\bf Example}. Write down for $[-\pi,2\pi]$ the periodic extension $f$ of period $\pi$ for the base function $$f_0(x) = \left\{ \begin{array}{rl} \sin(x) & 0 \le x \le \pi/2, \\ 2(\pi - x)/\pi & \pi/2 < x \le \pi. \end{array} \right. $$ {\bf Solution}: $$f(x) = \left\{ \begin{array}{ll} f_0(x+\pi) & -\pi \le x < 0, \\ f_0(x) & 0 \le x < \pi, \\ f_0(x-\pi) & \pi \le x < 2\pi \end{array} \right. = \left\{ \begin{array}{ll} \sin(x+\pi) & -\pi \le x < -\pi/2, \\ 2(- x)/\pi & -\pi/2 \le x < 0,\\ \sin(x) & 0 \le x \le \pi/2, \\ 2(\pi - x)/\pi & \pi/2 \le x < \pi. \\ \sin(x-\pi) & \pi \le x < 3\pi/2,\\ 2(2\pi - x)/\pi & 3\pi/2 \le x < 2\pi. \end{array} \right. $$ The second formula is more specific than the first but neither are readily usable for computing or graphing. There is more than one way to prepare an extension formula for plotting purposes. } \begin{exercises} \probx{1a.}{(Extension Formulas)} Let $f_1(x)=|\sin(x)|$ on $[-\pi,0]$. Define $g(x)$ to be the even $2\pi$--periodic extension of $f_1(x)$ to the whole real line. Define $h(x)$ to be the odd $2\pi$--periodic extension of $f_1(x)$ to the whole real line. Give simple algebraic formulas for $g(x)$ and $h(x)$. \SOL{1a.}{ Answers: $g(x) = |\sin(x)|$, $h(x)= -\sin(x)$. The answers were obtained graphically by identifying the base function on $[-\pi,\pi]$ for each of $g$ and $h$. Odd graphs are formed by drawing the graph on $x\ge 0$, then rotate that graph $90^\circ$ clockwise and flip the graph over. The graph first drawn will be on the back side of the paper, showing through on the half-axis $x\le 0$. } \global\edef\FIRSTPROB{\LASTPROB}% \LASTPROB==last problem number \prob{2.}{(Graphing of Extensions)} Let $f(x)$ be the $2\pi$--periodic extension of $f_0(x)$ to $(-\infty,\infty)$. Graph $f(x)$ for $-4\pi sin(u)*u; f := x -> periodicextension(x,f0,-Pi,Pi); # This plot code works in maple V.1, V.2 and V.3 plot(f,-Pi..3*Pi); # Test use in numeric integration. evalf(Int(f,x=-Pi..3*Pi)); \end{verbatim} \normalsize\rm The {\tt proc} command in {\tt maple} creates a function subroutine. A simple example is {\tt f:= x -> sin(x)*cos(x):} which is coded with {\tt proc} as {\tt f:=proc(x) sin(x)*cos(x): end:}, that is, one is equivalent to the other. The variable $x$ in this coding is a {\tt dummy variable} or {\tt formal argument} to the function {\tt f}. The command {\tt Int} differs from {\tt int} in one upper case letter. It is not a mistake. The upper case version does not apply any integration algorithms. That will save time since no symbolic answer is requested. The {\tt evalf} command causes the integration to be done numerically by a clever adaptive scheme kept secret from most persons who use {\tt maple}. To find out the secrets, read the {\tt maple} help file on {\tt userinfo}. {\bf Maple plot error messages}. The plot command {\tt plot(f(x),0..1)} is often tried by beginners. It usually returns an {\tt empty plot} or in some cases the silly message {\tt can't evaluate boolean}. The fix: use {\tt plot(f,0..1)} or {\tt plot(f(x),x=0..1)}. \medskip {\bf Example}. Let $f_0(x)=\sin(x)$ for $0\le x<\pi$, $f_0(x)= -\sin(2x)$ for $-\pi \le x < 0$. Graph its $2\pi$--periodic extension on $[-\pi,4\pi]$. {\bf Solution}: The maple code for procedure {\tt periodicextension} will be used by reading it from the maple library, instead of copying it into the source. The twist is in the definition of $f_0(x)$, which is solved by using the Heaviside Unit Step function: $$ H(t) \equiv \left\{ \begin{array}{rl} 1 & t \ge 0, \\ 0 & t < 0, \end{array}\right. ~~~~~ H(t-a)-H(t-b) = \left\{ \begin{array}{rl} 1 & a \le t < b, \\ 0 & \mbox{otherwise}. \end{array}\right. $$ The second formula above can be justified directly from the definition of $H(t)$. In {\tt maple} the Heaviside function $H(t)$ is coded as {\tt Heaviside(t)}. The method for writing $f_0(x)$ in terms of Heaviside functions is as follows (it applies to any piecewise defined function): $$ f_0(x) \equiv \left\{ \begin{array}{rl} \sin(x) & 0 \le x < \pi \\ -\sin(2x) & -\pi \le x < 0 \end{array}\right. ~=~ \sin(x) \left\{ \begin{array}{rl} 1 & 0 \le x < \pi \\ 0 & \mbox{otherwise} \end{array}\right. ~+~ (-\sin(2x)) \left\{ \begin{array}{rl} 1 & -\pi \le x < 0 \\ 0 & \mbox{otherwise} \end{array}\right. $$ Therefore, from the Heaviside formulas, $$ f_0(x) = \sin(x)(H(x-0)-H(x-\pi)) + (-\sin(2x))(H(x-(-\pi))-H(x-0)). $$ This justifies the mysterious definition of $f_0$ in the maple code below. It explains the mechanism for coding piecewise defined functions in preparation for graphing. \footnotesize\rm \begin{verbatim} # maple V.1, V.2 and V.3. See above for periodic extension procedure read`/u/cl/maple/periodext`: readlib(Heaviside): H:=x -> Heaviside(evalf(x)): f0:= u -> -sin(2*u)*(H(u+Pi)-H(u)) + sin(u)*(H(u)-H(u-Pi)): f := x -> periodicextension(x,f0,-Pi,Pi); plot(f,-Pi..4*Pi); \end{verbatim} \normalsize\rm } {\large\bf Fourier Series}. Let the function $f_0(x)$ be defined on $0 \le x < 2\pi$ by the formula $$ f_0(x)=\left\{ \begin{array}{rl} \displaystyle -1,&\text{~~for~~} \displaystyle 0 \le x < 2\pi/3\\ \displaystyle 0,&\text{~~for~~} \displaystyle 2\pi/3 \le x < 4\pi/3\\ \displaystyle 1,&\text{~~for~~} \displaystyle 4\pi/3 \le x < 2\pi. \end{array}\right. $$ \begin{exercises} \prob{3.}{(Graphing Periodic Extensions)} Define $f(x)$ to be the $2\pi$--periodic extension of $f_0(x)$ to the whole real line. Sketch by hand the graph of $f$ for $-2\pi Heaviside(evalf(t)): f0:= x -> 3*x*H(x)-x*H(x-1)-2*x*H(x+1): f:= x -> f0(x+2)*H(x+2)+(f0(x)-f0(x+2))*H(x+1)+ (f0(x-2)-f0(x))*H(x-1)-f0(x-2)*H(x-3): plot(f,-1..3); \end{verbatim} \normalsize\rm } \begin{exercises} \prob{4.}{(Fourier Series)} Let $f$ be the $2\pi$--periodic extension of $f_0$ to $(-\infty,\infty)$. Show that $f$ is {\em odd}, $f(x)=-f(-x)$, except possibly at jumps of $f$, and calculate the coefficients $a_n$ and $b_n$ for the Fourier series of $f$. \SOL{4.}{ The property of oddness can be shown on the symmetric interval $[-\pi,\pi]$ because $f$ is by definition $2\pi$--periodic. Since $f$ is zero in this symmetric interval except for $|x| \le 2\pi/3$, we only need to test oddness on $|x|\le 2\pi/3$. A quick graph of $f$ on this interval shows that $f$ is odd (rotate $90$ degrees clockwise and flip the graph over). The details can also be chased down with algebra, because $f(x) = 1$ on $-2\pi/3 \le x < 0$ and $f(x)=-1$ on $0\le x < 2\pi/3$. The Fourier series of an odd function consists entirely of sine terms, therefore $a_n=0$. The coefficients $b_n$ are given by the formula $b_n=\frac{1}{\pi}\int_{-\pi}^\pi f_0(x)\sin(nx)dx =\frac{2}{\pi}\int_{0}^{2\pi/3} (-1)\sin(nx)dx=\frac{2}{\pi}\frac{1}{n}(\cos(2n\pi/3)-1)$. The reduction to twice the integral over $[0,2\pi/3]$ results from the integrand being {\em even}. On this new interval, $f_0(x)=-1$, except for the points $x=2\pi/3$ to $x=\pi$, where $f_0(x)=0$. } \probx{4a.}{(Sum of a Fourier Series)} Let $f$ be the $2\pi$--periodic extension of $f_0$ to $(-\infty,\infty)$. Determine the sum of the Fourier series of $f$ for each value of $x$. \SOL{4a.}{ The series was calculated in the previous problem. Fourier's convergence theorem, if it applies, says that the series converges not to $f(x)$, but to $F(x)\equiv\frac12(f(x+0)+f(x-0))$ for each $x$. This answer agrees with $f$ except at $x=\pm \pi/3$ and $x=0$ (and $2\pi$ multiples of these). The result is that the sum of the Fourier series, denoted $F(x)$, satisfies $F(x)=f(x)$ except at these points, where $F(\pi/3)=-1/2$, $F(-\pi/3)=1/2$, $F(0)=0$. } \probx{4b.}{(Convergence Theorem)} Let $f$ be the $2\pi$--periodic extension of $f_0$ to $(-\infty,\infty)$. For which values of $x$ does the graph of $f(x)$ differ from the graph of the Fourier series of $f$? Reference a theorem. \SOL{4b.}{ Let $F$ be the Fourier series of $f$. Then Fourier's convergence theorem says that $F(x)=f(x)$ at points of continuity of $f$. Specific information is available at the points of discontinuity of $f$. At $x=0$, $F(x)=0 \ne f(x)$; at $x=\pi/3$, $F(x)=-1/2 \ne f(x)$ and at $x=-\pi/3$, $F(x)=1/2 \ne f(x)$. Therefore, the graphs differ at $x=0$, $x=\pm\pi/3$ and all $2\pi$--multiples of these points. } \end{exercises} \NOTES{ {\large\bf Evaluating Fourier coefficients}. The evaluation of $a_n$ and $b_n$ amounts to evaluation of integrals. {\bf Example}. Evaluate the Fourier coefficients $a_n$ and $b_n$ for the $2\pi$--periodic function $f$ determined by the base function $$f_0(x) = \left\{ \begin{array}{rl} 1 & -\pi \le x < 0 \\ -1 & 0 \le x < \pi \end{array} \right. $$ {\bf Calculus Solution}: The integrals can be split up as $\int_{-\pi}^0$ plus $\int_0^\pi$. In the first integral $f$ is $1$ and in the second integral $f$ is $-1$. For example, $$ b_n=\frac{1}{\pi}\int_{-\pi}^0 \cos(nx)(1)dx + \frac{1}{\pi}\int_0^{\pi} \cos(nx)(-1)dx =\frac{1-\cos n\pi}{n\pi}. $$ The answers are $a_0=0$, $a_n=0$, $b_n = (1-\cos(n\pi))/(n\pi)$. Therefore, the Fourier Series of $f$ is $$ \sum_{n=1}^\infty \frac{1-\cos(n\pi)}{n\pi}\sin(n\pi x)$$ {\bf Maple Solution}: The function $f$ can be written in terms of Heaviside's function $H(t)$ as $f(x)=(1)(H(x+\pi)-H(x-0)) + (-1)(H(x-0)-H(x-\pi)) = H(x+\pi) - H(x-\pi) -2H(x)$. The integrals to evaluate are: $$a_n = \frac{1}{\pi}\int_{-\pi}^{\pi} f(x)\cos(nx)dx,~~~ b_n = \frac{1}{\pi}\int_{-\pi}^{\pi} f(x)\sin(nx)dx. $$ Using {\tt maple} as the integral table: \footnotesize\rm \begin{verbatim} # Maple V.2 readlib(Heaviside): H:= t -> Heaviside(evalf(t)): f:= x -> H(x+Pi)-H(x-Pi)-2*H(x): p1:=int(H(x+Pi)*cos(n*x),x=-Pi..Pi); q1:=int(H(x+Pi)*sin(n*x),x=-Pi..Pi); p2:=int(H(x-Pi)*cos(n*x),x=-Pi..Pi); q2:=int(H(x-Pi)*sin(n*x),x=-Pi..Pi); p3:=int(H(x)*cos(n*x),x=-Pi..Pi); q3:=int(H(x)*sin(n*x),x=-Pi..Pi); a0:=int(f(x),x=-Pi..Pi); # assume n>0 hereafter an:=simplify((p1-p2-2*p3)/Pi); bn:=simplify((q1-q2-2*q3)/Pi); \end{verbatim} \normalsize\rm The reason for the strange disassembly of the function into terms with a single Heaviside function has to do with the integration routines in {\tt maple}. It works when coded as above. It fails for the composite {\tt int(f(x)*cos(n*x),x=-Pi..Pi)}. The answers agree with those reported for the hand calculation, above. {\bf Example}. Determine the values of $x$ for which the graph of $f(x)$ differs from the graph of its formal Fourier series. {\bf Solution}: The {\bf sum of the Fourier series} is the formal series given in fourier's convergence theorem. The content of the theorem is that the sum is {\em precisely} $(f(x+0)+f(x-0))/2$. Since $f$ is a step function, the sum is $f(x)$ except at the jumps. For example, at jump $x=0$, the series sums to $(1+ (-1))/2=0$. At $x=-1$, the series sums to $((-1)+(1))/2=0$. So the graphs of $f$ and the Fourier series differ at $x=n$, $n=0,\pm 1, \pm 2, \ldots$. } {\large\bf Series Identities}. Apply the theory of Fourier series to prove the following identities. You may cite any Fourier series formula without proof (a textbook and page number will suffice). Include all details in the application of the theory (theorem used, value of $x$ substituted). \NOTES{ {\bf Example}. Verify the formula $$\frac{\pi}{4}= 1-\frac13+\frac15-\frac17 +- \cdots.$$ {\bf Solution}: To verify it, use the Fourier series for the periodic extension $f$ of the {\em even} base function $f_0(x)=1$, $-\pi/2 < x < \pi/2$, $f_0(x)=0$, $\pi/2 < x < 3\pi/2$. Then substitute $x=0$ in the series formula $$ \frac{f(x+0)+f(x-0)}{2} = \frac12+\frac{2}{\pi}\left(\cos x - \frac13\cos 3x + \frac15\cos 5x - \cdots\right)$$ The result is $$ 1 = \frac12+\frac{2}{\pi}\left(\cos 0 - \frac13\cos 0 + \frac15\cos 0 - \cdots\right)$$ which reduces to the claimed formula. } \begin{exercises} \prob{5.}{(Series Identity)} $\dd \frac{\pi^2}{8} = \sum_{n=0}^\infty \frac{1}{(2n+1)^2}$. \SOL{5.}{ To prove the {\bf series identity} requires an explicit Fourier series expansion for the function $f(x)=|x|$ on $-\pi \le x \le \pi$, namely, $$|x| = \frac{\pi}{2} - \frac{4}{\pi}\sum_{n=0}^\infty \frac{1}{(2n+1)^2}.$$ The idea is to apply the {\bf Fourier Convergence Theorem} at some particular value of $x$ in $-\pi < x < \pi$, in order to evaluate the infinite series. Suitable values are those for the standard angles in trigonometry: $x=0,\pi/6,\pi/4,\pi/3,\pi/2$. The Fourier series of $f$ converges always to $$ \frac12\lim_{t \to x+} f(t) + \frac12\lim_{t \to x-} f(t),$$ which reduces to $f(x)$ at a point of continuity of $f$. The $f$ in this formula is the {\em sawtooth wave}, formed from the base function $|x|$ on the interval $[-\pi,\pi]$. } \probx{5a.}{(Series Identity)} $\dd \frac{\pi^2}{6} = \sum_{n=1}^\infty \frac{1}{n^2}$. \end{exercises} \SOL{5a.}{ This formula results by substitution of some value of $x$ into the Fourier series for $f(x)=x^2$ on $-\pi < x < \pi$. } {\large\bf Half Range Expansions}. \NOTES{ The formulas for {\bf Half-Range expansions} are identical to the standard Fourier coefficient formulas applied to the even or odd extension of the given function. In practise, this is hidden because of efficiency considerations, which dictate that coefficients be calculated from the Half-Range formulas. For an even function $f$ defined on $[0,T]$ the formulas are $$f = \frac{a_0}{2} + \sum_{n=1}^\infty a_n\cos(n\pi x/T),~~~ a_n = \frac{2}{T} \int_0^T f(x)\cos(n\pi x/T)dx$$ {\bf Example}. Find the first three coefficients of the even Half-Range expansion for $f(x)=\sin(x)$ on $[0,\pi]$ and graph on $[-\pi,\pi]$ the Fourier series truncated to three terms. {\bf Maple Solution}: Apply the previous formulas as follows: \footnotesize\rm \begin{verbatim} f:=x -> sin(x): T:=Pi: an:=(2/T)*int(f(x)*cos(n*PI*x/T),x=0..T); a0:=limit(an,n=0); a1:=limit(an,n=1); a2:=limit(an,n=2); g:=x -> a0/2+a1*cos(Pi*x/T)+a2*cos(2*Pi*x/T); plot(g,-Pi..Pi); \end{verbatim} \normalsize\rm The answers: $a_0=4/\pi$, $a_1=0$, $a_2= -4/(3\pi)$. The truncated series is $g(x)=4(1 -\cos(2x)/3)/\pi$. The {\tt limit} command in {\tt maple} is used here to resolve fractions that might otherwise cause a division by zero under direct substitution. } \medskip {\large\bf Function $f_1(x)$}. Let the function $f_1(x)$ be defined on $[0,\pi]$ by the formula $$ f_1(x)=\left\{ \begin{array}{rl} \displaystyle \sin(x),&\text{~~for~~} \displaystyle 0 \le x \le \pi/2,\\ \displaystyle 2(\pi-x)/\pi,&\text{~~for~~} \displaystyle \pi/2 \le x \le \pi. \end{array}\right. $$ \begin{exercises} \prob{6.}{(Half Range Expansions)} Define $g(x)$ to be the even $2\pi$--periodic extension of $f_1(x)$ to the whole real line. Define $h(x)$ to be the odd $2\pi$--periodic extension of $f_1(x)$ to the whole real line. Give formulas for $g(x)$ and $h(x)$ valid on $[2\pi,4\pi]$ and find the first three Fourier coefficients of the even half-range expansion for $f_1(x)$. \SOL{6.}{ The answers to the problem may be found by direct integration starting with the formula for $n\ge 1$: $$a_n = \frac{2}{\pi} \int_0^\pi f_1(x)\cos(n x)dx.$$ To integrate, split the integral into two pieces, one over $[0,\pi/2]$ and the other over $[\pi/2,\pi]$. This generates two calculus integrals to be evaluated, each of which offers respectable resistance. The first integrand, $\sin(x)\cos(nx)$, is handled by applying the trigonometric identity $$\sin(a+b)+\sin(a-b)=2\sin(a)\cos(b).$$ The second integrand, $2(\pi-x)\cos(nx)$, requires integration by parts, using $u=2(\pi-x)$ and $dv=\cos(nx)dx$. The answers are $\dd a_0=\frac{4+\pi}{4\pi}$, $\dd a_1=\frac{4-\pi}{\pi^2}$, $\dd a_2=\frac{-2}{3}\frac{3+\pi}{\pi^2}$. } \end{exercises} \NOTES{ {\large\bf Maple notes on \LASTPROB}. The tedium of integral evaluation can be lessened by using {\tt maple}: \footnotesize\rm {\tt a0:=(1/Pi)*(int(sin(x),x=0..Pi/2)+ int((2/Pi)*(Pi-x),x=Pi/2..Pi));}\\ {\tt a1:=(2/Pi)*(int(sin(x)*cos(x),x=0..Pi/2)+int((2/Pi)*(Pi-x)*cos(x),x=Pi/2..Pi));} \\ {\tt a2:=(2/Pi)*(int(sin(x)*cos(2*x),x=0..Pi/2)+int((2/Pi)*(Pi-x)*cos(2*x),x=Pi/2..Pi));} \normalsize\rm } \begin{exercises} \probx{6a.}{(Half Range Expansions)} Graph $g(x)$, the even $2\pi$--periodic extension of $f_1(x)$ to the whole real line, and on the same axes, graph the 3-termed even half-range expansion of $f_1(x)$, $-4\pi \le x \le 4\pi$. \SOL{6a.}{ This is an interesting computer graphics problem and it can waste a lot of time. The plot to do first is the 3-termed even half--range expansion $$h(x)= \frac{4+\pi}{4\pi} + \frac{4-\pi}{\pi^2}\cos(x) + \frac{-2}{3}\frac{3+\pi}{\pi^2}\cos(2x). $$ Definitely use a plot program or a graphing calculator. After getting this plot on paper, plot the function $g(x)$ by hand on the same axes, starting with the even extension of the base function $f_1(x)$ to $[-\pi,\pi]$, which is the sine curve $|\sin x|$ on $|x|\le \pi/2$ plus two straight lines. Four copies of this curve are added to the graph, on the intervals $[n\pi,n\pi+2\pi]$ for $n=1,3,-3,-5$. Clip the graph to give the final window $[-4\pi,4\pi]$. The curve $g(x)$ looks like a piece of decorative molding: it's entirely in the upper half plane, no negative cycles! } \end{exercises} {\large\bf Fourier Integral}. Consider the following table where it is assumed that $a>0$: \begin{center} \def\ds{\displaystyle} \begin{tabular}{|c|c|c|} \hline \small\bf $\ds F(x)$ & \small\bf $\ds A(t)$ & \small\bf $\ds B(t)$ \\ \hline \hline &&\\[-1em] $ \left\{ \begin{array}{ll} \ds 1 & |x| < a \\ \ds 0 & \mbox{otherwise} \end{array}\right. $ & $\ds \frac{2\sin(at)}{\pi t}$ & $\ds 0$ \\ &&\\[-1em] \hline &&\\[-1em] $ \left\{ \begin{array}{ll} \ds e^{-ax} & x \ge 0 \\ \ds 0 & \mbox{otherwise} \end{array}\right. $ & $\ds \frac{1}{\pi}\frac{a}{t^2+a^2}$ & $\ds \frac{1}{\pi}\frac{t}{t^2+a^2}$ \\ &&\\[-1em] \hline &&\\[-1em] $\ds \frac{a}{x^2+a^2}$ & $ \ds \exp(-at) $ & $\ds 0$ \\[1em] \hline $ \left\{ \begin{array}{ll} \ds\sin(x) & |x| < \pi \\ \ds 0 & \mbox{otherwise} \end{array}\right. $ & $\ds 0$ & $\ds \frac{2}{\pi}\frac{\sin(\pi t)}{1-t^2}$ \\[1ex] \hline \end{tabular} \end{center} The methods used to verify the table are as follows. \begin{itemize} \item[] {\bf Method 1.} Insert the definition of $F$ into the corresponding {\em sine} and {\em cosine} integrals $$A(t)=\frac{1}{\pi}\int_{-\infty}^{\infty} F(u)\cos(tu)du,$$ $$B(t)=\frac{1}{\pi}\int_{-\infty}^{\infty} F(u)\sin(tu)du$$ and perform the indicated integration over $(-\infty,\infty)$. \item[]{\bf Method 2.} Begin with the formulas for $A(t)$ and $B(t)$ in table columns 2 and 3, insert them into the Fourier Integral Formula $$ \int_0^{\infty} \left[A(t)\cos(tx)+B(t)\sin(tx)\right] dt$$ and perform the indicated integration over $(0,\infty)$ to obtain the $F$ in table column 1. \end{itemize} \NOTES{ The hand computations implied by the two methods of verification can be reduced to first year calculus problems. An example: {\bf Example}. Find the Fourier sine and cosine integrals for the function $f(x)=1$, $0 \pi/2. \end{array}\right. $$ Calculate the integrals $$\displaystyle A(t)=\frac{1}{\pi}\int_{-\infty}^{\infty} F(u)\cos(tu)du,$$ $$B(t)=\frac{1}{\pi}\int_{-\infty}^{\infty} F(u)\sin(tu)du.$$ For which values of $x$ is the Fourier integral of $F(x)$ equal to $F(x)$? Cite a theorem by textbook and page number to justify your claim. \SOL{8.}{ Since the function $F$ is odd, $A(t)=0$. Then the integrand in $B$ is even and $$B(t)=\frac{2}{\pi}\int_0^{\pi/2} \sin(x)\sin(tx)dx =\frac{2t}{\pi}\frac{\cos(\pi t/2)}{1-t^2}.$$ %% (1/Pi)*int(sin(x)*sin(t*x),x=0..Pi/2);expand("); The Fourier convergence theorem implies that $F$ and the Fourier integral of $F$ agree except possibly at discontinuities of $F$. Indeed, the integral evaluates to half the jump at $x=\pm\pi/2$, which is $\pm 1/2$. } \end{exercises} \medskip \NOTES{ {\large\bf Maple Notes on Fourier Integrals}. The complex integral $$ A(t)-iB(t)= \int_{-\infty}^{\infty} f(u)\cos(tu)du -i \int_{-\infty}^{\infty} f(u)\sin(tu)du, $$ because of Euler's Formula $\exp(-i\theta)=\cos(\theta)-i\sin(\theta)$, equals the {\bf Complex Fourier Transform} $$T\equiv \int_{-\infty}^{\infty} f(u)\exp(-itu)du$$ Maple computes $T$ in the {\tt fourier} package, so there is a way to automate the calculation of Sine and Cosine Fourier integrals. These improper integrals are uneasy calculations at best, so it is handy to have machinery for error-free evaluation. {\bf Example}. Compute $A(t)$ and $B(t)$, the Fourier Cosine and Sine integrals, for the function $$F(x) = \left\{ \begin{array}{rl} \cos(x) & |x| < \pi \\ 0 & \mbox{otherwise} \end{array}\right. $$ {\bf Solution}: First write $F$ in terms of Heaviside's function $H(t)$ as $F(x) = \cos(x)(H(x+Pi)-H(x-Pi))$. To evaluate the improper integrals $A$ and $B$ we compute the complex Fourier transform $T$ and use the identity $A=\mbox{Real($T$)}$, $B=-\mbox{Imaginary($T$)}$. In maple, we use the package {\tt fourier} and the complex evaluation function {\tt evalc} to simplify answers. \footnotesize\rm \begin{verbatim} # Maple V.2 readlib(fourier); # Read fourier library F:= u -> cos(u)*(Heaviside(u+Pi)-Heaviside(u-Pi)); # Define function F p:=simplify(fourier(F(u),u,t)); # Take Fourier integral and simplify At:=evalc(Re(p)); # Cosine integral is real part Bt:=evalc(-Im(p)); # Sine integral is imaginary part \end{verbatim} \normalsize\rm The answers: $A(t)=2t\sin(\pi t)/(t^2-1)$, $B(t)=0$. {\large\bf Technical Issues}. The technical assumptions of the {\bf Fourier integral theorem} must be emphasized. The mechanical part of the problem has been settled by the methods developed above. The assumption is that $F$ is piecewise continuous on $(-\infty,\infty)$ and the right and left hand derivative of $F$ exist at every point $x$. In practice, $F$ is the sum of terms of the form $ g(x)(H(x-a)-H(x-b))$ where $H$ is Heaviside's function, $(a,b)$ is a finite or infinite interval and $g$ is {\bf differentiable} on $(a,b)$. In order for $F$ to satisfy the technical assumption in Fourier's theorem it suffices to verify that $g$ has a one-sided derivative at each finite endpoint $x=a$ and/or $x=b$. In engineering problems the technical assumption is usually satisfied. \medskip {\bf Example}. Does Fourier's integral theorem apply to $F(x)=\sqrt{x}H(x)$? {\bf Solution}: No. The function $g(x)=\sqrt{x}$ is continuously differentiable except at $x=0$. However, $g$ doesn't have a left-hand derivative at $x=0$ (in fact $g'(0+)=\infty$). The conclusion of Fourier's integral theorem is that the Fourier integral equals the average of the right hand and left hand limits of $F$ at the point $x$. Since $F$ is normally represented by smooth functions times Heaviside functions, the conclusion is suited exactly to engineering applications! In summary, the {\em only} time it is obvious that $F(x)$ equals its Fourier integral is in the limited case when $F$ is a smooth function with no Heaviside functions in its definition. \medskip {\bf Example}. Assume $F(x)=\cos(x)(H(x) - H(x-\pi))$. What is the value of the Fourier integral of $F$ at $x=0$? {\bf Solution}: The answer is {\bf not} found by plugging $x=0$ into the equation for $F$. That only succeeds when $F$ is continuous at $x=0$. This function isn't continuous at $x=0$. The value in the discontinuous case is the average of the right and left hand limits of $F$ at $x=0$. Answer: $\bf 1/2$! } %% End of NOTES {\large\bf Gibb's Phenomenon}. The function $$F(x) \equiv H(x+1)-H(x-1)= \left\{ \begin{array}{ll} 1 & -1 \le x < 1 \\[1ex] 0 & x<-1,~~ x\ge 1, \end{array}\right. $$ where $H$ is Heaviside's function, has Fourier cosine and sine integrals $A(t)=2\sin(t)/(\pi t)$ and $B(t)=0$ (see the table above). Let $$ \begin{array}{l} \dd F_N(x) = \int_0^N \frac{2\sin t \cos tx}{\pi t}dt \\[2ex] \hidem{\dd F_N(x)} \dd =\frac{1}{\pi}\left[ Si(N+Nx)-Si(Nx-N) \right]. \end{array} $$ \begin{exercises} \prob{9.}{(Gibb's Phenomenon)} Evaluate $G(x)$ and also the limit $\dd \lim_{N\to\infty} F_N(x)$ by considering the cases $x<-1$, $x=-1$, $-1 1$. As might be guessed from this form, the LHS will also be expressed as a piecewise function over the same five intervals. This calculation is based upon the definitions of right and left hand limits plus the definition $H(u)=1$ for $u\ge 0$, $H(u)=0$ for $u<0$. The answer: % F(x)=H(x+1)-H(x-1) $$G(x) = \frac12[F(x+0)+F(x-0)] = \left\{ \begin{array}{cl} 0 & x < -1, \\ \frac12 & x=-1, \\ 1 & -11. \end{array}\right. $$ The LHS is the integral $$ \begin{array}{l} \dd\int_0^\infty 2\frac{\sin t}{\pi t}\cos(tx) + (0)\sin(tx)dt= \\[2ex] \dd\frac{2}{\pi}\int_0^\infty \frac{\sin t}{t}\cos(tx)dt= \\[2ex] \dd\frac{1}{\pi} \int_0^\infty\frac{\sin(t+tx) + \sin(t-tx)}{t}dt=\\[2ex] \dd\frac{1}{\pi}\lim_{N \to \infty} \left[\int_0^{(1+x)N} \frac{\sin u}{u}du-\int_0^{(x-1)N} \frac{\sin w}{w}dw \right] \end{array} $$ The third integral uses the trig identity $\sin(t+tx) + \sin(t-tx) = 2\sin(t)\cos(tx)$, which is recalled from the general trigonometric identity $\sin(a+b) + \sin(a-b) = 2\sin(a)\cos(b)$. The last equality is the result of two changes of variable: $u=t+tx$, $du=(1+x)dt$ and $w=t-tx$, $dw=(1-x)dt$. To compute the limits, break up the analysis into the five cases aforementioned and use the relation (see the textbook) $$\int_0^\infty \frac{\sin t}{t}dt = \frac{\pi}{2}.$$ The limits of integration $x+xN$ and $-N+xN$ tend to values $0$, $-\infty$ and $\infty$, depending upon the signs of $(1+x)$ and $(1-x)$. For example, if $1+x<0$, then $1-x>0$ and $F_N(x)$ limits to $$\frac{1}{\pi}\int_0^{-\infty} \frac{\sin(u)}{u}du -\frac{1}{\pi}\int_0^\infty \frac{\sin(u)}{u}du = -\frac{-\pi}{2\pi} - \frac{\pi}{2\pi}=0.$$ The other four cases are treated similarly. The answer for the LHS: $$ \int_0^\infty 2\frac{\sin t}{\pi t}\cos(tx) + (0)\sin(tx)dt= \lim_{N\to\infty} F_N(x) = \left\{ \begin{array}{cl} 0 & x < -1, \\ \frac12 & x=-1, \\ 1 & -11. \end{array}\right. $$ } \probx{9a.}{(Graphical Approximation)} Make 4 separate graphs of $F_N(x)$ on $(-\infty,\infty)$, for $N=8$, $N=16$, $N=32$, $N=64$. Explain in a short paragraph the connection with {\bf Gibb's Phenomenon}. \SOL{9a.}{ Gibb's Phenomenon for {\em Fourier series} is illustrated by the square wave example. In this example there is a characteristic cycle in the graph of the $n^{\mbox{th}}$ partial sum $s_n(x)$ that is of fixed amplitude and moves towards the points of discontinuity of $F$ as $n\to\infty$. Elsewhere in the graph of $s_n(x)$ the amplitudes decay to zero as $n\to\infty$ because of the approximation $F(x)=s_n(x)+\frac{M}{2n+1}$ with $M$ of bounded magnitude. In short, at each point $x$ where $F$ has two derivatives, the approximations $s_n(x)$ converge to the function value $F(x)$. It is known that the convergence is uniform on each closed interval consisting entirely of points of continuity of $F$ and $F'$. The only bizarre convergence possible is at the isolated points of discontinuity of $F$. Gibb's Phenomenon for {\em Fourier integrals} is explained in a similar way. A paragraph of explanation might include graphical evidence of the phenomenon at the points $x = \pm 1$. The number of plot points can grossly affect the graphical evidence because the {\em blip} at $x=\pm 1$ can `fall through the cracks' if the step size $h$ is too large. Suggestion: use 200 points for $n=64$. Show in the graphs that the {\em blips} get thinner but the height of each of the four blips stays around $0.13$. } \end{exercises} \NOTES{ {\large\bf Maple notes}. Some remarks will be given about the initial conversion of the formula using {\tt maple} and the method for graphing on an infinite interval. The function $Si(x)=\int_0^x t^{-1}\sin(t)dt$ can be used to reduce the partial Fourier integral $$\int_0^N \frac{\sin t \cos tx}{t}dt$$ by the {\tt maple} integration command \footnotesize \begin{verbatim} int((1/t)*sin(t)*cos(t*x),t=0..N); # 1/2 Si(N + x N) - 1/2 Si(- N + x N) \end{verbatim} \normalsize To make graphs it suffices to code as the function $f(x)$ either the integral or the answer in terms of $Si$. An example for the plot at $N=32$: \footnotesize \begin{verbatim} f:= x -> int(sin(t)*cos(t*x)/t,t=0..n); n:=32: m:=100: a:=-1.2: b:=1.2: h:=(b-a)/m: L:=[a,f(a)]: x:=a: for i from 1 to m do x:=x+h: L:=L,[x,f(x)]: od: plot([L]); \end{verbatim} \normalsize The idea of this plot method is to compute a list of points for the plot (called $L$) and then plot the list of points. Generally this method is faster than a normal function-based plot and it has the advantage of saving the plot points in a list for later replotting or review. } {\large\bf Convergence of Fourier Series}. The $n^{\mbox{th}}$ Fourier Partial Sum is $$ s_n(x) \equiv \displaystyle \frac{a_0}{2} + \sum_{k=1}^n \left( a_k\cos(kx) + b_k\sin(kx)\right) .$$ Fourier's Convergence theorem says $$ \begin{array}{l} \dd\lim_{n\to\infty} s_n(x)= \frac{\lim_{t\to x+} F(t) + \lim_{t\to x-} F(t)}{2} \\[1ex] \hidem{\dd\lim_{n\to\infty} s_n(x)} \dd=F(x) \end{array} $$ at points $x$ of continuity of $F$. The purpose of the problems below is to establish the {\bf convergence result} $$F(x) = s_n(x) + \frac{M}{2n+1}$$ where $M$ is a certain integral (see the fourth problem below). The formula says that the graphs of $F$ and its Fourier approximation $s_n(x)$ are nearly identical for large values of $n$. The derivation here assumes $F$ twice continuously differentiable. The sliding fixed-magnitude oscillations of Gibb's Phenomenon cannot satisfy $F(x)=s_n(x) + \frac{M}{2n+1}$ and therefore estimation formulas like this are impossible for discontinuous $F$. {\em Your task} is to establish the results by assuming the computations of the following theorem. \medskip \fbox{ \begin{minipage}{3.3in} {\large\bf Theorem}. Let $F$ be piecewise continuous and $2\pi$--periodic. Assume $x$ given and define $$G(u)\equiv (F(x+u)-F(x))/\sin(u/2).$$ The partial sums $s_n$ of the Fourier series of $F$ satisfy the following identities: $$ \begin{array}{l} % Fourier series Integrals with Dirichlet kernels \dd s_n(x)= \frac{1}{2\pi} \int_{-\pi}^{\pi} F(x+u) \frac{\sin\left(\left(2n+1\right)u/2\right)} {\sin\left(u/2\right)}du \\[2ex] \dd s_n(x)-F(x) = \dd \frac{1}{2\pi}\int_{-\pi}^{\pi} G(u)\sin((2n+1)u/2)du. \end{array} $$ \end{minipage} } \begin{exercises} \prob{10.}{(Convergence of Fourier Series)} Establish the trigonometric identity $$ D_n(u) \equiv \frac{\sin\left(\left(2n+1\right)u/2\right)} {2\sin\left(u/2\right)} = \frac12 + \sum_{k=1}^n \cos(ku) $$ The fraction is called {\bf Dirichlet's Kernel} in the literature. \SOL{10.}{ To prove the trigonometric identity, prove instead $$ \frac{\sin\left(\left(2n+1\right)u/2\right)} {\sin\left(u/2\right)} = -1 + 2\sum_{k=0}^n \cos(ku) $$ This formula arises by adding and subtracting the term $\cos(0)$ in the summation, so as to include index $0$ in the sum. The identity will be proved by showing that the left hand side (LHS) equals the right hand side (RHS). {\bf Trig Method}. This method requires only the trigonometric identity $2\sin a \cos b = \sin(a+b) - \sin(a-b)$; it establishes the equivalent formula $$2\sin(u/2)[\frac12 + \sum_{k=1}^N \cos(ku)] = \sin((2N+1)u/2)$$ The left hand side (LHS) can be written as $$LHS=\sin(u/2) + \sum_{k=1}^N 2\sin(u/2)\cos(ku)$$ By the trig identity $$LHS=\sin(u/2) + \sum_{k=1}^N \sin((k+\frac12)u)-\sin((k-\frac12)u)$$ or simply $$LHS=\sin(u/2) + [\sin(3u/2)-\sin(u/2)] + [\sin(5u/2)-\sin(3u/2)]+\cdots$$ The terms all cancel except one, which is the RHS of the above identity (calulus books study this idea as {\em telescoping series}). } \prob{11.}{(Convergence of Fourier Series)} Integrate the formula of the previous problem on $[-\pi,\pi]$ to obtain $\int_{-\pi}^{\pi} D_n(u)du = \pi$. \SOL{11.}{ To prove the integration identity, integrate over $[-\pi,\pi]$ across the equation for Dirichlet's kernel: $$ D_n(u) = \frac12 + \sum_{k=1}^n \cos(ku). $$ The value of the integral $\int_{-\pi}^\pi \cos(ku)du$ should be familiar. } \prob{12.}{(Convergence of Fourier Series)} Assume $x$ fixed and $G(u)\equiv[F(x+u)-F(x)]/\sin(u/2)$ continuously differentiable on $[-\pi,\pi]$. \footnote{ If $F$ is continuously differentiable on $[-\pi,\pi]$ and $F(x+u)-F(x)=uL(u)$ where $L$ is continuously differentiable at $u=0$, then $G$ can be used in the computation. By application of L'Hospital's Rule for the $0/0$ case it follows that $G$ is continuous at $u=0$ with value $G(0)=2F'(x)$. } Establish using integration by parts: $$ \frac{1}{2\pi}\int_{-\pi}^{\pi} G(u)\sin((2n+1)u/2)du $$ $$= \frac{\frac{1}{\pi}\int_{-\pi}^{\pi} G'(u)\cos((2n+1)u/2)du}{2n+1} $$ \SOL{12.}{ To obtain the identity, begin by writing down the statement of integration by parts learned in calculus: $\int udv = uv - \int vdu$. You have to identify $u$, $dv$, $du$, $v$. Finally, the free term $uv$ will be evaluated at $-\pi$ and $\pi$. Justify why the free term is zero. } \prob{13.}{(Convergence of Fourier Series)} Conclude from the preceding problems that in the convergence formula $F(x)=s_n(x)+\displaystyle\frac{M}{2n+1}$ the value of $M$ is the integral $$M=\frac{-1}{\pi}\int_{-\pi}^{\pi} G'(u)\cos((2n+1)u/2)du,$$ which is bounded by $2\max |G'(u)|$ on $|u|\le \pi$. \SOL{13.}{ The methods are those of calculus. Mathematical background will affect how quickly you can digest the details. The written solution collects the results of the preceding problems as follows. Start with the second equation in the theorem. Isolate $F(x)$ in this formula and replace the integral by its equivalent fraction $\frac{M}{2n+1}$ where $M$ is $(-1)$ times the integral whose integrand contains $G'$. Estimates of the size of $M$ are obtained from the calculus inequalities $|\cos\theta| \le 1$ and $\left|\int_a^b f(x)dx\right| \le \int_a^b |f(x)|dx$. } \end{exercises} \NOTES{ {\large\bf Notes on the convergence theory}. \footnotesize\rm {\bf Theorem}. Let $F$ be piecewise continuous and $2\pi$--periodic. Assume $x$ given and define $$G(u)\equiv (F(x+u)-F(x))/\sin(u/2).$$ The partial sums $s_n$ of the Fourier series of $F$ satisfy the following identities: $$ \begin{array}{l} % Fourier series Integrals with Dirichlet kernels \dd s_n(x)= \frac{1}{2\pi} \int_{-\pi}^{\pi} F(x+u) \frac{\sin\left(\left(2n+1\right)u/2\right)} {\sin\left(u/2\right)}du \\[2ex] \dd s_n(x)-F(x) = \dd \frac{1}{2\pi}\int_{-\pi}^{\pi} G(u)\sin((2n+1)u/2)du. \end{array} $$ {\bf Proof of the Theorem}: The first equality results from the formula for $s_n(x)$ by substitution of the Euler formulas for $a_n$ and $b_n$, as follows: \WORKUPaaMatrix{ $\dd s_n(x) \equiv \displaystyle a_0 + \sum_{k=1}^n \left( a_k\cos(kx) + b_k\sin(kx)\right) $}{\rr Definition.} \WORKUPaaMatrix{ \hidem{\dd s_n(x)} $\dd =\frac{1}{2\pi}\int_{-\pi}^{\pi} f(x)dx + \sum_{k=1}^n \frac{1}{\pi}\int_{-\pi}^{\pi} f(t)\cos(kt)dt\cos(kx)$\newline \hidem{\dd s_n(x)\equiv} $\dd +\sum_{k=1}^n \frac{1}{\pi}\int_{-\pi}^{\pi} f(t)\sin(kt)dt\sin(kx) $}{\rr Euler's formulas.} \WORKUPaaMatrix{ \hidem{\dd s_n(x)} $\dd =\frac{1}{2\pi}\int_{-\pi}^{\pi} f(x)\left[1 + 2\sum_{k=1}^n \cos(kt)\cos(kx) + \sin(kt)\sin(kx)\right]dt $}{\rr Each term has a common integral sign. Factors $\cos(kx)$ and $\sin(kx)$ are independent of $t$. } \WORKUPaaMatrix{ \hidem{\dd s_n(x)} $\dd =\frac{1}{2\pi}\int_{-\pi}^{\pi} f(x)\frac{\sin(\frac{2k+1}{2}(t-x))}{\sin\frac{t-x}{2}}dt $}{\rr Apply the identity of the first problem below.} \WORKUPaaMatrix{ \hidem{\dd s_n(x)} $\dd =\frac{1}{2\pi}\int_{-\pi}^{\pi} f(x+u)\frac{\sin(\frac{2k+1}{2}u)}{\sin\frac{u}{2}}dt $}{\rr Use $u=t-x$, $du=dt$.} The above identity depends upon the formula $\int_a^b g(t)dt=\int_{-\pi}^{\pi} g(t)dt$, which is valid for any $2\pi$-periodic function $g(t)$ on an interval $[a,b]$ of length $2\pi$. This rather tricky fact is best proved in isolation, within the scope of integration theory and periodic functions, remotely distant from the theory of Fourier series. The last equality of the theorem is established by expansion of the right side into two integrals. The first is $s_n(x)$ by the first equality of the theorem. The second is $$ -\int_{-\pi}^{\pi} \frac{1}{\pi} F(x)D_n(u)du = F(x)\int_{-\pi}^{\pi} \frac{-1}{\pi} D_n(u)du, ~~ D_n(u)\equiv \frac{\sin\frac{(2n+1)u}{2}}{\sin\frac{u}{2}} $$ because $F(x)$ is constant relative to integration on $u$. By problem 11, the integral equals $F(x)$ times $-1$, which proves the theorem. \normalsize\rm } \NOTES{ {\large\bf Complex Method for the Dirichlet Kernel Formula}. Substitute for $\cos(ku)$ the expression $\Re(\exp(iku))$. Therefore, the sum is $\sum_{k=0}^n r^k$ where $r=\exp(iu)$. Now apply the geometric series formula $$\sum_{k=0}^n r^k = (r^{n+1} - 1)/(r-1)$$ with $r=\exp)iu)$. To simplify this complex fraction of exponentials multiply the top and bottom of the fraction by the conjugate of the bottom: $\exp(-iu)-1$. The fraction becomes $$ \frac{2\exp(inu)-\exp(-iu)-2\exp(i(n+1)u)+\exp(iu)}{|\exp(iu)-1|^2} $$ The first reduction is to write $|\exp(iu)-1|^2=2(1-\cos(u))$. Now take the real part of the expression on the right hand side to obtain: $$\frac{\cos[nu]-\cos[(n+1)u]}{1-\cos(u)}$$ To the top of the fraction apply the trigonometric identity $$\cos(a) - \cos(b) = 2\sin\frac{a+b}{2} \sin \frac{b-a}{2}$$ and to the bottom of the fraction apply the half-angle identity $$1-\cos(u)=2\sin^2 \frac{u}{2}$$ The conclusion is that LHS=RHS. } \begin{exercises} \prob{14.}{(Boundary Value Problems)} Find all the {\bf eigenvalues} $\lambda$, and the corresponding {\bf eigenfunctions} $x(t)$, for the boundary value problem $$ \begin{array}{l} \dd x''+\lambda x=0,~~ 0 \le t \le 1,\\[1ex] \dd x(0)=0,~~ x'(1)=0. \end{array} $$ \SOL{14.}{ The details will be explained below by examples 1 and 2. The first is a solution done entirely with pencil and paper. The second uses a computer algebra assist for certain portions of the computation. The present problem can be solved by analyzing recipe cases 1 and 2, real roots to the characteristic equation, eventually showing that the only solution in these two cases is $x(t)=0$. For recipe case 3, the complex case, a nonzero solution $x(t)$ must be of the form $x(t)=c_2\sin(\sqrt{\lambda} t)$ and therefore the eigenvalues are $\sqrt{\lambda}=(2n+1)\pi/2$ for $n=0,1,2,\ldots$ with corresponding eigenfunctions $\sin((2n+1)\pi t/2)$. } \probx{14a.}{(Boundary Value Problems)} Find all the {\bf eigenvalues} $\lambda$, and the corresponding {\bf eigenfunctions} $x(t)$, for the boundary value problem $$ \begin{array}{l} \dd x''+\lambda x=0,~~ -1 \le t \le 1,\\[1ex] \dd x(-1)=0,~~ x'(1)=0. \end{array} $$ \SOL{14a.}{ All that has changed from the previous problem is the interval. The eigenpairs are given by $\sqrt{\lambda_n}=(2n+1)\pi/2$, $\sin(\sqrt{\lambda_n}(t+1)/2)$. } \probx{14b.}{(Sturm-Liouville Theory)} Find all the eigenvalues and corresponding eigenfunctions of the Sturm-Liouville problem $$ \begin{array}{l} \dd y''-12y'+(36+4\lambda)y=0, ~~ 0 \le x \le 5,\\[1ex] \dd y(0)=y(5)=0. \end{array} $$ \SOL{14b.}{ The recipe contains the variable $\lambda$, so use $\Lambda$ instead to avoid confusion. Then $\Lambda^2-12\Lambda+36+4\lambda=0$ is the characteristic equation. The boundary conditions $y(0)=y(5)=0$ lead to the trivial solution $y=0$, if the roots $\Lambda$ are real and distinct or real and equal (Cases 1 and 2 of the {\em recipe}). Only in Case 3 of the recipe is there a chance of producing a nonzero solution $y$. In that case, $\Lambda$ is complex, obtained by solving $(\Lambda-6)^2 + 4\lambda=0$, so $\lambda>0$ and $\Lambda=6+2\sqrt{\lambda}i$. Then $y = c_1e^{-6x}\cos(2\sqrt{\lambda}x)+ c_2e^{-6x}\sin(2\sqrt{\lambda}x)$, where $c_1$ and $c_2$ are to be determined such that $y(0)=0=y(5)$. Considering these two equations leads to $c_1=0$. But $c_2\ne 0$, provided $\sin(10\sqrt{\lambda})=0$. This relations gives the eigenvalues $10\sqrt{\lambda_n}=n\pi$ ($n=1,2,3,\ldots$) and eigenfunctions $y_n=e^{-6x}\sin(n\pi x/5)$. } \prob{15.}{(Orthogonal Functions)} Show that set of functions below is {\bf orthogonal} on the given interval and determine the corresponding {\bf orthonormal} set. $$ \sin\left((n+\frac{1}{2})\pi t\right),~~ n=0,1,2,\cdots,~~ 0 \le t \le 1.$$ \SOL{15.}{ The orthogonality issue is $\int_0^1 f(t)g(t)dt=0$, where $f(t)=\sin(at)$, $g(t)=\sin(bt)$, $a=(2n+1)/2$, $b=(2m+1)/2$ and $n\ne m$. The problem can be done with calculus methods, by subtracting the trig identities $$\cos(A+B)=\cos A \cos B - \sin A \sin B,~~ \cos(A-B)=\cos A \cos B + \sin A \sin B.$$ To make an orthogonal set $\{\phi_n(t)\}$ on $[a,b]$ into an orthonormal set, divide $\phi_n$ by $\| \phi_n \| \equiv \sqrt{\int_a^b |\phi_n(t)|^2dt}$. In this case, $\phi_n(t)=\sin\left((n+\frac{1}{2})\pi t\right)$ and the division is by $\sqrt{c_n}$ where $c_n=\int_0^1 |\sin\left((n+\frac{1}{2})\pi t\right)|^2dt=1/2$. The orthonormal set is $$ \sqrt{2}\sin\left((n+\frac{1}{2})\pi t\right),~~ n=0,1,2,\cdots,~~ 0 \le t \le 1.$$ } \probx{15a.}{(Orthogonal Functions)} Show that set of functions below is {\bf orthogonal} on the given interval and determine the corresponding {\bf orthonormal} set. $$ \cos\left(3nt\right),~~ n=0,1,2,\cdots,~~ |t| \le \pi.$$ \SOL{15a.}{ The orthogonality question reduces to the integration of products $\cos(at)\cos(bt)$ over $[-\pi/2,\pi/2]$, where $a=3n$ and $b=3m$ ($a\ne b$). Calculus textbooks suggest using the trig identity $2\cos(at)\cos(bt)=\cos(at+bt)-\cos(at-bt)$ in order to facilitate the integration. The orthonormal set is the original set, each term divided by the constant $\sqrt{3\pi/6}$. } \end{exercises} \NOTES{ The study will apply the {\em recipe} for second order linear differential equations with constant coefficients. {\large\bf Recipe}. Given a constant coefficient differential equation $ax'' + bx' + cx = 0$ with characteristic equation $a\lambda^2+b\lambda+c=0$, the general solution is given as $x(t)=c_1x_1(t)+c_2x_2(t)$ where $c_1$ and $c_2$ are arbitrary constants and the {\em basis solutions} $x_1$ and $x_2$ are given as follows in terms of the roots $\lambda_1$, $\lambda_2$ of the characteristic equation. \begin{topics}{\bf Case 3} \item[{\bf Case 1}.] Distinct Real Roots $\lambda_1 \ne \lambda_2$. Then $$x_1(t) = e^{\lambda_1 t},~~ x_2(t)=e^{\lambda_2 t}.$$ \item[{\bf Case 2}.] Double Roots $\lambda_1=\lambda_2$. Then $$x_1(t) = e^{\lambda_1 t},~~ x_2(t)=te^{\lambda_1 t}.$$ \item[{\bf Case 3}.] Complex Conjugate Roots $\lambda_1=A+iB$, $\lambda_2=A-iB$. Then $$x_1(t) = e^{A t}\cos(Bt),~~ x_2(t)=e^{A t}\sin(Bt).$$ \end{topics} {\large\bf Example 1}. Solve $x'' + \lambda x=0$ subject to the boundary conditions $x(0)=0$, $x(1)=0$. Show that the eigenvalues are $(n\pi)^2$ with corresponding eigenfunction $\sin(n\pi t)$ for integers $n > 0$. \medskip {\large\bf Solution}: The discriminant of the characteristic equation $\Lambda^2 + \lambda=0$ is $D\equiv 0^2 - 4(1)(\lambda)$. The sign of the discriminant selects one of the three cases of the classic recipe. {\bf Case 1. ($\lambda > 0$)} In this case let $\lambda=k^2$ with $k>0$. Then the recipe implies that the general solution of the differential equation is $x=c_1\cos(kt)+c_2\sin(kt)$. The boundary conditions put these demands on the solution: %\footnotesize $$ \left\{ \begin{array}{l} 0=x(0) = c_1\cos(0)+c_2\sin(0)\\ 0=x(1)= c_1\cos(k)+c_2\sin(k) \end{array} \right. \mbox{~~or~~} \left\{ \begin{array}{l} 0= c_1(1) + c_2(0) \\ 0= c_1\cos(k)+c_2\sin(k) \end{array} \right. $$ %\normalsize The constants $c_1$ and $c_2$ cannot both be zero. A homogeneous system of linear algebraic equations has a nontrivial solution if and only if the determinant of the matrix of coefficients vanishes: $$ 0 = \left| \begin{array}{ll} 1 & 0 \\ \cos(k) & \sin(k) \end{array} \right| = \sin(k)$$ The solution of this equation is $k=n\pi$. The eigenfunction $x(t)=c_1\cos(kt)+c_2\sin(kt)$ is determined by solving for $c_1$ and $c_2$: $$ \left[ \begin{array}{ll} 1 & 0 \\ \cos(k) & \sin(k) \end{array} \right] % \left[ \begin{array}{l} c_1 \\ c_2 \end{array} \right] = \left[ \begin{array}{l} 0 \\ 0 \end{array} \right] $$ The answer is $c_1=0$ and $c_2$ arbitrary. In summary, $k=n\pi$ for integral $n \ge 1$ and the eigenfunction is $\sin(n\pi t)$. {\bf Case 2. ($\lambda=0$)} The recipe implies that the general solution is $x=c_1+c_2t$. The boundary conditions put these demands on the solution: $$ \left\{ \begin{array}{l} 0=x(0) = c_1+c_2(0) \\ 0=x(1)= c_1+c_2(1) \end{array} \right. $$ The matrix of coefficients has determinant $\displaystyle \left| \begin{array}{ll} 1 & 0 \\ 1 & 1 \end{array} \right| = 1$. By Cramer's Rule, $c_1=c_2=0$. In summary, $\lambda=0$ is {\em not an eigenvalue}. {\bf Case 3. ($\lambda < 0$)} In this case let $\lambda= -k^2$ with $k>0$. Then the recipe implies that the general solution of the differential equation is $x=c_1\exp(kt)+c_2\exp(-kt)$. The boundary conditions put these demands on the solution: %\footnotesize $$ \left\{ \begin{array}{l} 0=x(0) = c_1\exp(0)+c_2\exp(0) \\ 0=x(1)= c_1\exp(k)+c_2\exp(-k) \end{array} \right. \mbox{~~or~~} \left\{ \begin{array}{l} 0= c_1(1) + c_2(1) \\ 0= c_1e^k+c_2e^{-k} \end{array} \right. $$ %\normalsize The determinant of the matrix of coefficients is $\displaystyle \left| \begin{array}{ll} 1 & 1 \\ e^k & e^{-k} \end{array} \right| = e^{-k}-e^k$. The determinant is negative, because $k>0$ implies $e^{2k} > 1$ and therefore $e^{-k} - e^{k} = e^{-k}(1-e^{2k}) < 0$. By Cramer's Rule, $c_1=c_2=0$. In summary, {\em $\lambda$ is not an eigenvalue} for $\lambda < 0$. \bigskip {\large\bf Example 2}. Solve $x'' + \lambda x=0$ subject to the boundary conditions $x(0)=0$, $x(1)=0$ and show that the eigenvalues are $(n\pi)^2$ with corresponding eigenfunction $\sin(n\pi t)$ for integers $n > 0$. \medskip {\large\bf Maple--Assisted Solution}: The code below uses {\tt maple} functions {\tt diff}, {\tt dsolve}, {\tt rhs}, {\tt unapply}, {\tt subs}. It is helpful to invoke {\tt maple}'s `?' help function on each new maple operation (eg, `? rhs'); {\bf Case 1. ($\lambda>0$)} \footnotesize\rm\begin{verbatim} # Case lambda = k^2 > 0. a:=0: b:=1: # Define interval eq1:=diff(x(t),t,t)+k^2*x(t)=0: # Define differential equation ic:= x(0)=c, D(x)(0)=d: # Define initial conditions at t=0 dsolve({eq1,ic},x(t)): # Symbolic solve of DE rhs("): # Strip right hand side of answer x1:=unapply(",t); # Make it into a function of t x1(a)=0; x1(b)=0; # Print boundary conditions subs({x1(a)=0},x1(b)=0); # Substitute first boundary condition into the second \end{verbatim}\normalsize\rm If $x(t)$ is the eigenfunction, then $k^2$ is the eigenvalue. We don't assume $c$ or $d$ nonzero. They get determined by the boundary conditions and the requirement that $x(t)$ be nonzero. The boundary conditions reduce to $c=0$, $d\sin(k)=0$. For an eigenfunction we need $c$ or $d$ nonzero. Choose $d=1$. The final equation for $k$: $\sin(k)=0$. Don't expect {\tt maple} to solve this equation. Answer: $k=n\pi$, $n=1,2,\ldots$ ($n=0$ can't happen with $k>0$). In summary, the eigenvalues are $(n\pi)^2$, the eigenfunctions are $\sin(n\pi t)$, $n=1,2,3,\ldots$. {\bf Case 2. ($\lambda=0$)} The case $\lambda=0$ can be decided by solving $x''=0$ subject to $x(0)=0$, $x(1)=0$. It can be done mechanically as below. Ultimately, the calculation shows $x\equiv 0$. Conclusion: $\lambda=0$ is not an eigenvalue. \footnotesize\rm \begin{verbatim} # Case lambda = 0. a:=0: b:=1: # Define interval eq2:=diff(x(t),t,t)=0: # Define differential equation ic:= x(0)=c, D(x)(0)=d: # Define initial conditions at t=0 dsolve({eq2,ic},x(t)): # Symbolic solve of DE rhs("): # Strip right hand side of answer x2:=unapply(",t); # Make it into a function of t x2(a)=0; x2(b)=0; # Print boundary conditions subs({x2(a)=0},x2(b)=0); # Substitute first boundary condition into the second \end{verbatim} \normalsize\rm {\bf Case 3. ($\lambda<0$)} The case $\lambda<0$ can be decided by solving $x''-k^2x=0$ subject to $x(0)=0$, $x(1)=0$. It is done below. Ultimately, the calculation shows $x\equiv 0$. Conclusion: No negative $\lambda$ is an eigenvalue. Used below are some new {\tt maple} functions: the linear algebra package {\tt linalg}, the function {\tt genmatrix}, which extracts the matrix of coeffients from a system of linear equations, {\tt numer} which extracts the numerator of a quotient, {\tt det} which computes the determinant of a matrix. \footnotesize\rm \begin{verbatim} # Case lambda = -k^2 < 0. a:=0: b:=1: # Define interval eq3:=diff(x(t),t,t)-k^2*x(t)=0: # Define differential equation ic:= x(0)=c, D(x)(0)=d: # Define initial conditions at t=0 dsolve({eq3,ic},x(t)): # Symbolic solve of DE rhs("): # Strip right hand side of answer x3:=unapply(",t); # Make it into a function of t with(linalg): # Load functions genmatrix, det bc:=[x3(a)=0, x3(b)=0]: # Define boundary conditions var:=[c,d]: # Define variables C:=genmatrix(bc,var); # get matrix of coefficients simplify(numer(det(C)))=0; # determinant, numerator, simplify \end{verbatim} \normalsize\rm The result is $e^k - e^{-k}$ for the numerator of the determinant. This is nonzero for $k>0$, therefore Cramer's Rule implies $c=d=0$. No eigenfunction exists for any $\lambda < 0$. \bigskip {\bf Boundary Conditions with Derivatives}. In problems with boundary conditions involving derivatives, for example, $x'(-1)=0$, the coding in {\tt maple} language changes, e.g., {\tt D(x)(-1)=0}. {\bf Symbol Conflicts}. The upper case {\tt D} is a reserved symbol in {\tt maple}. In particular, your maple program cannot use {\tt D} as a variable name. Since mathematics texts like to use $D$ for the discriminant in a quadratic equation, there is potential conflict in naming conventions. }% End NOTES \NOTES{ {\large\bf Maple notes on orthogonality}. The orthogonality relation $\int_0^1 y_n(t)y_m(t)dt=0$, where $y_n(t)=\sin((2n+1)t/2)$, $n\ne m$, can be coded in {\tt maple} as follows. \footnotesize\rm \begin{verbatim} yn:=x -> sin((2*n+1)*Pi*x/2): ym:=x -> sin((2*m+1)*Pi*x/2): int(yn(x)*ym(x),x=0..1); \end{verbatim} \normalsize\rm The integration will give answers involving terms of the form $\sin(k\pi)$ where $k$ is an integer. Without further use of {\tt maple}, decide that the terms of the answer must be zero for $n \ne m$. {\large\bf Additional notes on orthogonality}. In problem 15a, the functions $\cos(2nt)$ on $[-\pi,\pi]$ are claimed to be {\bf orthogonal functions}. This means that $$\int_{-\pi}^{\pi} \cos(2pt)\cos(2qt)dt = 0$$ for integers $p$ and $q$ selected from the list $n=0,1,2,3,\ldots$, $p \ne q$. The question of validity of the formula rests with a suitable integral table. Use of a table is troubled by the conversion of the problem at hand to the variables in the table. After all table applications the limits of integration must be inserted. This process has two potential sources of error. The application of maple to the problem solves the trouble with conversion of variables: {\em no conversions are required}. That eliminates one source of error. Secondly, substitution of the limits of integration is done {\em automatically} in maple, eliminating the other source of error. Conclusion: maple is a good replacement for an integral table! {\bf Example}. Compute $\displaystyle \int_{-\pi}^{\pi} \cos(2px)\cos(2qx)dx$ for integers $p \ne q$. {\bf Solution}: Employ {\tt maple} as below. \footnotesize\rm \begin{verbatim} a:=-Pi: b:=Pi: # Define interval f:=x -> cos(2*p*x): g:=x ->cos(2*q*x): # Define functions f and g int(f(x)*g(x),x=a..b): # Symbolic exact integration simplify("); # Simplify output of previous calculation \end{verbatim} \normalsize\rm The output of maple: $$ - ~{\frac {\sin(-2\,p\pi +2\,q\pi )p+\sin(-2\,p\pi +2\,q\pi )q+\sin(-2\, p\pi -2\,q\pi )p-\sin(-2\,p\pi -2\,q\pi )q}{2\,p^{2}-2\,q^{2}}} $$ To finish the example it is required to show that the result is zero. This depends upon elementary properties of trigonometric functions. The reason maple does not simplify the integral answer to zero is that $p$ and $q$ are not assumed to be integers. However, you know that, so you can simplify it. } % End NOTES \begin{exercises} \prob{16.}{(Sturm-Liouville Theory)} Show that the Sturm-Liouville problem $$ \begin{array}{l} \dd (x^2 y')'+\lambda x^{-2}y=0,~~ 1/2 \le x \le 1,\\[1ex] \dd y(1/2)=y(1)=0, \end{array} $$ has eigenvalues and eigenfunctions $$ \lambda_n=n^2 \pi^2, ~~~ y_n=\sin(n \pi / x), ~~~ n=1,2,3, \cdots. $$ Then establish by {\em direct integration} the {\bf orthogonality property} $$ \int_{1/2}^1 y_n(x)y_m(x)x^{-2}dx = 0,~~~ n \ne m.$$ \SOL{16}{ Write out the differential equation in standard form as $x^2y'' + 2xy' + \lambda x^{-2}y=0$, then clear fractions to obtain the new equation $x^4y'' +2x^3y' + \lambda y=0$. This equation has variable coefficients and it cannot be solved directly from the {\em Recipe} for constant equations nor from the theory of Euler equations (it's not an Euler equation). The information that $\sin(n\pi/x)$ is a solution, given in the problem, suggests to make the change of variables $t=1/x$. If $u(t)=y(1/t)$, then the differential equation in $u$ and $t$ is $u'' + \lambda u=0$ and the new boundary conditions are $u(2)=u(1)=0$. The solution for this problem is $u(t)=c_1\cos(\sqrt{\lambda}t) + c_2\sin(\sqrt{\lambda}t)$ and the boundary condition $u(2)=0$ implies $\sin(\sqrt{\lambda})=0$, giving $\sqrt{\lambda}=n\pi$. Then $u(t)=\sin(n\pi t)$. Transforming back gives the eigenfunction $y(x)=\sin(n\pi/x)$ corresponding to eigenvalue $\lambda_n=(n\pi)^2$. To verify the orthogonality relation, treat the problem as a calculus integration problem $$\int_{1/2}^1 \sin(n\pi/x)\sin(m\pi/x)x^{-2}dx$$ for $n\ne m$. The secret to evaluation of the integral is the change of variables used in the differential equation itself: $t=1/x$. } \end{exercises} \NOTES{ {\large\bf Eigenvalues and Eigenfunctions of Sturm-Liouville Systems}. To find the eigenvalues involves solving symbolically a differential equation of second order and then determining when the solution represents an eigenfunction. The maple solution is discussed here. {\bf Example}. Solve $y'' -3y'+(18+2\lambda)y=0$ in {\tt maple} subject to the boundary conditions $y(0)=0$, $y(1)=0$. {\bf Solution}: First we decide that monotone solutions to the differential equation cannot satisfy the boundary conditions. Therefore, the second order equation must have solutions that oscillate, hence the discriminant is negative. The discriminant condition is $63+8\lambda > 0$. Define $k^2=63+8\lambda$ and rewrite the problem as $y'' -3y'+(18+k^2/4-63/4))y=0$, $y(0)=0$, $y(1)=0$. \footnotesize\rm \begin{verbatim} # Case k > 0. with(linalg): # Load functions genmatrix, det eq:=diff(y(t),t,t)+(-3)*diff(y(t),t)+ (18+(k^2)/4-63/4)*y(t)=0: # Define differential equation a:=0: b:=1: ic:= y(a)=c, D(y)(a)=d: # Define initial conditions dsolve({eq,ic},y(t)): # find general solution rhs("): # isolate right hand side y1:=unapply(",t); # create function from equation bc:=[y1(a)=0, y1(b)=0]: var:=[c,d]: # boundary conditions, variables C:=genmatrix(bc,var); # get matrix of coefficients simplify(numer(det(C)))=0; # determinant, numerator, simplify \end{verbatim} \normalsize\rm After $y$ is inserted into the boundary conditions there results two equations in variables $c$, $d$ with unknown $k$. The function {\tt genmatrix} determines the matrix of coefficients of this linear system in variables $c$ and $d$. The system has a nontrivial solution $(c,d)$ exactly when the determinant of $C$ is zero: $\det(C)=0$ is the equation for the unknown $k$. The determinant condition on $k$ is $\sin(k/2)=0$. Don't expect {\tt maple} to solve this equation --- by inspection $k/2$ is an integer multiple of $\pi$. Answer: $k=2n\pi$, $n= 1, 2, 3,\ldots$ subject to $k^2=63+8\lambda$. Therefore, $\lambda=((2n\pi)^2-63)/8$, $n=1,2,3,\ldots$. The eigenfunction for value $k$ is of the form $y=cu_1 + du_2$ where $$ u_1(t)= -{\frac {3\,{e^{{\frac {3\,t}{2}}}}\sin({\frac {kt}{2}})}{k}}+{ e^{{\frac {3\,t}{2}}}}\cos({\frac {kt}{2}}) ,~~~ u_2(t)= {\frac {2{e^{{ \frac {3\,t}{2}}}}\sin({\frac {kt}{2}})}{k}} $$ We don't assume $c$ or $d$ nonzero because $c$ and $d$ are determined by the boundary conditions. The requirement on $y(t)$ is that it be nonzero. To find $(c,d)$, substitute the value of $k$ into matrix $C$ and solve the homogeneous system of equations. As an example, we solve below for $c$ and $d$ when $k=2\pi$: \footnotesize\rm \begin{verbatim} A:=subs(k=2*Pi,eval(C)); # Substitute 2*Pi for k in matrix C kernel(A); # Solve AX=0 for X=[c,d] # Answer: [0,1] \end{verbatim} \normalsize\rm Taking $c=0$ and $d=1$ in $cu_1+du_2$ gives eigenfunction $u_2$ (see formula above) for eigenvalue $\lambda=(k^2-63)/8$ where $k=2\pi$. {\large\bf Additional notes on problem 16}. The equation $(x^2y')' + \lambda x^{-2}y=0$ is a {\bf Self-Adjoint Sturm-Liouville Problem}. More generally, an equation of the form $(py')' + (q+\lambda r) y=0$ is an eigenvalue problem of self-adjoint form. The equation itself is coded in {\tt maple} as follows: \footnotesize\rm \begin{verbatim} p:= x -> x^2: q:= x -> 0: r:= x -> 1/x^2: # Define p, q, r eq:=diff((p(t)*diff(y(t),t)),t)+(q(t)+(k^2)*r(t))*y(t)=0: # Define equation \end{verbatim} \normalsize\rm The preceding example can be used as a model to solve the equation and to find the determining equation for the eigenvalues. To solve the orthogonality question, integrate with {\tt maple}: \footnotesize\rm \begin{verbatim} yn:=x -> sin(n*Pi/x): ym:=x -> sin(m*Pi/x): w:= x -> 1/x^2: int(yn(x)*ym(x)*w(x),x=1/2..1); \end{verbatim} \normalsize\rm Well, {\tt maple} is not very efficient for deciding that the answer is zero, but it is easily decided by inspection of the answer. } % End NOTES \begin{exercises} \prob{17.}{(Sturm-Liouville Theory, Numerical)} Find accurate to 3 decimal places the first 5 eigenvalues and eigenfunctions of the Sturm-Liouville problem $$ \begin{array}{l} \dd y''+\lambda y=0, ~~ 0 \le x \le 1, \\[1ex] \dd y(0)+y'(0)=0,~~ y(1)=0. \end{array} $$ Plot the five eigenfunctions and count the zeros of each in the open interval $0 < x < 1$. \SOL{17}{ For $\lambda=0$ the eigenfunction is $1-x$. For positive eigenvalues, solve $\tan(k)=k$ numerically where $\lambda=k^2$. The first of the infinitely many positive solutions $k$ is about 4.49. The eigenfunctions for positive $\lambda=k^2$ are $k \cos(k x) - \sin(k x)$. There are no eigenvalues $\lambda<0$. To graph the eigenfunctions see the maple example below. } \end{exercises} \NOTES{ {\large\bf Example}. Solve for the first five eigenvalues and eigenfunctions of the Sturm-Liouville problem $y'' + \lambda y=0$, $y'(0)-y(0)=0$, $y(1)=0$, accurate to 3 decimal places. Plot the five eigenfunctions on $[0,1]$ {\bf Solution}: First we solve for the determining equations. The discriminant of the characteristic equation is $-4\lambda$ and its sign determines the character of the solution. We consider three cases: $\lambda=k^2$, $\lambda=0$, $\lambda= -k^2$ for $k>0$. \footnotesize\rm \begin{verbatim} # Case lambda > 0. with(linalg): eq:=diff(y(t),t,t)+k^2*y(t)=0: # Define equation a:=0: b:=1: ic:= y(a)=c, D(y)(a)=d: # Define initial conditions y1:=unapply(rhs(dsolve({eq,ic},y(t))),t); # find general solution bc:=[D(y1)(a)-y1(a)=0, y1(b)=0]: # boundary conditions var:=[c,d]: # variables C:=genmatrix(bc,var); # get matrix of coefficients simplify(numer(det(C)))=0; # determinant, numerator, simplify \end{verbatim} \normalsize\rm The result is the determining equation $k\cos(k)+\sin(k)=0$ for $k>0$. If $k>0$ is a root, then the eigenfunction is $k\cos(kx)+\sin(kx)$. The ideas used apply to problem 3, in which the answers are similar (see the answer given in problem 3). To check the validity of the eigenfunction, whether obtained by {\tt maple} or by hand, use these ideas: \footnotesize\rm \begin{verbatim} p:=x -> k*cos(k*x)+sin(k*x): D(p)(0) - p(0); # Left boundary condition p(1); # Right boundary condition diff(p(x),x,x)+k^2*p(x): # Differential equation simplify("); # Answers: 0, k*cos(k)+sin(k), 0 \end{verbatim} \normalsize\rm The cases $\lambda=0$ and $\lambda<0$ can be handled by repeating the above analysis using substitutions: \footnotesize\rm \begin{verbatim} # Case lambda = 0. with(linalg): eq:=diff(y(t),t,t)=0: # Define equation a:=0: b:=1: ic:= y(a)=c, D(y)(a)=d: # Define initial conditions y1:=unapply(rhs(dsolve({eq,ic},y(t))),t); # find general solution bc:=[D(y1)(a)-y1(a)=0, y1(b)=0]: # boundary conditions var:=[c,d]: # variables C:=genmatrix(bc,var); # get matrix of coefficients simplify(numer(det(C)))=0; # determinant, numerator, simplify \end{verbatim} \normalsize\rm The calculation results in the {\em inconsistent} equation $-2=0$, therefore {\em there is no eigenvalue for $\lambda=0$}. For $\lambda < 0$: \footnotesize\rm \begin{verbatim} # Case lambda < 0. with(linalg): eq:=diff(y(t),t,t)-k^2*y(t)=0: # Define equation a:=0: b:=1: ic:= y(a)=c, D(y)(a)=d: # Define initial conditions y1:=unapply(rhs(dsolve({eq,ic},y(t))),t); # find general solution bc:=[D(y1)(a)-y1(a)=0, y1(b)=0]: # boundary conditions var:=[c,d]: # variables C:=genmatrix(bc,var); # get matrix of coefficients simplify(numer(det(C)))=0; # determinant, numerator, simplify \end{verbatim} \normalsize\rm The determining equation for $k$ is $-k{e^{-k}}+{e^{-k}}-k{e^{k}}-{e^{k}}=0 $ or $e^{2k}=(1-k)/(1+k)$. Graphing $y=\exp(2x)$ and $y=(1-x)/(1+x)$ on the same axes over $(-\infty,\infty)$ gives evidence that $k=0$ is the only root. \footnotesize\rm \begin{verbatim} plot({exp(2*k),(1-k)/(1+k)},k=-infinity..infinity); # only root is k=0 \end{verbatim} \normalsize\rm Therefore, $\det(C) \ne 0$ for $k>0$ and {\em there are no negative eigenvalues} $\lambda = -k^2 < 0$. The numeric values of the first five eigenvalues are determined by the following {\tt maple} code. The process is complicated by graphics issues, so a separate but more transparent example appears below, which covers the main ideas used to solve this example. \footnotesize\rm \begin{verbatim} plot({tan(x),-x},x=0..6*Pi,y=-30..1,numpoints=100,xtickmarks=15); # Click on intersections of curves to get approximate answers for x # Maple displays (x,y) of click in upper left corner of plot fsolve(tan(x)+x=0,x=2..3); fsolve(tan(x)+x=0,x=4.9..5.2); # Solving is tricky. fsolve(tan(x)+x=0,x=7.9..8.1); # Sometimes choose an interval fsolve(tan(x)+x=0,x=11..11.2); # more than once to get an answer fsolve(tan(x)+x=0,x=14.2..14.3); # Answers: 2.03, 4.91, 7.98, 11.09, 14.21 \end{verbatim} \normalsize\rm The answers for the example: \begin{center} \begin{tabular}{|l|l|} \hline \bf Eigenvalue & \bf Eigenfunction \\ \hline $4.115858365$ & $2.028757838\cos(2.028757838 x) + \sin(2.028757838 x)$\\ $24.13934203$ & $4.913180439\cos(4.913180439 x) + \sin(4.913180439 x)$\\ $63.65910654$ & $7.978665712\cos(7.978665712 x) + \sin(7.978665712 x)$\\ $122.8891618$ & $11.08553841\cos(11.08553841 x) + \sin(11.08553841 x)$\\ $201.8512584$ & $14.20743673\cos(14.20743673 x) + \sin(14.20743673 x)$\\ \hline \end{tabular} \end{center} To plot the five eigenfunctions use this {\tt maple} code: \footnotesize\rm \begin{verbatim} g:= (x,k) -> k*cos(k*x)+sin(k*x): k:=fsolve(tan(x)+x=0,x=2..3); plot(g(x,k),x=0..1,title=`lambda=4.11`); k:=fsolve(tan(x)+x=0,x=4.9..5.2); plot(g(x,k),x=0..1,title=`lambda=24.2`); k:=fsolve(tan(x)+x=0,x=7.9..8.1); plot(g(x,k),x=0..1,title=`lambda=63.7`); k:=fsolve(tan(x)+x=0,x=11..11.2); plot(g(x,k),x=0..1,title=`lambda=122.9`); k:=fsolve(tan(x)+x=0,x=14.2..14.3); plot(g(x,k),x=0..1,title=`lambda=201.9`); \end{verbatim} \normalsize\rm The number of zeros of each eigenfunction is easily counted from the graph. In the count, zeros at the endpoints are not included. The answer: $0$, $1$, $2$, $3$, $4$, respectively. {\bf Sturm-Liouville Numerics}. The equation $\det(C)=0$ which determines the eigenvalues of the Sturm-Liouville problem is often impossible to solve symbolically. The method for handling this equation, and finding the roots, depends heavily on graphical methods to localize the solutions. {\bf Example}. Solve numerically the equation $\tan(x)=0.6x$ for the first five roots in $(0,\infty)$. {\bf Solution}: First, graph the equations $y=\tan(x)$ and $y=0.6x$ on the same axes over a suitable interval and solve for the roots, approximately. \footnotesize\rm \begin{verbatim} plot({tan(x),0.6*x},x=0..6*Pi,y=-1..30,numpoints=100,xtickmarks=15); # Click on intersections of curves to get approximate answers for x # Maple displays (x,y) of left-click in upper left corner of plot # Warning: Vertical lines in a maple graph may not be on a curve! \end{verbatim} \normalsize\rm From the graph the roots are approximately $4.5$, $7.5$, $11$, $14$, $17.2$. The supposed roots are bracketed with endpoints $a$ and $b$ and then we apply the {\tt maple} equation solver called {\tt fsolve}: \footnotesize\rm \begin{verbatim} a:=4.0: b:=5.0: fsolve(tan(x)-0.6*x=0,x=a..b); # 4.346208044 is the answer for guess 4.5 \end{verbatim} \normalsize\rm The process is repeated many times to obtain the five roots: $4.34$, $7.639$, $10.84$, $14.01$, $17.18$. }% End NOTES \begin{exercises} \prob{18.}{(Periodic Eigenvalue Problem)} Consider the eigenvalue problem $$ \begin{array}{l} \dd y'' + \lambda y =0,~~ -\pi \le x \le \pi,\\[1ex] \dd y(-\pi)=y(\pi),~~~ y'(-\pi)=y'(\pi). \end{array} $$ Verify the table below and show by direct integration that the `eigenfunctions' are pairwise orthogonal. \begin{center} {\rr \begin{tabular}{|l|p{0.9in}|p{1.25in}|} \hline \bf Values & \bf Eigenvalue & \bf Eigenfunction \\ \hline \hline $\lambda < 0$ & No eigenvalue & \parbox{1.2in}{\rr No eigenfunction } \\ \hline $\lambda=0$ & $\lambda=0$ & \parbox{1.2in}{\rr $y=1$ is the only eigenfunction} \\ \hline $\lambda>0$ & \rr $\lambda=n^2$, $n=1,2,3,\ldots$ & \parbox{1.2in}{\rr Two linearly independent eigenfunctions $y_1=\cos(nx)$, $y_2=\sin(nx)$}\\ \hline \end{tabular} } \end{center} \end{exercises} \SOL{18}{ This example from Sturm-Liouville theory is known as the {\bf Periodic Boundary Value Problem}. The difficulty is in the analysis of the boundary conditions $y(-\pi)=y(\pi)$, $y'(-\pi)=y'(\pi)$, which are known as {\bf periodic boundary conditions} because they are valid for {\em any} $2\pi$-periodic differentiable function $y$. If attacked entirely by hand calculation, the three cases of the {\em Recipe} are treated through the discriminant of the characteristic equation, which is $0^2-4\lambda$. Because solutions of the differential equation are monotone in Case 1, the only periodic solution in that case is $y=0$. In case 2, solutions look like $y=c_1+c_2x$, and since constants are periodic, this case produces eigenvalue $\lambda=0$ and eigenfunction $y=1$. In case 3, solutions are indeed periodic, $y=c_1\cos bx + c_2\sin bx$, where $\lambda=b^2$. The boundary conditions are satisfied provided $b$ is a positive integer. This determines the eigenvalues $\lambda=n^2$, $n=1,2,3,\ldots$. The eigenfunctions use {\em both} constants $c_1$ and $c_2$, with the restriction that not both are zero. A more careful analysis of case 3 can be done using linear algebra techniques, as in the maple notes below. It is expected that any complete solution use such an argument, modeled after previous solutions given above. Orthogonality is proved by direct evaluation of integrals of the form $$\int_{-\pi}^\pi f(x)gx)dx$$ where $f(x)$ and $g(x)$ are selected from the list of functions $1$, $\cos x$, $\sin x$, $\cos 2x$, $\sin 2x$, \ldots, with the requirement that $f(x)\ne g(x)$. There are four cases to consider, but all are handled by calculus methods. } \NOTES{ {\large\bf Maple Notes}. The {\tt maple} attack begins by simplifying the problem to exclude the discriminant cases $\lambda <0$ and $\lambda=0$. To be absolutely convincing the cases $\lambda <0$ and $\lambda=0$ should be treated separately by code similar to what appears below. \footnotesize\rm \begin{verbatim} with(linalg): eq:=diff(y(t),t,t)+k^2*y(t)=0: # Define equation a:=-Pi: b:=Pi: ic:= y(a)=c, D(y)(a)=d: # Define initial conditions y1:=unapply(rhs(dsolve({eq,ic},y(t))),t); # find general solution bc:=[y1(a)-y1(b)=0, D(y1)(a)-D(y1)(b)=0]: # boundary conditions var:=[c,d]: # variables C:=genmatrix(bc,var); # get matrix of coefficients simplify(numer(det(C)))=0; # determinant, numerator, simplify \end{verbatim} \normalsize\rm Unfortunately, {\tt Maple} confuses issues by printing the answer as $-4 +4\cos(k\pi)^2=0$ instead of the correct $-4+4\cos^2(k \pi)=0$. Don't expect {\tt maple} to solve this equation for $k$. Use trigonometry. The surprising result is in solving for the variables $c$ and $d$. For example, suppose $k=1$. The following {\tt maple} code organizes the solving for $c$ and $d$: \footnotesize\rm \begin{verbatim} A:=subs(k=1,eval(C)); # Evaluate C at first eigenvalue map(simplify,"); # Simplify output (it's the zero matrix) kernel("); # Solve AX=0 for X # What? two answers [1,0] and [0,1]! \end{verbatim} \normalsize\rm The {\tt kernel} command in {\tt maple} gives two answers, effectively meaning that $[c,d]=c_1[1,0]+c_2[0,1]$ for {\em arbitrary} $c_1$ and $c_2$! This shows $c$ and $d$ are arbitrary! However, for an eigenfunction, $c$ and $d$ cannot both be zero. The two combinations $c=1$, $d=0$ and $c=0$, $d=1$ given by {\tt maple}'s {\tt kernel} command generate two linearly independent eigenfunctions for each eigenvalue. To prove orthogonality, verify $\int_a^b \phi(x)\psi(x)dx=0$ where $a=-\pi$, $b=\pi$, $\phi=c_1\cos(nx)+c_2\sin(nx)$ and $\psi$ is either $1$ or else $d_1\cos(mx)+d_2\sin(mx)$, $n\ne m$. The organization of this computation in {\tt maple} uses the functions {\tt int} and {\tt simplify}. It is necessary to examine the final answer and justify that it is zero (without using {\tt maple}). } %% End NOTES {\large\bf Wave Equation and D'Alemberts Solution}. Let the initial wave $f$ be given by the formula $$f(u)=\left\{ \begin{array}{rl} \sin(\pi u) & ~~~ -1\le u\le 1, \\ 0 & ~~~ \mbox{otherwise}. \end{array}\right.$$ The variable $x$ represents position along the string, $t$ is time measured from 0 and $y$ is the string deflection. \begin{exercises} \prob{19.}{(D'Alemberts Solution and Traveling Waves)} Consider the traveling wave $y=f(x-ct)$. For each of $c=1$, $c=3$ make four 2D ``snapshots'' on $-2 \le x \le 2$. Start at $t=0$ and end at some value of $t$ where the wave is about $75$\% off the graph. \SOL{19}{ Due to D'Alembert's solution with velocity zero, the wave $f(x-ct)$ can be viewed as the rigid motion to the right of the original shape function $f(x)$, moving at speed $c$. The full solution is the average of this wave and its mirror image $f(x+ct)$, which moves to the left with speed $c$. Such waves are called {\em traveling waves}, because they move in time. This graphing problem is easiest by hand, because the graphs requested are copies of the original, moved to a new origin. In short, if you can draw $f$, then you can trace the others. See below for computer graphing techniques. } \end{exercises} \NOTES{ The answer in maple with two snapshots omitted: \footnotesize\rm \begin{verbatim} readlib(Heaviside): f:=u -> sin(Pi*u)*(Heaviside(u+1)-Heaviside(u-1)): a:=-2: b:=2: c:=1: z:=(x,t) -> f(x-c*t): plot3d(z,a..b,0..2,axes=BOXED,orientation=[-150,60]); plot(z(x,0.0),x=a..b); plot(z(x,2.0),x=a..b); \end{verbatim}\normalsize\rm The Heaviside function is used automatically in {\tt maple V.2} and therefore the {\tt readlib} command can be omitted. } \begin{exercises} \probx{19a.}{(D'Alemberts Solution and Snapshots)} The solution for an infinite taut string with $y(x,0)=f(x)$, $y_t(x,0)=0$ is $y=\frac{1}{2}\left[ f(x-ct)+f(x+ct) \right]$. For each of $c=1$, $c=5$ make four 2D ``snapshots'' on $-5 \le x \le 5$. Start at $t=0$. Choose the other values of $t$ to represent interesting snapshots of the solution. \SOL{19a}{ Like the previous problem, there is an interpretation of the component pieces in the solution $y(x,t)$ as traveling waves. However, the graph is now to reflect both the left and right traveling waves, ultimately showing the constructive and destructive superpositions of the waves. While it is certainly possible to do the graphing by hand, please use instead computing machinery. Most hand calculators cannot handle this problem due to the inability to enter piecewise-defined functions. See the maple solution below for ideas. } \end{exercises} \NOTES{ The answer in maple with two snapshots omitted: \footnotesize\rm \begin{verbatim} readlib(Heaviside): f:=u -> sin(Pi*u)*(Heaviside(u+1)-Heaviside(u-1)): a:=-5: b:=5: c:=3: y:=(x,t) -> (f(x-c*t)+f(x+c*t))/2: plot3d(y,a..b,0..2,axes=BOXED,orientation=[-150,60]); plot(y(x,0.0),x=a..b); plot(y(x,1.8),x=a..b); \end{verbatim}\normalsize\rm An analysis of the snapshots can be done with {\tt maple V.2} command {\tt animate}. See below for examples. } {\large\bf The Plucked String}. A plucked string lies along the interval $0 \le x \le 1$. It has wave speed $c=1$, initial velocity zero and initial shape $$f(x)=\left\{ \begin{array}{rl} 2x & ~~~ 0\le x\le 1/2, \\ 2(1-x) & ~~~ 1/2 \le x \le 1. \end{array}\right.$$ Notation: $x$ represents position along the string, $t$ is time measured from $0$ and $y$ is the string deflection. \begin{exercises} \prob{20.}{(Plucked String)} The odd periodic extension of $f$ of period $2$ defined on $(-\infty,\infty)$ will be called $h$. The exact solution $y$ can be written in terms of $h$: $$y = \frac{1}{2}\left[ h(x-t)+h(x+t)\right]$$ Graph the solution surface $y$ on $0\le x \le 1$, $0\le t \le 1$ and also the ``snapshots'' of $y$ at times $t=0$, $t=0.2$, $t=0.4$, $t=0.6$, $t=0.8$, $t=1$. \SOL{20}{ The initial shape can be written for a computer graphing program in terms of the Heaviside function $H(x)$ as $$f(x)= 2xH(x) + (2-4x)H(x-1/2).$$ The odd extension $g$ of this function to $[-1,1]$ is given by $$g(u) = \mbox{signum}(u)f(|u|)$$ where $\mbox{signum}(0)=u/|u|$ for $u\ne 0$ (it has values $\pm 1$). There are various ways to program piecewise defined functions and their even or odd extensions. The above is only one of the possibilities. If you are trying to use a hand calculator to do the graphing, then most of these ideas won't help. The periodic extension of $g$ to $(-\infty,\infty)$, of period $2$, can be defined by a separate recursive function $h(u)$. We recommend using the maple version described below or its equivalent. } \end{exercises} \NOTES{ Shown below is the {\tt maple} code for plotting the true solution. Four of the snapshots are omitted. \footnotesize\rm \begin{verbatim} read`/u/cl/maple/periodext`: h:=periodicextension: # Library for periodic extension readlib(Heaviside): H:=x -> Heaviside(evalf(x)): # Library for Heaviside unit step f:= x -> 2*x*H(x) + (2-4*x)*H(x-1/2): # Initial shape g:= u -> signum(u)*f(abs(u)): # odd extension of f to [-1,1] a:=0:b:=1:y:=(x,t) -> (h(x-t,g,-1,1)+h(x+t,g,-1,1))/2: # D'Alembert solution plot3d(y,a..b,0..2,axes=BOXED,orientation=[-150,60]); # 3D plot plot('y(x,0.0)',x=a..b); # snapshot #1 plot('y(x,1.0)',x=a..b); # snapshot #6 with(plots): animate('y(x,t)',x=a..b,t=0..1,frames=7); # animation \end{verbatim}\normalsize\rm Surprise: {\tt plot('y(x,0.5)', x=a..b)} makes a big mistake. The plot has $y$-range about $10^{-10}$. The plot is not useful. The single quotes (not backquotes!) are needed to keep {\tt maple} from evaluating $y(x,t)$. As can be seen from {\tt maple}'s own {\tt sum} function, quotes are sometimes necessary and unavoidable constructs. An error message like {\em cannot evaluate boolean} sometimes signals the need for single quotes around functional expressions. An alternative but equally unpleasant fix for the quotes problem is to code like this: {\tt animate (y, a..b, 0..1, frames=7)}, However, the latter leaves the quote problem unexplained. } \begin{exercises} \probx{20a.}{(Plucked String)} Separation of variables leads to the series representation of the solution $$y(x,t)= \sum_{n=1}^\infty B_n\sin(n\pi x)\cos(n\pi t).$$ Determine the first five terms of the series. Let $z$ be the partial sum of the first five terms. Graph the solution surface $z$ on $0\le x \le 1$, $0\le t \le 1$ and also the ``snapshots'' of $z$ at times $t=0$, $t=0.2$, $t=0.4$, $t=0.6$, $t=0.8$, $t=1$. \SOL{20a}{ To compute the approximate solution $z$ of five terms, use orthogonality of the functions $\sin(n\pi x)$ plus the identity $y(x,0)=f(x)$ to obtain $$ \int_0^1 f(x)\sin(n\pi x) dx = B_n\int_0^1 (\sin(n\pi x))^2dx $$ The problem is to evaluate the integrals, hence computing $B_n$, for $n=1$ to $n=5$. Then plot the truncated series $$z(x,t) = \sum_{n=1}^5 B_n\sin(n\pi x)\cos(n\pi t)$$ in a series of snapshots, for $t=0$ to $t=1$ in steps of $0.2$. The integration steps and the plots are best done in a computer algebra system like {\tt maple}. } \end{exercises} \NOTES{ In {\tt maple} evaluate the integrals and solve for $B_n$: \footnotesize\rm \begin{verbatim} readlib(Heaviside): # Read in library function Heaviside H:=x -> Heaviside(evalf(x)): # Abbreviation f:= x -> 2*x*H(x) + (2-4*x)*H(x-1/2): # Initial shape p:=int(sin(n*Pi*x)*2*x*H(x),x=0..1) + # Have to break up Heaviside int(sin(n*Pi*x)*(2-4*x)*H(x-1/2),x=0..1); # integrands q:=int(sin(n*Pi*x)*sin(n*Pi*x),x=0..1); # But ordinary integrals Ok r:=simplify(p/q); # Simplify trig B:=m -> subs(n=m,eval(r)); # Define coefficient sequence for series \end{verbatim}\normalsize\rm The partial sum $z$ is obtained and graphed as follows: \footnotesize\rm \begin{verbatim} a:=0: b:=1: z:= (x,t) -> sum('B(n)*sin(n*Pi*x)*cos(n*Pi*t)','n'=1..5); plot3d(z,a..b,0..2,axes=BOXED,orientation=[-150,60]); # 3D plot plot(z(x,0.0),x=a..b); # snapshot #1 plot(z(x,1.0),x=a..b); # snapshot #6 with(plots): animate(z(x,t),x=a..b,t=0..1,frames=7); # animation \end{verbatim}\normalsize\rm The reason the quotes are absent in the plots (e.g., \verb|'z(x,0.0)'|) was explained earlier. For the explanation of quotes in the {\tt sum} expression see {\tt maple} help on {\tt sum}. } \begin{exercises} \probx{20b.}{(Sturm-Liouville Theory: Hermite Equation)} Let $H_n$ be the Hermite Polynomial solution of Hermite's equation $(e^{-x^2}y')' + 2ne^{-x^2}y=0$. Prove that the set of functions $\{H_n(x)\}$ is orthogonal with respect to weight function $e^{-x^2}$ on the interval $(-\infty,\infty)$. \SOL{20b}{ Mimic the proof for a finite interval. First, write down the two differential equations, one for $H_n$ and one for $H_m$ ($m\ne n$). Then multiply the first by $H_m$ and the second by $H_n$. Subtract the equations and integrate over $(-\infty,\infty)$. Finally, integrate by parts to obtain a formula for the integral of $(n-m)H_nH_me^{-x^2}$. Only boundary terms remain from the integration by parts, and these terms are readily seen to be zero. } \probx{20c.}{(Sturm-Liouville Theory: Hermite Equation)} Verify by direct integration, using formulas for the Hermite polynomials $H_{10}$ and $H_7$, the {\bf orthogonality} relation $$\int_{-\infty}^{\infty} H_{10}(x)H_7(x)e^{-x^2}dx=0.$$ \SOL{20c}{ Use integration by parts or apply maple's {\tt int} function. } \end{exercises} {\large\bf The Method of Separation of Variables}. Apply in each of the problems below the method of separation of variables to find the solution of the partial differential equation. \begin{exercises} \prob{21.}{(Separation of Variables and a Vibrating String)} Find the displacement $y(x,t)$ of the vibrating string problem on domain $t \ge 0$, $0 \le x \le 5\pi$: % (see Derrick-Grossman page 758, equation (12)) $$ \begin{array}{l} \dd \frac{\partial^2 y}{\partial t^2} = 4\frac{\partial^2 y}{\partial x^2},\\[2ex] \dd y(0,t)=y(5\pi,t)=0,~~ y_t(x,0)=0,\\[1ex] \dd y(x,0)=f(x)=11\sin 2x -4\sin 5x. \end{array} $$ Verify that your solution is correct. Graph five snapshots of the solution for values $t=0$ to $t=0.5$. \SOL{21}{ The eigenvalue problem is $X'' + \lambda X=0$, $X(0)=X(L)=0$ where $L=5\pi$ with eigenvalues $\lambda_n = (n\pi/L)^2$ and eigenfunctions $\sin(n\pi x/L)$. The associated problem in $t$ is $T'' + c^2\lambda T=0$ with general solution $c_1\cos(cn\pi t/L) + c_2\sin(cn\pi t/L)$. By the boundary condition $T'(0)=0$ it follows $c_2=0$. Therefore the solution to the problem is $$ \sum_{n=1}^\infty B_n \cos(n\pi ct/L)\sin(n\pi x/L)$$ where $B_n=(2/L)\int_0^L f(x)\sin(n\pi x/L)dx$. For the special case $f(x)=11\sin 2x -4\sin 5x$ we have $B_n= 0$ except for $n=10$ and $n=25$. Therefore, $y(x,t)=11\cos(4t)\sin(2x)-4\cos(10t)\sin(5x)$. It is important to observe that all graphs and animation can be obtained from the formula $y(x,t)=11\cos(4t)\sin(2x)-4\cos(10t)\sin(5x)$. The remarks below discuss a machine method for arriving at the formula for $y(x,t)$. The usual plotting methods apply to produce snapshots of $y(x,t)$. For example, $t=0.2$ appears below. Also included is animation code. } \end{exercises} \NOTES{ Plotting in {\tt maple}: \begin{verbatim} plot(y(x,0.2),x=0..L,title=`t=0.2`); with(plots): animate(y(x,t),x=0..L,t=0..Pi/2,frames=50); \end{verbatim} The {\tt maple} solution for $y(x,t)$ amounts to evaluation of integrals and sums. The syntax is tricky so study the example and consult the {\tt maple} help file on {\tt sum} to explain why the single quotes appear. \footnotesize\rm \begin{verbatim} f:=x -> 11*sin (2*x) -4*sin (5*x); L:=5*Pi: c:=2: g:=(x,m) -> (2/L)*sin(m*Pi*x/L)*f(x); B:=n -> int(g(x,n),x=0..L); y:= (x,t) -> sum('B(k)*cos(k*Pi*c*t/L)*sin(k*Pi*x/L)', 'k'=1..30 ); \end{verbatim}\normalsize\rm To verify the solution in the case of the special $f$: \footnotesize\rm \begin{verbatim} y(0,t); y(L,t); simplify(subs(t=0,diff(y(x,t),t))); y(x,0); diff(y(x,t),t,t)-4*diff(y(x,t),x,x); \end{verbatim}\normalsize\rm The answers should be $0$, $0$, $0$, $11\sin 2x -4\sin 5x$, $0$. To decide on the snapshots, do an animation as follows and selectively choose the values of $t$. \footnotesize\rm \begin{verbatim} with(plots): animate(y(x,t),x=0..L,t=0..Pi/2,frames=50); \end{verbatim}\normalsize\rm } \begin{exercises} \prob{22.}{(Separation of Variables and the Heat Equation)} Find the temperature $T(x,t)$ in an insulated rod whose ends are kept at $0$ degrees. The domain is $t\ge 0$ and $0 \le x \le 3\pi $. The equation is % (see Derrick-Grossman page 763, equation (12)) $$\frac{\partial T}{\partial t} = 2\frac{\partial^2 T}{\partial x^2},$$ and the boundary conditions are $$ \begin{array}{l} \dd T(x,0)=f(x) = 10\sin 2x-3\sin 5x,\\[1ex] \dd T(0,t)=T(3\pi,t)=0. \end{array} $$ After obtaining an explicit formula for $T$, graph $5$ snapshots on $0 \le x \le 3\pi$ for $t=0$, $t=0.1$, $t=0.2$, $t=0.3$, $t=0.4$. \SOL{22}{ %% %% 10 sin(2 x) exp(- 8 t) - 3 sin(5 x) exp(- 50 t) %% The solution of the differential equation with boundary conditions by the method of separation of variables leads to a Sturm-Liouville boundary value problem $X''+k^2X=0$, $X(0)=X(L)=0$ where $L=3\pi$. The problem has eigenvalues $(n\pi/L)^2$ with eigenfunctions $B_n\sin(n\pi x/L)$. The problem for $t$ results in factor $\exp(-(n\pi/L)^2(2t))$. By analyzing the terms of the series for $f$ in the eigenfunctions it is found that $B_6$ and $B_{15}$ contribute but the others are zero. Finally, the answer is $T=10\,\sin(2\,x){e^{-8\,t}}-3\,\ \sin(5\,x){e^{-50\,t}}$. This answer is checked in {\tt maple} by the command \begin{quote}\tt T:=(x,t) -> 10*sin(2*x)*exp(-36*(2*t)/9) \\ -3*sin(5*x)*exp(-225*(2*t)/9); \\ diff(T(x,t),t) - 2*diff(T(x,t),x,x); \end{quote} To do the snapshots it is useful to invoke the animation panel in {\tt maple} using these commands: \begin{quote}\tt with(plots): \\ animate( T(x,t),x=0..3*Pi,t=0..0.5,numpoints=100,frames=50); \end{quote} The control panels allows selection of the appropriate snapshot frames for printing. With some experiments the frame count can be adjusted so that the exact values of $t$ are known for the snapshots. } \end{exercises} \begin{exercises} \probx{22a.}{(Separation of Variables and the Heat Equation)} Find a formula for the temperature $u(x,t)$ in a bar of length $1$ satisfies the conditions $$ \begin{array}{l} \dd u_t=u_{xx} \quad \text{for} \quad 0 \le x \le 1,\ \ t \ge 0, \\[1ex] \dd u(0,t)=0 \quad \text{and} \quad u_x(1,t)=0 \quad \text{for} \quad t \ge 0, \\[1ex] \dd u(x,0)=9\sin(3\pi x / 2)-4\sin(11\pi x / 2) \end{array} $$ for $0 \le x \le 1$. Plot $5$ snapshots of the solution on $0 \le x \le 1$ for $t$ in the domain $t=0$ to $t=0.08$ that show interesting behavior of the solution. %The relevant %Sturm-Liouville problem is worked in the text as Example 11.6.1 (pp. % Derrick-Grossman 703-704). \SOL{22a}{ This is a classical textbook application of the method of separation of variables. The method leads to the Sturm-Liouville problem $X'' + k^2X=0$, $X(0)=0$, $X'(1)=0$. The eigenvalues are $((n+1/2)\pi)^2$ and the eigenfunctions are $\sin((n+1/2)\pi x)$. The terms of the series solution $\sum_{n=1}^\infty B_n\sin(\lambda_n x)\exp(-\lambda_n t)$ that contribute are $n=2$ and $n=6$. The answer is $u=9\sin(3\pi x / 2)\exp(-(3\pi/2)^2t)-4\sin(11\pi x / 2)\exp(-(11\pi/2)^2t)$. To check the boundary conditions: \begin{quote}\tt u(0,t); \\ u(x,0); \\ subs(x=1,diff(u(x,t),x)); \end{quote} The answers are supposed to be $0$, $9\,\sin(\frac {3\pi x}{2})-4\,\sin({\frac {11\,\pi \,x}{2}})$, $0$. To do the snapshots it is useful to invoke the animation panel in {\tt maple} using these commands: \begin{quote}\tt with(plots): \\ animate( u(x,t),x=0..1,t=0..0.05,numpoints=100,frames=12); \end{quote} The control panel allow selection of the appropriate snapshot frames for printing. The exact values of $t$ should be reported. } \end{exercises} {\large\bf Laplace Equation in a Rectangle}. The steady state temperature $T=T(x,y)$ in a rectangular plate satisfies the Laplace equation $$ T_{xx}+T_{yy}=0, ~~~ 0 \le x \le 2, ~~ 0 \le y \le 1. $$ \begin{exercises} \prob{23.}{(Laplace Equation in a Rectangle)} Assume boundary conditions $T(x,0)=T(x,1)=0$, $0 \le x \le 2$, $T(0,y)=f(y)$, $T(2,y)=0$, $0 \le y \le 1$. Apply the method of separation of variables to construct $T=T(x,y)$. Find the exact solution for the case that $f(y)=4\sin(6 \pi y)-3\sin(8 \pi y)$. Plot the exact solution $T$ over the rectangular plate and mark the surface points where the temperature is locally at a maximum. %Answer: %T:=(x,y) -> (4/sinh(12*Pi))*sinh(6*Pi*(2-x))*sin(6*Pi*y) - % (3/sinh(16*Pi))*sinh(8*Pi*(2-x))*sin(8*Pi*y); % From Wilcox 353 final exam % \SOL{23}{ The answer to the problem is $$ T= {\frac {4\,\sinh(6\,\pi \,\left (2-x\right ))\sin(6\, \pi \,y)}{\sinh(12\,\pi )}}-{\frac {3\,\sinh(8\ \,\pi \,\left (2-x \right ))\sin(8\,\pi \,y)}{\sinh(16\,\pi )}} $$ } \end{exercises} \NOTES{ To check the answer: \footnotesize\rm \begin{verbatim} T:=(x,y) -> (4/sinh(12*Pi))*sinh(6*Pi*(2-x))*sin(6*Pi*y) - (3/sinh(16*Pi))*sinh(8*Pi*(2-x))*sin(8*Pi*y); T(x,0); T(x,1); T(2,y); T(0,y); \end{verbatim}\normalsize\rm The answers should be $0$, $0$, $0$, $4\,\sin(6\,\pi \,y)-3\,\sin(8\,\pi \,y)$. The plot in 3D is done with the following command: \footnotesize\rm \begin{verbatim} plot3d(T,0..2,0..1,axes=BOXED,style=PATCH,orientation=[-150,60]); \end{verbatim}\normalsize\rm From this plot, analyze the surface for points of local maximum temperature. To adjust the perspective of the plot hold down \MOUSELEFT{} and drag the wire frame into the new viewing position, then click \MOUSEMIDDLE{} to redisplay the plot. A mathematical analysis of extrema can be done using calculus. Either the extrema are on the boundary or else interior to the rectangle at a point $(x,y)$ where the gradient vanishes. By the graph, no tangent plane is horizontal, therefore the gradient cannot vanish interior to the rectangle. Conclusion: the extrema are on the boundary of the rectangle. To evaluate the extrema it is enough to consider the curves $T(x,0)$, $T(x,1)$, $T(0,y)$, $T(2,y)$ and find their extrema. Four plots will help with the logic. To actually find the extrema requires calculus: the extrema are at the endpoints or at a point where the derivative vanishes. Approximate answers for local maxima are: $(0.14,0,3.1)$, $(0.4,0,6.8)$, $(0.7,0,5.5)$. Find the exact answers to several digits. For example, to find the exact answer for $0.7$: \footnotesize\rm \begin{verbatim} p:=fsolve(diff(T(0,y),y)=0,y=.6..0.8): x1:=evalf(p,8): y1:=0: z1:=evalf(T(0,p),8): print(x1,y1,z1); \end{verbatim}\normalsize\rm \begin{quote} \footnotesize\rm The {\em method of contour plots} uses a plane graph consisting of curves of equal temperature. The same method is used in geological survey maps to engrave curves of constant altitude and in weather maps to show lines of contant barometric pressure. A useful contour plot can be made by selecting {\tt contour} from the {\tt Style} menu on the upper portion of the 3D plot window. This information is useful for shrinking the plot domain in order to find extrema. In this example, $0 \le x \le 0.1$ and $0 \le y \le 1$ is the region of interest. In {\tt maple} there is the function {\tt contourplot} in the {\tt plots} package that is made for production of contour plots. It is somewhat different from the one described above. An example: \begin{verbatim} T:=(x,y) -> (4/sinh(12*Pi))*sinh(6*Pi*(2-x))*sin(6*Pi*y) - (3/sinh(16*Pi))*sinh(8*Pi*(2-x))*sin(8*Pi*y); with(plots): contourplot(T,0..0.1,0..1,grid=[49,49],axes=BOXED,tickmarks=[5,20,1]); \end{verbatim} It is not claimed that extrema can be determined solely from a contour plot. In {\tt maple}'s contour plots the expected labels for the temperature $T$ along a contour are missing. However, with knowledge of the location of the extrema, the contour map can determine the extrema up to the fineness of the tick marks. \end{quote} } \begin{exercises} \probx{23a.}{(Laplace Equation in a Rectangle)} Assume boundary conditions $T(x,0)=0$, $T(x,1)=g(x)$, $0 \le x \le 1$, $T(0,y)=T(1,y)=0$, $0 \le y \le 1$. Apply the method of separation of variables to construct $T=T(x,y)$. Find the exact solution as an infinite series for the case that $g(x)=x(1-x)$. Let $w$ be the approximation obtained by the first five terms of the infinite series. Plot $w$ over the rectangular plate and mark the surface point of maximum temperature. % Page 768, problem 2, with numbers inserted. \SOL{23a}{ Variables separate and produce two problems. The first is $u'' + k^2u=0$, $u(0)=u(1)=0$ with eigenvalues $(n\pi)^2$ and eigenfunctions $\sin(n\pi x)$. The second problem $v'' - k^2v=0$, $v(0)=0$ has solution $\sinh(n\pi y)$. The temperature is $$T=\sum_{n=1}^\infty B_n\sin(n\pi x)\sinh(n\pi y).$$ To determine $B_n$ requires a quick orthogonality calculation giving $$B_n=(4-4\cos(n\pi))/(n^3\pi^3\sinh(n\pi)).$$ Let $w$ be the series truncated to five terms. As a check on computations evaluate $w(0,y)$, $w(1,y)$, $w(x,0)$, $w(x,1)$. The first 3 should be $0$ and the last a function that approximates $x(1-x)$. The graphs of $w(x,1)$ and $x(1-x)$ should be nearly identical. All plots are done in {\tt maple}; see below. } \end{exercises} \NOTES{ The implementation in {\tt maple}: \footnotesize\rm \begin{verbatim} q:=int((sin(n*Pi*x))^2,x=0..1)*sinh(n*Pi); # ==sinh(n*Pi)/2 p:=int(x*(1-x)*sin(n*Pi*x),x=0..1); # p=B(n)*q B:= m -> subs(n=m,eval(simplify(p/q))); w:=(x,y) -> sum(B(n)*sin(n*Pi*x)*sinh(n*Pi*y),n=1..5); w(0,y); w(1,y); w(x,0); w(x,1); plot(w(x,1),x=0..1); plot(x*(1-x),x=0..1); plot3d(w,0..1,0..1,axes=BOXED,style=PATCH); \end{verbatim}\normalsize\rm The 3D plot takes more than 30 seconds. } {\large\bf Heat Diffusion with Newton's Law of Cooling}. The method of separation of variables applies to determine the temperature $T(x,t)$ in an insulated rod, given the left end $x=0$ radiates heat into a medium at zero temperature ({\bf Newton's Law of Cooling}) while the the right end $x=1$ is held at zero temperature. The model for $T(x,t)$ is $$ \begin{array}{l} \dd \frac{\partial T}{\partial t} - \frac{\partial^2 T}{\partial x^2} =0,\\[2ex] \dd T_x(0,t) - T(0,t) = 0, \\[1ex] \dd T(1,t)=0, \\[1ex] \dd T(x,0)=f(x),~~ 0 \le x \le 1,~~ t\ge 0. \end{array} $$ \begin{exercises} \prob{24.}{(Newton's Cooling Law and Diffusion)} Solve for $T(x,t)$ explicitly for the special case $$f(x) = 1 + x - 2x^2$$ using just the first three terms of the series. Plot the first three eigenfunctions on $[0,1]$, the initial data $T(x,0)$ given by the first three terms, the actual initial data $f(x)$ and a series of five snapshots of $T(x,t)$ for $t=0$ to $t=0.8$. \SOL{24}{ The problem will be discussed for initial data $f(x) = 1 + \sin(x) - x^2(1+\sin(1))$. Much of what is said below applies to an arbitrary function $f$. Computations and graphs require the formula for $f$ to be used explicitly. The eigenvalue problem is $y'' + \lambda y=0$, $y'(0)-y(0)=0$, $y(1)=0$ with eigenvalues $\lambda=k^2 > 0$ and $k$ is the solution of $\tan(k)+k=0$ for $k>0$ (asymptotic to $(2n-1)\pi/2$ as $n\to\infty$). The eigenfunctions are $k\cos(kx)+\sin(kx)$. The problem in $t$ has solution $\exp(-k^2t)$. The eigenvalue analysis can be carried out in {\tt maple}. Presented here is the case for $\lambda > 0$; a complete solution would consider also $\lambda=0$ and $\lambda < 0$. } \end{exercises} \NOTES{ \footnotesize\rm \begin{verbatim} # Case lambda > 0. with(linalg): eq:=diff(y(t),t,t)+(k^2)*y(t)=0: # Define differential equation a:=0: b:=1: ic:= y(a)=c, D(y)(a)=d: # Define initial conditions y1:=unapply(rhs(dsolve({eq,ic},y(t))),t); # find general solution bc:=[D(y1)(a)-y1(a)=0, y1(b)=0]: # Define boundary conditions var:=[c,d]: # Define variables C:=genmatrix(bc,var); # get matrix of coefficients simplify(numer(det(C)))=0; # determinant, numerator, simplify \end{verbatim} \normalsize\rm To solve for the first three eigenvalues graph $\tan(x)$ and $-x$ and record where they cross, approximately: $x=2$, $x=5$, $x=8$. A {\tt maple} command which gives a rough graph is \begin{verbatim} plot({tan(x),-x},x=0..4*Pi,y=-15..0,numpoints=100); \end{verbatim} The answers can be generated by choosing intervals about the guesses, for example \footnotesize\rm \begin{verbatim} fsolve(tan(x)+x=0,x,1.9..2.2); # Answer: 2.03 approximately \end{verbatim}\normalsize\rm The other answers are $4.91$, $7.98$ approximately. Clicking on the left mouse button over the graph intersections gave numbers very close to these values. The truncated series solution to three terms is extracted from $$ T(x,t)=\sum_{n=1}^\infty B_n(\sqrt{\lambda_n}\cos(\sqrt{\lambda_n}x)+ \sin(\sqrt{\lambda_n}x))\exp(-\lambda_n t) $$ An approximate answer is $$ \begin{array}{lcl} T(x,t)&=& 0.517\left ( 2.03\cos( 2.03x)+\sin( 2.03x)\right ){e^{- 4.12t}}+\\ && -0.013\left ( 4.91\cos( 4.91x)+\sin( 4.91x)\right ){e^{- 24.1t}}+\\ && 0.002\left ( 7.98\cos( 7.98x)+\sin( 7.98x)\right ){e^{- 63.7t}} \end{array} $$ which is obtained by evaluating $B_n$ from the equation $$ \int_0^1 f(x)(\sqrt{\lambda_n}\cos(\sqrt{\lambda_n}x)+ \sin(\sqrt{\lambda_n}x))dx = B_n \int_0^1 (\sqrt{\lambda_n}\cos(\sqrt{\lambda_n}x)+ \sin(\sqrt{\lambda_n}x))^2dx $$ If you choose to use {\tt maple} to help with the computations, then here are some useful ideas for the first three terms: \footnotesize\rm \begin{verbatim} plot({tan(x),-x},x=0..4*Pi,y=-15..0,numpoints=100); k1:=fsolve(tan(x)+x=0,x,1.9..2.2); # First eigenvalue k1*k1=4.1 k2:=fsolve(tan(x)+x=0,x,4.9..5.1); # Second eigenvalue k2*k2=24.1 k3:=fsolve(tan(x)+x=0,x,7.9..8.0); # Third eigenvalue k3*k3=63.7 f:=x -> 1+sin(x)-x^2*(1+sin(1)): # Initial data g:=(x,k) -> k*cos(k*x)+sin(k*x): # Define Eigenfunction p:=int(g(x,k)*f(x),x=0..1): # Find B coefficients using q:=int((g(x,k))^2,x=0..1): # orthogonality+series formula B:= n -> subs(k=n,eval(p/q)); # Coefficients B b1:=B(k1); b2:=B(k2); b3:=B(k3); # 0.517, -0.013, 0.002 T:=(x,t) -> b1*g(x,k1)*exp(-k1*k1*t) # Insert coefficients in +b2*g(x,k2)*exp(-k2*k2*t) # formal series to make the +b3*g(x,k3)*exp(-k3*k3*t): # three-termed series T plot3d(T,0..1,0..1,axes=BOXED,style=PATCH); with(plots): animate(T(x,t),x=0..1,t=0..1,frames=25);# Animation of surface slices plot(T(x,0),x=0..1); plot(f,0..1); # plots of initial data plot(g(x,k1),x=0..1); # Plot first eigenfunction \end{verbatim}\normalsize\rm Some answers will print out unevaluated with sine and cosine terms. Use {\tt evalf} to convert the answers to decimals. } % Need problems on regular points, RSP, recursions, bessel functions, % legendre polynomials. \begin{exercises} \prob{25.}{(Elementary Power Series)} Prove that $\displaystyle\frac{1}{1+x} = \sum_{n=0}^\infty (-1)^nx^n$, valid for $|x| < 1$. Use this formula and methods of integration and differentiation of series to develop series formulas for $\ln(1+x)$, $\arctan(x)$, $1/(1+x)^2$, $1/(1+x)^3$. \SOL{25}{ Let $f(x)=1/(1+x)$. Then $f$ has a Taylor series expansion about $x_0=0$ given by Taylor's theorem on power series: $$f(x) = \sum_{n=0}^\infty \frac{f^{(n)}(x_0)}{n!} (x-x_0)^n.$$ To do the proof, apply Taylor's theorem to produce the claimed identity. The remaining results can be established mathematically by formally integrating and differentiating this series identity. } \end{exercises} \NOTES{ Since $f$ has a Taylor series expansion about $x_0=0$, the coefficient $f^{(n)}(0)/n!$ can be computed in {\tt maple} as follows: \footnotesize\rm \begin{verbatim} y:=1/(1+x): subs(x=0,y/0!); y:=diff(y,x): subs(x=0,y/1!); y:=diff(y,x): subs(x=0,y/2!); \end{verbatim}\normalsize\rm The formal series can be manipulated by {\tt maple} commands: \footnotesize\rm \begin{verbatim} sum('(-x)^n','n'=0..infinity); y:=unapply(",x); diff(y(x),x); int(y(x),x); \end{verbatim}\normalsize\rm The action of {\tt maple} is to find closed-form answers for series problems. Therefore, the answers printed are related to the claimed identities but {\tt maple} will not produce for you the expected {\em proof}. } \begin{exercises} \prob{26.}{(Regular-Point Series Solutions)} The equation $y'=y+x^3$, $y(0)=1$ has a power series solution of the form $y(x)=\sum_{n=0}^\infty c_nx^n$. Find a recursion relation for the coefficients $\{ c_n\}$ and determine the first five terms explicitly. \SOL{26}{ The subject of series solutions of differential equations has been criticised for its tedious calculations and resulting unpleasant recursion relations. Computer algebra systems do away with the tedious calculations. The approach is elegant and with the help of the computer algebra system, also the solution is elegant. Also settled below is the mystery of multiple forms of the answer, which are sometimes an infinite series and other times a closed-form expression in terms of elementary functions. The example below is selected to illustrate the use of the {\bf Linear differential equations series package} called {\tt ODEseries}. The functions to be used: \begin{center} \small\rm\begin{tabular}{|l|l|l|} \hline \bf Function &\bf Package &\bf Description\\ \hline \tt dsolve &maple library &solve a differential equation\\ \tt dsolve[series] &maple library &solve a differential equation\\ & &by the method of series\\ \tt diffeqtorec &ODEseries &find recursion relation for a linear\\ & &differential equation \\ \tt rectoproc &ODEseries &convert recursion to a function of n\\ & & It solves for $c_n$ in $\sum c_nx^n$.\\ \tt sum &maple library &sum of a series\\ \tt plot &maple library &plot a function\\ \hline \end{tabular}\normalsize\rm \end{center} The basic mathematical algorithm for solving a linear differential equation by the method of series is: \begin{itemize} \item[(a)] Assume $y = \sum_{n=0}^\infty c_nx^n$. Formally insert $y$ and its derivatives into the differential equation. Any terms in the differential equation that are not already power series are themselves expanded into series. \item[(b)] Collect terms on powers of $x$. \item[(c)] Set all coefficients of $x$ equal to zero to generate initial conditions for the coefficients $c_n$ and also a recursion relation for the $c_n$'s. \item[(d)] Solve the recursion. Substitute the answer back into the series for $y$. \end{itemize} {\large\bf Example}. Solve the differential equation $y' - y =x^2 - x$, $y(0)=1$ by two methods: (1) Find a closed form solution, (2) Find a series solution. Plot each on $[0,1]$ and compare the graphs. {\large\bf Solution}: To learn about series methods and gain intuition, we carry out steps (a) to (c) by hand, then return to the computer algebra system to verify the work and fill in the gaps. Assume $y$ is the series given by the algorithm. The derivative is: $$y' = \sum_{n=1}^\infty nc_nx^{n-1}$$ This equation and the original for $y$ are inserted into the differential equation to get $$\sum_{n=1}^\infty nc_nx^{n-1} - \sum_{n=0}^\infty c_nx^{n}=x^2-x$$ The terms on the right are a finite Taylor series. Collect all terms left, isolating the terms for $x^0$, $x^1$, $x^2$ from the rest: $$(c_1-c_0)x^0 + (2c_2-c_1+1)x^1 + (3c_3-c_2-1)x^2 + \sum_{n=4}^\infty nc_nx^{n-1} - \sum_{n=3}^\infty c_nx^{n}=0$$ The two infinite series can be added by changing index in the first by the rule $k=n-1$ (that makes the powers of $x$ match in the two series). $$(c_1-c_0)x^0 + (2c_2-c_1+1)x^1 + (3c_3-c_2-1)x^2 + \sum_{k=3}^\infty (k+1)c_{k+1}x^{k} - \sum_{n=3}^\infty c_nx^{n}=0$$ Set the coefficients to zero to obtain the initial conditions on the coefficients and also the recursion relations: $$(c_1-c_0)=0,~~ (2c_2-c_1+1)=0,~~ (3c_3-c_2-1)=0,~~ (k+1)c_{k+1} - c_k=0~~(k\ge 3)$$ The final recursion is solved to give $c_{k+1}=6c_3/(k+1)!$ for $k\ge 3$. The value of $c_0$ is $y(0)$ or $1$ from the series representation of $y$. Therefore, $c_0=1$, $c_1=1$, $c_2=0$, $c_3=1/3$. Finally, $$y = \sum_{n=0}^\infty c_nx^n = 1 + x + x^3/3 + x^4/12 + x^5/60 +\cdots$$ } \end{exercises} \NOTES{ The solution in {\tt maple} appears below. The first part is a standard example of how to solve a first order differential equation symbolically, then plot the solution. The second part does the main steps of the series solution algorithm above. In particular, steps (a) to (c) are done by {\tt diffeqtorec}. Step (d) is done by {\tt rectoproc} and {\tt sum}. It is emphasized that your {\em intuition} about series solutions originates in your own hand calculations. Use {\tt maple} to check your answer, to shorten tedious steps and to correct calculation errors. In particular, use {\tt maple} sparingly on series problems. \footnotesize\rm \begin{verbatim} # PART I. dsolve exact solution de:=diff(y(x),x)-y(x)=x^2-x: ic:=y(0)=1: dsolve({de,ic},y(x)); Y:=unapply(rhs("),x): plot(Y,0..1,title=`Y:=Exact Solution`); # PART II. classical series solution with recursion relations read(`/u/cl/maple/ODEseries`): p:=diffeqtorec({de,ic},y(x),c(n)); C:=rectoproc(p,c(n)): sum('C(n)*x^n','n'=0..6); S:=unapply(",x): plot(S,0..1,title=`S=Series Solution, n=0..6`); \end{verbatim}\normalsize\rm } \begin{exercises} \probx{26a.}{(Regular-Point Series Basis)} The general solution $y = Ay_1(x) + By_2(x)$ with arbitrary constants $A$ and $B$ of the differential equation $y'' + 2xy' + y=0$ is expressed in terms of power series solutions $y_1$ and $y_2$. Assume the initial conditions $y_1(0)=y_2'(0)=1$, $y_1'(0)=y_2(0)=0$. Find the first five terms of the power series for $y_1$ and $y_2$. \SOL{26a}{ To solve $y'' + 2xy' + y=0$ requires something besides the normal engineering `recipe' for linear differential equations. The recipe suggests using a characteristic equation $a\lambda^2 + b\lambda + c=0$ where $a$, $b$, $c$ are the coefficients of the second order linear differential equation. However, for this equation $b=2x$, which is {\em not constant}, therefore the recipe {\em fails to apply}! There is no elementary closed--form solution! Here is {\tt maple} code for solving the differential equation in closed--form: } \end{exercises} \NOTES{ \footnotesize\rm \begin{verbatim} eq:=diff(y(x),x,x) + 2*x*diff(y(x),x) + y(x)=0: dsolve(de,y(x)); \end{verbatim}\normalsize\rm There is no answer printed by {\tt maple}. The interpretation of the {\tt maple} answer is this: {\em no solution was found}! The series method applies to find a basis of solutions: {\large\bf Example}. Solve the initial value problem $y'' + 3xy' + y=0$, $y(0)=0$, $y'(0)=1$ by the series method for the first five terms. Plot the 5-term solution. {\large\bf Solution}: Assume $y = \sum_{n=0}^\infty c_nx^n$, then $y' = \sum_{n=0}^\infty nc_nx^{n-1}$, $y'' = \sum_{n=0}^\infty n(n-1)c_nx^{n-2}$. Inserting into the differential equation gives \footnotesize $$\sum_{n=0}^\infty n(n-1)c_nx^{n-2} +3x\sum_{n=0}^\infty nc_nx^{n-1} + \sum_{n=0}^\infty c_nx^n=0$$ \normalsize After re-indexing and term collection the series sums take this form: \footnotesize $$\sum_{n=0}^\infty (k+2)(k+1)c_{k+2}x^{k} +\sum_{n=0}^\infty 3nc_nx^{n} + \sum_{n=0}^\infty c_nx^n=0$$ \normalsize Therefore the recursion relation is $$c_0=y(0)=0,~~ c_1=y'(0)=1,~~ (n+2)(n+1)c_{n+2} + (3n+1)c_n=0$$ The recursions are solved for the first six terms. The answers: $$c_0=0,~~ c_1=1,~~ c_2=1/2,~~ c_3= -2/3,~~ c_4=0,~~ c_5=1/3$$ The answer to the differential equation, to six terms, is: $$y=x-{\frac {2\,x^{3}}{3}}+{\frac {\ x^{5}}{3}}$$ {\bf Checking the Answer in Maple}. A check in {\tt maple} reveals that correct answers were obtained in the major steps of the solution: \footnotesize\rm \begin{verbatim} read`/u/cl/maple/ODEseries`: eq:=diff(y(x),x,x) + 3*x*diff(y(x),x) + y(x)=0: ic:=y(0)=0,D(y)(0)=1: diffeqtorec({eq,ic},y(x),c(n)); # Determine recurrence relations C:=rectoproc(",c(n)): # Function for coefficients sum('C(n)*x^n','n'=0..6); # Series solution to 7 terms Y:=unapply(",x): plot(Y,0..1); # Plot approximate solution \end{verbatim}\normalsize\rm It is possible to check only the steps from the recursion to the end, as follows: \footnotesize\rm \begin{verbatim} r:={(n+2)*(n+1)*c(n+2) + (3*n+1)*c(n),c(0)=0,c(1)=1}: C:=rectoproc(r,c(n)): sum('C(n)*x^n','n'=0..6); \end{verbatim}\normalsize\rm Shown below is a method for obtaining the answer with {\tt dsolve}. The method obtains only the final answer, and not the steps between. \footnotesize\rm \begin{verbatim} eq:=diff(y(x),x,x) + 3*x*diff(y(x),x) + y(x)=0: # Define Differential Equation ic:=y(0)=0,D(y)(0)=1: # Define initial conditions dsolve({eq,ic},y(x),series); # Solve the differential equation q:=rhs("); # Grab right hand side of answer w:=subs(O(1)=0,q); # Drop higher order terms Y:=unapply(w,x); # Make a function plot(Y,0..1); # Plot the solution approximation \end{verbatim}\normalsize\rm } \begin{exercises} \probx{26b.}{(Airy's Equation)} Compute the recursion relations and the first five terms of the series solution $y(x)=\sum_{n=0}^\infty c_nx^n$ for Airy's differential equation $y'' - xy=0$, $y(0)=1$, $y'(0)=0$. \SOL{26b}{ It is known that Airy's differential equation has solution expressed in terms of Bessel functions. If you apply {\tt maple}'s function {\tt dsolve} to the Airy equation, then these special functions will be displayed. The expected solution is to be done by hand. Substitute a Taylor series into the differential equation, collect terms on powers of $x$ and set the coefficients to zero in order to obtain the recursion relations and initial conditions. Solve the recursion for the first six terms. The answer: $$ c_0=1,~~ c_1=c_2=0,~~ c_3=1/6,~~, c_4=c_5=0, ~~ y = 1 + x^3/6$$ } \end{exercises} \NOTES{ Some {\tt maple} source appears below to show how to check the series answer against the hand computation. \footnotesize\rm \begin{verbatim} # Airy's differential equation de:=diff(y(x),x,x)-x*y(x)=0: # Define Differential equation ic:=y(0)=1,D(y)(0)=0: # Define initial conditions Order:=6: # Set series remainder O(x^6) dsolve({de,ic},y(x),series); # Solve the DE by series method op(1,"): # Extract truncated answer rhs("): # Isolate terms on right side Y:=unapply(",x): # Make a real function of x Y(x); # Print answer for first 6 terms plot(Y,0..1); # Plot the approximate solution \end{verbatim}\normalsize\rm The {\tt ODEseries} package can be used to check the recursions and also the series solution. \footnotesize\rm \begin{verbatim} # Recursion relations and series expansion read`/u/cl/maple/ODEseries`: de:=diff(y(x),x,x)-x*y(x)=0: # Define Differential equation ic:=y(0)=1,D(y)(0)=0: # Define initial conditions diffeqtorec({de,ic},y(x),c(n)); # Determine recurrence relations C:=rectoproc(",c(n)); # Function for coefficients sum('C(n)*x^n','n'=0..5); # Series solution to 6 terms \end{verbatim}\normalsize\rm } \begin{exercises} \prob{27.}{(Euler Differential Equation Recipe)} An Euler equation $ax^2y'' + bxy' + cy=0$ can be solved by a {\bf recipe} similar to the one used for an equation with constant coefficients: \begin{quote}\small\rm Solve for the roots $\lambda$ in the Euler auxiliary equation $a\lambda(\lambda-1) + b\lambda + c =0$ and apply the recipe for the solution of constant equations to obtain a solution $z(t)$. Replace $t$ by $\ln(x)$ to obtain the general solution $y(x)$ of the Euler differential equation $ax^2y'' + bxy' + cy=0$. \end{quote} Apply the recipe to solve \begin{quote} (a) $x^2y'' + 4xy' + 2y=0$, (b) $x^2y'' + 3xy' + y=0$, \\ (c) $x^2y'' + 3xy' + 2y=0$. \end{quote} \SOL{27}{ The basic ideas can be seen in the example to follow. } \end{exercises} \NOTES{ {\large\bf Example}. Solve the Euler equation $x^2y'' + 4xy' + 2y=0$. {\large\bf Solution}: The Euler auxiliary equation is $\lambda(\lambda-1)+4\lambda + 2=0$ which factors as $(\lambda+1)(\lambda+2)=0$, giving roots of $-1$ and $-2$. By the recipe for constant equations $z(t)=c_1\exp(-t)+c_2\exp(-2t)$. Setting $t=\ln(x)$ or $x=\exp(t)$ gives $z(t)=c_1(1/x) + c_2(1/x^2)$. Therefore, the Euler equation general solution is $y(x)=c_1(1/x) + c_2(1/x^2)$. } \begin{exercises} \probx{27a.}{(Euler Equation Change of Variables)} Apply the change of variables $x=e^t$, $z(t)=y(e^t)$ to the Euler differential equation $ax^2y'' + bxy' + cy=0$ and obtain the constant equation $a(z'' - z') + bz' + cz=0$ whose characteristic equation is $a\lambda(\lambda-1) + b\lambda + c =0$. This calculation justifies the recipe for Euler equations claimed above. \SOL{27a}{ The chain from calculus is applied repeatedly to compute $dz/dt$ and $d^2z/dt^2$: $$\frac{dz}{dt} = \frac{dy}{dx}\frac{dx}{dt} = y'(x)e^t$$ $$\frac{d^2z}{dt^2} = \frac{d}{dt}\left( y'(x)e^t\right) = y'(x)\frac{d}{dt}\left( e^t\right) + \frac{d}{dt}\left( y'(x)\right)e^t$$ The last equality is from the product rule. To simplify further the last equality, use $\frac{d}{dt}(e^t) = e^t$ and also the chain rule once again: $$\frac{d}{dt}\left( y'(x)\right) = \frac{d}{dx}\left(y'(x)\right)\frac{dx}{dt} = y''(x)e^t$$ A final substitution on the right of $e^t$ by $x$ will connect all these calculations to the original differential equation, which contains terms of the form $x^2y''(x)$, $xy'(x)$ and $y(x)$. Substitute the formulas obtained in the original Euler differential equation to get the new equation in variables $z$ and $t$. } \end{exercises} \begin{exercises} \prob{28.}{(Regular Singular Points and Series Solutions)} Let $ax^2y'' + bxy' + cy=0$ be a second order Euler differential equation with arbitrary constant coefficients $a$, $b$, $c$, and $a\ne 0$. Prove that Euler's equation has a regular singular point at $x=0$. Observe that the classic Euler Auxiliary Equation $a\lambda(\lambda-1) + b\lambda + c =0$ is identical to the Indicial Equation of the Frobenius Theory. \SOL{28}{ The proof amounts to arranging the equation in the form $x^2p(x)y'' + xq(x)y' + r(x)y=0$ where $p$, $q$, $r$ are power series about $x=0$ and $p(0) \ne 0$. The classic Euler auxiliary equation is found formally by substitution of $\exp(\lambda x)$ into the differential equation. The Indicial equation of Frobenius theory is $p(0)\lambda(\lambda-1) + q(0)\lambda + r(0)=0$. } \end{exercises} \begin{exercises} \probx{28a.}{(Perturbations of Euler Equations)} Let $Ax^2y'' + B(x)xy' + C(x)y=0$ be a second order differential equation with regular singular point at $x=0$, $A\ne 0$ a constant and $B(x)$ and $C(x)$ power series convergent for $|x| x^4*BesselJ(2,x) - 2*x^3*BesselJ(3,x):|~~~\verb|evalf(f(4)-f(0));| \end{center} } \begin{exercises} \probx{30a.}{(Bessel Functions Integral Formulas)} Prove the formula $$\int_0^x J_1(x)\sin(x)dx = xJ_1(x)\sin(x)$$ $$ + (x\cos(x)-\sin(x))J_0(x)$$ \SOL{30a}{ Both sides of the formula have the same value at $x=0$. It suffices to show that the left and right side of the identity have identical derivatives. The derivative of the first term on the right is $$(x\sin(x)J_1(x))' = \sin(x)J_1(x) + x\cos(x)J_1(x) + x\sin(x)J_1'(x)$$ Apply the identity $xJ_1'(x)=xJ_0(x)-J_1(x)$ to simplify. The derivative of the second term on the right is $$(x\cos(x)-\sin(x))J_0(x))' = -x\sin(x)J_0(x) + x\cos(x)J_0'(x) - \sin(x)J_0'(x)$$ Apply the identity $J_0'(x)=-J_1(x)$ to simplify. } \end{exercises} %%\item[(e)] Solve for the general solution to the Airy differential %%equation $y'' + xy=0$. Transform it to a Bessel equation in $u$ and %%$z$ %%of order $p=1/3$ using $u=y/\sqrt(x)$, $z=(2/3)x^{3/2}$ and then apply %%the theory of Bessel's equation. %% %%\begin{solution}{3e} %% %%The transformation requires the chain rule to compute $y''$ in terms %%of %%$u$ and $z$. Here are the main steps: %%$$y' = \frac{dy}{dx} = \frac{d(\sqrt(x)u(x)}{dz}\frac{dz}{dx} = %%\left((1/2)x^{-1/2}\frac{dx}{dz}u(x) %%+\sqrt{x}\frac{du}{dz}\frac{dz}{dx}\right) \frac{dz}{dx}$$ %%This can be simplified because $dz/dx = \sqrt{x}$ into %%$y' = u(x)/(2\sqrt{x}) + x(du/dz)$. This is repeated to find a formula %%for $y''$: %%$$y'' = \frac{d}{dz}\left(u(x)/(2\sqrt{x})\right)\frac{dz}{dx} + %%\frac{du}{dz} + x\frac{d^2u}{dz^2}\frac{dz}{dx}$$ %% %% %%\end{solution} %% %%\item[(f)] Give three Bessel equations that illustrate cases 1, 2, 3 %%of %%the Frobenius theory. For each example, apply the Frobenius theory and %%display the general solution, using only the first three terms of each %%power series. \begin{exercises} \prob{31.}{(Legendre Polynomials)} Calculate the first six Legendre polynomials and display them in a table. \SOL{31}{ One way to compute the Legendre polynomials is by hand, using Rodrigue's formula. This method requires repeated differentiation to obtain the answer. The other way to do it by hand is to use the {\em Bonnet Recursion Formula} for Legendre polynomials, namely, $$(n+1)P_{n+1}(x)=(2n+1)xP_n(x) - nP_{n-1}(x).$$ Start with $P_0(x)=1$ and $P_1(x)=x$, then using $n=1$ in Bonnet's formula gives an equation for $P_2$: $2P_2=3xP_1-P_0$ or $2P_2=3x^3-1$. The polynomials $P_3$, $P_4$, \ldots are found similarly. } \end{exercises} \NOTES{ The {\tt maple} package {\tt orthopoly} contains routines for all the major sets of orthogonal polynomials. To display the Legendre polynomials: \footnotesize\rm \begin{verbatim} with(orthopoly): P(0,x);P(1,x);P(2,x);P(3,x);P(4,x);P(5,x); \end{verbatim}\normalsize\rm The Legendre polynomial of degree 5 is: $${\frac {63\,x^{5}}{8}}-{\frac {\ 35\,x^{3}}{4}}+{\frac {15\,x}{8\ }}$$ } \begin{exercises} \probx{31a.}{(Graphs of Legendre Polynomials)} Present five separate graphs showing the first five Legendre polynomials on $[-1,1]$, $P_0$ to $P_4$. \SOL{31a}{ All plots will be done in {\tt maple}. } \end{exercises} \NOTES{ The single plots can be done one at a time using the {\tt maple} command {\tt plot}, e.g., to plot $P_2$, use the command \begin{center} \verb|plot((3*x^2-1)/2,x=-1..1);| \end{center} Once the five plots are on the screen, shrink them to size, and print the screen. The {\tt maple} command {\tt plot} is able to plot several functions on a single graph. The syntax is {\t plot(\{ f,g\},a..b)} to plot two functions $f$ and $g$ on $[a,b]$. It works for any number of functions. A typical plot command for Legendre polynomials: \begin{center} \tt plot(\{P(2,x),P(3,x)\},x=-1..1); \end{center} } \begin{exercises} \prob{32.}{(Roots of Legendre Polynomials)} Produce a table which displays the roots of each of the first six (6) Legendre polynomials. \SOL{32}{ The {\tt orthopoly} package in {\tt maple} together with the {\tt solve} routine will find roots of Legendre polynomials. There are numerical routines for the roots, a classic reference being {\em Numerical Recipes}. The {\tt maple} commands for the $5^{\mbox{th}}$ Legendre polynomial: \begin{quote}\tt with(orthopoly): \\ evalf(solve(P(5,x)=0,x)); \\ \end{quote} Without the {\tt evalf} function, radicals will print, which is not too useful in a table. } \end{exercises} \begin{exercises} \probx{32a.}{(Legendre Polynomial Integral formula)} Let $P_n$ be the $n^{\mbox{th}}$ Legendre polynomial given by Rodrigues' formula, $n=1,2,\ldots$. Assume that $\int_{-1}^1 P_n^2(x)dx = 2/(2n+1)$ and $\int_{-1}^1 P_n(x)P_m(x)dx=0$ for $n\ne m$. Prove that $$\displaystyle\int_{-1}^1 xP_n(x)P_{n-1}(x)dx = \frac{2n}{4n^2-1}$$ \SOL{32a}{ The key identity to be used in the solution is the {\em Bonnet Recursion Formula} $$xP_n(x) = \frac{n+1}{2n+1}P_{n+1}(x) + \frac{n}{2n+1}P_{n-1}(x)$$ The idea of the derivation is to multiply the above equation by something to produce on the left the integrand of the claimed identity. The result $\int_{-1}^1 P_n(x)P_m(x)dx=0$ for $n\ne m$ is proved in a number of textbooks on engineering mathematics. The proof is based upon the fact that $P_n(x)$ is a solution of the Legendre differential equation $(1-x^2)y'' - 2xy' + n(n+1)y=0$. This equation is multiplied by $P_m(x)$. Then $n$ and $m$ are switched to obtain a second equation. Subtract the two equations and integrate. The method is used in Sturm-Liouville theory to prove that eigenfunctions are orthogonal. The result $\int_{-1}^1 P_n^2(x)dx = 2/(2n+1)$ is also proved in engineering mathematics textbooks. It follows from rather deep properties of orthogonality and the {\em Bonnet recursion} mentioned above. } \end{exercises} \begin{exercises} \prob{33.}{(Steady-State Temperature in a Plate)} The steady state temperature $T=T(x,y)$ in a rectangular plate satisfies the Laplace equation with boundary conditions $$ \begin{array}{l} \dd T_{xx}+T_{yy}=0,~~ 0 \le x \le 2\pi,~~ 0 \le y \le \pi,\\[1ex] \dd T(x,0)=0, \quad T(x,\pi)=0,~~ 0 \le x \le 2\pi,\\[1ex] \dd T(0,y)=f(y),~~ T(2\pi,y)=0,~~ 0 \le y \le \pi. \end{array} $$ Show how the method of separation of variables can be used to construct $T=T(x,y)$ for {\it any} function $f(y)$. Then find the solution for the special case that $f(y)=y$ for $0 \le y \le \pi$. \SOL{33}{ The boundary value problem for the $y$--variable after separation of variables is $Y'' + \lambda Y=0$, $Y(0)=Y(\pi)=0$. The problem for the $x$--variable is $X'' - \lambda X=0$, $X(2\pi)=0$. The series solution is $T=\sum_{n=1}^\infty c_n\sinh(n(x-2\pi))\sin(ny)$ where $c_n$ is given by the formula $$ c_n\sinh(-2n\pi) = \frac{2}{\pi}\int_0^\pi f(y)\sin(ny)dy$$ For the function $f(y)=y$ the answer is $c_n=\frac{2(-1)^{n}}{n\sinh(2n\pi)}$. The use of {\tt maple} in this problem is to check the answer and perform some standard integration. } \end{exercises} \NOTES{ \footnotesize\rm \begin{verbatim} m:=10: T:=(x,y) -> sum('c(n)*sinh(n*(x-2*Pi))*sin(n*y)','n'=1..m); diff(T(x,y),x,x)+diff(T(x,y),y,y); # Check T a solution of the pde int(sin(n*y)^2,y=0..Pi); # coefficient calculation int(y*sin(n*y),y=0..Pi); # coefficient calculation \end{verbatim}\normalsize\rm } \begin{exercises} \prob{34.}{(Steady-State Temperature in a Sector)} The steady state temperature $T=T(r,\theta)$ in the sector $0 \le r \le 1$, $0 \le \theta \le \pi/3$ (polar coordinates) satisfies the Laplace equation and boundary conditions $$ \begin{array}{l} \dd T_{rr}+{1 \over r}T_{r}+{1 \over r^2}T_{\theta \theta}=0,~~ 0 < r \le 1,~~ 0 \le \theta \le \frac{\pi}{3}, \\[1.5ex] \dd T(r,0)=0,~~ T(r,\pi/3)=0, ~~ 0 \le r \le 1,\\[1ex] \dd T(1,\theta)=f(\theta), ~~~ 0 \le \theta \le \pi/3. \end{array} $$ Show how the method of separation of variables can be used to construct $T=T(r,\theta)$, then find the solution for the case that $f(\theta)=5\sin(6 \theta)-11\sin(27 \theta)$. \SOL{34}{ The temperature $T$ is separated as the product $R(r)\Theta(\theta)$. The differential equations are $r^2R'' + rR' -\lambda R=0$ and $\Theta'' + \lambda \Theta=0$. The boundary conditions for the $\theta$--variable are $\Theta(0)=\Theta(\pi/3)=0$. This Sturm-Liouville problem has solution $\lambda=9n^2$, $\Theta(\theta)=\sin(3n\theta)$. The Euler auxiliary equation for $r^2R'' + rR' - 9n^2R=0$ has two real roots $\pm 3n$. The positive root $3n$ produces the continuous solution $R=r^{3n}$ while the negative root produces a singular solution. The temperature $T$ is given by the series formula $$ T = \sum_{n=1}^\infty c_n r^{3n}\sin(3n\theta)$$ The coefficients are determined by setting $r=1$ and $T(r,\theta)=f(\theta)$ in the formula. The result is a trigonometric series expansion for $f$ which can be solved by using orthogonality properties of the functions $\sin(3n\theta)$ on $[0,\pi/3]$. The answer: $c_n= \frac{6}{\pi}\int_0^{\pi/3} f(\theta)\sin(3n\theta)d\theta$. } \end{exercises} \NOTES{ The evaluation of the coefficients $c_n$ for the function $f(\theta)=5\sin(6 \theta)-11\sin(27 \theta)$ can be performed in {\tt maple} or evaluated with an integral table. Here is an example for a slightly different function $f$: \footnotesize\rm \begin{verbatim} f:= x -> 2*sin(6*x)-7*sin(33*x): # A different f! int(sin(3*n*x)^2,x=0..Pi/3); # Pi/6 p:=int(sin(3*n*x)*f(x),x=0..Pi/3); # coefficient calculation # 3*sin(n*Pi)*(-88+9*n^2)/((n-2)*(n+2)*(n-11)*(n+11)) # This is zero except for n=2 and n=11 limit(6*p/Pi,n=2); # Answer=2 limit(6*p/Pi,n=11); # Answer=-7 \end{verbatim}\normalsize\rm The answer for this $f$: $c_n=0$ except for $n=2$ and $n=11$, $c_2=2$, $c_{11}=-7$. } \begin{exercises} \prob{35.}{(Vibrating Drum)} Assume the displacement $u(r,t)$ of a vibrating drum of radius $R$ and wave speed $c$ is given by the Bessel Series $$\sum_{n=1}^\infty \left(a_n\cos\left(\frac{s_nct}{R}\right) +b_n\sin\left(\frac{s_nct}{R}\right)\right) J_0\left(\frac{s_nr}{R}\right)$$ where $\{s_n\}_{n=1}^\infty$ is the sequence of positive zeros of $J_0(x)$. It is known that the Bessel functions $\{J_0(s_nr/R)\}$ are orthogonal on $[0,R]$. %developed %in Problems 13-19, pp.778-779 of Derrick-Grossman %Use equation (23), p.779 Take $c=1$ and $R=1$. Assume initial values $$ u(r,0)=0, ~~ u_t(r,0)=f(r),~~ 0 \le r \le 1. $$ Find $a_n$ and $b_n$ in terms of integrals of $J_0$ and the sequence of roots $\{s_n\}_{n=1}^\infty$. Verify that for the case $f(r)=J_0(s_2\,r)-0.5J_0(s_4\,r)$ that $$u(r,t)=\frac{\sin(s_2t)}{s_2}J_0(s_2r) - \frac12\frac{\sin(s_4t)}{s_4}J_0(s_4r)$$ Finally, produce a snapshot graph of $u$ at $t=0.2$ and discuss its relationship to the drum vibration. \SOL{35}{ The boundary condition $u(r,0)=0$ gives rise to an equation with $0$ on the left and on the right a Bessel series with coefficients $a_n$. Multiply the equation by $J_0(s_nr)$. Integrate over $[0,1]$ to determine $a_n$. The boundary condition $u_t(r,0)=f(r)$ in a similar way gives rise to an equation with a Bessel series on the right. Repeat the ideas already outlined to determine $b_n$. When $f$ is specialized to $f(r)=J_0(s_2\,r)-0.5J_0(s_4\,r)$ most of the integrals in the formulas are zero. To justify claims, apply orthogonality properties. There are only two nonzero coefficients in the Bessel series, namely $b_2$ and $b_4$. } \end{exercises} \NOTES{ To find the roots $s_n$ in {\tt maple}: \footnotesize\rm \begin{verbatim} plot(BesselJ(0,x),x=0..20,y=-1..1); # Plot Bessel function J0 # Click near root: 2.4 approximately s1:=fsolve(BesselJ(0,x)=0,x=2.3..2.6); # find root 2.4 on [2.3,2.6] \end{verbatim}\normalsize\rm The answers for $s_1$ to $s_4$: 2.404825558, 5.520078110, 8.653727913, 11.79153444. The plot of snapshots of the solution can be accomplished by finding the second and fourth roots of $J_0(x)$ and forming the solution equation: \footnotesize\rm \begin{verbatim} s2:=fsolve(BesselJ(0,x)=0,x=5.5..5.6); # find root near x=5.5 s4:=fsolve(BesselJ(0,x)=0,x=11.7..11.9); # find root near x=11.8 u:=(r,t) -> sin(s2*t)*BesselJ(0,s2*r)/s2 -(1/2)*sin(s4*t)*BesselJ(0,s4*r)/s4: # Define solution plot(u(r,0.2),r=0..1); # Plot on [0,1] \end{verbatim}\normalsize\rm The calculations produce the deflection height of the drum head, which is a signed number between $-0.2$ and $0.2$. To visualize the motion of the drum head in time, apply animation to coordinates $x$, $y$, $z$ given by $x=r\cos(\theta)$, $y=r\sin(\theta)$, $z=u(r,t)$. The $t$ variable is sampled at time $t=0$ to $t=2$ in order to show the full deflections of the drum head. The implementation below for {\tt maple} is recommended only for high speed work stations. \footnotesize\rm \begin{verbatim} s2:=fsolve(BesselJ(0,x)=0,x=5.5..5.6); # find root near x=5.5 s4:=fsolve(BesselJ(0,x)=0,x=11.7..11.9); # find root near x=11.8 u:=(r,t) -> sin(s2*t)*BesselJ(0,s2*r)/s2 -(1/2)*sin(s4*t)*BesselJ(0,s4*r)/s4: # Define solution with(plots): drum:=[r*cos(y),r*sin(y),u(r,t)]: # Drum head surface, cylindrical coordinates animate3d(drum,r=0..1,y=0..2*Pi,t=0..2,frames=20,style=PATCH,axes=BOXED); \end{verbatim}\normalsize\rm Printing of the plot can be done with the unix shell command {\tt paw} discussed earlier. There is no support for printing animation, but it is possible to produce a 3D plot for one of the frames and then print that plot: \footnotesize\rm \begin{verbatim} # Assume u(r,t) already defined as above drum1:=[r*cos(y),r*sin(y),u(r,0.15)]: # Drum head surface, t=0.15 plot3d(drum1,r=0..1,y=0..2*Pi,style=PATCH,axes=BOXED); \end{verbatim}\normalsize\rm } \begin{exercises} \prob{36.}{(Steady-State Temperature in a Sphere)} The steady state temperature $T=T(r,\phi)$ in the sphere $0 \le r \le 1$, $0 \le \phi \le \pi$ (spherical coordinates) satisfies the Laplace equation $$ T_{rr}+{2 \over r}T_{r}+{1 \over r^2\sin\phi}(\sin\phi T_{\phi})_{\phi}=0 $$ for $ 0 < r \le 1,\quad 0 < \phi < \pi$, as well as the boundary condition $$ T(1,\phi)=f(\phi), \quad \text{for} \quad 0 \le \phi \le \pi. $$ Let $P_n(x)$ be the $nth$ Legendre polynomial. The product solutions of the Laplace equation are known to be $$ T_n(r,\phi)=r^nP_n(\cos\phi) %(see p. 287 for $P_0,\cdots,P_4$). $$ Verify by a direct argument that the product solutions are indeed solutions of the Laplace equation. Construct from the products the solution $T=T(r,\phi)$ for {\it any} function $f(\phi)$. Then find $T=T(r,\phi)$ for the special case that $f(\phi)= 4\cos^3 \phi - 3\cos \phi$. \SOL{36}{ The solution is on slides. Maple notes below. } \end{exercises} \NOTES{ The solution begins by verifying that the given product solution is indeed a solution of the Laplace equation. \footnotesize\rm \begin{verbatim} with(orthopoly): # Load package for Legendre polynomials u:= (r,phi) -> r^n*p(cos(phi)): # Define product solution e1:=diff(u(r,phi),r,r); # First term of pde e2:=(2/r)*diff(u(r,phi),r); # Second and third terms of pde e3:=(1/r^2)*(1/sin(phi))*diff(sin(phi)*diff(u(r,phi),phi),phi); simplify(e1+e2+e3=0); # simplify pde simplify("/r^(n-2)); # divide out common factor subs(cos(phi)=x,"); # replace cosine terms by x \end{verbatim}\normalsize\rm The final result is the differential equation $$ n^2p(x) + np(x) - 2xp'(x) + p''(x) -x^2p''(x) = 0$$ This is recognized as Legendre's differential equation of order $n$, therefore $p=P_n(x)$ is a solution. The differential equation is verified. Assume that an abritrary $f$ is given. To construct $T$ the orthogonality of the Legendre polynomials on $[-1,1]$ will be used. The equation to be satisfied is $$ f(\phi) = T(1,\phi) = \sum_{n=1}^\infty c_nP_n(\cos(\phi))$$ The basic idea is to multiply the above equation by $P_n(\cos(\phi))\sin(\phi)$ and integrate over $[0,\pi]$ to determine $c_n$. If this is carried out for $f(\phi)=4\cos^3 \phi - 3\cos \phi$, then $$c_1=3/5,~~ c_3=2/5,~~ c_n = 0,~~ \mbox{otherwise} $$ The process of evaluation of the $c_n$'s involves knowledge of Legendre polynomials. In particular, by a change of variables, $$\int_0^\pi \cos^3(\phi)P_n(\cos(\phi))\sin(\phi)d\phi = \int_{-1}^{1} -x^3P_n(x)dx$$ To evaluate the integral exactly use the formulas $P_3(x)=(5x^3-3x)/2$ and $P_1(x)=x$ to write $-x^3 = (-2/5)P_3(x)+(-3/5)P_1(x)$. Then $$\int_0^\pi \cos^3(\phi)P_n(\cos(\phi))\sin(\phi)d\phi = (-2/5)\int_{-1}^{1} P_3(x)P_n(x)dx + (-3/5)\int_{-1}^{1} P_1(x)P_n(x)dx$$ By orthogonality the only contributions to the sum for $T$ are for indices $n=1$ and $n=3$. } \ifsolutions\else\end{multicols}\fi \eject {\Large\bf Notes on Printing Multiple Graphs} The problems suggest making several graphs. While easy on a video monitor it generates too much paper in a report. Explained here for our local unix system is {\em a method for printing several graphs on a single sheet of paper}. The idea is to {\bf display} the graphs on the screen, {\bf resize} and {\bf move} each until the screen is correct, then {\bf print the screen}. Up to 8 graphs can be displayed by this mechanism. All will appear on a single sheet of paper. The {\tt resize} and {\tt move} commands are available on mouse menu 1. Press \MOUSELEFT{} to activate the menu. With these two commands it is possible to shrink a graphic and reposition it anywhere on the screen. Once the screen is loaded with the graphics it can be printed to the {\tt b129lab1} laser printer by the following unix shell command. The command is to be issued in a console or terminal window. \begin{quote}\tt alias paw "xdpr -root -device ps -portrait -Pb129lab1" \\ paw \end{quote} One issue of the {\tt alias} command above installs the command {\tt paw} and thereafter, until logout, you may use {\tt paw} to print the screen. The command name {\tt paw} stands for {\large\bf P}rint {\large\bf A}ll {\large\bf W}indows (sometimes called print screen). Unix wizards can install the command {\tt paw} permanently by copying the {\tt alias} line above to the unix resource file {\tt .cshrc} in their root directory. The file {\tt .cshrc} is read at login and therefore the command {\tt paw} will always be available. Another possibility is to copy the following lines into the file {\tt \$HOME/.twmrc} (the first line should already be in the source!): \begin{quote} \verb|"Print Window" !"xdpr -device ps &"| \verb|"Print Screen" !"xdpr -root -device ps -portrait &"| \end{quote} Newer accounts probably already have these features installed. If you have an old account, then the file ".twmrc" can be updated by the following steps in a local terminal window: \begin{center} \verb|mv $HOME/.twmrc $HOME/.twmrc.old| \\ \verb|cp -p /usr/skel/.twmrc $HOME/.twmrc| \end{center} \bigskip \begin{quote}\small\sf \noindent {\large\bf Instructions for Math 353, Applied Partial Differential Equations}. Below is a list of thirty-six (36) problems required for 100\% credit. A solution must be attempted on each of the 36 problems. Collection of the take-home exam will be made in weekly increments. \end{quote} \SCHEDULE \end{document}