Math 353 Applied Partial Differential Equations Problem Notes on Take-Home Exams

Periodic Functions. Let the function f_{0}(x) be defined on [−π, π] by the formula

f_{0}(x) = \{ \begin{matrix}0, & for − π ≤ x < 0, \\ sin(x), & for 0 ≤ x ≤ π. \end{matrix}

1. (Periodic Extension) Define f(x) to be the 2π–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π, 5π]. The formula should depend only upon the symbols f_{0} and 2π and not upon the actual piecewise definition of the function f_{0}.

Solution to 1. The formula is a piecewise defined function which equals f_{0}(x) on [−π, π] and on the intervals [−3π, −π], [π, 3π], [3π, 5π] takes the function values f_{0}(x + 2π), f_{0}(x − 2π), f_{0}(x − 4π), respectively. All of this is verified graphically.

An expected and essential part of the solution is to chase down a reference to translations f_{0}(x + c) in a calculus textbook. According to the theory of rigid motions, each 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.

Periodic Extension. The idea of a 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 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 ≤ t ≤ T. Example. Write down for [−π, 2π] the periodic extension f of period π for the base function

f_{0}(x) = \{ \begin{matrix}sin(x) & 0 ≤ x ≤ π/2, \\ 2(π − x)/π & π/2 < x ≤ π. \end{matrix}

Solution:

f(x) = \{ \begin{matrix}f_{0}(x + π) & −π ≤ x < 0, \\ f_{0}(x) & 0 ≤ x < π, \\ f_{0}(x − π) & π ≤ x < 2π \end{matrix} = \{ \begin{matrix}sin(x + π) & −π ≤ x < −π/2, \\ 2(−x)/π & −π/2 ≤ x < 0, \\ sin(x) & 0 ≤ x ≤ π/2, \\ 2(π − x)/π & π/2 ≤ x < π. \\ sin(x − π) & π ≤ x < 3π/2, \\ 2(2π − x)/π & 3π/2 ≤ x < 2π. \end{matrix}

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.

1a. (Extension Formulas) Let f_{1}(x) = | sin(x)| on [−π, 0]. Define g(x) to be the even 2π–periodic extension of f_{1}(x) to the whole real line. Define h(x) to be the odd 2π–periodic extension of f_{1}(x) to the whole real line. Give simple algebraic formulas for g(x) and h(x).

Solution to 1a. Answers: g(x) = | sin(x)|, h(x) = − sin(x). The answers were obtained graphically by identifying the base function on [−π, π] for each of g and h.

Odd graphs are formed by drawing the graph on x ≥ 0, then rotate that graph 90^{∘} 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 ≤ 0.

2. (Graphing of Extensions) Let f(x) be the 2π–periodic extension of f_{0}(x) to (−∞, ∞). Graph f(x) for −4π < x < 2π.

Solution to 2. The graph can be drawn manually by regarding the graph of the base function f_{0} on [−π, π] as the xerox master. Make three xerox copies and paste the four duplicates side–by–side to make an initial graph of f. The exact size is obtained by snipping off the left and right edge of the figure, to make the interval [−4π, 2π].

Computer graphing can be done using a trick. The equation 2f_{0}(x) = | sin x| + sin x can be guessed by examining this expression on the interval [−π, π]. Therefore, feed the function \frac{1}{2}(| sin x| + sin x) and the domain [−4π, 2π] to your favorite computer graphing program.

This problem is not about formulas and there is no algebra expected. However, that is because all the groundwork was done in the previous problem. See in particular therein the remarks about translations and the relationship of formulas to the panels constructed here as xerox copies.

2a. (Even and Odd Periodic Functions) Let f(x) be the 2π–periodic extension of f_{0}(x) to (−∞, ∞). An approximation to 2f(x) is the trigonometric polynomial 2/π+sin(x)−(4/3) cos(2x)/π−(4/15) cos(4x)/π. Which terms of the expression will approximate an even function and which terms will approximate an odd function? Explain.

Solution to 2a. The constant term and the cosine terms represent an even function. The other terms make up an odd function. To verify the statements about even and odd functions, use the definitions: F (−x) = F (x) (F even), F (−x) = −F (x) (F odd).

2b. (Graphing Finite Trigonometric Series) Graph on [−π, π] the approximation to f_{0}(x) given by the trigono-

metric series

\frac{1}{π} + \frac{1}{2} sin(x) − \frac{2}{3π} cos(2x) − \frac{2}{15π} cos(4x).

Solution to 2b. Use a computer graphing program or a graphing calculator to get the approximate shape. The graph should look something like the base function f_{0}.

Maple notes on 1a–2b.

Example. Let f_{0}(x) = x sin(x) on [−π, π]. Plot its 2π–periodic extension f on [−π, 3π]. Solution: The maple code below is succinct and deserving of study since it is recursive. The interpretation of the graph: This is the sum of the Fourier series of f_{0}(x).

periodicextension:=proc(x,f0,a,b)
if 0<evalf(a-x) then
RETURN(periodicextension(x+b-a,f0,a,b)):
fi:
if 0<evalf(x-b) then
RETURN(periodicextension(x-b+a,f0,a,b)):
fi:
RETURN(eval(f0(x))):
end:
#
f0:= u -> 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));

The proc command in maple creates a function subroutine. A simple example is f:= x -> sin(x)*cos(x): which is coded with proc as f:=proc(x) sin(x)*cos(x): end:, that is, one is equivalent to the other. The variable x in this coding is a dummy variable or formal argument to the function f. The command Int differs from 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 evalf command causes the integration to be done numerically by a clever adaptive scheme kept secret from most persons who use maple. To find out the secrets, read the maple help file on userinfo. Maple plot error messages. The plot command plot(f(x),0..1) is often tried by beginners. It usually returns an empty plot or in some cases the silly message can't evaluate boolean. The fix: use

plot(f,0..1) or plot(f(x),x=0..1).

Example. Let f_{0}(x) = sin(x) for 0 ≤ x < π, f_{0}(x) = − sin(2x) for −π ≤ x < 0. Graph its 2π–periodic extension on [−π, 4π]. Solution: The maple code for procedure 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) ≡ \{ \begin{matrix}1 & t ≥ 0, \\ 0 & t < 0, \end{matrix} H(t − a) − H(t − b) = \{ \begin{matrix}1 & a ≤ t < b, \\ 0 & otherwise. \end{matrix}

The second formula above can be justified directly from the definition of H(t). In maple the Heaviside function H(t) is coded as 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) ≡ \{ \begin{matrix}sin(x) & 0 ≤ x < π \\ − sin(2x) & −π ≤ x < 0 \end{matrix} = sin(x) \{ \begin{matrix}1 & 0 ≤ x < π \\ 0 & otherwise \end{matrix} + (− sin(2x)) \{ \begin{matrix}1 & −π ≤ x < 0 \\ 0 & otherwise \end{matrix}

Therefore, from the Heaviside formulas,

f_{0}(x) = sin(x)(H(x − 0) − H(x − π)) + (− sin(2x))(H(x − (−π)) − 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.

# 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);

Fourier Series. Let the function f_{0}(x) be defined on 0 ≤ x < 2π by the formula

f_{0}(x) = \{ \begin{matrix}−1, & for 0 ≤ x < 2π/3 \\ 0, & for 2π/3 ≤ x < 4π/3 \\ 1, & for 4π/3 ≤ x < 2π. \end{matrix}

3. (Graphing Periodic Extensions) Define f(x) to be the 2π–periodic extension of f_{0}(x) to the whole real line. Sketch by hand the graph of f for −2π < x < 4π.

Solution to 3. Graph f_{0} on [0, 2π], which is a step function with values −1, 0 and 1. This is the base function (i.e., the xerox master). Assemble three copies and paste them together horizontally. The resulting figure is the answer.

3a. (Heaviside Formula) Give an explicit formula for f(x) when −2π ≤ x < 4π in terms of the Heaviside function H(t). The answer is to be a piecewise constant function with values −1, 1, 0 containing no reference to f_{0}.

Solution to 3a. Represent the base function f_{0} by the Heaviside expression f_{0}(x) = (−1)[H(t − a) − H(t)] + (1)[H(t − c) − H(t − b)] with c = 2π, b = 2c/3 and a = c/3. Then f_{0}(x) = 0 unless 0 ≤ x ≤ 2π. Copies of the base function shifted left and right are f_{0}(x + 2π) and f_{0}(x − 2π). Using the graph of f as a guide, the answer is f(x) = f_{0}(x + 2π) + f_{0}(x) + f_{0}(x + 2π). See the examples below for the fine detail and method.

Heaviside’s method for piecewise–defined functions. Example. Let f_{0}(x) = x for 0 < x < 1, f_{0}(x) = −2x for −1 < x < 0 and define ω = 2. Graph the ω–periodic extension of f_{0} on [−2, 3]. Solution: This solution uses only the Heaviside function to prepare the function f for plotting. We use:

H(t) ≡ \{ \begin{matrix}1 & t ≥ 0, \\ 0 & t < 0, \end{matrix} H(t − a) − H(t − b) = \{ \begin{matrix}1 & a ≤ t < b, \\ 0 & otherwise. \end{matrix}

First, we must write f_{0} in terms of H:

f_{0}(x) ≡ \{ \begin{matrix}x & 0 ≤ x < 1 \\ −2x & −1 ≤ x < 0 \end{matrix} = x \{ \begin{matrix}1 & 0 ≤ x < 1 \\ 0 & otherwise \end{matrix} + (−2x) \{ \begin{matrix}1 & −1 ≤ x < 0 \\ 0 & otherwise \end{matrix}

Therefore, from the Heaviside formulas,

f_{0}(x) = x(H(x − 0) − H(x − 1)) + (−2x)(H(x − (−1)) − H(x − 0)) = 3xH(x) − xH(x − 1) − 2xH(x + 1).

Finally, we write f in terms of H on [−2, 3]:

f(x) ≡ \{ \begin{matrix}f_{0}(x + ω) & −2 ≤ x < −1 \\ f_{0}(x) & −1 ≤ x < 1 \\ f_{0}(x − ω) & 1 ≤ x < 3 \end{matrix}

f(x) = f_{0}(x + 2) \{ \begin{matrix}1 & −2 ≤ x < −1 \\ 0 & otherwise \end{matrix} + f_{0}(x) \{ \begin{matrix}1 & −1 ≤ x < 1 \\ 0 & otherwise \end{matrix} + f_{0}(x − 2) \{ \begin{matrix}1 & 1 ≤ x < 3 \\ 0 & otherwise \end{matrix}

f(x) = f_{0}(x + 2)(H(x + 2) − H(x + 1)) + f_{0}(x)(H(x + 1) − H(x − 1))
+f_{0}(x − 2)(H(x − 1) − H(x − 3))
= f_{0}(x + 2)H(x + 2) + (f_{0}(x) − f_{0}(x + 2))H(x + 1)
+(f_{0}(x − 2) − f_{0}(x))H(x − 1) − f_{0}(x − 2)H(x − 3)

This form of f is suitable for input into a graph program that has a Heaviside function intrinsic. For example:

# maple V.1, V.2, V.3
readlib(Heaviside): H:= t -> 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);

4. (Fourier Series) Let f be the 2π–periodic extension of f_{0} to (−∞, ∞). Show that f is 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.

Solution to 4. The property of oddness can be shown on the symmetric interval [−π, π] because f is by definition 2π–periodic. Since f is zero in this symmetric interval except for |x| ≤ 2π/3, we only need to test oddness on |x| ≤ 2π/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π/3 ≤ x < 0 and f(x) = −1 on 0 ≤ x < 2π/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}{π} ∫_{−π}^{π} f_{0}(x) sin(nx)dx = \frac{2}{π} ∫_{0}^{2π/3}(−1) sin(nx)dx = \frac{2}{π} \frac{1}{n}(cos(2nπ/3) − 1). The reduction to twice the integral over [0, 2π/3] results from the integrand being even. On this new interval, f_{0}(x) = −1, except for the points x = 2π/3 to x = π, where f_{0}(x) = 0.

4a. (Sum of a Fourier Series) Let f be the 2π–periodic extension of f_{0} to (−∞, ∞). Determine the sum of the Fourier series of f for each value of x.

Solution to 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) ≡ \frac{1}{2}(f(x + 0) + f(x − 0)) for each x. This answer agrees with f except at x = ±π/3 and x = 0 (and 2π 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 (π/3) = −1/2, F (−π/3) = 1/2, F (0) = 0.

4b. (Convergence Theorem) Let f be the 2π–periodic extension of f_{0} to (−∞, ∞). For which values of x does the graph of f(x) differ from the graph of the Fourier series of f? Reference a theorem.

Solution to 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 ≠ f(x); at x = π/3, F (x) = −1/2 ≠ f(x) and at x = −π/3, F (x) = 1/2 ≠ f(x). Therefore, the graphs differ at x = 0, x = ±π/3 and all 2π–multiples of these points.

Evaluating Fourier coefficients. The evaluation of a_{n} and b_{n} amounts to evaluation of integrals. Example. Evaluate the Fourier coefficients a_{n} and b_{n} for the 2π–periodic function f determined by the base function

f_{0}(x) = \{ \begin{matrix}1 & −π ≤ x < 0 \\ −1 & 0 ≤ x < π \end{matrix}

Calculus Solution: The integrals can be split up as ∫_{−π}^{0} plus ∫_{0}^{π}. In the first integral f is 1 and in the second integral f is −1. For example,

b_{n} = \frac{1}{π} ∫_{−π}^{0} cos(nx)(1)dx + \frac{1}{π} ∫_{0}^{π} cos(nx)(−1)dx = \frac{1 − cos nπ}{nπ} .

The answers are a_{0} = 0, a_{n} = 0, b_{n} = (1 − cos(nπ))/(nπ). Therefore, the Fourier Series of f is

∑_{n=1}^{∞} \frac{1 − cos(nπ)}{nπ} sin(nπx)

Maple Solution: The function f can be written in terms of Heaviside’s function H(t) as f(x) = (1)(H(x + π) − H(x − 0)) + (−1)(H(x − 0) − H(x − π)) = H(x + π) − H(x − π) − 2H(x). The integrals to evaluate are:

a_{n} = \frac{1}{π} ∫_{−π}^{π} f(x) cos(nx)dx, b_{n} = \frac{1}{π} ∫_{−π}^{π} f(x) sin(nx)dx.

Using maple as the integral table:

# 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);

The reason for the strange disassembly of the function into terms with a single Heaviside function has to do with the integration routines in maple. It works when coded as above. It fails for the composite

int(f(x)*cos(n*x),x=-Pi..Pi).

The answers agree with those reported for the hand calculation, above.

Example. Determine the values of x for which the graph of f(x) differs from the graph of its formal Fourier series. Solution: The 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 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, ±1, ±2, . . .. 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). Example. Verify the formula

\frac{π}{4} = 1 − \frac{1}{3} + \frac{1}{5} − \frac{1}{7} + − · · · .

Solution: To verify it, use the Fourier series for the periodic extension f of the even base function f_{0}(x) = 1, −π/2 < x < π/2, f_{0}(x) = 0, π/2 < x < 3π/2. Then substitute x = 0 in the series formula

\frac{f(x + 0) + f(x − 0)}{2} = \frac{1}{2} + \frac{2}{π} (cos x − \frac{1}{3} cos 3x + \frac{1}{5} cos 5x − ⋯)

The result is

1 = \frac{1}{2} + \frac{2}{π} (cos 0 − \frac{1}{3} cos 0 + \frac{1}{5} cos 0 − ⋯)

which reduces to the claimed formula.

5. (Series Identity) \frac{π^{2}}{8} = ∑_{n=0}^{∞} \frac{1}{(2n + 1)^{2}} .

Solution to 5. To prove the series identity requires an explicit Fourier series expansion for the function f(x) = |x| on −π ≤ x ≤ π, namely,

|x| = \frac{π}{2} − \frac{4}{π} ∑_{n=0}^{∞} \frac{1}{(2n + 1)^{2}} .

The idea is to apply the Fourier Convergence Theorem at some particular value of x in −π < x < π, in order to evaluate the infinite series. Suitable values are those for the standard angles in trigonometry: x = 0, π/6, π/4, π/3, π/2. The Fourier series of f converges always to

\frac{1}{2} lim_{t→x+} f(t) + \frac{1}{2} lim_{t→x−} f(t),

which reduces to f(x) at a point of continuity of f. The f in this formula is the sawtooth wave, formed from the base function |x| on the interval [−π, π].

5a. (Series Identity) \frac{π^{2}}{6} = ∑_{n=1}^{∞} \frac{1}{n^{2}} .

Solution to 5a. This formula results by substitution of some value of x into the Fourier series for f(x) = x^{2} on −π < x < π. Half Range Expansions. The formulas for 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} + ∑_{n=1}^{∞} a_{n} cos(nπx/T), a_{n} = \frac{2}{T} ∫_{0}^{T} f(x) cos(nπx/T)dx

Example. Find the first three coefficients of the even Half-Range expansion for f(x) = sin(x) on [0, π] and graph on [−π, π] the Fourier series truncated to three terms. Maple Solution: Apply the previous formulas as follows:

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);

The answers: a_{0} = 4/π, a_{1} = 0, a_{2} = −4/(3π). The truncated series is g(x) = 4(1 − cos(2x)/3)/π. The limit command in maple is used here to resolve fractions that might otherwise cause a division by zero under direct substitution.

Function f_{1}(x). Let the function f_{1}(x) be defined on [0, π] by the formula

f_{1}(x) = \{ \begin{matrix}sin(x), & for 0 ≤ x ≤ π/2, \\ 2(π − x)/π, & for π/2 ≤ x ≤ π. \end{matrix}

6. (Half Range Expansions) Define g(x) to be the even 2π–periodic extension of f_{1}(x) to the whole real line. Define h(x) to be the odd 2π–periodic extension of f_{1}(x) to the whole real line. Give formulas for g(x) and h(x) valid on [2π, 4π] and find the first three Fourier coefficients of the even half-range expansion for f_{1}(x).

Solution to 6. The answers to the problem may be found by direct integration starting with the formula for n ≥ 1:

a_{n} = \frac{2}{π} ∫_{0}^{π} f_{1}(x) cos(nx)dx.

To integrate, split the integral into two pieces, one over [0, π/2] and the other over [π/2, π]. 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(π − x) cos(nx), requires integration by parts, using u = 2(π − x) and dv = cos(nx)dx.

The answers are a_{0} = \frac{4 + π}{4π}, a_{1} = \frac{4 − π}{π^{2}} , a_{2} = \frac{−2}{3} \frac{3 + π}{π^{2}} .

Maple notes on 6. The tedium of integral evaluation can be lessened by using maple:

a0:=(1/Pi)*(int(sin(x),x=0..Pi/2)+ int((2/Pi)*(Pi-x),x=Pi/2..Pi));
a1:=(2/Pi)*(int(sin(x)*cos(x),x=0..Pi/2)+int((2/Pi)*(Pi-x)*cos(x),x=Pi/2..Pi));
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));

6a. (Half Range Expansions) Graph g(x), the even 2π–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π ≤ x ≤ 4π.

Solution to 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 + π}{4π} + \frac{4 − π}{π^{2}} cos(x) + \frac{−2}{3} \frac{3 + π}{π^{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 [−π, π], which is the sine curve | sin x| on |x| ≤ π/2 plus two straight lines. Four copies of this curve are added to the graph, on the intervals [nπ, nπ + 2π] for n = 1, 3, −3, −5. Clip the graph to give the final window [−4π, 4π]. The curve g(x) looks like a piece of decorative molding: it’s entirely in the upper half plane, no negative cycles!

Fourier Integral. Consider the following table where it is assumed that a > 0:

F (x)  A(t)  B(t)

\{ \begin{matrix}1 & |x| < a \\ 0 & otherwise \end{matrix}

\frac{2 sin(at)}{πt} 0

\{ \begin{matrix}e^{−ax} & x ≥ 0 \\ 0 & otherwise \end{matrix}

\frac{1}{π} \frac{a}{t^{2} + a^{2}} \frac{1}{π} \frac{t}{t^{2} + a^{2}}
\frac{a}{x^{2} + a^{2}}  exp(−at)  0

\{ \begin{matrix}sin(x) & |x| < π \\ 0 & otherwise \end{matrix}

0 \frac{2}{π} \frac{sin(πt)}{1 − t^{2}}

The methods used to verify the table are as follows.

Method 1. Insert the definition of F into the corresponding sine and cosine integrals

A(t) = \frac{1}{π} ∫_{−∞}^{∞} F (u) cos(tu)du,

B(t) = \frac{1}{π} ∫_{−∞}^{∞} F (u) sin(tu)du

and perform the indicated integration over (−∞, ∞).

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

∫_{0}^{∞} [A(t) cos(tx) + B(t) sin(tx)] dt

and perform the indicated integration over (0, ∞) to obtain the F in table column 1.

The hand computations implied by the two methods of verification can be reduced to first year calculus problems. An example: Example. Find the Fourier sine and cosine integrals for the function f(x) = 1, 0 < x < 2, f(x) = 0, otherwise. Solution: The integrals and their evaluation:

A(t) = \frac{1}{π} ∫_{−∞}^{∞} f(x) cos(tx)dx = \frac{1}{π} ∫_{0}^{2} cos(tx)dx = \frac{sin 2t}{πt}.

B(t) = \frac{1}{π} ∫_{−∞}^{∞} f(x) sin(tx)dx = \frac{1}{π} ∫_{0}^{2} sin(tx)dx = \frac{1 − cos 2t}{πt}.

Example. Insert the Fourier sine and cosine integrals formulas

A(t) = \frac{sin 2t}{πt}, B(t) = \frac{1 − cos 2t}{πt}

into the Fourier integral formula

f(x) = ∫_{0}^{∞}[A(t) cos tx + B(t) sin tx]dt

and integrate. The answer should be the function f(x) = 1, 0 < x < 2, f(x) = 0, otherwise, except at jumps of f. Solution: The integration problem is

∫_{0}^{∞} [\frac{sin 2t}{πt} cos tx + \frac{1 − cos 2t}{πt} sin tx] dt.

It reduces to

\frac{1}{π} ∫_{0}^{∞} [\frac{sin 2t cos tx + sin tx − cos 2t sin tx}{t}] dt.

Integrating over [0, N] instead of [0, ∞) gives the answer

\frac{1}{π} (Si(xN) − Si(xN − 2N)) .

Let N → ∞, then the result is sign(x) − sign(x − 2), which agrees with f except possibly at the jumps.

7. (Fourier Integral) Verify lines 1 and 2 of the table.

Solution to 7. Method 1 is applied to establish line 1 of the table. Then B(t) = 0 because F is even and

A(t) = \frac{1}{π} ∫_{−∞}^{∞} f(x) cos(tx)dx = \frac{1}{π} ∫_{−a}^{a} f(x) cos(tx)dx = \frac{2}{π} \frac{sin at}{t}.

Method 1 is applied to obtain line 2 of the table. The objective is to insert the formula for F into the formulas for A and B to derive the two table entries \frac{1}{π} \frac{a}{a^{2}+t^{2}} and \frac{1}{π} \frac{t}{a^{2}+t^{2}} . The first integral is

A(t) = \frac{1}{π} ∫_{−∞}^{∞} F (x) cos txdx = \frac{1}{π} ∫_{0}^{∞} e^{−ax} cos txdx.

The integral on the right can be viewed as the direct Laplace transform of the function x → cos(tx) evaluated at s = a. From transform tables the Laplace transform is s/(s^{2} + t^{2}), hence after substitution it reduces to \frac{1}{π} \frac{a}{a^{2}+t^{2}} . The second integral B(t) is done the same way.

7a. (Fourier Integral) Verify lines 3 and 4 of the table.

Solution to 7a. Method 2 is applied to obtain line 3 of the table. In the Fourier integral replace A(t) by e^{−a|t|} and B(t) by 0 to obtain the integration problem

∫_{0}^{∞}[A(t) cos tx + B(t) sin tx]dt = ∫_{0}^{∞} e^{−at} cos txdt.

The integral on the right is the direct Laplace transform of the function t → cos(tx) (x fixed) evaluated at s = a. Therefore the integral evaluates to \frac{s}{s^{2}+x^{2}} with s = a, or \frac{a}{a^{2}+x^{2}} .

Method 1 is applied to get line 4 of the table. The given F (x) is odd, therefore A(t) = 0. The integration required is

B(t) = \frac{1}{π} ∫_{−∞}^{∞} F (x) sin txdx = \frac{1}{π} ∫_{−π}^{π} sin x sin txdx = \frac{2}{π} \frac{sin(πt)}{1 − t^{2}}.

8. (Fourier Sine and Cosine Integrals) Let the function F (x) be defined by

F (x) = \{ \begin{matrix}sin x, & for |x| ≤ π/2, \\ 0, & for |x| > π/2. \end{matrix}

Calculate the integrals

A(t) = \frac{1}{π} ∫_{−∞}^{∞} F (u) cos(tu)du,

B(t) = \frac{1}{π} ∫_{−∞}^{∞} 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.

Solution to 8. Since the function F is odd, A(t) = 0. Then the integrand in B is even and

B(t) = \frac{2}{π} ∫_{0}^{π/2} sin(x) sin(tx)dx = \frac{2t}{π} \frac{cos(πt/2)}{1 − t^{2}}.

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 = ±π/2, which is ±1/2.

Maple Notes on Fourier Integrals. The complex integral

A(t) − iB(t) = ∫_{−∞}^{∞} f(u) cos(tu)du − i ∫_{−∞}^{∞} f(u) sin(tu)du,

because of Euler’s Formula exp(−iθ) = cos(θ) − i sin(θ), equals the Complex Fourier Transform

T ≡ ∫_{−∞}^{∞} f(u) exp(−itu)du

Maple computes T in the 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. Example. Compute A(t) and B(t), the Fourier Cosine and Sine integrals, for the function

F (x) = \{ \begin{matrix}cos(x) & |x| < π \\ 0 & otherwise \end{matrix}

Solution: First write F in terms of Heaviside’s function H(t) as F (x) = cos(x)(H(x + P i) − H(x − P i)). To evaluate the improper integrals A and B we compute the complex Fourier transform T and use the identity A = Real(T), B = −Imaginary(T). In maple, we use the package fourier and the complex evaluation function evalc to simplify answers.

# 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

The answers: A(t) = 2t sin(πt)/(t^{2} − 1), B(t) = 0. Technical Issues. The technical assumptions of the 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 (−∞, ∞) 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 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.

Example. Does Fourier’s integral theorem apply to F (x) = \sqrt{x}H(x)? 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+) = ∞). 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 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.

Example. Assume F (x) = cos(x)(H(x) − H(x − π)). What is the value of the Fourier integral of F at x = 0? Solution: The answer is 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: 1/2!

Gibb’s Phenomenon. The function

F (x) ≡ H(x + 1) − H(x − 1) = \{ \begin{matrix}1 & −1 ≤ x < 1 \\ 0 & x < −1, x ≥ 1, \end{matrix}

where H is Heaviside’s function, has Fourier cosine and sine integrals A(t) = 2 sin(t)/(πt) and B(t) = 0 (see the table above). Let

\begin{matrix}F_{N}(x) = ∫_{0}^{N} \frac{2 sin t cos tx}{πt}dt \\ = \frac{1}{π} [Si(N + Nx) − Si(Nx − N)] . \end{matrix}

9. (Gibb’s Phenomenon) Evaluate G(x) and also the limit lim_{N→∞} F_{N}(x) by considering the cases x < −1, x = −1, −1 < x < 1, x = 1, 1 < x separately. Use the result of the computation to verify the integral equality of the Fourier integral theorem:

∫_{0}^{∞} [A(t) cos(tx) + B(t) sin(tx)] dt = G(x)

where G(x) = \frac{1}{2}[F (x + 0) + F (x − 0)].

Solution to 9. The idea of the verification for the identity

∫_{0}^{∞} [A(t) cos(tx) + B(t) sin(tx)] dt = G(x) ≡ \frac{1}{2}[F (x + 0) + F (x − 0)]

is to show by direct calculation that the left hand side (LHS) and the right hand side (RHS) are the same. This involves two computations. The first produces an explicit piecewise expression for G(x). The second is a limit evaluation, because the left side is an improper integral, defined in calculus to be the limit as N → ∞ of the partial sum F_{N}(x).

The RHS G(x) = (F (x + 0) + F (x − 0))/2 can be written out as a piecewise defined function over the five sets x < −1, x = −1, −1 < x < 1, x = 1, x > 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 ≥ 0, H(u) = 0 for u < 0. The answer:

G(x) = \frac{1}{2}[F (x + 0) + F (x − 0)] = \{ \begin{matrix}0 & x < −1, \\ \frac{1}{2} & x = −1, \\ 1 & −1 < x < 1, \\ \frac{1}{2} & x = 1, \\ 0 & x > 1. \end{matrix}

The LHS is the integral

∫_{0}^{∞} 2\frac{sin t}{πt} cos(tx) + (0) sin(tx)dt =

\frac{2}{π} ∫_{0}^{∞} \frac{sin t}{t} cos(tx)dt =

\frac{1}{π} ∫_{0}^{∞} \frac{sin(t + tx) + sin(t − tx)}{t}dt =

\frac{1}{π} lim_{N→∞} [∫_{0}^{(1+x)N} \frac{sin u}{u}du − ∫_{0}^{(x−1)N} \frac{sin w}{w} dw]

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)

∫_{0}^{∞} \frac{sin t}{t}dt = \frac{π}{2} .

The limits of integration x + xN and −N + xN tend to values 0, −∞ and ∞, 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}{π} ∫_{0}^{−∞} \frac{sin(u)}{u}du − \frac{1}{π} ∫_{0}^{∞} \frac{sin(u)}{u}du = −\frac{−π}{2π} − \frac{π}{2π} = 0.

The other four cases are treated similarly. The answer for the LHS:

∫_{0}^{∞} 2\frac{sin t}{πt} cos(tx) + (0) sin(tx)dt = lim_{N→∞} F_{N}(x) = \{ \begin{matrix}0 & x < −1, \\ \frac{1}{2} & x = −1, \\ 1 & −1 < x < 1, \\ \frac{1}{2} & x = 1, \\ 0 & x > 1. \end{matrix}

9a. (Graphical Approximation) Make 4 separate graphs of F_{N}(x) on (−∞, ∞), for N = 8, N = 16, N = 32, N = 64. Explain in a short paragraph the connection with Gibb’s Phenomenon.

Solution to 9a. Gibb’s Phenomenon for Fourier series is illustrated by the square wave example. In this example there is a characteristic cycle in the graph of the nth partial sum s_{n}(x) that is of fixed amplitude and moves towards the points of discontinuity of F as n → ∞. Elsewhere in the graph of s_{n}(x) the amplitudes decay to zero as n → ∞ 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 Fourier integrals is explained in a similar way. A paragraph of explanation might include graphical evidence of the phenomenon at the points x = ±1. The number of plot points can grossly affect the graphical evidence because the blip at x = ±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 blips get thinner but the height of each of the four blips stays around 0.13.

Maple notes. Some remarks will be given about the initial conversion of the formula using maple and the method for graphing on an infinite interval. The function Si(x) = ∫_{0}^{x} t^{−1} sin(t)dt can be used to reduce the partial Fourier integral

∫_{0}^{N} \frac{sin t cos tx}{t}dt

by the maple integration command

int((1/t)*sin(t)*cos(t*x),t=0..N);
#  1/2 Si(N + x N) - 1/2 Si(- N + x N)

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:

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]);

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.

Convergence of Fourier Series. The nth Fourier Partial Sum is

s_{n}(x) ≡ \frac{a_{0}}{2} + ∑_{k=1}^{n} (a_{k} cos(kx) + b_{k} sin(kx)) .

Fourier’s Convergence theorem says

lim_{n→∞} s_{n}(x) = \frac{lim_{t→x+} F (t) + lim_{t→x−} F (t)}{2}

= F (x)

at points x of continuity of F. The purpose of the problems below is to establish the 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. Your task is to establish the results by assuming the computations of the following theorem.

Theorem. Let F be piecewise continuous and 2π–periodic. Assume x given and define

G(u) ≡ (F (x + u) − F (x))/ sin(u/2).

The partial sums s_{n} of the Fourier series of F satisfy the following identities:

s_{n}(x) = \frac{1}{2π} ∫_{−π}^{π} F (x + u)\begin{matrix}sin & ((2n + 1) u/2) \\ sin & (u/2) \end{matrix}du

s_{n}(x) − F (x) = \frac{1}{2π} ∫_{−π}^{π} G(u) sin((2n + 1)u/2)du.

10. (Convergence of Fourier Series) Establish the trigonometric identity

D_{n}(u) ≡ \frac{sin ((2n + 1) u/2)}{2 sin (u/2)} = \frac{1}{2} + ∑_{k=1}^{n} cos(ku)

The fraction is called Dirichlet’s Kernel in the literature. Solution to 10. To prove the trigonometric identity, prove instead

\begin{matrix}sin & ((2n + 1) u/2) \\ sin & (u/2) \end{matrix} = −1 + 2 ∑_{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). 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)[\frac{1}{2} + ∑_{k=1}^{N} cos(ku)] = sin((2N + 1)u/2)

The left hand side (LHS) can be written as

LHS = sin(u/2) + ∑_{k=1}^{N} 2 sin(u/2) cos(ku)

By the trig identity

LHS = sin(u/2) + ∑_{k=1}^{N} sin((k + \frac{1}{2})u) − sin((k − \frac{1}{2})u)

or simply

LHS = sin(u/2) + [sin(3u/2) − sin(u/2)] + [sin(5u/2) − sin(3u/2)] + ⋯

The terms all cancel except one, which is the RHS of the above identity (calulus books study this idea as telescoping series).

11. (Convergence of Fourier Series) Integrate the formula of the previous problem on [−π, π] to obtain ∫_{−π}^{π} D_{n}(u)du = π. Solution to 11. To prove the integration identity, integrate over [−π, π] across the equation for Dirichlet’s kernel:

D_{n}(u) = \frac{1}{2} + ∑_{k=1}^{n} cos(ku).

The value of the integral ∫_{−π}^{π} cos(ku)du should be familiar.

12. (Convergence of Fourier Series) Assume x fixed and G(u) ≡ [F (x + u) − F (x)]/ sin(u/2) continuously differentiable on [−π, π]. ^{1} Establish using integration by parts:

\frac{1}{2π} ∫_{−π}^{π} G(u) sin((2n + 1)u/2)du

= \frac{\frac{1}{π} ∫_{−π}^{π} G′(u) cos((2n + 1)u/2)du}{2n + 1}

Solution to 12. To obtain the identity, begin by writing down the statement of integration by parts learned in calculus: ∫ udv = uv − ∫ vdu. You have to identify u, dv, du, v. Finally, the free term uv will be evaluated at −π and π. Justify why the free term is zero.

13. (Convergence of Fourier Series) Conclude from the preceding problems that in the convergence formula

F (x) = s_{n}(x) + \frac{M}{2n + 1} the value of M is the integral

M = \frac{−1}{π} ∫_{−π}^{π} G′(u) cos((2n + 1)u/2)du,

which is bounded by 2 max |G′(u)| on |u| ≤ π.

Solution to 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 θ| ≤ 1 and |∫_{a}^{b} f(x)dx| ≤ ∫_{a}^{b} |f(x)|dx.

Notes on the convergence theory.

Theorem. Let F be piecewise continuous and 2π–periodic. Assume x given and define

G(u) ≡ (F (x + u) − F (x))/ sin(u/2).

The partial sums s_{n} of the Fourier series of F satisfy the following identities:

s_{n}(x) = \frac{1}{2π} ∫_{−π}^{π} F (x + u) \begin{matrix}sin & ((2n + 1) u/2) \\ sin & (u/2) \end{matrix} du

s_{n}(x) − F (x) = \frac{1}{2π} ∫_{−π}^{π} G(u) sin((2n + 1)u/2)du.

b_{n}, as follows:

= \frac{1}{2π} ∫_{−π}^{π} f(x)dx + ∑_{k=1}^{n} \frac{1}{π} ∫_{−π}^{π} f(t) cos(kt)dt cos(kx)

+ ∑_{k=1}^{n} \frac{1}{π} ∫_{−π}^{π} f(t) sin(kt)dt sin(kx)

Euler’s formulas.

Each term has a common integral sign. Factors cos(kx) and sin(kx) are independent of t.

u = 0 with value G(0) = 2F′(x).

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 s_{n}(x) ≡ a_{0} + ∑_{k=1}^{n} (a_{k} cos(kx) + b_{k} sin(kx))  Definition.

= \frac{1}{2π} ∫_{−π}^{π} f(x) [1 + 2 ∑_{k=1}^{n} cos(kt) cos(kx) + sin(kt) sin(kx)] dt 1If F is continuously differentiable on [−π, π] 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

= \frac{1}{2π} ∫_{−π}^{π} f(x) \frac{sin( \frac{2k+1}{2} (t − x))}{sin \frac{t−x}{2}} dt Apply the identity of

the first problem below.

= \frac{1}{2π} ∫_{−π}^{π} f(x + u) \frac{sin( \frac{2k+1}{2} u)}{sin \frac{u}{2}} dt  {Use u = t − x, \atop du = dt.}

The above identity depends upon the formula ∫_{a}^{b} g(t)dt = ∫_{−π}^{π} g(t)dt, which is valid for any 2π-periodic function g(t) on an interval [a, b] of length 2π. 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

− ∫_{−π}^{π} \frac{1}{π} F (x)D_{n}(u)du = F (x) ∫_{−π}^{π} \frac{−1}{π} D_{n}(u)du, D_{n}(u) ≡ \begin{matrix}sin & \frac{(2n+1)u}{2} \\ sin & \frac{u}{2} \end{matrix}

because F (x) is constant relative to integration on u. By problem 11, the integral equals F (x) times −1, which proves the theorem.

Complex Method for the Dirichlet Kernel Formula. Substitute for cos(ku) the expression ℜ(exp(iku)).

Therefore, the sum is ∑_{k=0}^{n} r^{k} where r = exp(iu). Now apply the geometric series formula

∑_{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.

14. (Boundary Value Problems) Find all the eigenvalues λ, and the corresponding eigenfunctions x(t), for the boundary value problem

\begin{matrix}x′′ + λx = 0, 0 ≤ t ≤ 1, \\ x(0) = 0, x′(1) = 0. \end{matrix}

Solution to 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{λ}t) and therefore the eigenvalues are \sqrt{λ} = (2n+1)π/2 for n = 0, 1, 2, … with corresponding eigenfunctions sin((2n+1)πt/2).

14a. (Boundary Value Problems) Find all the eigenvalues λ, and the corresponding eigenfunctions x(t), for the boundary value problem

\begin{matrix}x′′ + λx = 0, −1 ≤ t ≤ 1, \\ x(−1) = 0, x′(1) = 0. \end{matrix}

Solution to 14a. All that has changed from the previous problem is the interval. The eigenpairs are given by \sqrt{λ_{n}} = (2n + 1)π/2, sin(\sqrt{λ_{n}}(t + 1)/2).

14b. (Sturm-Liouville Theory) Find all the eigenvalues and corresponding eigenfunctions of the Sturm-Liouville problem

\begin{matrix}y′′ − 12y′ + (36 + 4λ)y = 0, 0 ≤ x ≤ 5, \\ y(0) = y(5) = 0. \end{matrix}

Solution to 14b. The recipe contains the variable λ, so use Λ instead to avoid confusion. Then Λ^{2} − 12Λ + 36 + 4λ = 0 is the characteristic equation. The boundary conditions y(0) = y(5) = 0 lead to the trivial solution y = 0, if the roots Λ are real and distinct or real and equal (Cases 1 and 2 of the recipe). Only in Case 3 of the recipe is there a chance of producing a nonzero solution y. In that case, Λ is complex, obtained by solving (Λ − 6)^{2} + 4λ = 0, so λ > 0 and Λ = 6 + 2\sqrt{λ}i. Then y = c_{1}e^{−6x} cos(2\sqrt{λ}x) + c_{2}e^{−6x} sin(2\sqrt{λ}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} ≠ 0, provided sin(10\sqrt{λ}) = 0. This relations gives the eigenvalues 10\sqrt{λ_{n}} = nπ (n = 1, 2, 3, …) and eigenfunctions y_{n} = e^{−6x} sin(nπx/5).

15. (Orthogonal Functions) Show that set of functions below is orthogonal on the given interval and determine the corresponding orthonormal set.

sin ((n + \frac{1}{2})πt) , n = 0, 1, 2, ⋯, 0 ≤ t ≤ 1.

Solution to 15. The orthogonality issue is ∫_{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 ≠ 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 \{ϕ_{n}(t)\} on [a, b] into an orthonormal set, divide ϕ_{n} by ∥ϕ_{n}∥ ≡ \sqrt{∫_{a}^{b} |ϕ_{n}(t)|^{2}dt}.

In this case, ϕ_{n}(t) = sin ((n + \frac{1}{2})πt) and the division is by \sqrt{c_{n}} where c_{n} = ∫_{0}^{1} | sin ((n + \frac{1}{2})πt) |^{2}dt = 1/2. The orthonormal set is

\sqrt{2} sin ((n + \frac{1}{2})πt) , n = 0, 1, 2, ⋯, 0 ≤ t ≤ 1.

15a. (Orthogonal Functions) Show that set of functions below is orthogonal on the given interval and determine the corresponding orthonormal set.

cos (3nt) , n = 0, 1, 2, ⋯, |t| ≤ π.

Solution to 15a. The orthogonality question reduces to the integration of products cos(at) cos(bt) over [−π/2, π/2], where a = 3n and b = 3m (a ≠ 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π/6}.

The study will apply the recipe for second order linear differential equations with constant coefficients. Recipe. Given a constant coefficient differential equation ax′′ + bx′ + cx = 0 with characteristic equation aλ^{2} + bλ + c = 0, the general solution is given as x(t) = c_{1}x_{1}(t) + c_{2}x_{2}(t) where c_{1} and c_{2} are arbitrary constants and the basis solutions x_{1} and x_{2} are given as follows in terms of the roots λ_{1}, λ_{2} of the characteristic equation.

Case 1. Distinct Real Roots λ_{1} ≠ λ_{2}. Then

x_{1}(t) = e^{λ_{1}t}, x_{2}(t) = e^{λ_{2}t}.

Case 2. Double Roots λ_{1} = λ_{2}. Then

x_{1}(t) = e^{λ_{1}t}, x_{2}(t) = te^{λ_{1}t}.

Case 3. Complex Conjugate Roots λ_{1} = A + iB, λ_{2} = A − iB. Then

x_{1}(t) = e^{At} cos(Bt), x_{2}(t) = e^{At} sin(Bt).

Example 1. Solve x′′ + λx = 0 subject to the boundary conditions x(0) = 0, x(1) = 0. Show that the eigenvalues are (nπ)^{2} with corresponding eigenfunction sin(nπt) for integers n > 0.

Solution: The discriminant of the characteristic equation Λ^{2} + λ = 0 is D ≡ 0^{2} − 4(1)(λ). The sign of the discriminant selects one of the three cases of the classic recipe. Case 1. (λ > 0) In this case let λ = 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:

\{ \begin{matrix}0 = x(0) = c_{1} cos(0) + c_{2} sin(0) \\ 0 = x(1) = c_{1} cos(k) + c_{2} sin(k) \end{matrix} or \{ \begin{matrix}0 = c_{1}(1) + c_{2}(0) \\ 0 = c_{1} cos(k) + c_{2} sin(k) \end{matrix}

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 = | \begin{matrix}1 & 0 \\ cos(k) & sin(k) \end{matrix} | = sin(k)

The solution of this equation is k = nπ. The eigenfunction x(t) = c_{1} cos(kt) + c_{2} sin(kt) is determined by solving for c_{1} and c_{2}:

[ \begin{matrix}1 & 0 \\ cos(k) & sin(k) \end{matrix} ] [ \begin{matrix}c_{1} \\ c_{2} \end{matrix} ] = [ \begin{matrix}0 \\ 0 \end{matrix} ]

The answer is c_{1} = 0 and c_{2} arbitrary. In summary, k = nπ for integral n ≥ 1 and the eigenfunction is sin(nπt). Case 2. (λ = 0) The recipe implies that the general solution is x = c_{1} + c_{2}t. The boundary conditions put these demands on the solution:

\{ \begin{matrix}0 = x(0) = c_{1} + c_{2}(0) \\ 0 = x(1) = c_{1} + c_{2}(1) \end{matrix}

The matrix of coefficients has determinant | \begin{matrix}1 & 0 \\ 1 & 1 \end{matrix} | = 1.

By Cramer’s Rule, c_{1} = c_{2} = 0. In summary, λ = 0 is not an eigenvalue. Case 3. (λ < 0)

In this case let λ = −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:

\{ \begin{matrix}0 = x(0) = c_{1} exp(0) + c_{2} exp(0) \\ 0 = x(1) = c_{1} exp(k) + c_{2} exp(−k) \end{matrix} or \{ \begin{matrix}0 = c_{1}(1) + c_{2}(1) \\ 0 = c_{1}e^{k} + c_{2}e^{−k} \end{matrix}

The determinant of the matrix of coefficients is | \begin{matrix}1 & 1 \\ e^{k} & e^{−k} \end{matrix} | = 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, λ is not an eigenvalue for λ < 0.

Example 2. Solve x′′ + λx = 0 subject to the boundary conditions x(0) = 0, x(1) = 0 and show that the eigenvalues are (nπ)^{2} with corresponding eigenfunction sin(nπt) for integers n > 0.

Maple–Assisted Solution: The code below uses maple functions diff, dsolve, rhs, unapply, subs. It is helpful to invoke maple’s ‘?’ help function on each new maple operation (eg, ‘? rhs’); Case 1. (λ > 0)

# 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

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 maple to solve this equation. Answer: k = nπ, n = 1, 2, … (n = 0 can’t happen with k > 0). In summary, the eigenvalues are (nπ)^{2}, the eigenfunctions are sin(nπt), n = 1, 2, 3, . . .. Case 2. (λ = 0) The case λ = 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 ≡ 0. Conclusion: λ = 0 is not an eigenvalue.

# 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

Case 3. (λ < 0) The case λ < 0 can be decided by solving x′′ − k^{2}x = 0 subject to x(0) = 0, x(1) = 0. It is done below. Ultimately, the calculation shows x ≡ 0. Conclusion: No negative λ is an eigenvalue. Used below are some new maple functions: the linear algebra package linalg, the function genmatrix, which extracts the matrix of coeffients from a system of linear equations, numer which extracts the numerator of a quotient, det which computes the determinant of a matrix.

# 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

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 λ < 0.

Boundary Conditions with Derivatives. In problems with boundary conditions involving derivatives, for example, x′(−1) = 0, the coding in maple language changes, e.g., D(x)(-1)=0. Symbol Conflicts. The upper case D is a reserved symbol in maple. In particular, your maple program cannot use 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. Maple notes on orthogonality. The orthogonality relation ∫_{0}^{1} y_{n}(t)y_{m}(t)dt = 0, where y_{n}(t) = sin((2n+ 1)t/2), n ≠ m, can be coded in maple as follows.

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);

The integration will give answers involving terms of the form sin(kπ) where k is an integer. Without further use of maple, decide that the terms of the answer must be zero for n ≠ m. Additional notes on orthogonality. In problem 15a, the functions cos(2nt) on [−π, π] are claimed to be orthogonal functions. This means that

∫_{−π}^{π} cos(2pt) cos(2qt)dt = 0

for integers p and q selected from the list n = 0, 1, 2, 3, …, p ≠ 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: no conversions are required. That eliminates one source of error. Secondly, substitution of the limits of integration is done automatically in maple, eliminating the other source of error. Conclusion: maple is a good replacement for an integral table!

Example. Compute ∫_{−π}^{π} cos(2px) cos(2qx)dx for integers p ≠ q.

Solution: Employ maple as below.

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

The output of maple:

− \frac{sin(−2 pπ + 2 qπ)p + sin(−2 pπ + 2 qπ)q + sin(−2 pπ − 2 qπ)p − sin(−2 pπ − 2 qπ)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.

16. (Sturm-Liouville Theory) Show that the Sturm-Liouville problem

\begin{matrix}(x^{2}y′)′ + λx^{−2}y = 0, 1/2 ≤ x ≤ 1, \\ y(1/2) = y(1) = 0, \end{matrix}

has eigenvalues and eigenfunctions

λ_{n} = n^{2}π^{2}, y_{n} = sin(nπ/x), n = 1, 2, 3, · · · .

Then establish by direct integration the orthogonality property

∫_{1/2}^{1} y_{n}(x)y_{m}(x)x^{−2}dx = 0, n ≠ m.

Solution to 16. Write out the differential equation in standard form as x^{2}y′′ + 2xy′ + λx^{−2}y = 0, then clear fractions to obtain the new equation x^{4}y′′ + 2x^{3}y′ + λy = 0. This equation has variable coefficients and it cannot be solved directly from the Recipe for constant equations nor from the theory of Euler equations (it’s not an Euler equation). The information that sin(nπ/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′′ + λ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{λ}t) + c_{2} sin(\sqrt{λ}t) and the boundary condition u(2) = 0 implies sin(\sqrt{λ}) = 0, giving \sqrt{λ} = nπ. Then u(t) = sin(nπt). Transforming back gives the eigenfunction y(x) = sin(nπ/x) corresponding to eigenvalue λ_{n} = (nπ)^{2}. To verify the orthogonality relation, treat the problem as a calculus integration problem

∫_{1/2}^{1} sin(nπ/x) sin(mπ/x)x^{−2}dx

for n ≠ m. The secret to evaluation of the integral is the change of variables used in the differential equation itself: t = 1/x.

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. Example. Solve y′′ − 3y′ + (18 + 2λ)y = 0 in maple subject to the boundary conditions y(0) = 0, y(1) = 0. 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λ > 0. Define k^{2} = 63 + 8λ and rewrite the problem as y′′ − 3y′ + (18 + k^{2}/4 − 63/4))y = 0, y(0) = 0, y(1) = 0.

# 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

After y is inserted into the boundary conditions there results two equations in variables c, d with unknown k. The function 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 maple to solve this equation — by inspection k/2 is an integer multiple of π. Answer: k = 2nπ, n = 1, 2, 3, … subject to k^{2} = 63 + 8λ. Therefore, λ = ((2nπ)^{2} − 63)/8, n = 1, 2, 3, . . .. 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{2e^{\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π:

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]

Taking c = 0 and d = 1 in cu_{1} + du_{2} gives eigenfunction u_{2} (see formula above) for eigenvalue λ = (k^{2} − 63)/8 where k = 2π. Additional notes on problem 16. The equation (x^{2}y′)′ + λx^{−2}y = 0 is a Self-Adjoint Sturm-Liouville Problem. More generally, an equation of the form (py′)′ + (q + λr)y = 0 is an eigenvalue problem of self-adjoint form. The equation itself is coded in maple as follows:

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

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 maple:

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);

Well, maple is not very efficient for deciding that the answer is zero, but it is easily decided by inspection of the answer.

17. (Sturm-Liouville Theory, Numerical) Find accurate to 3 decimal places the first 5 eigenvalues and eigenfunctions of the Sturm-Liouville problem

\begin{matrix}y′′ + λy = 0, 0 ≤ x ≤ 1, \\ y(0) + y′(0) = 0, y(1) = 0. \end{matrix}

Plot the five eigenfunctions and count the zeros of each in the open interval 0 < x < 1.

Solution to 17. For λ = 0 the eigenfunction is 1 − x. For positive eigenvalues, solve tan(k) = k numerically where λ = k^{2}. The first of the infinitely many positive solutions k is about 4.49. The eigenfunctions for positive λ = k^{2} are k cos(kx) − sin(kx). There are no eigenvalues λ < 0. To graph the eigenfunctions see the maple example below.

Example. Solve for the first five eigenvalues and eigenfunctions of the Sturm-Liouville problem y′′ + λy = 0, y′(0) − y(0) = 0, y(1) = 0, accurate to 3 decimal places. Plot the five eigenfunctions on [0, 1] Solution: First we solve for the determining equations. The discriminant of the characteristic equation is −4λ and its sign determines the character of the solution. We consider three cases: λ = k^{2}, λ = 0, λ = −k^{2} for k > 0.

# 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

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 maple or by hand, use these ideas:

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

The cases λ = 0 and λ < 0 can be handled by repeating the above analysis using substitutions:

# 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

The calculation results in the inconsistent equation −2 = 0, therefore there is no eigenvalue for λ = 0. For λ < 0:

# 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

The determining equation for k is −ke^{−k} + e^{−k} − ke^{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 (−∞, ∞) gives evidence that k = 0 is the only root.

plot({exp(2*k),(1-k)/(1+k)},k=-infinity..infinity);
# only root is k=0

Therefore, det(C) ≠ 0 for k > 0 and there are no negative eigenvalues λ = −k^{2} < 0. The numeric values of the first five eigenvalues are determined by the following 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.

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

The answers for the example:

Eigenvalue Eigenfunction
4.115858365 2.028757838 cos(2.028757838x) + sin(2.028757838x)
24.13934203 4.913180439 cos(4.913180439x) + sin(4.913180439x)
63.65910654 7.978665712 cos(7.978665712x) + sin(7.978665712x)
122.8891618 11.08553841 cos(11.08553841x) + sin(11.08553841x)
201.8512584 14.20743673 cos(14.20743673x) + sin(14.20743673x)

To plot the five eigenfunctions use this maple code:

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`);

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. 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. Example. Solve numerically the equation tan(x) = 0.6x for the first five roots in (0, ∞). 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.

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!

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 maple equation solver called fsolve:

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

The process is repeated many times to obtain the five roots: 4.34, 7.639, 10.84, 14.01, 17.18.

18. (Periodic Eigenvalue Problem) Consider the eigenvalue problem

\begin{matrix}y′′ + λy = 0, −π ≤ x ≤ π, \\ y(−π) = y(π), y′(−π) = y′(π). \end{matrix}

Verify the table below and show by direct integration that the ‘eigenfunctions’ are pairwise orthogonal.

Values Eigenvalue

Eigenfunction

λ < 0 No eigenvalue No eigenfunction

λ = 0

λ = 0

y = 1 is the only eigenfunction

λ > 0 {λ = n^{2}, \atop n = 1, 2, 3, …}

Two linearly independent eigenfunctions y_{1} = cos(nx), y_{2} = sin(nx)

Solution to 18. This example from Sturm-Liouville theory is known as the Periodic Boundary Value Problem. The difficulty is in the analysis of the boundary conditions y(−π) = y(π), y′(−π) = y′(π), which are known as periodic boundary conditions because they are valid for any 2π-periodic differentiable function y. If attacked entirely by hand calculation, the three cases of the Recipe are treated through the discriminant of the characteristic equation, which is 0^{2} − 4λ. 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_{2}x, and since constants are periodic, this case produces eigenvalue λ = 0 and eigenfunction y = 1. In case 3, solutions are indeed periodic, y = c_{1} cos bx + c_{2} sin bx, where λ = b^{2}. The boundary conditions are satisfied provided b is a positive integer. This determines the eigenvalues λ = n^{2}, n = 1, 2, 3, . . .. The eigenfunctions use 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

∫_{−π}^{π} 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, …, with the requirement that f(x) ≠ g(x). There are four cases to consider, but all are handled by calculus methods. Maple Notes. The maple attack begins by simplifying the problem to exclude the discriminant cases λ < 0 and λ = 0. To be absolutely convincing the cases λ < 0 and λ = 0 should be treated separately by code similar to what appears below.

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

Unfortunately, Maple confuses issues by printing the answer as −4 + 4 cos(kπ)^{2} = 0 instead of the correct −4 + 4 cos^{2}(kπ) = 0. Don’t expect 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 maple code organizes the solving for c and d:

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]!

The kernel command in maple gives two answers, effectively meaning that [c, d] = c_{1}[1, 0] + c_{2}[0, 1] for 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 maple’s kernel command generate two linearly independent eigenfunctions for each eigenvalue. To prove orthogonality, verify ∫_{a}^{b} ϕ(x)ψ(x)dx = 0 where a = −π, b = π, ϕ = c_{1} cos(nx) + c_{2} sin(nx) and ψ is either 1 or else d_{1} cos(mx) + d_{2} sin(mx), n ≠ m. The organization of this computation in maple uses the functions int and simplify. It is necessary to examine the final answer and justify that it is zero (without using maple). Wave Equation and D’Alemberts Solution. Let the initial wave f be given by the formula

f(u) = \{ \begin{matrix}sin(πu) & −1 ≤ u ≤ 1, \\ 0 & otherwise. \end{matrix}

The variable x represents position along the string, t is time measured from 0 and y is the string deflection.

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 ≤ x ≤ 2. Start at t = 0 and end at some value of t where the wave is about 75% off the graph.

Solution to 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 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.

The answer in maple with two snapshots omitted:

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);

The Heaviside function is used automatically in maple V.2 and therefore the readlib command can be omitted.

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} [f(x − ct) + f(x + ct)]. For each of c = 1, c = 5 make four 2D “snapshots” on −5 ≤ x ≤ 5. Start at t = 0. Choose the other values of t to represent interesting snapshots of the solution.

Solution to 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.

The answer in maple with two snapshots omitted:

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);

An analysis of the snapshots can be done with maple V.2 command animate. See below for examples. The Plucked String. A plucked string lies along the interval 0 ≤ x ≤ 1. It has wave speed c = 1, initial velocity zero and initial shape

f(x) = \{ \begin{matrix}2x & 0 ≤ x ≤ 1/2, \\ 2(1 − x) & 1/2 ≤ x ≤ 1. \end{matrix}

Notation: x represents position along the string, t is time measured from 0 and y is the string deflection.

20. (Plucked String) The odd periodic extension of f of period 2 defined on (−∞, ∞) will be called h. The exact solution y can be written in terms of h:

y = \frac{1}{2} [h(x − t) + h(x + t)]

Graph the solution surface y on 0 ≤ x ≤ 1, 0 ≤ t ≤ 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.

Solution to 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) = signum(u)f(|u|)

where signum(0) = u/|u| for u ≠ 0 (it has values ±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 (−∞, ∞), of period 2, can be defined by a separate recursive function h(u). We recommend using the maple version described below or its equivalent.

Shown below is the maple code for plotting the true solution. Four of the snapshots are omitted.

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

Surprise: 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 maple from evaluating y(x, t). As can be seen from maple’s own sum function, quotes are sometimes necessary and unavoidable constructs. An error message like 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: animate (y, a..b, 0..1, frames=7), However, the latter leaves the quote problem unexplained.

20a. (Plucked String) Separation of variables leads to the series representation of the solution

y(x, t) = ∑_{n=1}^{∞} B_{n} sin(nπx) cos(nπ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 ≤ x ≤ 1, 0 ≤ t ≤ 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.

Solution to 20a. To compute the approximate solution z of five terms, use orthogonality of the functions sin(nπx) plus the identity y(x, 0) = f(x) to obtain

∫_{0}^{1} f(x) sin(nπx)dx = B_{n} ∫_{0}^{1}(sin(nπx))^{2}dx

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) = ∑_{n=1}^{5} B_{n} sin(nπx) cos(nπ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 maple.

In maple evaluate the integrals and solve for B_{n}:

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

The partial sum z is obtained and graphed as follows:

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

The reason the quotes are absent in the plots (e.g., 'z(x,0.0)') was explained earlier. For the explanation of quotes in the sum expression see maple help on sum.

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 (−∞, ∞).

Solution to 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 ≠ n). Then multiply the first by H_{m} and the second by H_{n}. Subtract the equations and integrate over (−∞, ∞). Finally, integrate by parts to obtain a formula for the integral of (n − m)H_{n}H_{m}e^{−x^{2}}. Only boundary terms remain from the integration by parts, and these terms are readily seen to be zero.

20c. (Sturm-Liouville Theory: Hermite Equation) Verify by direct integration, using formulas for the Hermite polynomials H_{10} and H_{7}, the orthogonality relation

∫_{−∞}^{∞} H_{10}(x)H_{7}(x)e^{−x^{2}}dx = 0.

Solution to 20c. Use integration by parts or apply maple’s int function.

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.

21. (Separation of Variables and a Vibrating String) Find the displacement y(x, t) of the vibrating string problem on domain t ≥ 0, 0 ≤ x ≤ 5π:

\frac{∂^{2}y}{∂t^{2}} = 4\frac{∂^{2}y}{∂x^{2}} ,

y(0, t) = y(5π, t) = 0, y_{t}(x, 0) = 0,

y(x, 0) = f(x) = 11 sin 2x − 4 sin 5x.

Verify that your solution is correct. Graph five snapshots of the solution for values t = 0 to t = 0.5.

Solution to 21. The eigenvalue problem is X′′ + λX = 0, X(0) = X(L) = 0 where L = 5π with eigenvalues λ_{n} = (nπ/L)^{2} and eigenfunctions sin(nπx/L). The associated problem in t is T′′ +c^{2}λT = 0 with general solution c_{1} cos(cnπt/L) + c_{2} sin(cnπt/L). By the boundary condition T′(0) = 0 it follows c_{2} = 0. Therefore the solution to the problem is

∑_{n=1}^{∞} B_{n} cos(nπct/L) sin(nπx/L)

where B_{n} = (2/L) ∫_{0}^{L} f(x) sin(nπ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.

Plotting in maple:

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);

The maple solution for y(x, t) amounts to evaluation of integrals and sums. The syntax is tricky so study the example and consult the maple help file on sum to explain why the single quotes appear.

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 );

To verify the solution in the case of the special f:

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);

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.

with(plots): animate(y(x,t),x=0..L,t=0..Pi/2,frames=50);

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 ≥ 0 and 0 ≤ x ≤ 3π. The equation is

\frac{∂T}{∂t} = 2\frac{∂^{2}T}{∂x^{2}} ,

and the boundary conditions are

\begin{matrix}T (x, 0) = f(x) = 10 sin 2x − 3 sin 5x, \\ T (0, t) = T (3π, t) = 0. \end{matrix}

After obtaining an explicit formula for T, graph 5 snapshots on 0 ≤ x ≤ 3π for t = 0, t = 0.1, t = 0.2, t = 0.3, t = 0.4.

Solution to 22. 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^{2}X = 0, X(0) = X(L) = 0 where L = 3π. The problem has eigenvalues (nπ/L)^{2} with eigenfunctions B_{n} sin(nπx/L). The problem for t results in factor exp(−(nπ/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 maple by the command

    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);

To do the snapshots it is useful to invoke the animation panel in maple using these commands:

    with(plots):
    animate( T(x,t),x=0..3*Pi,t=0..0.5,numpoints=100,frames=50);

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.

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

u_{t} = u_{xx} for 0 ≤ x ≤ 1, t ≥ 0,
u(0, t) = 0 and u_{x}(1, t) = 0 for t ≥ 0,
u(x, 0) = 9 sin(3πx/2) − 4 sin(11πx/2)

for 0 ≤ x ≤ 1. Plot 5 snapshots of the solution on 0 ≤ x ≤ 1 for t in the domain t = 0 to t = 0.08 that show interesting behavior of the solution.

Solution to 22a. This is a classical textbook application of the method of separation of variables. The method leads to the Sturm-Liouville problem X′′ + k^{2}X = 0, X(0) = 0, X′(1) = 0. The eigenvalues are ((n + 1/2)π)^{2} and the eigenfunctions are sin((n + 1/2)πx). The terms of the series solution ∑_{n=1}^{∞} B_{n} sin(λ_{n}x) exp(−λ_{n}t) that contribute are n = 2 and n = 6. The answer is u = 9 sin(3πx/2) exp(−(3π/2)^{2}t) − 4 sin(11πx/2) exp(−(11π/2)^{2}t).

To check the boundary conditions:

    u(0,t);
    u(x,0);
    subs(x=1,diff(u(x,t),x));

The answers are supposed to be 0, 9 sin(\frac{3πx}{2}) − 4 sin(\frac{11 π x}{2}), 0.

To do the snapshots it is useful to invoke the animation panel in maple using these commands:

         with(plots):
         animate( u(x,t),x=0..1,t=0..0.05,numpoints=100,frames=12);

The control panel allow selection of the appropriate snapshot frames for printing. The exact values of t should be reported.

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 ≤ x ≤ 2, 0 ≤ y ≤ 1.

23. (Laplace Equation in a Rectangle) Assume boundary conditions T (x, 0) = T (x, 1) = 0, 0 ≤ x ≤ 2, T (0, y) = f(y), T (2, y) = 0, 0 ≤ y ≤ 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πy) − 3 sin(8πy). Plot the exact solution T over the rectangular plate and mark the surface points where the temperature is locally at a maximum.

Solution to 23. The answer to the problem is

T = \frac{4 sinh(6 π (2 − x)) sin(6 π y)}{sinh(12 π)} − \frac{3 sinh(8 π (2 − x)) sin(8 π y)}{sinh(16 π)}

To check the 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);
T(x,0); T(x,1); T(2,y); T(0,y);

The answers should be 0, 0, 0, 4 sin(6 π y) − 3 sin(8 π y). The plot in 3D is done with the following command:

plot3d(T,0..2,0..1,axes=BOXED,style=PATCH,orientation=[-150,60]);

From this plot, analyze the surface for points of local maximum temperature. To adjust the perspective of the plot hold down Mouse Left Button and drag the wire frame into the new viewing position, then click Mouse Middle Button 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:

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);

The 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 contour from the 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 ≤ x ≤ 0.1 and 0 ≤ y ≤ 1 is the region of interest.

In maple there is the function contourplot in the plots package that is made for production of contour plots. It is somewhat different from the one described above. An example:

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]);

It is not claimed that extrema can be determined solely from a contour plot. In 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.

23a. (Laplace Equation in a Rectangle) Assume boundary conditions T (x, 0) = 0, T (x, 1) = g(x), 0 ≤ x ≤ 1, T (0, y) = T (1, y) = 0, 0 ≤ y ≤ 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.

Solution to 23a. Variables separate and produce two problems. The first is u′′ + k^{2}u = 0, u(0) = u(1) = 0 with eigenvalues (nπ)^{2} and eigenfunctions sin(nπx). The second problem v′′−k^{2}v = 0, v(0) = 0 has solution sinh(nπy). The temperature is

T = ∑_{n=1}^{∞} B_{n} sin(nπx) sinh(nπy).

To determine B_{n} requires a quick orthogonality calculation giving

B_{n} = (4 − 4 cos(nπ))/(n^{3}π^{3} sinh(nπ)).

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 maple; see below.

The implementation in maple:

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);

The 3D plot takes more than 30 seconds. 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 (Newton’s Law of Cooling) while the the right end x = 1 is held at zero temperature. The model for T (x, t) is

\frac{∂T}{∂t} − \frac{∂^{2}T}{∂x^{2}} = 0,

T_{x}(0, t) − T (0, t) = 0,

T (1, t) = 0,

T (x, 0) = f(x), 0 ≤ x ≤ 1, t ≥ 0.

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.

Solution to 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′′ + λy = 0, y′(0) − y(0) = 0, y(1) = 0 with eigenvalues λ = k^{2} > 0 and k is the solution of tan(k) + k = 0 for k > 0 (asymptotic to (2n − 1)π/2 as n → ∞). The eigenfunctions are k cos(kx) + sin(kx). The problem in t has solution exp(−k^{2}t).

The eigenvalue analysis can be carried out in maple. Presented here is the case for λ > 0; a complete solution would consider also λ = 0 and λ < 0.

# 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

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 maple command which gives a rough graph is

plot({tan(x),-x},x=0..4*Pi,y=-15..0,numpoints=100);

The answers can be generated by choosing intervals about the guesses, for example

fsolve(tan(x)+x=0,x,1.9..2.2);
# Answer: 2.03 approximately

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) = ∑_{n=1}^{∞} B_{n}(\sqrt{λ_{n}} cos(\sqrt{λ_{n}}x) + sin(\sqrt{λ_{n}}x)) exp(−λ_{n}t)

An approximate answer is

T (x, t) = 0.517 (2.03 cos(2.03x) + sin(2.03x)) e^{−4.12t}+
−0.013 (4.91 cos(4.91x) + sin(4.91x)) e^{−24.1t}+
0.002 (7.98 cos(7.98x) + sin(7.98x)) e^{−63.7t}

which is obtained by evaluating B_{n} from the equation

∫_{0}^{1} f(x)(\sqrt{λ_{n}} cos(\sqrt{λ_{n}}x) + sin(\sqrt{λ_{n}}x))dx = B_{n} ∫_{0}^{1}(\sqrt{λ_{n}} cos(\sqrt{λ_{n}}x) + sin(\sqrt{λ_{n}}x))^{2}dx

If you choose to use maple to help with the computations, then here are some useful ideas for the first three terms:

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

Some answers will print out unevaluated with sine and cosine terms. Use evalf to convert the answers to decimals.

25. (Elementary Power Series) Prove that \frac{1}{1 + x} = ∑_{n=0}^{∞}(−1)^{n}x^{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}.

Solution to 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) = ∑_{n=0}^{∞} \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.

Since f has a Taylor series expansion about x_{0} = 0, the coefficient f^{(n)}(0)/n! can be computed in maple as follows:

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!);

The formal series can be manipulated by maple commands:

sum('(-x)^n','n'=0..infinity);
y:=unapply(",x);
diff(y(x),x);
int(y(x),x);

The action of maple is to find closed-form answers for series problems. Therefore, the answers printed are related to the claimed identities but maple will not produce for you the expected proof.

26. (Regular-Point Series Solutions) The equation y′ = y + x^{3}, y(0) = 1 has a power series solution of the form y(x) = ∑_{n=0}^{∞} c_{n}x^{n}. Find a recursion relation for the coefficients {c_{n}} and determine the first five terms explicitly.

Solution to 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 Linear differential equations series package called ODEseries. The functions to be used:

Function  Package  Description
dsolve  maple library solve a differential equation
dsolve[series] maple library solve a differential equation

by the method of series

diffeqtorec ODEseries find recursion relation for a linear

differential equation

rectoproc ODEseries convert recursion to a function of n

It solves for c_{n} in ∑ c_{n}x^{n}.

sum  maple library sum of a series
plot  maple library plot a function

The basic mathematical algorithm for solving a linear differential equation by the method of series is:

(a) Assume y = ∑_{n=0}^{∞} c_{n}x^{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.

(b) Collect terms on powers of x.

(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.

(d) Solve the recursion. Substitute the answer back into the series for y.

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.

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′ = ∑_{n=1}^{∞} nc_{n}x^{n−1}

This equation and the original for y are inserted into the differential equation to get

∑_{n=1}^{∞} nc_{n}x^{n−1} − ∑_{n=0}^{∞} c_{n}x^{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} + ∑_{n=4}^{∞} nc_{n}x^{n−1} − ∑_{n=3}^{∞} c_{n}x^{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} + ∑_{k=3}^{∞}(k + 1)c_{k+1}x^{k} − ∑_{n=3}^{∞} c_{n}x^{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 ≥ 3)

The final recursion is solved to give c_{k+1} = 6c_{3}/(k + 1)! for k ≥ 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 = ∑_{n=0}^{∞} c_{n}x^{n} = 1 + x + x^{3}/3 + x^{4}/12 + x^{5}/60 + ⋯

The solution in 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 diffeqtorec. Step (d) is done by

rectoproc and sum.

It is emphasized that your intuition about series solutions originates in your own hand calculations. Use maple to check your answer, to shorten tedious steps and to correct calculation errors. In particular, use maple sparingly on series problems.

# 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`);

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}.

Solution to 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λ^{2} + bλ + c = 0 where a, b, c are the coefficients of the second order linear differential equation. However, for this equation b = 2x, which is not constant, therefore the recipe fails to apply!

There is no elementary closed–form solution! Here is maple code for solving the differential equation in closed–form:

eq:=diff(y(x),x,x) + 2*x*diff(y(x),x) + y(x)=0:
dsolve(de,y(x));

There is no answer printed by maple. The interpretation of the maple answer is this: no solution was found! The series method applies to find a basis of solutions: 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.

Solution: Assume y = ∑_{n=0}^{∞} c_{n}x^{n}, then y′ = ∑_{n=0}^{∞} nc_{n}x^{n−1}, y′′ = ∑_{n=0}^{∞} n(n − 1)c_{n}x^{n−2}. Inserting into the differential equation gives

∑_{n=0}^{∞} n(n − 1)c_{n}x^{n−2} + 3x ∑_{n=0}^{∞} nc_{n}x^{n−1} + ∑_{n=0}^{∞} c_{n}x^{n} = 0

After re-indexing and term collection the series sums take this form:

∑_{n=0}^{∞}(k + 2)(k + 1)c_{k+2}x^{k} + ∑_{n=0}^{∞} 3nc_{n}x^{n} + ∑_{n=0}^{∞} c_{n}x^{n} = 0

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}

Checking the Answer in Maple. A check in maple reveals that correct answers were obtained in the major steps of the solution:

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

It is possible to check only the steps from the recursion to the end, as follows:

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);

Shown below is a method for obtaining the answer with dsolve. The method obtains only the final answer, and not the steps between.

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

26b. (Airy’s Equation) Compute the recursion relations and the first five terms of the series solution y(x) = ∑_{n=0}^{∞} c_{n}x^{n} for Airy’s differential equation y′′ − xy = 0, y(0) = 1, y′(0) = 0.

Solution to 26b. It is known that Airy’s differential equation has solution expressed in terms of Bessel functions. If you apply maple’s function 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

Some maple source appears below to show how to check the series answer against the hand computation.

# 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

The ODEseries package can be used to check the recursions and also the series solution.

# 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

27. (Euler Differential Equation Recipe) An Euler equation ax^{2}y′′ + bxy′ + cy = 0 can be solved by a recipe similar to the one used for an equation with constant coefficients:

Solve for the roots λ in the Euler auxiliary equation aλ(λ − 1) + bλ + 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^{2}y′′ + bxy′ + cy = 0.

Apply the recipe to solve

(a) x^{2}y′′ + 4xy′ + 2y = 0, (b) x^{2}y′′ + 3xy′ + y = 0,
(c) x^{2}y′′ + 3xy′ + 2y = 0.

Solution to 27. The basic ideas can be seen in the example to follow.

Example. Solve the Euler equation x^{2}y′′ + 4xy′ + 2y = 0. Solution: The Euler auxiliary equation is λ(λ − 1) + 4λ + 2 = 0 which factors as (λ + 1)(λ + 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}).

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^{2}y′′ + bxy′ + cy = 0 and obtain the constant equation a(z′′ − z′) + bz′ + cz = 0 whose characteristic equation is aλ(λ − 1) + bλ + c = 0. This calculation justifies the recipe for Euler equations claimed above.

Solution to 27a. The chain from calculus is applied repeatedly to compute dz/dt and d^{2}z/dt^{2}:

\frac{dz}{dt} = \frac{dy}{dx} \frac{dx}{dt} = y′(x)e^{t}

\frac{d^{2}z}{dt^{2}} = \frac{d}{dt} (y′(x)e^{t}) = y′(x)\frac{d}{dt} (e^{t}) + \frac{d}{dt} (y′(x)) 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} (y′(x)) = \frac{d}{dx} (y′(x)) \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^{2}y′′(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.

28. (Regular Singular Points and Series Solutions) Let ax^{2}y′′+bxy′+cy = 0 be a second order Euler differential equation with arbitrary constant coefficients a, b, c, and a ≠ 0. Prove that Euler’s equation has a regular singular point at x = 0. Observe that the classic Euler Auxiliary Equation aλ(λ − 1) + bλ + c = 0 is identical to the Indicial Equation of the Frobenius Theory.

Solution to 28. The proof amounts to arranging the equation in the form x^{2}p(x)y′′+xq(x)y′+r(x)y = 0 where p, q, r are power series about x = 0 and p(0) ≠ 0.

The classic Euler auxiliary equation is found formally by substitution of exp(λx) into the differential equation. The Indicial equation of Frobenius theory is p(0)λ(λ − 1) + q(0)λ + r(0) = 0.

28a. (Perturbations of Euler Equations) Let Ax^{2}y′′ + B(x)xy′ + C(x)y = 0 be a second order differential equation with regular singular point at x = 0, A ≠ 0 a constant and B(x) and C(x) power series convergent for |x| < R. Prove that the equation can be rearranged into the equivalent form

ax^{2}y′′ + bxy′ + cy = f(x)y′ + g(x)y

where a = A ≠ 0, b and c are constants and f(x) and g(x) are power series with f(0) = f′(0) = g(0) = 0.

In particular, every equation with regular singular point is a perturbation of an Euler differential equation whose auxiliary (characteristic) equation is identical to the Indicial Equation. The solutions in the Frobenius theory can be viewed as power series perturbations of Euler equation solutions.

Solution to 28a. The method of the conversion can be extracted from the example below.

Example. Show that 2x^{2}y′′ + (3x + 6x^{2})y′ + (1 + exp(x))y = 0 has a regular singular point at x = 0 and write the equation in the form ax^{2}y′′ + bxy′ + cy = f(x)y′ + g(x)y. Solution: Let p(x) = 2x^{2}/x^{2}, q(x) = (3x + 6x^{2})/x, r(x) = (1 + exp(x))/1. Then p, q, r are power series at x = 0 and the given equation has a regular singular point at x = 0. To obtain the form claimed, write r(x) = 2 + x + x^{2}/2! + · · ·. This is not necessary for p and q since they are finite power series already. Insert the power series into the equation and move all but the leading terms to the right side of the equation:

2x^{2}y′′ + (3x)y′ + (2)y = 0 − 6x^{2}y′ − (x + x^{2}/2! + x^{3}/3! + ⋯)y

This form matches ax^{2}y′′+bxy′+cy = f(x)y′+g(x)y. The series f(x) = −6x^{2} and g(x) = −(x+x^{2}/2!+x^{3}/3!+ ⋯) satisfy the requirement f(0) = f′(0) = g(0) = 0. In particular, the Euler part is 2x^{2}y′′ + (3x)y′ + (2)y = 0 and the Euler auxiliary equation for the Euler part is exactly the Indicial equation, 2λ(λ − 1) + 3λ + 2 = 0.

29. (Bessel Functions J_{n}) Graph the Bessel functions J_{0}, J_{1}, J_{2} on the interval [0, 20]. Display the number of zeros of each function in [0, 20] in a brief table. In a second table display the equation which most closely approximates each function near x = 0, chosen from the list y = 1, y = x, y = x^{2}, y = x^{3}, y = x^{4}.

Solution to 29. It is expected that you will use a graphing program on this problem. Suitable ones are maple, gnuplot and microcomputer plotting programs.

The simplest method is to use maple because the Bessel functions are intrinsic to the library. The names:

Mathematical name Maple name

                                  J_{0}(x)  BesselJ(0,x)
                                  J_{1}(x)  BesselJ(1,x)
                                  J_{2}(x)  BesselJ(2,x)

A usual plot command for maple is

                         plot(BesselJ(1,x),x=0..20,y=-1..1);

The zeros of J_{1}(x) can be read off the left upper corner of the plot by clicking the left mouse button at a zero. The answers: 0.0, 3.8, 7.0, 10.2, 13.3, 16.5, 19.7. There are six zeros in the open interval (0, 20).

To understand the shape of the curve y = J_{1}(x) near x = 0, shorten the plot domain to [0, 0.1]. Plot also the sample curves y = 1, y = x, y = x^{2}, y = x^{3}, y = x^{4} on [0, 0.1]. Compare shapes to select the sample curve that matches J_{1}(x) near x = 0.

29a. (Bessel Functions Y_{n}) Graph the Bessel functions of the second kind Y_{1}, Y_{2}, Y_{3} on the interval [0.5, 20]. Display the number of zeros of each function in [0.5, 20] in a table. In a second table, display one equation of the form y = ax^{−n} which most closely approximates each function near x = 0, n = 0, …, 5. The function Y_{0} is approximately 2 ln(x)/π near x = 0 and therefore this analysis does not apply to Y_{0}.

Solution to 29a. The best method is maple because the Bessel functions of the Second kind are intrinsic to the library. The names:

Mathematical name Maple name

                                  Y_{1}(x)  BesselY(1,x)
                                  Y_{2}(x)  BesselY(2,x)
                                  Y_{3}(x)  BesselY(3,x)

A normal plot command for maple is

                            plot(BesselY(1,x),x=0.5..20);

The zeros of Y_{1}(x) can be read off the left upper corner of the plot by clicking the left mouse button at a zero. The answers: 2.5, 5.5, 8.6, 11.8, 14.9, 18.2. There are six zeros in the open interval (0.5, 20).

To decide upon y = ax^{−n} perform experiments like this:

                       plot(BesselY(2,x)*x*x,x=0..0.1,y=-2..0)

From the graph it is apparent that in a graphical sense Y_{2}(x)x^{2} = a for some constant a. The value of a is about −1.275. So the equation that most closely approximates Y_{2}(x) near x = 0 is y = −1.275x^{−2}.

30. (Integrals of Bessel Functions) Evaluate the integral ∫_{0}^{4} x^{4}J_{1}(x)dx in two ways: numerically, and by reduction to evaluations of Bessel functions.

Solution to 30. This problem requires numerical integration, best carried out in maple or matlab.

The numerical integration is done in maple:

                       evalf(Int(BesselJ(1,x)*x^4,x=0..4))

The reason for the upper case “I” on “Int” is to return the integral without trying to apply integration rules. The function evalf causes the integral to be evaluated numerically. The answer: 38.15485668. The method for evaluating the integral by Bessel function identities depends upon the formulas

(x^{2}J_{2}(x))′ = x^{2}J_{1}(x), (x^{3}J_{3}(x))′ = x^{3}J_{2}(x)

and integration by parts. The final evaluation can be done in maple:

            f:=x -> x^4*BesselJ(2,x) - 2*x^3*BesselJ(3,x):  evalf(f(4)-f(0));

30a. (Bessel Functions Integral Formulas) Prove the formula

∫_{0}^{x} J_{1}(x) sin(x)dx = xJ_{1}(x) sin(x)
+(x cos(x) − sin(x))J_{0}(x)

Solution to 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.

31. (Legendre Polynomials) Calculate the first six Legendre polynomials and display them in a table.

Solution to 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 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}, … are found similarly.

The maple package orthopoly contains routines for all the major sets of orthogonal polynomials. To display the Legendre polynomials:

with(orthopoly):
P(0,x);P(1,x);P(2,x);P(3,x);P(4,x);P(5,x);

The Legendre polynomial of degree 5 is:

\frac{63 x^{5}}{8} − \frac{35 x^{3}}{4} + \frac{15 x}{8}

31a. (Graphs of Legendre Polynomials) Present five separate graphs showing the first five Legendre polynomials on [−1, 1], P_{0} to P_{4}.

Solution to 31a. All plots will be done in maple.

The single plots can be done one at a time using the maple command plot, e.g., to plot P_{2}, use the command

                           plot((3*x^2-1)/2,x=-1..1);

Once the five plots are on the screen, shrink them to size, and print the screen. The maple command plot is able to plot several functions on a single graph. The syntax is 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:

                          plot({P(2,x),P(3,x)},x=-1..1);

32. (Roots of Legendre Polynomials) Produce a table which displays the roots of each of the first six (6) Legendre polynomials.

Solution to 32. The orthopoly package in maple together with the solve routine will find roots of Legendre polynomials. There are numerical routines for the roots, a classic reference being Numerical Recipes. The maple commands for the 5th Legendre polynomial:

    with(orthopoly):
    evalf(solve(P(5,x)=0,x));

Without the evalf function, radicals will print, which is not too useful in a table.

32a. (Legendre Polynomial Integral formula) Let P_{n} be the nth Legendre polynomial given by Rodrigues’ formula, n = 1, 2, . . .. Assume that ∫_{−1}^{1} P_{n}^{2}(x)dx = 2/(2n + 1) and ∫_{−1}^{1} P_{n}(x)P_{m}(x)dx = 0 for n ≠ m. Prove that

∫_{−1}^{1} xP_{n}(x)P_{n−1}(x)dx = \frac{2n}{4n^{2} − 1}

Solution to 32a. The key identity to be used in the solution is the 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 ∫_{−1}^{1} P_{n}(x)P_{m}(x)dx = 0 for n ≠ 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 ∫_{−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 Bonnet recursion mentioned above.

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

T_{xx} + T_{yy} = 0, 0 ≤ x ≤ 2π, 0 ≤ y ≤ π,
T (x, 0) = 0, T (x, π) = 0, 0 ≤ x ≤ 2π,
T (0, y) = f(y), T (2π, y) = 0, 0 ≤ y ≤ π.

Show how the method of separation of variables can be used to construct T = T (x, y) for any function f(y). Then find the solution for the special case that f(y) = y for 0 ≤ y ≤ π.

Solution to 33. The boundary value problem for the y–variable after separation of variables is Y′′ + λY = 0, Y (0) = Y (π) = 0. The problem for the x–variable is X′′ − λX = 0, X(2π) = 0.

The series solution is T = ∑_{n=1}^{∞} c_{n} sinh(n(x − 2π)) sin(ny) where c_{n} is given by the formula

c_{n} sinh(−2nπ) = \frac{2}{π} ∫_{0}^{π} f(y) sin(ny)dy

For the function f(y) = y the answer is c_{n} = \frac{2(−1)^{n}}{n sinh(2nπ)}.

The use of maple in this problem is to check the answer and perform some standard integration.

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

34. (Steady-State Temperature in a Sector) The steady state temperature T = T (r, θ) in the sector 0 ≤ r ≤ 1, 0 ≤ θ ≤ π/3 (polar coordinates) satisfies the Laplace equation and boundary conditions

T_{rr} + \frac{1}{r}T_{r} + \frac{1}{r^{2}} T_{θθ} = 0, 0 < r ≤ 1, 0 ≤ θ ≤ \frac{π}{3} ,

T (r, 0) = 0, T (r, π/3) = 0, 0 ≤ r ≤ 1,

T (1, θ) = f(θ), 0 ≤ θ ≤ π/3.

Show how the method of separation of variables can be used to construct T = T (r, θ), then find the solution for the case that f(θ) = 5 sin(6θ) − 11 sin(27θ).

Solution to 34. The temperature T is separated as the product R(r)Θ(θ). The differential equations are r^{2}R′′ + rR′ − λR = 0 and Θ′′ + λΘ = 0.

The boundary conditions for the θ–variable are Θ(0) = Θ(π/3) = 0. This Sturm-Liouville problem has solution λ = 9n^{2}, Θ(θ) = sin(3nθ).

The Euler auxiliary equation for r^{2}R′′ + rR′ − 9n^{2}R = 0 has two real roots ±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 = ∑_{n=1}^{∞} c_{n}r^{3n} sin(3nθ)

The coefficients are determined by setting r = 1 and T (r, θ) = f(θ) 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θ) on [0, π/3].

The answer: c_{n} = \frac{6}{π} ∫_{0}^{π/3} f(θ) sin(3nθ)dθ.

The evaluation of the coefficients c_{n} for the function f(θ) = 5 sin(6θ) − 11 sin(27θ) can be performed in maple or evaluated with an integral table. Here is an example for a slightly different function f:

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

The answer for this f: c_{n} = 0 except for n = 2 and n = 11, c_{2} = 2, c_{11} = −7.

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

∑_{n=1}^{∞} (a_{n} cos (\frac{s_{n}ct}{R}) + b_{n} sin (\frac{s_{n}ct}{R})) J_{0} (\frac{s_{n}r}{R} )

where {s_{n}}_{n=1}^{∞} is the sequence of positive zeros of J_{0}(x). It is known that the Bessel functions {J_{0}(s_{n}r/R)} are orthogonal on [0, R].

Take c = 1 and R = 1. Assume initial values

u(r, 0) = 0, u_{t}(r, 0) = f(r), 0 ≤ r ≤ 1.

Find a_{n} and b_{n} in terms of integrals of J_{0} and the sequence of roots {s_{n}}_{n=1}^{∞}. 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_{2}t)}{s_{2}}J_{0}(s_{2}r) − \frac{1}{2} \frac{sin(s_{4}t)}{s_{4}}J_{0}(s_{4}r)

Finally, produce a snapshot graph of u at t = 0.2 and discuss its relationship to the drum vibration.

Solution to 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_{n}r). 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}.

To find the roots s_{n} in maple:

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]

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:

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]

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(θ), y = r sin(θ), 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 maple is recommended only for high speed work stations.

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);

Printing of the plot can be done with the unix shell command 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:

# 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);

36. (Steady-State Temperature in a Sphere) The steady state temperature T = T (r, ϕ) in the sphere 0 ≤ r ≤ 1, 0 ≤ ϕ ≤ π (spherical coordinates) satisfies the Laplace equation

T_{rr} + \frac{2}{r}T_{r} + \frac{1}{r^{2} sin ϕ}(sin ϕT_{ϕ})_{ϕ} = 0

for 0 < r ≤ 1, 0 < ϕ < π, as well as the boundary condition

T (1, ϕ) = f(ϕ), for 0 ≤ ϕ ≤ π.

Let P_{n}(x) be the nth Legendre polynomial. The product solutions of the Laplace equation are known to be

T_{n}(r, ϕ) = r^{n}P_{n}(cos ϕ)

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, ϕ) for any function f(ϕ). Then find T = T (r, ϕ) for the special case that f(ϕ) = 4 cos^{3} ϕ − 3 cos ϕ.

Solution to 36. The solution is on slides. Maple notes below.

The solution begins by verifying that the given product solution is indeed a solution of the Laplace equation.

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

The final result is the differential equation

n^{2}p(x) + np(x) − 2xp′(x) + p′′(x) − x^{2}p′′(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(ϕ) = T (1, ϕ) = ∑_{n=1}^{∞} c_{n}P_{n}(cos(ϕ))

The basic idea is to multiply the above equation by P_{n}(cos(ϕ)) sin(ϕ) and integrate over [0, π] to determine c_{n}. If this is carried out for f(ϕ) = 4 cos^{3} ϕ − 3 cos ϕ, then

c_{1} = 3/5, c_{3} = 2/5, c_{n} = 0, otherwise

The process of evaluation of the c_{n}’s involves knowledge of Legendre polynomials. In particular, by a change of variables,

∫_{0}^{π} cos^{3}(ϕ)P_{n}(cos(ϕ)) sin(ϕ)dϕ = ∫_{−1}^{1} −x^{3}P_{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

∫_{0}^{π} cos^{3}(ϕ)P_{n}(cos(ϕ)) sin(ϕ)dϕ = (−2/5) ∫_{−1}^{1} P_{3}(x)P_{n}(x)dx + (−3/5) ∫_{−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.

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 a method for printing several graphs on a single sheet of paper. The idea is to display the graphs on the screen, resize and move each until the screen is correct, then print the screen. Up to 8 graphs can be displayed by this mechanism. All will appear on a single sheet of paper. The resize and move commands are available on mouse menu 1. Press Mouse Left Button 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 b129lab1 laser printer by the following unix shell command. The command is to be issued in a console or terminal window.

     alias paw "xdpr -root -device ps -portrait -Pb129lab1"
     paw

One issue of the alias command above installs the command paw and thereafter, until logout, you may use paw to print the screen. The command name paw stands for Print All Windows (sometimes called print screen). Unix wizards can install the command paw permanently by copying the alias line above to the unix resource file .cshrc in their root directory. The file .cshrc is read at login and therefore the command paw will always be available. Another possibility is to copy the following lines into the file $HOME/.twmrc (the first line should already be in the source!):

     "Print Window"          !"xdpr -device ps &"
     "Print Screen"          !"xdpr -root -device ps -portrait &"

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:

                            mv  $HOME/.twmrc  $HOME/.twmrc.old
                          cp  -p  /usr/skel/.twmrc  $HOME/.twmrc

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.

Take-home I [1,2,3,4]: Tuesday of week 2

Take-home II [5,6,7,8]: Tuesday of week 3

Take-home III [9,10,11,12] Tuesday of week 4

Take-home IV [14,15,16,17]: Tuesday of week 5

Take-home V [18,19,20,21,22]: Tuesday of week 6

Take-home VI [23,24,25,26,27] Tuesday of week 7

Take-home VII [28,29,30,31,32]: Tuesday of week 8

Take-home VIII [33,34,35,36]: Tuesday of week 9