Chapter 4
Separation of Variables and
the Sturm-Liouville Problem

4.2 Orthogonality of Functions

Example 4.1

Integrate can be used to evaluate symbolically the overlap integrals of the functions Sin[k_n x], which were determined in the preceding section to be solutions of the wave equation with the boundary conditions (4.11).  With k_n=n(Pi)/L, where n is an integer, we obtain

<<Algebra`Trigonometry`

(Integrate[TrigToComplex[Sin[(n Pi x)/L]*
     Sin[(m Pi x)/L]],x]/.x->L)-
(Integrate[TrigToComplex[Sin[(n Pi x)/L]*
     Sin[(m Pi x)/L]],x]/.x->0)

Apart[ComplexToTrig[%]]

If m≠n, both terms in this expression vanish, which shows that Sin[k_n x] and Sin[k_m x] are orthogonal.  If m=n, the corresponding commands yield

Integrate[TrigToComplex[(Sin[(n Pi x)/L])^2],{x,0,L}]

Simplify[ComplexToTrig[%]]

The second term vanishes, so the integral takes the value L/2, independent of n.  Thus, the orthogonality relation of the sine functions over the interval (0,L) can be summarized as in Equation (4.26).

4.3 Fourier Series

Example 4.2 (Figures 4.1 and 4.2)

To illustrate some of the points just discussed, we consider the Fourier series representation of two different types of function.  Consider first the function in Equation (4.33). This function is continuous everywhere, but has discontinuities in its first derivative at x+Pi and x=0.  To generate the Fourier series of this function within Mathematica, we first construct the function in (4.33),

f[x_]:=If[x<=0,(x+Pi)/Pi,(Pi-x)/Pi]

The Fourier series is now constructed by first defining the expansion coefficients as in (4.28) and (4.29).  When we do so, we find that the b[n] all vanish for reasons that will be explained below.  Thus, we need only consider the a[n], which are calculated using Integrate:

a[n_]:=(1/Pi)Integrate[((Pi+t)/Pi) Cos[n t],
     {t,-Pi,0}]+(1/Pi)Integrate[((Pi-t)/Pi) Cos[n t],
     {t,0,Pi}]

The partial sums of the Fourier series are now constructed using Sum.  The
distinction between a[0] and the a[n] for n≠0 as shown in (4.28) and (4.29) can be taken into account by defining new coefficients A[k] using If:

A[k_]:=If[k>0,a[k],(1/2)a[0]]

F[x_,n_]:=Sum[A[k]Cos[k x],{k,0,n}]

The partial sums of the Fourier series can now be generated up to any desired
order.  For example,

f1[x_]=F[x,1]

f5[x_]=F[x,5]

Plot[{f1[x],f[x]},{x,-3,3},PlotRange->{0,1},
     AxesLabel->{x,f}]

Plot[{f5[x],f[x]},{x,-3,3},PlotRange->{0,1},
     AxesLabel->{x,f}]

The Fourier representations obtained from F[x,1] and F[x,5] are  shown
in Figure 4.1.  Even the first term reproduces the profile of the function quite
closely.  The only marked differences between the partial sum and the function are
at the points of discontinuity of the first derivative.  Such features cannot be
reproduced with a finite series of smooth functions, since any finite sum of functions with continuous first derivatives will also have a continuous first
derivative.

Consider now the function in Equation (4.34),which has discontinuities at x=Pi/8 and x=7Pi/8.  The Fourier series for this function is generated by following similar steps to those used for obtaining the series in (4.33) and the partial sums including the first 5, the first 10 and the first 20 terms are shown in Figure 4.2.

f[x_]:=If[Abs[x-Pi/2]<3Pi/8,1,0]

a[k_]:=
     ((1/Pi)Integrate[Cos[k t],t]/.t->7Pi/8)-
     ((1/Pi)Integrate[Cos[k t],t]/.t->Pi/8);
A[k_]:=If[k>0,a[k],(1/2)a[0]];
B[k_]:=
     ((1/Pi)Integrate[Sin[k t],t]/.t->7Pi/8)-
     ((1/Pi)Integrate[Sin[k t],t]/.t->Pi/8);

F[x_,n_]:=A[0]+Sum[A[k]Cos[k x]+B[k]Sin[k x],{k,1,n}]

f5[x_]=F[x,5];

Plot[{f5[x],f[x]},{x,-Pi,Pi},PlotRange->{-0.1,1.2},
     AxesLabel->{x,f}]

f10[x_]=F[x,10];

Plot[{f10[x],f[x]},{x,-Pi,Pi},PlotRange->{-0.1,1.2},
     AxesLabel->{x,f}]

f20[x_]=F[x,20];

Plot[{f20[x],f[x]},{x,-Pi,Pi},PlotRange->{-0.1,1.2},
     AxesLabel->{x,f}]

Although the general profile of the function is reproduced even by including only
first five terms, appreciable differences are evident between the Fourier representation and the function even after the first twenty terms have been
included in the partial sum.  The slow convergence is due to the rapid variation
of the function near the discontinuities, which requires the higher-order trigonometric functions to be included in the partial sum because of their short periods.

However, including more terms in the Fourier expansion of this function does not lead to an improved representation in all respects.  While the value of the function is better reproduced at points near the discontinuities, the series representation at the discontinuities shows marked deviations from the profile of the function regardless of how many terms are included in the Fourier series.  This effect, which is known as the Gibbs phenomenon and is a general feature of Fourier series of discontinuous functions, is discussed in Problem 8.

4.4 Fourier Series Solutions of Differential Equations

Example 4.4 (Figure 4.4)

The form of the solution in (4.64) provides some interesting insights into the way solutions of the diffusion equation develop in time.  The important point is that the rate of decay of the nth term in the sum in (4.64) is proportional to n^2.  As we have seen in Example 4.2, the higher-order terms in a Fourier series represent the sharp features of a function.  Thus, this solution shows that the sharpest features of the initial function have the fastest rate of decay. In particular, any sharp features such as discontinuities in the initial function are quickly smoothed out.

To see this behavior in a specific example consider the solution (4.64) with
the initial function f[x] given by a form similar to that in (4.33), i.e. f[x]=x if 0<x<1/2 and f[x]=1-x if 1/2<x<1.  This function is continuous but has a discontinuous first derivative at x=1/2. The Fourier series in Equation (4.65) is constructed by following steps analogous to those in Example 4.3.  The time development of the solution from the initial profile in (4.66) is shown in Figure 4.4 for the times t=0, t=0.025 and t=0.1.  The first twenty-five nonvanishing terms were included in the Fourier sine series.  As expected, the sharply-peaked
structure in f[x] near x=1/2 is lost quickly, since the higher-order Fourier components that are required to reproduce this behavior decay most rapidly.   As t->Infinity, the solution is seen to approach u=0, as all of the heat contained in the interior region has diffused to the boundaries, where it has been absorbed.

f[x_]:=If[Abs[x-0.5]<0.5,2(0.5-Abs[0.5-x]),0]

b[k_]:=2 NIntegrate[f[s]Sin[k Pi s],{s,0,1}]

Do[B[k]=b[2k+1],{k,0,12}]

u[x_,t_,n_]:=Sum[B[k]Exp[-((2k+1)Pi)^2 t]*
     Sin[(2k+1) Pi x],{k,0,n}]

Plot[{u[x,0,12],u[x,0.025,12],u[x,0.1,12]},{x,0,1},
     PlotRange->{0,1},AxesLabel->{x,u}]