First Order Differential Equations

Typically, the first differential equations encountered are first order equations. A first order differential equation takes the form

$$F\left(y^{\prime}, y, x\right)=0 \nonumber$$

There are two general forms for which one can formally obtain a solution. The first is the separable case and the second is a first order equation. We indicate that we can formally obtain solutions, as one can display the needed integration that leads to a solution. However, the resulting integrals are not always reducible to elementary functions nor does one obtain explicit solutions when the integrals are doable.

A first order equation is separable if it can be written the form

$$\dfrac{d y}{d x}=f(x) g(y) \nonumber$$

Special cases result when either $f(x)=1$ or $g(y)=1$. In the first case the equation is said to be autonomous.

The general solution to equation (1.5) is obtained in terms of two integrals:

$$\int \dfrac{d y}{g(y)}=\int f(x) d x+C, \nonumber$$

where $C$ is an integration constant. This yields a 1-parameter family of solutions to the differential equation corresponding to different values of $C$. If one can solve (1.6) for $y(x)$, then one obtains an explicit solution. Otherwise, one has a family of implicit solutions. If an initial condition is given as well, then one might be able to find a member of the family that satisfies this condition, which is often called a particular solution.

Example 1.1. $y^{\prime}=2 x y, y(0)=2$.

Applying (1.6), one has

$$\int \dfrac{d y}{y}=\int 2 x d x+C \nonumber$$

Integrating yields

$$\ln |y|=x^{2}+C . \nonumber$$

Exponentiating, one obtains the general solution,

$$y(x)=\pm e^{x^{2}+C}=A e^{x^{2}} . \nonumber$$

Here we have defined $A=\pm e^{C}$. Since $C$ is an arbitrary constant, $A$ is an arbitrary constant. Several solutions in this 1-parameter family are shown in Figure 1.1.

Next, one seeks a particular solution satisfying the initial condition. For $y(0)=2$, one finds that $A=2$. So, the particular solution satisfying the initial conditions is $y(x)=2 e^{x^{2}}$.

Example 1.2. $y y^{\prime}=-x$.

Following the same procedure as in the last example, one obtains:

$$\int y d y=-\int x d x+C \Rightarrow y^{2}=-x^{2}+A, \quad \text { where } \quad A=2 C . \nonumber$$

Thus, we obtain an implicit solution. Writing the solution as $x^{2}+y^{2}=A$, we see that this is a family of circles for $A>0$ and the origin for $A=0$. Plots of some solutions in this family are shown in Figure 1.2.

In this case one seeks an integrating factor, $\mu(x)$, which is a function that one can multiply through the equation making the left side a perfect derivative. Thus, obtaining,

$$\dfrac{d}{d x}[\mu(x) y(x)]=\mu(x) q(x) . \nonumber$$

The integrating factor that works is $\mu(x)=\exp \left(\int^{x} p(\xi) d \xi\right)$. One can show this by expanding the derivative in Equation (1.8),

$$\mu(x) y^{\prime}(x)+\mu^{\prime}(x) y(x)=\mu(x) q(x) \nonumber$$

and comparing this equation to the one obtained from multiplying (1.7) by $\mu(x)$ :

$$\mu(x) y^{\prime}(x)+\mu(x) p(x) y(x)=\mu(x) q(x) . \nonumber$$

Note that these last two equations would be the same if

$$\dfrac{d \mu(x)}{d x}=\mu(x) p(x) . \nonumber$$

This is a separable first order equation whose solution is the above given form for the integrating factor,

$$\mu(x)=\exp \left(\int^{x} p(\xi) d \xi\right) \nonumber$$

Equation (1.8) is easily integrated to obtain

$$y(x)=\dfrac{1}{\mu(x)}\left[\int^{x} \mu(\xi) q(\xi) d \xi+C\right] . \nonumber$$

Example 1.3. $x y^{\prime}+y=x, \quad x>0, y(1)=0$.

One first notes that this is a linear first order differential equation. Solving for $y^{\prime}$, one can see that the original equation is not separable. However, it is not in the standard form. So, we first rewrite the equation as

$$\dfrac{d y}{d x}+\dfrac{1}{x} y=1 . \nonumber$$

Noting that $p(x)=\dfrac{1}{x}$, we determine the integrating factor

$$\mu(x)=\exp \left[\int^{x} \dfrac{d \xi}{\xi}\right]=e^{\ln x}=x . \nonumber$$

Multiplying equation (1.13) by $\mu(x)=x$, we actually get back the original equation! In this case we have found that $x y^{\prime}+y$ must have been the derivative of something to start. In fact, $(x y)^{\prime}=x y^{\prime}+x$. Therefore, equation (1.8) becomes

$$(x y)^{\prime}=x \nonumber$$

Integrating one obtains

$$x y=\dfrac{1}{2} x^{2}+C, \nonumber$$

or

$$y(x)=\dfrac{1}{2} x+\dfrac{C}{x} . \nonumber$$

Inserting the initial condition into this solution, we have $0=\dfrac{1}{2}+C$. Therefore, $C=-\dfrac{1}{2}$. Thus, the solution of the initial value problem is $y(x)=\dfrac{1}{2}\left(x-\dfrac{1}{x}\right)$

Example 1.4. $(\sin x) y^{\prime}+(\cos x) y=x^{2} \sin x$.

Actually, this problem is easy if you realize that

$$\dfrac{d}{d x}((\sin x) y)=(\sin x) y^{\prime}+(\cos x) y \nonumber$$

But, we will go through the process of finding the integrating factor for practice.

First, rewrite the original differential equation in standard form:

$$y^{\prime}+(\cot x) y=x^{2} \nonumber$$

Then, compute the integrating factor as

$$\mu(x)=\exp \left(\int^{x} \cot \xi d \xi\right)=e^{-\ln (\sin x)}=\dfrac{1}{\sin x} \nonumber$$

Using the integrating factor, the original equation becomes

$$\dfrac{d}{d x}((\sin x) y)=x^{2} . \nonumber$$

Integrating, we have

$$y \sin x=\dfrac{1}{3} x^{3}+C \nonumber$$

So, the solution is

$$y=\left(\dfrac{1}{3} x^{3}+C\right) \csc x \nonumber$$

There are other first order equations that one can solve for closed form solutions. However, many equations are not solvable, or one is simply interested in the behavior of solutions. In such cases one turns to direction fields. We will return to a discussion of the qualitative behavior of differential equations later in the course.

Second Order Linear Differential Equations

Second order differential equations are typically harder than first order. In most cases students are only exposed to second order linear differential equations. A general form for a second order linear differential equation is given by

$$a(x) y^{\prime \prime}(x)+b(x) y^{\prime}(x)+c(x) y(x)=f(x) . \nonumber$$

One can rewrite this equation using operator terminology. Namely, one first defines the differential operator $L=a(x) D^{2}+b(x) D+c(x)$, where $D=\dfrac{d}{d x}$. Then equation (1.14) becomes

$$L y=f \nonumber$$

The solutions of linear differential equations are found by making use of the linearity of $L$. Namely, we consider the vector space ${ }^{1}$ consisting of realvalued functions over some domain. Let $f$ and $g$ be vectors in this function space. $L$ is a linear operator if for two vectors $f$ and $g$ and scalar $a$, we have that

a. $L(f+g)=L f+L g$

b. $L(a f)=a L f$.

One typically solves (1.14) by finding the general solution of the homogeneous problem,

$$L y_{h}=0 \nonumber$$

and a particular solution of the nonhomogeneous problem,

$$L y_{p}=f \text {. } \nonumber$$

Then the general solution of (1.14) is simply given as $y=y_{h}+y_{p}$. This is true because of the linearity of $L$. Namely,

$$\begin{aligned} L y &=L\left(y_{h}+y_{p}\right) \\[4pt] &=L y_{h}+L y_{p} \\[4pt] &=0+f=f \end{aligned} \nonumber$$

There are methods for finding a particular solution of a differential equation. These range from pure guessing to the Method of Undetermined Coefficients, or by making use of the Method of Variation of Parameters. We will review some of these methods later.

Determining solutions to the homogeneous problem, $L y_{h}=0$, is not always easy. However, others have studied a variety of second order linear equations

${ }^{1}$ We assume that the reader has been introduced to concepts in linear algebra. Late in the text we will recall the definition of a vector space and see that linear algebra is in the background of the study of many concepts in the solution of differential equations. and have saved us the trouble for some of the differential equations that often appear in applications.

Again, linearity is useful in producing the general solution of a homogeneous linear differential equation. If $y_{1}$ and $y_{2}$ are solutions of the homogeneous equation, then the linear combination $y=c_{1} y_{1}+c_{2} y_{2}$ is also a solution of the homogeneous equation. In fact, if $y_{1}$ and $y_{2}$ are linearly independent, ${ }^{2}$ then $y=c_{1} y_{1}+c_{2} y_{2}$ is the general solution of the homogeneous problem. As you may recall, linear independence is established if the Wronskian of the solutions in not zero. In this case, we have

$$W\left(y_{1}, y_{2}\right)=y_{1}(x) y_{2}^{\prime}(x)-y_{1}^{\prime}(x) y_{2}(x) \neq 0 \nonumber$$

Constant Coefficient Equations

The simplest and most seen second order differential equations are those with constant coefficients. The general form for a homogeneous constant coefficient second order linear differential equation is given as

$$a y^{\prime \prime}(x)+b y^{\prime}(x)+c y(x)=0, \nonumber$$

where $a, b$, and $c$ are constants.

Solutions to (1.18) are obtained by making a guess of $y(x)=e^{r x}$. Inserting this guess into (1.18) leads to the characteristic equation

$$a r^{2}+b r+c=0 . \nonumber$$

The roots of this equation in turn lead to three types of solution depending upon the nature of the roots as shown below.

Example 1.5. $y^{\prime \prime}-y^{\prime}-6 y=0 y(0)=2, y^{\prime}(0)=0$.

The characteristic equation for this problem is $r^{2}-r-6=0$. The roots of this equation are found as $r=-2,3$. Therefore, the general solution can be quickly written down:

$$y(x)=c_{1} e^{-2 x}+c_{2} e^{3 x} . \nonumber$$

Note that there are two arbitrary constants in the general solution. Therefore, one needs two pieces of information to find a particular solution. Of course, we have the needed information in the form of the initial conditions.

One also needs to evaluate the first derivative

$$y^{\prime}(x)=-2 c_{1} e^{-2 x}+3 c_{2} e^{3 x} \nonumber$$

${ }^{2}$ Recall, a set of functions $\left\{y_{i}(x)\right\}_{i=1}^{n}$ is a linearly independent set if and only if

$$c_{1} y\left(1(x)+\ldots+c_{n} y_{n}(x)=0\right. \nonumber$$

implies $c_{i}=0$, for $i=1, \ldots, n$. in order to attempt to satisfy the initial conditions. Evaluating $y$ and $y^{\prime}$ at $x=0$ yields

$$\begin{aligned} &2=c_{1}+c_{2} \\[4pt] &0=-2 c_{1}+3 c_{2} \end{aligned} \nonumber$$

These two equations in two unknowns can readily be solved to give $c_{1}=6 / 5$ and $c_{2}=4 / 5$. Therefore, the solution of the initial value problem is obtained as $y(x)=\dfrac{6}{5} e^{-2 x}+\dfrac{4}{5} e^{3 x}$.

Classification of Roots of the Characteristic Equation for Second Order Constant Coefficient ODEs

1. Real, distinct roots $r_{1}, r_{2}$. In this case the solutions corresponding to each root are linearly independent. Therefore, the general solution is simply $y(x)=c_{1} e^{r_{1} x}+c_{2} e^{r_{2} x}$.
2. Real, equal roots $r_{1}=r_{2}=r$. In this case the solutions corresponding to each root are linearly dependent. To find a second linearly independent solution, one uses the Method of Reduction of Order. This gives the second solution as $x e^{r x}$. Therefore, the general solution is found as $y(x)=\left(c_{1}+c_{2} x\right) e^{r x}$. [This is covered in the appendix to this chapter.]
3. Complex conjugate roots $r_{1}, r_{2}=\alpha \pm i \beta$. In this case the solutions corresponding to each root are linearly independent. Making use of Euler’s identity, $e^{i \theta}=\cos (\theta)+i \sin (\theta)$, these complex exponentials can be rewritten in terms of trigonometric functions. Namely, one has that $e^{\alpha x} \cos (\beta x)$ and $e^{\alpha x} \sin (\beta x)$ are two linearly independent solutions. Therefore, the general solution becomes $y(x)=e^{\alpha x}\left(c_{1} \cos (\beta x)+c_{2} \sin (\beta x)\right)$. [This is covered in the appendix to this chapter.]

Example 1.6. $y^{\prime \prime}+6 y^{\prime}+9 y=0$.

In this example we have $r^{2}+6 r+9=0$. There is only one root, $r=-3$. Again, the solution is easily obtained as $y(x)=\left(c_{1}+c_{2} x\right) e^{-3 x}$.

Example 1.7. $y^{\prime \prime}+4 y=0$.

The characteristic equation in this case is $r^{2}+4=0$. The roots are pure imaginary roots, $r=\pm 2 i$ and the general solution consists purely of sinusoidal functions: $y(x)=c_{1} \cos (2 x)+c_{2} \sin (2 x)$

Example 1.8. $y^{\prime \prime}+2 y^{\prime}+4 y=0$.

The characteristic equation in this case is $r^{2}+2 r+4=0$. The roots are complex, $r=-1 \pm \sqrt{3} i$ and the general solution can be written as $y(x)=$ $\left[c_{1} \cos (\sqrt{3} x)+c_{2} \sin (\sqrt{3} x)\right] e^{-x}$ One of the most important applications of the equations in the last two examples is in the study of oscillations. Typical systems are a mass on a spring, or a simple pendulum. For a mass $m$ on a spring with spring constant $k>0$, one has from Hooke’s law that the position as a function of time, $x(t)$, satisfies the equation

$$m x^{\prime \prime}+k x=0 . \nonumber$$

This constant coefficient equation has pure imaginary roots $(\alpha=0)$ and the solutions are pure sines and cosines. Such motion is called simple harmonic motion.

Adding a damping term and periodic forcing complicates the dynamics, but is nonetheless solvable. The next example shows a forced harmonic oscillator.

Example 1.9. $y^{\prime \prime}+4 y=\sin x$.

This is an example of a nonhomogeneous problem. The homogeneous problem was actually solved in Example 1.7. According to the theory, we need only seek a particular solution to the nonhomogeneous problem and add it to the solution of the last example to get the general solution.

The particular solution can be obtained by purely guessing, making an educated guess, or using the Method of Variation of Parameters. We will not review all of these techniques at this time. Due to the simple form of the driving term, we will make an intelligent guess of $y_{p}(x)=A \sin x$ and determine what $A$ needs to be. Recall, this is the Method of Undetermined Coefficients which we review in the next section. Inserting our guess in the equation gives $(-A+4 A) \sin x=\sin x$. So, we see that $A=1 / 3$ works. The general solution of the nonhomogeneous problem is therefore $y(x)=$ $c_{1} \cos (2 x)+c_{2} \sin (2 x)+\dfrac{1}{3} \sin x$

Method of Undetermined Coefficients

To date, we only know how to solve constant coefficient, homogeneous equations. How does one solve a nonhomogeneous equation like that in Equation $(1.14)$

$$a(x) y^{\prime \prime}(x)+b(x) y^{\prime}(x)+c(x) y(x)=f(x) \nonumber$$

Recall, that one solves this equation by finding the general solution of the homogeneous problem,

$$L y_{h}=0 \nonumber$$

and a particular solution of the nonhomogeneous problem,

$$L y_{p}=f \nonumber$$

Then the general solution of (1.14) is simply given as $y=y_{h}+y_{p}$. So, how do we find the particular solution? You could guess a solution, but that is not usually possible without a little bit of experience. So we need some other methods. There are two main methods. In the first case, the Method of Undetermined Coefficients, one makes an intelligent guess based on the form of $f(x)$. In the second method, one can systematically develop the particular solution. We will come back to this method the Method of Variation of Parameters, later in the book.

Let’s solve a simple differential equation highlighting how we can handle nonhomogeneous equations.

Example 1.10. Consider the equation

$$y^{\prime \prime}+2 y^{\prime}-3 y=4 \nonumber$$

The first step is to determine the solution of the homogeneous equation. Thus, we solve

$$y_{h}^{\prime \prime}+2 y_{h}^{\prime}-3 y_{h}=0 . \nonumber$$

The characteristic equation is $r^{2}+2 r-3=0$. The roots are $r=1,-3$. So, we can immediately write the solution

$$y_{h}(x)=c_{1} e^{x}+c_{2} e^{-3 x} . \nonumber$$

The second step is to find a particular solution of (1.22). What possible function can we insert into this equation such that only a 4 remains? If we try something proportional to $x$, then we are left with a linear function after inserting $x$ and its derivatives. Perhaps a constant function you might think. $y=4$ does not work. But, we could try an arbitrary constant, $y=A$.

Let’s see. Inserting $y=A$ into $(1.22)$, we obtain

$$-3 A=4 . \nonumber$$

Ah ha! We see that we can choose $A=-\dfrac{4}{3}$ and this works. So, we have a particular solution, $y_{p}(x)=-\dfrac{4}{3}$. This step is done.

Combining our two solutions, we have the general solution to the original nonhomogeneous equation (1.22). Namely,

$$y(x)=y_{h}(x)+y_{p}(x)=c_{1} e^{x}+c_{2} e^{-3 x}-\dfrac{4}{3} \nonumber$$

Insert this solution into the equation and verify that it is indeed a solution. If we had been given initial conditions, we could now use them to determine our arbitrary constants.

What if we had a different source term? Consider the equation

$$y^{\prime \prime}+2 y^{\prime}-3 y=4 x . \nonumber$$

The only thing that would change is our particular solution. So, we need a guess. We know a constant function does not work by the last example. So, let’s try $y_{p}=A x$. Inserting this function into Equation (??), we obtain

$$2 A-3 A x=4 x . \nonumber$$

Picking $A=-4 / 3$ would get rid of the $x$ terms, but will not cancel everything. We still have a constant left. So, we need something more general.

Let’s try a linear function, $y_{p}(x)=A x+B$. Then we get after substitution $\operatorname{into}(1.24)$

$$2 A-3(A x+B)=4 x . \nonumber$$

Equating the coefficients of the different powers of $x$ on both sides, we find a system of equations for the undetermined coefficients:

$$\begin{array}{r} 2 A-3 B=0 \\[4pt] -3 A=4 . \end{array} \nonumber$$

These are easily solved to obtain

$$\begin{aligned} &A=-\dfrac{4}{3} \\[4pt] &B=\dfrac{2}{3} A=-\dfrac{8}{9} \end{aligned} \nonumber$$

So, our particular solution is

$$y_{p}(x)=-\dfrac{4}{3} x-\dfrac{8}{9} . \nonumber$$

This gives the general solution to the nonhomogeneous problem as

$$y(x)=y_{h}(x)+y_{p}(x)=c_{1} e^{x}+c_{2} e^{-3 x}-\dfrac{4}{3} x-\dfrac{8}{9} \nonumber$$

There are general forms that you can guess based upon the form of the driving term, $f(x)$. Some examples are given in Table 1.1.4. More general applications are covered in a standard text on differential equations. However, the procedure is simple. Given $f(x)$ in a particular form, you make an appropriate guess up to some unknown parameters, or coefficients. Inserting the guess leads to a system of equations for the unknown coefficients. Solve the system and you have your solution. This solution is then added to the general solution of the homogeneous differential equation.

Example 1.11. As a final example, let’s consider the equation

$$y^{\prime \prime}+2 y^{\prime}-3 y=2 e^{-3 x} \nonumber$$

According to the above, we would guess a solution of the form $y_{p}=A e^{-3 x}$. Inserting our guess, we find

$$0=2 e^{-3 x} \nonumber$$

Oops! The coefficient, $A$, disappeared! We cannot solve for it. What went wrong?

The answer lies in the general solution of the homogeneous problem. Note that $e^{x}$ and $e^{-3 x}$ are solutions to the homogeneous problem. So, a multiple of $e^{-3 x}$ will not get us anywhere. It turns out that there is one further modification of the method. If our driving term contains terms that are solutions of the homogeneous problem, then we need to make a guess consisting of the smallest possible power of $x$ times the function which is no longer a solution of the homogeneous problem. Namely, we guess $y_{p}(x)=A x e^{-3 x}$. We compute the derivative of our guess, $y_{p}^{\prime}=A(1-3 x) e^{-3 x}$ and $y_{p}^{\prime \prime}=A(9 x-6) e^{-3 x}$. Inserting these into the equation, we obtain

$$[(9 x-6)+2(1-3 x)-3 x] A e^{-3 x}=2 e^{-3 x} \nonumber$$

$$-4 A=2 . \nonumber$$

So, $A=-1 / 2$ and $y_{p}(x)=-\dfrac{1}{2} x e^{-3 x}$.

Cauchy-Euler Equations

Another class of solvable linear differential equations that is of interest are the Cauchy-Euler type of equations. These are given by

$$a x^{2} y^{\prime \prime}(x)+b x y^{\prime}(x)+c y(x)=0 . \nonumber$$

Note that in such equations the power of $x$ in each of the coefficients matches the order of the derivative in that term. These equations are solved in a manner similar to the constant coefficient equations.

One begins by making the guess $y(x)=x^{r}$. Inserting this function and its derivatives,

$$y^{\prime}(x)=r x^{r-1}, \quad y^{\prime \prime}(x)=r(r-1) x^{r-2}, \nonumber$$

into Equation (1.28), we have

$$[a r(r-1)+b r+c] x^{r}=0 \nonumber$$

Since this has to be true for all $x$ in the problem domain, we obtain the characteristic equation

$$a r(r-1)+b r+c=0 \nonumber$$

Just like the constant coefficient differential equation, we have a quadratic equation and the nature of the roots again leads to three classes of solutions. These are shown below. Some of the details are provided in the next section.

Classification of Roots of the Characteristic Equation for Cauchy-Euler Differential Equations

1. Real, distinct roots $r_{1}, r_{2}$. In this case the solutions corresponding to each root are linearly independent. Therefore, the general solution is simply $y(x)=c_{1} x^{r_{1}}+c_{2} x^{r_{2}}$.
2. Real, equal roots $r_{1}=r_{2}=r$. In this case the solutions corresponding to each root are linearly dependent. To find a second linearly independent solution, one uses the Method of Reduction of Order. This gives the second solution as $x^{r} \ln |x|$. Therefore, the general solution is found as $y(x)=\left(c_{1}+c_{2} \ln |x|\right) x^{r}$.
3. Complex conjugate roots $r_{1}, r_{2}=\alpha \pm i \beta$. In this case the solutions corresponding to each root are linearly independent. These complex exponentials can be rewritten in terms of trigonometric functions. Namely, one has that $x^{\alpha} \cos (\beta \ln |x|)$ and $x^{\alpha} \sin (\beta \ln |x|)$ are two linearly independent solutions. Therefore, the general solution becomes $y(x)=$ $x^{\alpha}\left(c_{1} \cos (\beta \ln |x|)+c_{2} \sin (\beta \ln |x|)\right)$.

Example 1.12. $x^{2} y^{\prime \prime}+5 x y^{\prime}+12 y=0$

As with the constant coefficient equations, we begin by writing down the characteristic equation. Doing a simple computation,

$$\begin{aligned} 0 &=r(r-1)+5 r+12 \\[4pt] &=r^{2}+4 r+12 \\[4pt] &=(r+2)^{2}+8 \\[4pt] -8 &=(r+2)^{2} \end{aligned} \nonumber$$

one determines the roots are $r=-2 \pm 2 \sqrt{2} i$. Therefore, the general solution is $y(x)=\left[c_{1} \cos (2 \sqrt{2} \ln |x|)+c_{2} \sin (2 \sqrt{2} \ln |x|)\right] x^{-2}$ Example 1.13. $t^{2} y^{\prime \prime}+3 t y^{\prime}+y=0, \quad y(1)=0, y^{\prime}(1)=1$.

For this example the characteristic equation takes the form

$$r(r-1)+3 r+1=0, \nonumber$$

$$r^{2}+2 r+1=0 . \nonumber$$

There is only one real root, $r=-1$. Therefore, the general solution is

$$y(t)=\left(c_{1}+c_{2} \ln |t|\right) t^{-1} . \nonumber$$

However, this problem is an initial value problem. At $t=1$ we know the values of $y$ and $y^{\prime}$. Using the general solution, we first have that

$$0=y(1)=c_{1} \nonumber$$

Thus, we have so far that $y(t)=c_{2} \ln |t| t^{-1}$. Now, using the second condition and

$$y^{\prime}(t)=c_{2}(1-\ln |t|) t^{-2}, \nonumber$$

we have

$$1=y(1)=c_{2} \text {. } \nonumber$$

Therefore, the solution of the initial value problem is $y(t)=\ln |t| t^{-1}$.

Nonhomogeneous Cauchy-Euler Equations We can also solve some nonhomogeneous Cauchy-Euler equations using the Method of Undetermined Coefficients. We will demonstrate this with a couple of examples.

Example 1.14. Find the solution of $x^{2} y^{\prime \prime}-x y^{\prime}-3 y=2 x^{2}$.

First we find the solution of the homogeneous equation. The characteristic equation is $r^{2}-2 r-3=0$. So, the roots are $r=-1,3$ and the solution is $y_{h}(x)=c_{1} x^{-1}+c_{2} x^{3}$

We next need a particular solution. Let’s guess $y_{p}(x)=A x^{2}$. Inserting the guess into the nonhomogeneous differential equation, we have

$$\begin{aligned} 2 x^{2} &=x^{2} y^{\prime \prime}-x y^{\prime}-3 y=2 x^{2} \\[4pt] &=2 A x^{2}-2 A x^{2}-3 A x^{2} \\[4pt] &=-3 A x^{2} \end{aligned} \nonumber$$

So, $A=-2 / 3$. Therefore, the general solution of the problem is

$$y(x)=c_{1} x^{-1}+c_{2} x^{3}-\dfrac{2}{3} x^{2} . \nonumber$$

Example 1.15. Find the solution of $x^{2} y^{\prime \prime}-x y^{\prime}-3 y=2 x^{3}$.

In this case the nonhomogeneous term is a solution of the homogeneous problem, which we solved in the last example. So, we will need a modification of the method. We have a problem of the form

$$a x^{2} y^{\prime \prime}+b x y^{\prime}+c y=d x^{r} \nonumber$$

where $r$ is a solution of $a r(r-1)+b r+c=0$. Let’s guess a solution of the form $y=A x^{r} \ln x$. Then one finds that the differential equation reduces to $A x^{r}(2 a r-a+b)=d x^{r}$. [You should verify this for yourself.]

With this in mind, we can now solve the problem at hand. Let $y_{p}=$ $A x^{3} \ln x$. Inserting into the equation, we obtain $4 A x^{3}=2 x^{3}$, or $A=1 / 2$. The general solution of the problem can now be written as

$$y(x)=c_{1} x^{-1}+c_{2} x^{3}+\dfrac{1}{2} x^{3} \ln x . \nonumber$$

Overview of the Course

For the most part, your first course in differential equations was about solving initial value problems. When second order equations did not fall into the above cases, then you might have learned how to obtain approximate solutions using power series methods, or even finding new functions from these methods. In this course we will explore two broad topics: systems of differential equations and boundary value problems.

We will see that there are interesting initial value problems when studying systems of differential equations. In fact, many of the second order equations that you have seen in the past can be written as a system of two first order equations. For example, the equation for simple harmonic motion,

$$x^{\prime \prime}+\omega^{2} x=0 \nonumber$$

can be written as the system

$$\begin{gathered} x^{\prime}=y \\[4pt] y^{\prime}=-\omega^{2} x \end{gathered} . \nonumber$$

Just note that $x^{\prime \prime}=y^{\prime}=-\omega^{2} x$. Of course, one can generalize this to systems with more complicated right hand sides. The behavior of such systems can be fairly interesting and these systems result from a variety of physical models.

In the second part of the course we will explore boundary value problems. Often these problems evolve from the study of partial differential equations. Such examples stem from vibrating strings, temperature distributions, bending beams, etc. Boundary conditions are conditions that are imposed at more than one point, while for initial value problems the conditions are specified at one point. For example, we could take the oscillation equation above and ask when solutions of the equation would satisfy the conditions $x(0)=0$ and $x(1)=0$. The general solution, as we have determined earlier, is

$$x(t)=c_{1} \cos \omega t+c_{2} \sin \omega t \nonumber$$

Requiring $x(0)=0$, we find that $c_{1}=0$, leaving $x(t)=c_{2} \sin \omega t$. Also imposing that $0=x(1)=c_{2} \sin \omega$, we are forced to make $\omega=n \pi$, for $n=1,2, \ldots$. (Making $c_{2}=0$ would not give a nonzero solution of the problem.) Thus, there are an infinite number of solutions possible, if we have the freedom to choose our $\omega$. In the second half of the course we will investigate techniques for solving boundary value problems and look at several applications, including seeing the connections with partial differential equations and Fourier series.

Appendix: Reduction of Order and Complex Roots

In this section we provide some of the details leading to the general forms for the constant coefficient and Cauchy-Euler differential equations. In the first subsection we review how the Method of Reduction of Order is used to obtain the second linearly independent solutions for the case of one repeated root. In the second subsection we review how the complex solutions can be used to produce two linearly independent real solutions.

Introduction

In this chapter we will begin our study of systems of differential equations. After defining first order systems, we will look at constant coefficient systems and the behavior of solutions for these systems. Also, most of the discussion will focus on planar, or two dimensional, systems. For such systems we will be able to look at a variety of graphical representations of the family of solutions and discuss the qualitative features of systems we can solve in preparation for the study of systems whose solutions cannot be found in an algebraic form.

A general form for first order systems in the plane is given by a system of two equations for unknowns $x(t)$ and $y(t)$ :

$$\begin{aligned} &x^{\prime}(t)=P(x, y, t) \\[4pt] &y^{\prime}(t)=Q(x, y, t) \end{aligned} \nonumber$$

An autonomous system is one in which there is no explicit time dependence:

$$\begin{aligned} &x^{\prime}(t)=P(x, y) \\[4pt] &y^{\prime}(t)=Q(x, y) \end{aligned} \nonumber$$

Otherwise the system is called nonautonomous.

A linear system takes the form

$$\begin{aligned} x^{\prime} &=a(t) x+b(t) y+e(t) \\[4pt] y^{\prime} &=c(t) x+d(t) y+f(t) \end{aligned} \nonumber$$

A homogeneous linear system results when $e(t)=0$ and $f(t)=0$.

A linear, constant coefficient system of first order differential equations is given by

$$\begin{aligned} &x^{\prime}=a x+b y+e \\[4pt] &y^{\prime}=c x+d y+f \end{aligned} \nonumber$$

We will focus on linear, homogeneous systems of constant coefficient first order differential equations:

As we will see later, such systems can result by a simple translation of the unknown functions. These equations are said to be coupled if either $b \neq 0$ or $c \neq 0$.

We begin by noting that the system (2.5) can be rewritten as a second order constant coefficient linear differential equation, which we already know how to solve. We differentiate the first equation in system system (2.5) and systematically replace occurrences of $y$ and $y^{\prime}$, since we also know from the first equation that $y=\dfrac{1}{b}\left(x^{\prime}-a x\right)$. Thus, we have

$$\begin{aligned} x^{\prime \prime} &=a x^{\prime}+b y^{\prime} \\[4pt] &=a x^{\prime}+b(c x+d y) \\[4pt] &=a x^{\prime}+b c x+d\left(x^{\prime}-a x\right) . \end{aligned} \nonumber$$

Rewriting the last line, we have

$x^{\prime \prime}-(a+d) x^{\prime}+(a d-b c) x=0 .$

This is a linear, homogeneous, constant coefficient ordinary differential equation. We know that we can solve this by first looking at the roots of the characteristic equation

$$r^{2}-(a+d) r+a d-b c=0 \nonumber$$

and writing down the appropriate general solution for $x(t)$. Then we can find $y(t)$ using Equation (2.5):

$$y=\dfrac{1}{b}\left(x^{\prime}-a x\right) . \nonumber$$

We now demonstrate this for a specific example.

Example 2.1. Consider the system of differential equations

$$\begin{aligned} &x^{\prime}=-x+6 y \\[4pt] &y^{\prime}=x-2 y . \end{aligned} \nonumber$$

Carrying out the above outlined steps, we have that $x^{\prime \prime}+3 x^{\prime}-4 x=0$. This can be shown as follows:

$$\begin{aligned} x^{\prime \prime} &=-x^{\prime}+6 y^{\prime} \\[4pt] &=-x^{\prime}+6(x-2 y) \\[4pt] &=-x^{\prime}+6 x-12\left(\dfrac{x^{\prime}+x}{6}\right) \\[4pt] &=-3 x^{\prime}+4 x \end{aligned} \nonumber$$

The resulting differential equation has a characteristic equation of $r^{2}+3 r-$ $4=0$. The roots of this equation are $r=1,-4$. Therefore, $x(t)=c_{1} e^{t}+c_{2} e^{-4 t} .$ But, we still need $y(t)$. From the first equation of the system we have

$$y(t)=\dfrac{1}{6}\left(x^{\prime}+x\right)=\dfrac{1}{6}\left(2 c_{1} e^{t}-3 c_{2} e^{-4 t}\right) . \nonumber$$

Thus, the solution to our system is

$$\begin{aligned} &x(t)=c_{1} e^{t}+c_{2} e^{-4 t}, \\[4pt] &y(t)=\dfrac{1}{3} c_{1} e^{t}-\dfrac{1}{2} c_{2} e^{-4 t} \end{aligned} \nonumber$$

Sometimes one needs initial conditions. For these systems we would specify conditions like $x(0)=x_{0}$ and $y(0)=y_{0}$. These would allow the determination of the arbitrary constants as before.

Example 2.2. Solve

$$\begin{aligned} &x^{\prime}=-x+6 y \\[4pt] &y^{\prime}=x-2 y \end{aligned} \nonumber$$

given $x(0)=2, y(0)=0$.

We already have the general solution of this system in (2.11). Inserting the initial conditions, we have

$$\begin{aligned} &2=c_{1}+c_{2}, \\[4pt] &0=\dfrac{1}{3} c_{1}-\dfrac{1}{2} c_{2} . \end{aligned} \nonumber$$

Solving for $c_{1}$ and $c_{2}$ gives $c_{1}=6 / 5$ and $c_{2}=4 / 5$. Therefore, the solution of the initial value problem is

$$\begin{aligned} &x(t)=\dfrac{2}{5}\left(3 e^{t}+2 e^{-4 t}\right) \\[4pt] &y(t)=\dfrac{2}{5}\left(e^{t}-e^{-4 t}\right) \end{aligned} \nonumber$$

Polar Representation of Spirals

In the examples with a center or a spiral, one might be able to write the solutions in polar coordinates. Recall that a point in the plane can be described by either Cartesian $(x, y)$ or polar $(r, \theta)$ coordinates. Given the polar form, one can find the Cartesian components using

$$x=r \cos \theta \text { and } y=r \sin \theta . \nonumber$$

Given the Cartesian coordinates, one can find the polar coordinates using

$$r^{2}=x^{2}+y^{2} \text { and } \tan \theta=\dfrac{y}{x} . \nonumber$$

Since $x$ and $y$ are functions of $t$, then naturally we can think of $r$ and $\theta$ as functions of $t$. The equations that they satisfy are obtained by differentiating the above relations with respect to $t$.

Differentiating the first equation in $(2.27)$ gives

$$r r^{\prime}=x x^{\prime}+y y^{\prime} . \nonumber$$

Inserting the expressions for $x^{\prime}$ and $y^{\prime}$ from system $2.5$, we have

$$r r^{\prime}=x(a x+b y)+y(c x+d y) . \nonumber$$

In some cases this may be written entirely in terms of $r$ ’s. Similarly, we have that

$$\theta^{\prime}=\dfrac{x y^{\prime}-y x^{\prime}}{r^{2}} \nonumber$$

which the reader can prove for homework.

In summary, when converting first order equations from rectangular to polar form, one needs the relations below.

Matrix Formulation

We have investigated several linear systems in the plane and in the next chapter we will use some of these ideas to investigate nonlinear systems. We need a deeper insight into the solutions of planar systems. So, in this section we will recast the first order linear systems into matrix form. This will lead to a better understanding of first order systems and allow for extensions to higher dimensions and the solution of nonhomogeneous equations later in this chapter.

We start with the usual homogeneous system in Equation (2.5). Let the unknowns be represented by the vector

$$\mathbf{x}(t)=\left(\begin{array}{c} x(t) \\[4pt] y(t) \end{array}\right) \nonumber$$

Then we have that

$$\mathbf{x}^{\prime}=\left(\begin{array}{l} x^{\prime} \\[4pt] y^{\prime} \end{array}\right)=\left(\begin{array}{c} a x+b y \\[4pt] c x+d y \end{array}\right)=\left(\begin{array}{ll} a & b \\[4pt] c & d \end{array}\right)\left(\begin{array}{l} x \\[4pt] y \end{array}\right) \equiv A \mathbf{x} \nonumber$$

Here we have introduced the coefficient matrix $A$. This is a first order vector differential equation,

$$\mathbf{x}^{\prime}=A \mathbf{x} \nonumber$$

Formerly, we can write the solution as

$$\mathbf{x}=\mathbf{x}_{0} e^{A t} \nonumber$$

$\overline{1}$ The exponential of a matrix is defined using the Maclaurin series expansion We would like to investigate the solution of our system. Our investigations will lead to new techniques for solving linear systems using matrix methods.

We begin by recalling the solution to the specific problem (2.12). We obtained the solution to this system as

$$\begin{gathered} x(t)=c_{1} e^{t}+c_{2} e^{-4 t}, \\[4pt] y(t)=\dfrac{1}{3} c_{1} e^{t}-\dfrac{1}{2} c_{2} e^{-4 t} \end{gathered} \nonumber$$

This can be rewritten using matrix operations. Namely, we first write the solution in vector form.

$$\begin{aligned} \mathbf{x} &=\left(\begin{array}{c} x(t) \\[4pt] y(t) \end{array}\right) \\[4pt] &=\left(\begin{array}{c} c_{1} e^{t}+c_{2} e^{-4 t} \\[4pt] \dfrac{1}{3} c_{1} e^{t}-\dfrac{1}{2} c_{2} e^{-4 t} \end{array}\right) \\[4pt] &=\left(\begin{array}{c} c_{1} e^{t} \\[4pt] \dfrac{1}{3} c_{1} e^{t} \end{array}\right)+\left(\begin{array}{c} c_{2} e^{-4 t} \\[4pt] -\dfrac{1}{2} c_{2} e^{-4 t} \end{array}\right) \\[4pt] &=c_{1}\left(\begin{array}{c} 1 \\[4pt] \dfrac{1}{3} \end{array}\right) e^{t}+c_{2}\left(\begin{array}{c} 1 \\[4pt] -\dfrac{1}{2} \end{array}\right) e^{-4 t} \end{aligned} \nonumber$$

We see that our solution is in the form of a linear combination of vectors of the form

$$\mathbf{x}=\mathbf{v} e^{\lambda t} \nonumber$$

with $\mathbf{v}$ a constant vector and $\lambda$ a constant number. This is similar to how we began to find solutions to second order constant coefficient equations. So, for the general problem (2.3) we insert this guess. Thus,

$$\begin{aligned} \mathbf{x}^{\prime} &=A \mathbf{x} \Rightarrow \\[4pt] \lambda \mathbf{v} e^{\lambda t} &=A \mathbf{v} e^{\lambda t} . \end{aligned} \nonumber$$

For this to be true for all $t$, we have that

$$A \mathbf{v}=\lambda \mathbf{v} . \nonumber$$

This is an eigenvalue problem. $A$ is a $2 \times 2$ matrix for our problem, but could easily be generalized to a system of $n$ first order differential equations. We will confine our remarks for now to planar systems. However, we need to recall how to solve eigenvalue problems and then see how solutions of eigenvalue problems can be used to obtain solutions to our systems of differential equations..

$$e^{x}=\sum_{k=0}^{\infty}=1+x+\dfrac{x^{2}}{2 !}+\dfrac{x^{3}}{3 !}+\cdots \nonumber$$

So, we define

$$e^{A}=\sum_{k=0}^{\infty}=I+A+\dfrac{A^{2}}{2 !}+\dfrac{A^{3}}{3 !}+\cdots \nonumber$$

In general, it is difficult computing $e^{A}$ unless $A$ is diagonal.

Solving Constant Coefficient Systems in $2 \mathrm{D}$

Before proceeding to examples, we first indicate the types of solutions that could result from the solution of a homogeneous, constant coefficient system of first order differential equations. We begin with the linear system of differential equations in matrix form.

$$\dfrac{d \mathbf{x}}{d t}=\left(\begin{array}{ll} a & b \\[4pt] c & d \end{array}\right) \mathbf{x}=A \mathbf{x} . \nonumber$$

The type of behavior depends upon the eigenvalues of matrix $A$. The procedure is to determine the eigenvalues and eigenvectors and use them to construct the general solution.

If we have an initial condition, $\mathbf{x}\left(t_{0}\right)=\mathbf{x}_{0}$, we can determine the two arbitrary constants in the general solution in order to obtain the particular solution. Thus, if $\mathbf{x}_{1}(t)$ and $\mathbf{x}_{2}(t)$ are two linearly independent solutions ${ }^{2}$, then the general solution is given as

$$\mathbf{x}(t)=c_{1} \mathbf{x}_{1}(t)+c_{2} \mathbf{x}_{2}(t) \nonumber$$

Then, setting $t=0$, we get two linear equations for $c_{1}$ and $c_{2}$ :

$$c_{1} \mathbf{x}_{1}(0)+c_{2} \mathbf{x}_{2}(0)=\mathbf{x}_{0} \nonumber$$

The major work is in finding the linearly independent solutions. This depends upon the different types of eigenvalues that one obtains from solving the eigenvalue equation, $\operatorname{det}(A-\lambda I)=0$. The nature of these roots indicate the form of the general solution. On the next page we summarize the classification of solutions in terms of the eigenvalues of the coefficient matrix. We first make some general remarks about the plausibility of these solutions and then provide examples in the following section to clarify the matrix methods for our two dimensional systems.

The construction of the general solution in Case I is straight forward. However, the other two cases need a little explanation.

${ }^{2}$ Recall that linear independence means $c_{1} \mathbf{x}_{1}(t)+c_{2} \mathbf{x}_{2}(t)=\mathbf{0}$ if and only if $c_{1}, c_{2}=$ 0 . The reader should derive the condition on the $\mathbf{x}_{i}$ for linear independence.

1. Case I: Two real, distinct roots.

Solve the eigenvalue problem $A \mathbf{v}=\lambda \mathbf{v}$ for each eigenvalue obtaining two eigenvectors $\mathbf{v}_{1}, \mathbf{v}_{2}$. Then write the general solution as a linear combination $\mathbf{x}(t)=c_{1} e^{\lambda_{1} t} \mathbf{v}_{1}+c_{2} e^{\lambda_{2} t} \mathbf{v}_{2}$

Examples of the Matrix Method

Here we will give some examples for typical systems for the three cases mentioned in the last section.

Example 2.14. $A=\left(\begin{array}{ll}4 & 2 \\[4pt] 3 & 3\end{array}\right)$.

Eigenvalues: We first determine the eigenvalues.

$$0=\left|\begin{array}{cc} 4-\lambda & 2 \\[4pt] 3 & 3-\lambda \end{array}\right| \nonumber$$

Therefore,

$$\begin{aligned} &0=(4-\lambda)(3-\lambda)-6 \\[4pt] &0=\lambda^{2}-7 \lambda+6 \\[4pt] &0=(\lambda-1)(\lambda-6) \end{aligned} \nonumber$$

The eigenvalues are then $\lambda=1,6$. This is an example of Case I.

Eigenvectors: Next we determine the eigenvectors associated with each of these eigenvalues. We have to solve the system $A \mathbf{v}=\lambda \mathbf{v}$ in each case. Case $\lambda=1$

$$\begin{gathered} \left(\begin{array}{ll} 4 & 2 \\[4pt] 3 & 3 \end{array}\right)\left(\begin{array}{l} v_{1} \\[4pt] v_{2} \end{array}\right)=\left(\begin{array}{l} v_{1} \\[4pt] v_{2} \end{array}\right) \\[4pt] \left(\begin{array}{ll} 3 & 2 \\[4pt] 3 & 2 \end{array}\right)\left(\begin{array}{l} v_{1} \\[4pt] v_{2} \end{array}\right)=\left(\begin{array}{l} 0 \\[4pt] 0 \end{array}\right) \end{gathered} \nonumber$$

This gives $3 v_{1}+2 v_{2}=0$. One possible solution yields an eigenvector of

$$\left(\begin{array}{l} v_{1} \\[4pt] v_{2} \end{array}\right)=\left(\begin{array}{c} 2 \\[4pt] -3 \end{array}\right) \nonumber$$

Case $\lambda=6$

$$\begin{aligned} &\left(\begin{array}{cc} 4 & 2 \\[4pt] 3 & 3 \end{array}\right)\left(\begin{array}{l} v_{1} \\[4pt] v_{2} \end{array}\right)=6\left(\begin{array}{l} v_{1} \\[4pt] v_{2} \end{array}\right) \\[4pt] &\left(\begin{array}{cc} -2 & 2 \\[4pt] 3 & -3 \end{array}\right)\left(\begin{array}{l} v_{1} \\[4pt] v_{2} \end{array}\right)=\left(\begin{array}{l} 0 \\[4pt] 0 \end{array}\right) \end{aligned} \nonumber$$

For this case we need to solve $-2 v_{1}+2 v_{2}=0$. This yields

$$\left(\begin{array}{l} v_{1} \\[4pt] v_{2} \end{array}\right)=\left(\begin{array}{l} 1 \\[4pt] 1 \end{array}\right) \nonumber$$

General Solution: We can now construct the general solution.

$$\begin{aligned} \mathbf{x}(t) &=c_{1} e^{\lambda_{1} t} \mathbf{v}_{1}+c_{2} e^{\lambda_{2} t} \mathbf{v}_{2} \\[4pt] &=c_{1} e^{t}\left(\begin{array}{c} 2 \\[4pt] -3 \end{array}\right)+c_{2} e^{6 t}\left(\begin{array}{l} 1 \\[4pt] 1 \end{array}\right) \\[4pt] &=\left(\begin{array}{c} 2 c_{1} e^{t}+c_{2} e^{6 t} \\[4pt] -3 c_{1} e^{t}+c_{2} e^{6 t} \end{array}\right) . \end{aligned} \nonumber$$

Example 2.15. $A=\left(\begin{array}{ll}3 & -5 \\[4pt] 1 & -1\end{array}\right)$.

Eigenvalues: Again, one solves the eigenvalue equation.

$$0=\left|\begin{array}{cc} 3-\lambda & -5 \\[4pt] 1 & -1-\lambda \end{array}\right| \nonumber$$

Therefore,

$$\begin{aligned} &0=(3-\lambda)(-1-\lambda)+5 \\[4pt] &0=\lambda^{2}-2 \lambda+2 \\[4pt] &\lambda=\dfrac{-(-2) \pm \sqrt{4-4(1)(2)}}{2}=1 \pm i \end{aligned} \nonumber$$

The eigenvalues are then $\lambda=1+i, 1-i$. This is an example of Case III.

Eigenvectors: In order to find the general solution, we need only find the eigenvector associated with $1+i$.

$$\begin{gathered} \left(\begin{array}{l} 3-5 \\[4pt] 1-1 \end{array}\right)\left(\begin{array}{l} v_{1} \\[4pt] v_{2} \end{array}\right)=(1+i)\left(\begin{array}{l} v_{1} \\[4pt] v_{2} \end{array}\right) \\[4pt] \left(\begin{array}{cc} 2-i & -5 \\[4pt] 1 & -2-i \end{array}\right)\left(\begin{array}{l} v_{1} \\[4pt] v_{2} \end{array}\right)=\left(\begin{array}{l} 0 \\[4pt] 0 \end{array}\right) \end{gathered} \nonumber$$

We need to solve $(2-i) v_{1}-5 v_{2}=0$. Thus,

$$\left(\begin{array}{l} v_{1} \\[4pt] v_{2} \end{array}\right)=\left(\begin{array}{c} 2+i \\[4pt] 1 \end{array}\right) \nonumber$$

Complex Solution: In order to get the two real linearly independent solutions, we need to compute the real and imaginary parts of $\mathbf{v} e^{\lambda t}$.

$$\begin{aligned} e^{\lambda t}\left(\begin{array}{c} 2+i \\[4pt] 1 \end{array}\right) &=e^{(1+i) t}\left(\begin{array}{c} 2+i \\[4pt] 1 \end{array}\right) \\[4pt] &=e^{t}(\cos t+i \sin t)\left(\begin{array}{c} 2+i \\[4pt] 1 \end{array}\right) \\[4pt] &=e^{t}\left(\begin{array}{c} (2+i)(\cos t+i \sin t) \\[4pt] \cos t+i \sin t \end{array}\right) \\[4pt] &=e^{t}\left(\begin{array}{c} (2 \cos t-\sin t)+i(\cos t+2 \sin t) \\[4pt] \cos t+i \sin t \end{array}\right)+i e^{t}\left(\begin{array}{c} \cos t+2 \sin t \\[4pt] \sin t \end{array}\right) \\[4pt] &=e^{t}\left(\begin{array}{c} 2 \cos t-\sin t \\[4pt] \cos t \end{array}\right) \end{aligned} \nonumber$$

General Solution: Now we can construct the general solution.

$$\begin{aligned} \mathbf{x}(t) &=c_{1} e^{t}\left(\begin{array}{l} 2 \cos t-\sin t \\[4pt] \cos t \end{array}\right)+c_{2} e^{t}\left(\begin{array}{c} \cos t+2 \sin t \\[4pt] \sin t \end{array}\right) \\[4pt] &=e^{t}\left(\begin{array}{c} c_{1}(2 \cos t-\sin t)+c_{2}(\cos t+2 \sin t) \\[4pt] c_{1} \cos t+c_{2} \sin t \end{array}\right) \end{aligned} \nonumber$$

Note: This can be rewritten as

$$\mathbf{x}(t)=e^{t} \cos t\left(\begin{array}{c} 2 c_{1}+c_{2} \\[4pt] c_{1} \end{array}\right)+e^{t} \sin t\left(\begin{array}{c} 2 c_{2}-c_{1} \\[4pt] c_{2} \end{array}\right) \nonumber$$

Example 2.16. $A=\left(\begin{array}{cc}7 & -1 \\[4pt] 9 & 1\end{array}\right)$.

Eigenvalues:

$$0=\left|\begin{array}{cc} 7-\lambda & -1 \\[4pt] 9 & 1-\lambda \end{array}\right| \nonumber$$

Therefore,

$$\begin{aligned} &0=(7-\lambda)(1-\lambda)+9 \\[4pt] &0=\lambda^{2}-8 \lambda+16 \\[4pt] &0=(\lambda-4)^{2} \end{aligned} \nonumber$$

There is only one real eigenvalue, $\lambda=4$. This is an example of Case II.

Eigenvectors: In this case we first solve for $\mathbf{v}_{1}$ and then get the second linearly independent vector.

$$\begin{aligned} &\left(\begin{array}{cc} 7 & -1 \\[4pt] 9 & 1 \end{array}\right)\left(\begin{array}{l} v_{1} \\[4pt] v_{2} \end{array}\right)=4\left(\begin{array}{l} v_{1} \\[4pt] v_{2} \end{array}\right) \\[4pt] &\left(\begin{array}{ll} 3 & -1 \\[4pt] 9 & -3 \end{array}\right)\left(\begin{array}{l} v_{1} \\[4pt] v_{2} \end{array}\right)=\left(\begin{array}{l} 0 \\[4pt] 0 \end{array}\right) \end{aligned} \nonumber$$

Therefore, we have

$$3 v_{1}-v_{2}=0, \quad \Rightarrow \quad\left(\begin{array}{l} v_{1} \\[4pt] v_{2} \end{array}\right)=\left(\begin{array}{l} 1 \\[4pt] 3 \end{array}\right) \text {. } \nonumber$$

Second Linearly Independent Solution:

Now we need to solve $A \mathbf{v}_{2}-\lambda \mathbf{v}_{2}=\mathbf{v}_{1}$.

Expanding the matrix product, we obtain the system of equations

$$\begin{array}{r} 3 u_{1}-u_{2}=1 \\[4pt] 9 u_{1}-3 u_{2}=3 . \end{array} \nonumber$$

The solution of this system is $\left(\begin{array}{l}u_{1} \\[4pt] u_{2}\end{array}\right)=\left(\begin{array}{l}1 \\[4pt] 2\end{array}\right)$.

General Solution: We construct the general solution as

$$\begin{aligned} \mathbf{y}(t) &=c_{1} e^{\lambda t} \mathbf{v}_{1}+c_{2} e^{\lambda t}\left(\mathbf{v}_{2}+t \mathbf{v}_{1}\right) \\[4pt] &=c_{1} e^{4 t}\left(\begin{array}{l} 1 \\[4pt] 3 \end{array}\right)+c_{2} e^{4 t}\left[\left(\begin{array}{l} 1 \\[4pt] 2 \end{array}\right)+t\left(\begin{array}{l} 1 \\[4pt] 3 \end{array}\right)\right] \\[4pt] &=e^{4 t}\left(\begin{array}{c} c_{1}+c_{2}(1+t) \\[4pt] 3 c_{1}+c_{2}(2+3 t) \end{array}\right) \end{aligned} \nonumber$$

$$\begin{aligned} & \left(\begin{array}{cc}7 & -1 \\[4pt]9 & 1\end{array}\right)\left(\begin{array}{l}u_{1} \\[4pt]u_{2}\end{array}\right)-4\left(\begin{array}{l}u_{1} \\[4pt]u_{2}\end{array}\right)=\left(\begin{array}{l}1 \\[4pt]3\end{array}\right) \\[4pt] & \left(\begin{array}{ll}3 & -1 \\[4pt]9 & -3\end{array}\right)\left(\begin{array}{l}u_{1} \\[4pt]u_{2}\end{array}\right)=\left(\begin{array}{l}1 \\[4pt]3\end{array}\right) \text {. } \end{aligned} \nonumber$$

Theory of Homogeneous Constant Coefficient Systems

There is a general theory for solving homogeneous, constant coefficient systems of first order differential equations. We begin by once again recalling the specific problem (2.12). We obtained the solution to this system as

$$\begin{gathered} x(t)=c_{1} e^{t}+c_{2} e^{-4 t}, \\[4pt] y(t)=\dfrac{1}{3} c_{1} e^{t}-\dfrac{1}{2} c_{2} e^{-4 t} \end{gathered} \nonumber$$

Table 2.2. Three examples of systems with a repeated root of $\lambda=2$.

This time we rewrite the solution as

$$\begin{aligned} \mathbf{x} &=\left(\begin{array}{c} c_{1} e^{t}+c_{2} e^{-4 t} \\[4pt] \dfrac{1}{3} c_{1} e^{t}-\dfrac{1}{2} c_{2} e^{-4 t} \end{array}\right) \\[4pt] &=\left(\begin{array}{cc} e^{t} & e^{-4 t} \\[4pt] \dfrac{1}{3} e^{t}-\dfrac{1}{2} e^{-4 t} \end{array}\right)\left(\begin{array}{c} c_{1} \\[4pt] c_{2} \end{array}\right) \\[4pt] & \equiv \Phi(t) \mathbf{C} \end{aligned} \nonumber$$

Thus, we can write the general solution as a $2 \times 2$ matrix $\Phi$ times an arbitrary constant vector. The matrix $\Phi$ consists of two columns that are linearly independent solutions of the original system. This matrix is an example of what we will define as the Fundamental Matrix of solutions of the system. So, determining the Fundamental Matrix will allow us to find the general solution of the system upon multiplication by a constant matrix. In fact, we will see that it will also lead to a simple representation of the solution of the initial value problem for our system. We will outline the general theory.

Consider the homogeneous, constant coefficient system of first order differential equations

$$\begin{aligned} \dfrac{d x_{1}}{d t} &=a_{11} x_{1}+a_{12} x_{2}+\ldots+a_{1 n} x_{n} \\[4pt] \dfrac{d x_{2}}{d t} &=a_{21} x_{1}+a_{22} x_{2}+\ldots+a_{2 n} x_{n} \\[4pt] & \vdots \\[4pt] \dfrac{d x_{n}}{d t} &=a_{n 1} x_{1}+a_{n 2} x_{2}+\ldots+a_{n n} x_{n} \end{aligned} \nonumber$$

As we have seen, this can be written in the matrix form $\mathbf{x}^{\prime}=A \mathbf{x}$, where

$$\mathbf{x}=\left(\begin{array}{c} x_{1} \\[4pt] x_{2} \\[4pt] \vdots \\[4pt] x_{n} \end{array}\right) \nonumber$$

and

$$A=\left(\begin{array}{cccc} a_{11} & a_{12} & \cdots & a_{1 n} \\[4pt] a_{21} & a_{22} & \cdots & a_{2 n} \\[4pt] \vdots & \vdots & \ddots & \vdots \\[4pt] a_{n 1} & a_{n 2} & \cdots & a_{n n} \end{array}\right) \nonumber$$

Now, consider $m$ vector solutions of this system: $\phi_{1}(t), \phi_{2}(t), \ldots \phi_{m}(t)$. These solutions are said to be linearly independent on some domain if

$$c_{1} \phi_{1}(t)+c_{2} \phi_{2}(t)+\ldots+c_{m} \phi_{m}(t)=0 \nonumber$$

for all $t$ in the domain implies that $c_{1}=c_{2}=\ldots=c_{m}=0$.

Let $\phi_{1}(t), \phi_{2}(t), \ldots \phi_{n}(t)$ be a set of $n$ linearly independent set of solutions of our system, called a fundamental set of solutions. We construct a matrix from these solutions using these solutions as the column of that matrix. We define this matrix to be the fundamental matrix solution. This matrix takes the form

$$\Phi=\left(\begin{array}{lll} \phi_{1} & \ldots & \phi_{n} \end{array}\right)=\left(\begin{array}{cccc} \phi_{11} & \phi_{12} & \cdots & \phi_{1 n} \\[4pt] \phi_{21} & \phi_{22} & \cdots & \phi_{2 n} \\[4pt] \vdots & \vdots & \ddots & \vdots \\[4pt] \phi_{n 1} & \phi_{n 2} & \cdots & \phi_{n n} \end{array}\right) \nonumber$$

What do we mean by a "matrix" solution? We have assumed that each $\phi_{k}$ is a solution of our system. Therefore, we have that $\phi_{k}^{\prime}=A \phi_{k}$, for $k=1, \ldots, n$. We say that $\Phi$ is a matrix solution because we can show that $\Phi$ also satisfies the matrix formulation of the system of differential equations. We can show this using the properties of matrices.

$$\begin{aligned} \dfrac{d}{d t} \Phi &=\left(\phi_{1}^{\prime} \ldots \phi_{n}^{\prime}\right) \\[4pt] &=\left(A \phi_{1} \ldots A \phi_{n}\right) \\[4pt] &=A\left(\phi_{1} \ldots \phi_{n}\right) \\[4pt] &=A \Phi \end{aligned} \nonumber$$

Given a set of vector solutions of the system, when are they linearly independent? We consider a matrix solution $\Omega(t)$ of the system in which we have $n$ vector solutions. Then, we define the Wronskian of $\Omega(t)$ to be

$$W=\operatorname{det} \Omega(t) \nonumber$$

If $W(t) \neq 0$, then $\Omega(t)$ is a fundamental matrix solution. Before continuing, we list the fundamental matrix solutions for the set of examples in the last section. (Refer to the solutions from those examples.) Furthermore, note that the fundamental matrix solutions are not unique as one can multiply any column by a nonzero constant and still have a fundamental matrix solution.

Example 2.14 $A=\left(\begin{array}{ll}4 & 2 \\[4pt] 3 & 3\end{array}\right)$.

$$\Phi(t)=\left(\begin{array}{cc} 2 e^{t} & e^{6 t} \\[4pt] -3 e^{t} & e^{6 t} \end{array}\right) \nonumber$$

We should note in this case that the Wronskian is found as

$$\begin{aligned} W &=\operatorname{det} \Phi(t) \\[4pt] &=\left|\begin{array}{cc} 2 e^{t} & e^{6 t} \\[4pt] -3 e^{t} & e^{6 t} \end{array}\right| \\[4pt] &=5 e^{7 t} \neq 0 . \end{aligned} \nonumber$$

Example 2.15 $A=\left(\begin{array}{ll}3 & -5 \\[4pt] 1 & -1\end{array}\right)$.

$$\Phi(t)=\left(\begin{array}{cc} e^{t}(2 \cos t-\sin t) & e^{t}(\cos t+2 \sin t) \\[4pt] e^{t} \cos t & e^{t} \sin t \end{array}\right) \nonumber$$

Example $2.16 A=\left(\begin{array}{cc}7 & -1 \\[4pt] 9 & 1\end{array}\right)$.

$$\Phi(t)=\left(\begin{array}{cc} e^{4 t} & e^{4 t}(1+t) \\[4pt] 3 e^{4 t} & e^{4 t}(2+3 t) \end{array}\right) \nonumber$$

So far we have only determined the general solution. This is done by the following steps:

Nonhomogeneous Systems

Before leaving the theory of systems of linear, constant coefficient systems, we will discuss nonhomogeneous systems. We would like to solve systems of the form

$$\begin{aligned} & \mathbf{x}(t)=\left(\begin{array}{cc}e^{-t} & e^{2 t} \\[4pt]-2 e^{-t} & -e^{2 t}\end{array}\right)\left(\begin{array}{cc}-1 & -1 \\[4pt]2 & 1\end{array}\right)\left(\begin{array}{l}1 \\[4pt]2\end{array}\right) \\[4pt] & =\left(\begin{array}{cc}e^{-t} & e^{2 t} \\[4pt]-2 e^{-t} & -e^{2 t}\end{array}\right)\left(\begin{array}{c}-3 \\[4pt]4\end{array}\right) \\[4pt] & =\left(\begin{array}{c}-3 e^{-t}+4 e^{2 t} \\[4pt]6 e^{-t}-4 e^{2 t}\end{array}\right) \end{aligned} \nonumber$$

$$\mathbf{x}^{\prime}=A(t) \mathbf{x}+\mathbf{f}(t) \nonumber$$

We will assume that we have found the fundamental matrix solution of the homogeneous equation. Furthermore, we will assume that $A(t)$ and $\mathbf{f}(t)$ are continuous on some common domain.

As with second order equations, we can look for solutions that are a sum of the general solution to the homogeneous problem plus a particular solution of the nonhomogeneous problem. Namely, we can write the general solution as

$$\mathbf{x}(t)=\Phi(t) \mathbf{C}+\mathbf{x}_{p}(t), \nonumber$$

where $\mathbf{C}$ is an arbitrary constant vector, $\Phi(t)$ is the fundamental matrix solution of $\mathbf{x}^{\prime}=A(t) \mathbf{x}$, and

$$\mathbf{x}_{p}^{\prime}=A(t) \mathbf{x}_{p}+\mathbf{f}(t) . \nonumber$$

Such a representation is easily verified.

We need to find the particular solution, $\mathbf{x}_{p}(t)$. We can do this by applying The Method of Variation of Parameters for Systems. We consider a solution in the form of the solution of the homogeneous problem, but replace the constant vector by unknown parameter functions. Namely, we assume that

$$\mathbf{x}_{p}(t)=\Phi(t) \mathbf{c}(t) . \nonumber$$

Differentiating, we have that

$$\mathbf{x}_{p}^{\prime}=\Phi^{\prime} \mathbf{c}+\Phi \mathbf{c}^{\prime}=A \Phi \mathbf{c}+\Phi \mathbf{c}^{\prime} \nonumber$$

$$\mathbf{x}_{p}^{\prime}-A \mathbf{x}_{p}=\Phi \mathbf{c}^{\prime} . \nonumber$$

But the left side is $\mathbf{f}$. So, we have that,

$$\Phi \mathbf{c}^{\prime}=\mathbf{f} \nonumber$$

or, since $\Phi$ is invertible (why?),

$$\mathbf{c}^{\prime}=\Phi^{-1} \mathbf{f} \nonumber$$

In principle, this can be integrated to give c. Therefore, the particular solution can be written as

$$\mathbf{x}_{p}(t)=\Phi(t) \int^{t} \Phi^{-1}(s) \mathbf{f}(s) d s . \nonumber$$

This is the variation of parameters formula.

The general solution of Equation (2.70) has been found as

$$\mathbf{x}(t)=\Phi(t) \mathbf{C}+\Phi(t) \int^{t} \Phi^{-1}(s) \mathbf{f}(s) d s . \nonumber$$

We can use the general solution to find the particular solution of an initial value problem consisting of Equation (2.70) and the initial condition $\mathbf{x}\left(t_{0}\right)=$ $\mathbf{x}_{0}$. This condition is satisfied for a solution of the form

$$\mathbf{x}(t)=\Phi(t) \mathbf{C}+\Phi(t) \int_{t_{0}}^{t} \Phi^{-1}(s) \mathbf{f}(s) d s \nonumber$$

provided

$$\mathbf{x}_{0}=\mathbf{x}\left(t_{0}\right)=\Phi\left(t_{0}\right) \mathbf{C} . \nonumber$$

This can be solved for $\mathbf{C}$ as in the last section. Inserting the solution back into the general solution $(2.73)$, we have

$$\mathbf{x}(t)=\Phi(t) \Phi^{-1}\left(t_{0}\right) \mathbf{x}_{0}+\Phi(t) \int_{t_{0}}^{t} \Phi^{-1}(s) \mathbf{f}(s) d s \nonumber$$

This solution can be written a little neater in terms of the principal matrix solution, $\Psi(t)=\Phi(t) \Phi^{-1}\left(t_{0}\right)$ :

$$\mathbf{x}(t)=\Psi(t) \mathbf{x}_{0}+\Psi(t) \int_{t_{0}}^{t} \Psi^{-1}(s) \mathbf{f}(s) d s \nonumber$$

Finally, one further simplification occurs when $A$ is a constant matrix, which are the only types of problems we have solved in this chapter. In this case, we have that $\Psi^{-1}(t)=\Psi(-t)$. So, computing $\Psi^{-1}(t)$ is relatively easy.

Example 2.18. $x^{\prime \prime}+x=2 \cos t, x(0)=4, x^{\prime}(0)=0$. This example can be solved using the Method of Undetermined Coefficients. However, we will use the matrix method described in this section.

First, we write the problem in matrix form. The system can be written as

$$\begin{gathered} x^{\prime}=y \\[4pt] y^{\prime}=-x+2 \cos t \end{gathered} \nonumber$$

Thus, we have a nonhomogeneous system of the form

$$\mathbf{x}^{\prime}=A \mathbf{x}+\mathbf{f}=\left(\begin{array}{cc} 0 & 1 \\[4pt] -1 & 0 \end{array}\right)\left(\begin{array}{l} x \\[4pt] y \end{array}\right)+\left(\begin{array}{c} 0 \\[4pt] 2 \cos t \end{array}\right) \nonumber$$

Next we need the fundamental matrix of solutions of the homogeneous problem. We have that

$$A=\left(\begin{array}{cc} 0 & 1 \\[4pt] -1 & 0 \end{array}\right) \text {. } \nonumber$$

The eigenvalues of this matrix are $\lambda=\pm i$. An eigenvector associated with $\lambda=i$ is easily found as $\left(\begin{array}{l}1 \\[4pt] i\end{array}\right)$. This leads to a complex solution

$$\left(\begin{array}{l} 1 \\[4pt] i \end{array}\right) e^{i t}=\left(\begin{array}{c} \cos t+i \sin t \\[4pt] i \cos t-\sin t \end{array}\right) . \nonumber$$

From this solution we can construct the fundamental solution matrix

$$\Phi(t)=\left(\begin{array}{cc} \cos t & \sin t \\[4pt] -\sin t & \cos t \end{array}\right) \nonumber$$

So, the general solution to the homogeneous problem is

$$\mathbf{x}_{h}=\Phi(t) \mathbf{C}=\left(\begin{array}{c} c_{1} \cos t+c_{2} \sin t \\[4pt] -c_{1} \sin t+c_{2} \cos t \end{array}\right) \nonumber$$

Next we seek a particular solution to the nonhomogeneous problem. From Equation (2.73) we see that we need $\Phi^{-1}(s) \mathbf{f}(s)$. Thus, we have

$$\begin{aligned} \Phi^{-1}(s) \mathbf{f}(s) &=\left(\begin{array}{cc} \cos s & -\sin s \\[4pt] \sin s & \cos s \end{array}\right)\left(\begin{array}{c} 0 \\[4pt] 2 \cos s \end{array}\right) \\[4pt] &=\left(\begin{array}{c} -2 \sin s \cos s \\[4pt] 2 \cos ^{2} s \end{array}\right) \end{aligned} \nonumber$$

We now compute

therefore, the general solution is

$$\mathbf{x}=\left(\begin{array}{c} c_{1} \cos t+c_{2} \sin t \\[4pt] -c_{1} \sin t+c_{2} \cos t \end{array}\right)+\left(\begin{array}{c} t \sin t \\[4pt] \sin t+t \cos t \end{array}\right) \nonumber$$

The solution to the initial value problem is

$$\mathbf{x}=\left(\begin{array}{cc} \cos t & \sin t \\[4pt] -\sin t & \cos t \end{array}\right)\left(\begin{array}{l} 4 \\[4pt] 0 \end{array}\right)+\left(\begin{array}{c} t \sin t \\[4pt] \sin t+t \cos t \end{array}\right) \nonumber$$

$$\mathbf{x}=\left(\begin{array}{c} 4 \cos t+t \sin t \\[4pt] -3 \sin t+t \cos t \end{array}\right) \nonumber$$

Spring-Mass Systems

There are many problems in physics that result in systems of equations. This is because the most basic law of physics is given by Newton’s Second Law, which states that if a body experiences a net force, it will accelerate. In particular, the net force is proportional to the acceleration with a proportionality constant called the mass, $m$. This is summarized as

$$\sum \mathbf{F}=m \mathbf{a} \nonumber$$

Since $\mathbf{a}=\ddot{\mathbf{x}}$, Newton’s Second Law is mathematically a system of second order differential equations for three dimensional problems, or one second order differential equation for one dimensional problems. If there are several masses, then we would naturally end up with systems no matter how many dimensions are involved.

A standard problem encountered in a first course in differential equations is that of a single block on a spring as shown in Figure 2.18. The net force in this case is the restoring force of the spring given by Hooke’s Law,

$$F_{s}=-k x \nonumber$$

where $k>0$ is the spring constant. Here $x$ is the elongation of the spring, or the displacement of the block from equilibrium. When $x$ is positive, the spring force is negative and when $x$ is negative the spring force is positive. We have depicted a horizontal system sitting on a frictionless surface.

A similar model can be provided for vertically oriented springs. Place the block on a vertically hanging spring. It comes to equilibrium, stretching the spring by $\ell_{0}$. Newton’s Second Law gives

$$-m g+k \ell_{0}=0 . \nonumber$$

Now, pulling the mass further by $x_{0}$, and releasing it, the mass begins to oscillate. Letting $x$ be the displacement from the new equilibrium, Newton’s Second Law now gives $m \ddot{x}=-m g+k\left(\ell_{0}-x\right)=-k x$.

In both examples (a horizontally or vetically oscillating mass) Newton’s Second Law of motion reults in the differential equation

$$m \ddot{x}+k x=0 . \nonumber$$

This is the equation for simple harmonic motion which we have already encountered in Chapter $1 .$

This second order equation can be written as a system of two first order equations.

$$\begin{aligned} &x^{\prime}=y \\[4pt] &y^{\prime}=-\dfrac{k}{m} x \end{aligned} \nonumber$$

The coefficient matrix for this system is

$$A=\left(\begin{array}{cc} 0 & 1 \\[4pt] -\omega^{2} & 0 \end{array}\right) \nonumber$$

where $\omega^{2}=\dfrac{k}{m}$. The eigenvalues of this system are $\lambda=\pm i \omega$ and the solutions are simple sines and cosines,

$$\begin{aligned} &x(t)=c_{1} \cos \omega t+c_{2} \sin \omega t, \\[4pt] &y(t)=\omega\left(-c_{1} \sin \omega t+c_{2} \cos \omega t\right) . \end{aligned} \nonumber$$

We further note that $\omega$ is called the angular frequency of oscillation and is given in $\mathrm{rad} / \mathrm{s}$. The frequency of oscillation is

$$f=\dfrac{\omega}{2 \pi} \nonumber$$

It typically has units of $\mathrm{s}^{-1}$, cps, or Hz. The multiplicative inverse has units of time and is called the period,

$$T=\dfrac{1}{f} . \nonumber$$

Thus, the period of oscillation for a mass $m$ on a spring with spring constant $k$ is given by

$$T=2 \pi \sqrt{\dfrac{m}{k}} \nonumber$$

Of course, we did not need to convert the last problem into a system. In fact, we had seen this equation in Chapter 1 . However, when one considers

more complicated spring-mass systems, systems of differential equations occur naturally. Consider two blocks attached with two springs as shown in Figure 2.19. In this case we apply Newton’s second law for each block.

First, consider the forces acting on the first block. The first spring is stretched by $x_{1}$. This gives a force of $F_{1}=-k_{1} x_{1}$. The second spring may also exert a force on the block depending if it is stretched, or not. If both blocks are displaced by the same amount, then the spring is not displaced. So, the amount by which the spring is displaced depends on the relative displacements of the two masses. This results in a second force of $F_{2}=k_{2}\left(x_{2}-x_{1}\right)$.

There is only one spring connected to mass two. Again the force depends on the relative displacement of the masses. It is just oppositely directed to the force which mass one feels from this spring.

Combining these forces and using Newton’s Second Law for both masses, we have the system of second order differential equations

$$\begin{aligned} &m_{1} \ddot{x}_{1}=-k_{1} x_{1}+k_{2}\left(x_{2}-x_{1}\right) \\[4pt] &m_{2} \ddot{x}_{2}=-k_{2}\left(x_{2}-x_{1}\right) \end{aligned} \nonumber$$

One can rewrite this system of two second order equations as a system of four first order equations. This is done by introducing two new variables $x_{3}=\dot{x}_{1}$ and $x_{4}=\dot{x}_{2}$. Note that these physically are the velocities of the two blocks.

The resulting system of first order equations is given as

$$\begin{aligned} \dot{x}_{1} &=x_{3} \\[4pt] \dot{x}_{2} &=x_{4} \\[4pt] \dot{x}_{3} &=-\dfrac{k_{1}}{m_{1}} x_{1}+\dfrac{k_{2}}{m_{1}}\left(x_{2}-x_{1}\right) \\[4pt] \dot{x}_{4} &=-\dfrac{k_{2}}{m_{2}}\left(x_{2}-x_{1}\right) \end{aligned} \nonumber$$

We can write our new system in matrix form as

$$\left(\begin{array}{c} \dot{x}_{1} \\[4pt] \dot{x}_{2} \\[4pt] \dot{x}_{3} \\[4pt] \dot{x}_{4} \end{array}\right)=\left(\begin{array}{cccc} 0 & 0 & 1 & 0 \\[4pt] 0 & 0 & 0 & 1 \\[4pt] -\dfrac{k_{1}+k_{2}}{m_{1}} & \dfrac{k_{2}}{m_{1}} & 0 & 0 \\[4pt] \dfrac{k_{2}}{m_{2}} & -\dfrac{k_{2}}{m_{2}} & 0 & 0 \end{array}\right)\left(\begin{array}{l} x_{1} \\[4pt] x_{2} \\[4pt] x_{3} \\[4pt] x_{4} \end{array}\right) \nonumber$$

Electrical Circuits

Another problem often encountered in a first year physics class is that of an LRC series circuit. This circuit is pictured in Figure $2.20$. The resistor is a circuit element satisfying Ohm’s Law. The capacitor is a device that stores electrical energy and an inductor, or coil, stores magnetic energy.

The physics for this problem stems from Kirchoff’s Rules for circuits. Since there is only one loop, we will only need Kirchoff’s Loop Rule. Namely, the sum of the drops in electric potential are set equal to the rises in electric potential. The potential drops across each circuit element are given by

1. Resistor: $V_{R}=I R$.
2. Capacitor: $V_{C}=\dfrac{q}{C}$.
3. Inductor: $V_{L}=L \dfrac{d I}{d t}$.

Adding these potential drops and setting the sum equal to the voltage supplied by the voltage source, $V(t)$, we obtain

$$I R+\dfrac{q}{C}+L \dfrac{d I}{d t}=V(t) \nonumber$$

Furthermore, we recall that the current is defined as $I=\dfrac{d q}{d t}$. where $q$ is the charge in the circuit. Since both $q$ and $I$ are unknown, we can replace the current by its expression in terms of the charge to obtain

$$L \ddot{q}+R \dot{q}+\dfrac{1}{C} q=V(t) . \nonumber$$

This is a second order differential equation for $q(t)$. One can set up a system of equations and proceed to solve them. However, this is a constant coefficient differential equation and can also be solved using the methods in Chapter $1 .$

In the next examples we will look at special cases that arise for the series LRC circuit equation. These include $R C$ circuits, solvable by first order methods and $L C$ circuits, leading to oscillatory behavior.

Love Affairs

The next application is one that has been studied by several authors as a cute system involving relationships. One considers what happens to the affections that two people have for each other over time. Let $R$ denote the affection that Romeo has for Juliet and $J$ be the affection that Juliet has for Romeo. positive values indicate love and negative values indicate dislike.

One possible model is given by

$$\begin{aligned} &\dfrac{d R}{d t}=b J \\[4pt] &\dfrac{d J}{d t}=c R \end{aligned} \nonumber$$

with $b>0$ and $c<0$. In this case Romeo loves Juliet the more she likes him. But Juliet backs away when she finds his love for her increasing.

A typical system relating the combined changes in affection can be modeled

$$\begin{aligned} &\dfrac{d R}{d t}=a R+b J \\[4pt] &\dfrac{d J}{d t}=c R+d J \end{aligned} \nonumber$$

Several scenarios are possible for various choices of the constants. For example, if $a>0$ and $b>0$, Romeo gets more and more excited by Juliet’s love for him. If $c>0$ and $d<0$, Juliet is being cautious about her relationship with Romeo. For specific values of the parameters and initial conditions, one can explore this match of an overly zealous lover with a cautious lover.

Predator Prey Models

Another common model studied is that of competing species. For example, we could consider a population of rabbits and foxes. Left to themselves, rabbits would tend to multiply, thus

$$\dfrac{d R}{d t}=a R, \nonumber$$

with $a>0$. In such a model the rabbit population would grow exponentially. Similarly, a population of foxes would decay without the rabbits to feed on. So, we have that

$$\dfrac{d F}{d t}=-b F \nonumber$$

for $b>0$.

Now, if we put these populations together on a deserted island, they would interact. The more foxes, the rabbit population would decrease. However, the more rabbits, the foxes would have plenty to eat and the population would thrive. Thus, we could model the competing populations as

$$\begin{gathered} \dfrac{d R}{d t}=a R-c F, \\[4pt] \dfrac{d F}{d t}=-b F+d R, \end{gathered} \nonumber$$

where all of the constants are positive numbers. Studying this coupled system would lead to as study of the dynamics of these populations. We will discuss other (nonlinear) systems in the next chapter.

Mixture Problems

There are many types of mixture problems. Such problems are standard in a first course on differential equations as examples of first order differential equations. Typically these examples consist of a tank of brine, water containing a specific amount of salt with pure water entering and the mixture leaving, or the flow of a pollutant into, or out of, a lake.

In general one has a rate of flow of some concentration of mixture entering a region and a mixture leaving the region. The goal is to determine how much stuff is in the region at a given time. This is governed by the equation

$$\text { Rate of change of substance }=\text { Rate In }-\text { Rate Out. } \nonumber$$

This can be generalized to the case of two interconnected tanks. We provide some examples.

Chemical Kinetics

There are many problems that come from studying chemical reactions. The simplest reaction is when a chemical $A$ turns into chemical $B$. This happens at a certain rate, $k>0$. This can be represented by the chemical formula

In this case we have that the rates of change of the concentrations of $A,[A]$, and $B,[B]$, are given by

$$\begin{aligned} &\dfrac{d[A]}{d t}=-k[A] \\[4pt] &\dfrac{d[B]}{d t}=k[A] \end{aligned} \nonumber$$

Think about this as it is a key to understanding the next reactions.

A more complicated reaction is given by

$$A \underset{k_{1}}{\longrightarrow} B \underset{k_{2}}{\longrightarrow} C \text {. } \nonumber$$

In this case we can add to the above equation the rates of change of concentrations $[B]$ and $[C]$. The resulting system of equations is

$$\begin{aligned} \dfrac{d[A]}{d t} &=-k_{1}[A] \\[4pt] \dfrac{d[B]}{d t} &=k_{1}[A]-k_{2}[B] \\[4pt] \dfrac{d[C]}{d t} &=k_{2}[B] \end{aligned} \nonumber$$

One can further consider reactions in which a reverse reaction is possible. Thus, a further generalization occurs for the reaction

$$A \underset{k_{1}}{\stackrel{k_{3}}{\longrightarrow}} B \underset{k_{2}}{\longrightarrow} C \text {. } \nonumber$$

The resulting system of equations is

$$\begin{aligned} &\dfrac{d[A]}{d t}=-k_{1}[A]+k_{3}[B] \\[4pt] &\dfrac{d[B]}{d t}=k_{1}[A]-k_{2}[B]-k_{3}[B] \\[4pt] &\dfrac{d[C]}{d t}=k_{2}[B] \end{aligned} \nonumber$$

More complicated chemical reactions will be discussed at a later time.

Epidemics

Another interesting area of application of differential equation is in predicting the spread of disease. Typically, one has a population of susceptible people or animals. Several infected individuals are introduced into the population and one is interested in how the infection spreads and if the number of infected people drastically increases or dies off. Such models are typically nonlinear and we will look at what is called the SIR model in the next chapter. In this section we will model a simple linear model.

Let break the population into three classes. First, $S(t)$ are the healthy people, who are susceptible to infection. Let $I(t)$ be the number of infected people. Of these infected people, some will die from the infection and others recover. Let’s assume that initially there in one infected person and the rest, say $N$, are obviously healthy. Can we predict how many deaths have occurred by time t?

Let’s try and model this problem using the compartmental analysis we had seen in the mixing problems. The total rate of change of any population would be due to those entering the group less those leaving the group. For example, the number of healthy people decreases due infection and can increase when some of the infected group recovers. Let’s assume that the rate of infection is proportional to the number of healthy people, $a S$. Also, we assume that the number who recover is proportional to the number of infected, $r I$. Thus, the rate of change of the healthy people is found as

$$\dfrac{d S}{d t}=-a S+r I . \nonumber$$

Let the number of deaths be $D(t)$. Then, the death rate could be taken to be proportional to the number of infected people. So,

$$\dfrac{d D}{d t}=d I \nonumber$$

Finally, the rate of change of infectives is due to healthy people getting infected and the infectives who either recover or die. Using the corresponding terms in the other equations, we can write

$$\dfrac{d I}{d t}=a S-r I-d I . \nonumber$$

This linear system can be written in matrix form.

$$\dfrac{d}{d t}\left(\begin{array}{c} S \\[4pt] I \\[4pt] D \end{array}\right)=\left(\begin{array}{ccc} -a & r & 0 \\[4pt] a & -d-r & 0 \\[4pt] 0 & d & 0 \end{array}\right)\left(\begin{array}{c} S \\[4pt] I \\[4pt] D \end{array}\right) \nonumber$$

The eigenvalue equation for this system is

$$\lambda\left[\lambda^{2}+(a+r+d) \lambda+a d\right]=0 . \nonumber$$

The reader can find the solutions of this system and determine if this is a realistic model.

Appendix: Diagonalization and Linear Systems

As we have seen, the matrix formulation for linear systems can be powerful, especially for $n$ differential equations involving $n$ unknown functions. Our ability to proceed towards solutions depended upon the solution of eigenvalue problems. However, in the case of repeated eigenvalues we saw some additional complications. This all depends deeply on the background linear algebra. Namely, we relied on being able to diagonalize the given coefficient matrix. In this section we will discuss the limitations of diagonalization and introduce the Jordan canonical form.

We begin with the notion of similarity. Matrix $A$ is similar to matrix $B$ if and only if there exists a nonsingular matrix $P$ such that

$$B=P^{-1} A P . \nonumber$$

Recall that a nonsingular matrix has a nonzero determinant and is invertible.

We note that the similarity relation is an equivalence relation. Namely, it satisfies the following

1. $A$ is similar to itself.
2. If $A$ is similar to $B$, then $B$ is similar to $A$.
3. If $A$ is similar to $B$ and $B$ is similar to $C$, the $A$ is similar to $C$.

Also, if $A$ is similar to $B$, then they have the same eigenvalues. This follows from a simple computation of the eigenvalue equation. Namely,

$$\begin{aligned} 0 &=\operatorname{det}(B-\lambda I) \\[4pt] &=\operatorname{det}\left(P^{-1} A P-\lambda P^{-1} I P\right) \\[4pt] &=\operatorname{det}(P)^{-1} \operatorname{det}(A-\lambda I) \operatorname{det}(P) \\[4pt] &=\operatorname{det}(A-\lambda I) \end{aligned} \nonumber$$

Therefore, $\operatorname{det}(A-\lambda I)=0$ and $\lambda$ is an eigenvalue of both $A$ and $B$.

An $n \times n$ matrix $A$ is diagonalizable if and only if $A$ is similar to a diagonal matrix $D$; i.e., there exists a nonsingular matrix $P$ such that

$$D=P^{-1} A P . \nonumber$$

One of the most important theorems in linear algebra is the Spectral Theorem. This theorem tells us when a matrix can be diagonalized. In fact, it goes beyond matrices to the diagonalization of linear operators. We learn in linear algebra that linear operators can be represented by matrices once we pick a particular representation basis. Diagonalization is simplest for finite dimensional vector spaces and requires some generalization for infinite dimensional vectors spaces. Examples of operators to which the spectral theorem applies are self-adjoint operators (more generally normal operators on Hilbert spaces). We will explore some of these ideas later in the course. The spectral theorem provides a canonical decomposition, called the spectral decomposition, or eigendecomposition, of the underlying vector space on which it acts.

The next theorem tells us how to diagonalize a matrix:

Theorem 2.23. Let $A$ be an $n \times n$ matrix. Then $A$ is diagonalizable if and only if $A$ has $n$ linearly independent eigenvectors. If so, then

$$D=P^{-1} A P . \nonumber$$

If $\left\{v_{1}, \ldots, v_{n}\right\}$ are the eigenvectors of $A$ and $\left\{\lambda_{1}, \ldots, \lambda_{n}\right\}$ are the corresponding eigenvalues, then $v_{j}$ is the jth column of $P$ and $D_{j j}=\lambda_{j}$.

A simpler determination results by noting

Theorem 2.24. Let $A$ be an $n \times n$ matrix with $n$ real and distinct eigenvalues. Then $A$ is diagonalizable.

Therefore, we need only look at the eigenvalues and determine diagonalizability. In fact, one also has from linear algebra the following result.

Theorem 2.25. Let $A$ be an $n \times n$ real symmetric matrix. Then $A$ is diagonalizable.

Recall that a symmetric matrix is one whose transpose is the same as the matrix, or $A_{i j}=A_{j i}$.

Example 2.26. Consider the matrix

$$A=\left(\begin{array}{lll} 1 & 2 & 2 \\[4pt] 2 & 3 & 0 \\[4pt] 2 & 0 & 3 \end{array}\right) \nonumber$$

This is a real symmetric matrix. The characteristic polynomial is found to be

$$\operatorname{det}(A-\lambda I)=-(\lambda-5)(\lambda-3)(\lambda+1)=0 \nonumber$$

As before, we can determine the corresponding eigenvectors (for $\lambda=-1,3,5$, respectively) as

$$\left(\begin{array}{c} -2 \\[4pt] 1 \\[4pt] 1 \end{array}\right), \quad\left(\begin{array}{c} 0 \\[4pt] -1 \\[4pt] 1 \end{array}\right), \quad\left(\begin{array}{l} 1 \\[4pt] 1 \\[4pt] 1 \end{array}\right) \text {. } \nonumber$$

We can use these to construct the diagonalizing matrix $P$. Namely, we have

$$P^{-1} A P=\left(\begin{array}{ccc} -2 & 0 & 1 \\[4pt] 1 & -1 & 1 \\[4pt] 1 & 1 & 1 \end{array}\right)^{-1}\left(\begin{array}{lll} 1 & 2 & 2 \\[4pt] 2 & 3 & 0 \\[4pt] 2 & 0 & 3 \end{array}\right)\left(\begin{array}{ccc} -2 & 0 & 1 \\[4pt] 1 & -1 & 1 \\[4pt] 1 & 1 & 1 \end{array}\right)=\left(\begin{array}{ccc} -1 & 0 & 0 \\[4pt] 0 & 3 & 0 \\[4pt] 0 & 0 & 5 \end{array}\right) \nonumber$$

Now diagonalization is an important idea in solving linear systems of first order equations, as we have seen for simple systems. If our system is originally diagonal, that means our equations are completely uncoupled. Let our system take the form

$$\dfrac{d \mathbf{y}}{d t}=D \mathbf{y} \nonumber$$

where $D$ is diagonal with entries $\lambda_{i}, i=1, \ldots, n$. The system of equations, $y_{i}^{\prime}=\lambda_{i} y_{i}$, has solutions

$$y_{i}(t)=c_{c} e^{\lambda_{i} t} \nonumber$$

Thus, it is easy to solve a diagonal system.

Let $A$ be similar to this diagonal matrix. Then

$$\dfrac{d \mathbf{y}}{d t}=P^{-1} A P \mathbf{y} \nonumber$$

This can be rewritten as

$$\dfrac{d P \mathbf{y}}{d t}=A P \mathbf{y} \nonumber$$

Defining $\mathbf{x}=P \mathbf{y}$, we have

$$\dfrac{d \mathbf{x}}{d t}=A \mathbf{x} \nonumber$$

This simple derivation shows that if $A$ is diagonalizable, then a transformation of the original system in $\mathbf{x}$ to new coordinates, or a new basis, results in a simpler system in $\mathbf{y}$.

However, it is not always possible to diagonalize a given square matrix. This is because some matrices do not have enough linearly independent vectors, or we have repeated eigenvalues. However, we have the following theorem:

Theorem 2.27. Every $n \times n$ matrix $A$ is similar to a matrix of the form

$$J=\operatorname{diag}\left[J_{1}, J_{2}, \ldots, J_{n}\right] \nonumber$$

where

$$J_{i}=\left(\begin{array}{ccccc} \lambda_{i} & 1 & 0 & \cdots & 0 \\[4pt] 0 & \lambda_{i} & 1 & \cdots & 0 \\[4pt] \vdots & \ddots & \ddots & \ddots & \vdots \\[4pt] 0 & \cdots & 0 & \lambda_{i} & 1 \\[4pt] 0 & 0 & \cdots & 0 & \lambda_{i} \end{array}\right) \nonumber$$

We will not go into the details of how one finds this Jordan Canonical Form or proving the theorem. In practice you can use a computer algebra system to determine this and the similarity matrix. However, we would still need to know how to use it to solve our system of differential equations. Example 2.28. Let’s consider a simple system with the $3 \times 3$ Jordan block

$$A=\left(\begin{array}{lll} 2 & 1 & 0 \\[4pt] 0 & 2 & 1 \\[4pt] 0 & 0 & 2 \end{array}\right) \nonumber$$

The corresponding system of coupled first order differential equations takes the form

$$\begin{aligned} &\dfrac{d x_{1}}{d t}=2 x_{1}+x_{2}, \\[4pt] &\dfrac{d x_{2}}{d t}=2 x_{2}+x_{3}, \\[4pt] &\dfrac{d x_{3}}{d t}=2 x_{3} . \end{aligned} \nonumber$$

The last equation is simple to solve, giving $x_{3}(t)=c_{3} e^{2 t}$. Inserting into the second equation, you have a

$$\dfrac{d x_{2}}{d t}=2 x_{2}+c_{3} e^{2 t} \nonumber$$

Using the integrating factor, $e^{-2 t}$, one can solve this equation to get $x_{2}(t)=$ $\left(c_{2}+c_{3} t\right) e^{2 t}$. Similarly, one can solve the first equation to obtain $x_{1}(t)=$ $\left(c_{1}+c_{2} t+\dfrac{1}{2} c_{3} t^{2}\right) e^{2 t}$

This should remind you of a problem we had solved earlier leading to the generalized eigenvalue problem in (2.43). This suggests that there is a more general theory when there are multiple eigenvalues and relating to Jordan canonical forms.

Let’s write the solution we just obtained in vector form. We have

$$\mathbf{x}(t)=\left[c_{1}\left(\begin{array}{l} 1 \\[4pt] 0 \\[4pt] 0 \end{array}\right)+c_{2}\left(\begin{array}{l} t \\[4pt] 1 \\[4pt] 0 \end{array}\right)+c_{3}\left(\begin{array}{c} \dfrac{1}{2} t^{2} \\[4pt] t \\[4pt] 1 \end{array}\right)\right] e^{2 t} \nonumber$$

It looks like this solution is a linear combination of three linearly independent solutions,

$$\begin{aligned} &\mathbf{x}=\mathbf{v}_{1} e^{2 \lambda t} \\[4pt] &\mathbf{x}=\left(t \mathbf{v}_{1}+\mathbf{v}_{2}\right) e^{\lambda t} \\[4pt] &\mathbf{x}=\left(\dfrac{1}{2} t^{2} \mathbf{v}_{1}+t \mathbf{v}_{2}+\mathbf{v}_{3}\right) e^{\lambda t} \end{aligned} \nonumber$$

where $\lambda=2$ and the vectors satisfy the equations

$$\begin{aligned} &(A-\lambda I) \mathbf{v}_{1}=0 \\[4pt] &(A-\lambda I) \mathbf{v}_{2}=\mathbf{v}_{1} \\[4pt] &(A-\lambda I) \mathbf{v}_{3}=\mathbf{v}_{2} \end{aligned} \nonumber$$

and

$$\begin{aligned} (A-\lambda I) \mathbf{v}_{1} &=0 \\[4pt] (A-\lambda I)^{2} \mathbf{v}_{2} &=0 \\[4pt] (A-\lambda I)^{3} \mathbf{v}_{3} &=0 \end{aligned} \nonumber$$

It is easy to generalize this result to build linearly independent solutions corresponding to multiple roots (eigenvalues) of the characteristic equation.

Introduction

Most of your studies of differential equations to date have been the study linear differential equations and common methods for solving them. However, the real world is very nonlinear. So, why study linear equations? Because they are more readily solved. As you may recall, we can use the property of linear superposition of solutions of linear differential equations to obtain general solutions. We will see that we can sometimes approximate the solutions of nonlinear systems with linear systems in small regions of phase space.

In general, nonlinear equations cannot be solved obtaining general solutions. However, we can often investigate the behavior of the solutions without actually being able to find simple expressions in terms of elementary functions. When we want to follow the evolution of these solutions, we resort to numerically solving our differential equations. Such numerical methods need to be executed with care and there are many techniques that can be used. We will not go into these techniques in this course. However, we can make use of computer algebra systems, or computer programs, already developed for obtaining such solutions.

Nonlinear problems occur naturally. We will see problems from many of the same fields we explored in Section 2.9. One example is that of population dynamics. Typically, we have a certain population, $y(t)$, and the differential equation governing the growth behavior of this population is developed in a manner similar to that used previously for mixing problems. We note that the rate of change of the population is given by the Rate In minus the Rate Out. The Rate In is given by the number of the species born per unit time. The Rate Out is given by the number that die per unit time.

A simple population model can be obtained if one assumes that these rates are linear in the population. Thus, we assume that the Rate In $=b y$ and the Rate Out $=m y$. Here we have denoted the birth rate as $b$ and the mortality rate as $m$, . This gives the rate of change of population as

$$\dfrac{d y}{d t}=b y-m y \nonumber$$

Generally, these rates could depend upon time. In the case that they are both constant rates, we can define $k=b-m$ and we obtain the familiar exponential model:

$$\dfrac{d y}{d t}=k y . \nonumber$$

This is easily solved and one obtains exponential growth $(k>0)$ or decay $(k<$ $0)$. This model has been named after Malthus 1 , a clergyman who used this model to warn of the impending doom of the human race if its reproductive practices continued.

However, when populations get large enough, there is competition for resources, such as space and food, which can lead to a higher mortality rate. Thus, the mortality rate may be a function of the population size, $m=m(y)$. The simplest model would be a linear dependence, $m=\tilde{m}+c y$. Then, the previous exponential model takes the form

$$\dfrac{d y}{d t}=k y-c y^{2} \nonumber$$

This is known as the logistic model of population growth. Typically, $c$ is small and the added nonlinear term does not really kick in until the population gets large enough.

While one can solve this particular equation, it is instructive to study the qualitative behavior of the solutions without actually writing down the explicit solutions. Such methods are useful for more difficult nonlinear equations. We will investigate some simple first order equations in the next section. In the following section we present the analytic solution for completeness.

We will resume our studies of systems of equations and various applications throughout the rest of this chapter. We will see that we can get quite a bit of information about the behavior of solutions by using some of our earlier methods for linear systems.

Autonomous First Order Equations

In this section we will review the techniques for studying the stability of nonlinear first order autonomous equations. We will then extend this study to looking at families of first order equations which are connected through a parameter.

Recall that a first order autonomous equation is given in the form

${ }^{1}$ Malthus, Thomas Robert. An Essay on the Principle of Population. Library of Economics and Liberty. Retrieved August 2, 2007 from the World Wide Web: http://www.econlib.org/library/Malthus/malPop1.html

$$\dfrac{d y}{d t}=f(y) . \nonumber$$

We will assume that $f$ and $\dfrac{\partial f}{\partial y}$ are continuous functions of $y$, so that we know that solutions of initial value problems exist and are unique.

We will recall the qualitative methods for studying autonomous equations by considering the example

$$\dfrac{d y}{d t}=y-y^{2} . \nonumber$$

This is just an example of a logistic equation.

First, one determines the equilibrium, or constant, solutions given by $y^{\prime}=$ 0 . For this case, we have $y-y^{2}=0$. So, the equilibrium solutions are $y=0$ and $y=1$. Sketching these solutions, we divide the $t y$-plane into three regions. Solutions that originate in one of these regions at $t=t_{0}$ will remain in that region for all $t>t_{0}$ since solutions cannot intersect. [Note that if two solutions intersect then they have common values $y_{1}$ at time $t_{1}$. Using this information, we could set up an initial value problem for which the initial condition is $y\left(t_{1}\right)=y_{1}$. Since the two different solutions intersect at this point in the phase plane, we would have an initial value problem with two different solutions corresponding to the same initial condition. This contradicts the uniqueness assumption stated above. We will leave the reader to explore this further in the homework.]

Next, we determine the behavior of solutions in the three regions. Noting that $d y / d t$ gives the slope of any solution in the plane, then we find that the solutions are monotonic in each region. Namely, in regions where $d y / d t>0$, we have monotonically increasing functions. We determine this from the right side of our equation.

For example, in this problem $y-y^{2}>0$ only for the middle region and $y-y^{2}<0$ for the other two regions. Thus, the slope is positive in the middle region, giving a rising solution as shown in Figure 3.1. Note that this solution does not cross the equilibrium solutions. Similar statements can be made about the solutions in the other regions.

We further note that the solutions on either side of $y=1$ tend to approach this equilibrium solution for large values of $t$. In fact, no matter how close one is to $y=1$, eventually one will approach this solution as $t \rightarrow \infty$. So, the equilibrium solution is a stable solution. Similarly, we see that $y=0$ is an unstable equilibrium solution.

If we are only interested in the behavior of the equilibrium solutions, we could just construct a phase line. In Figure $3.2$ we place a vertical line to the right of the $t y$-plane plot. On this line one first places dots at the corresponding equilibrium solutions and labels the solutions. These points at the equilibrium solutions are end points for three intervals. In each interval one then places arrows pointing upward (downward) indicating solutions with positive (negative) slopes. Looking at the phase line one can now determine if a given equilibrium is stable (arrows pointing towards the point) or unstable

(arrows pointing away from the point). In Figure $3.3$ we draw the final phase line by itself.

Solution of the Logistic Equation

We have seen that one does not need an explicit solution of the logistic equation (3.2) in order to study the behavior of its solutions. However, the logistic equation is an example of a nonlinear first order equation that is solvable. It is an example of a Riccati equation.

The general form of the Riccati equation is

$$\dfrac{d y}{d t}=a(t)+b(t) y+c(t) y^{2} \nonumber$$

As long as $c(t) \neq 0$, this equation can be reduced to a second order linear differential equation through the transformation

$$y(t)=-\dfrac{1}{c(t)} \dfrac{\dot{x}(t)}{x(t)} . \nonumber$$

We will demonstrate this using the simple case of the logistic equation,

$$\dfrac{d y}{d t}=k y-c y^{2} . \nonumber$$

We let

$$y(t)=\dfrac{1}{c} \dfrac{\dot{x}}{x} \nonumber$$

Then

$$\begin{aligned} \dfrac{d y}{d t} &=\dfrac{1}{c}\left[\dfrac{\ddot{x}}{x}-\left(\dfrac{\dot{x}}{x}\right)^{2}\right] \\[4pt] &=\dfrac{1}{c}\left[\dfrac{\ddot{x}}{x}-(c y)^{2}\right] \\[4pt] &=\dfrac{1}{c} \dfrac{\ddot{x}}{x}-c y^{2} \end{aligned} \nonumber$$

Inserting this into the logistic equation (3.5), we have

$$\dfrac{1}{c} \dfrac{\ddot{x}}{x}-c y^{2}=k \dfrac{1}{c}\left(\dfrac{\dot{x}}{x}\right)-c y^{2}, \nonumber$$

or

$$\ddot{x}=k \dot{x} . \nonumber$$

This equation is readily solved to give

$$x(t)=A+B e^{k t} . \nonumber$$

Therefore, we have the solution to the logistic equation is

$$y(t)=\dfrac{1}{c} \dfrac{\dot{x}}{x}=\dfrac{k B e^{k t}}{c\left(A+B e^{k t}\right)} \nonumber$$

It appears that we have two arbitrary constants. But, we started out with a first order differential equation and expect only one arbitrary constant. However, we can resolve this by dividing the numerator and denominator by $k B e^{k t}$ and defining $C=\dfrac{A}{B}$. Then we have

$$y(t)=\dfrac{k / c}{1+C e^{-k t}}, \nonumber$$

showing that there really is only one arbitrary constant in the solution.

We should note that this is not the only way to obtain the solution to the logistic equation, though it does provide an introduction to Riccati equations. A more direct approach would be to use separation of variables on the logistic equation. The reader should verify this.

Nonlinear Pendulum

In this section we will introduce the nonlinear pendulum as our first example of periodic motion in a nonlinear system. Oscillations are important in many areas of physics. We have already seen the motion of a mass on a spring, leading to simple, damped, and forced harmonic motions. Later we will explore these effects on a simple nonlinear system. In this section we will introduce the nonlinear pendulum and determine its period of oscillation.

We begin by deriving the pendulum equation. The simple pendulum consists of a point mass $m$ hanging on a string of length $L$ from some support. [See Figure 3.11.] One pulls the mass back to some starting angle, $\theta_{0}$, and releases it. The goal is to find the angular position as a function of time, $\theta(t)$.

There are a couple of derivations possible. We could either use Newton’s Second Law of Motion, $F=m a$, or its rotational analogue in terms of torque. We will use the former only to limit the amount of physics background needed.

There are two forces acting on the point mass, the weight and the tension in the string. The weight points downward and has a magnitude of $m g$, where $g$ is the standard symbol for the acceleration due to gravity. At the surface of the earth we can take this to be $9.8 \mathrm{~m} / \mathrm{s}^{2}$ or $32.2 \mathrm{ft} / \mathrm{s}^{2}$. In Figure $3.12$ we show both the weight and the tension acting on the mass. The net force is also shown.

The tension balances the projection of the weight vector, leaving an unbalanced component of the weight in the direction of the motion. Thus, the magnitude of the sum of the forces is easily found from this unbalanced component as $F=m g \sin \theta$.

Newton’s Second Law of Motion tells us that the net force is the mass times the acceleration. So, we can write

$$m \ddot{x}=-m g \sin \theta . \nonumber$$

Next, we need to relate $x$ and $\theta . x$ is the distance traveled, which is the length of the arc traced out by our point mass. The arclength is related to the angle, provided the angle is measured in radians. Namely, $x=r \theta$ for $r=L$. Thus, we can write

$$m L \ddot{\theta}=-m g \sin \theta \nonumber$$

Canceling the masses, leads to the nonlinear pendulum equation

$$L \ddot{\theta}+g \sin \theta=0 . \nonumber$$

There are several variations of Equation (3.8) which will be used in this text. The first one is the linear pendulum. This is obtained by making a small angle approximation. For small angles we know that $\sin \theta \approx \theta$. Under this approximation (3.8) becomes

$$L \ddot{\theta}+g \theta=0 . \nonumber$$

We can also make the system more realistic by adding damping. This could be due to energy loss in the way the string is attached to the support or due to the drag on the mass, etc. Assuming that the damping is proportional to the angular velocity, we have equations for the damped nonlinear and damped linear pendula:

$$\begin{gathered} L \ddot{\theta}+b \dot{\theta}+g \sin \theta=0 . \\[4pt] L \ddot{\theta}+b \dot{\theta}+g \theta=0 . \end{gathered} \nonumber$$

Finally, we can add forcing. Imagine that the support is attached to a device to make the system oscillate horizontally at some frequency. Then we could have equations such as

$$L \ddot{\theta}+b \dot{\theta}+g \sin \theta=F \cos \omega t . \nonumber$$

We will look at these and other oscillation problems later in the exercises. These are summarized in the table below.

In Search of Solutions

Before returning to studying the equilibrium solutions of the nonlinear pendulum, we will look at how far we can get at obtaining analytical solutions. First, we investigate the simple linear pendulum.

The linear pendulum equation (3.9) is a constant coefficient second order linear differential equation. The roots of the characteristic equations are $r=$ $\pm \sqrt{\dfrac{g}{L}} i$. Thus, the general solution takes the form

$$\theta(t)=c_{1} \cos \left(\sqrt{\dfrac{g}{L}} t\right)+c_{2} \sin \left(\sqrt{\dfrac{g}{L}} t\right) \nonumber$$

We note that this is usually simplified by introducing the angular frequency

$$\omega \equiv \sqrt{\dfrac{g}{L}} . \nonumber$$

One consequence of this solution, which is used often in introductory physics, is an expression for the period of oscillation of a simple pendulum. REcall that the period is the time it takes to complete one cycle of the oscillation. The period is found to be

$$T=\dfrac{2 \pi}{\omega}=2 \pi \sqrt{\dfrac{L}{g}} \nonumber$$

This value for the period of a simple pendulum is based on the linear pendulum equation, which was derived assuming a small angle approximation. How good is this approximation? What is meant by a small angle? We recall the Taylor series approximation of $\sin \theta$ about $\theta=0$ :

$$\sin \theta=\theta-\dfrac{\theta^{3}}{3 !}+\dfrac{\theta^{5}}{5 !}+\ldots \nonumber$$

One can obtain a bound on the error when truncating this series to one term after taking a numerical analysis course. But we can just simply plot the relative error, which is defined as Relative Error $=\left|\dfrac{\sin \theta-\theta}{\sin \theta}\right| \times 100 \%$.

A plot of the relative error is given in Figure 3.13. We note that a one percent relative error corresponds to about $0.24$ radians, which is less that fourteen degrees. Further discussion on this is provided at the end of this section.

We now turn to the nonlinear pendulum. We first rewrite Equation (3.8) in the simpler form

$$\ddot{\theta}+\omega^{2} \sin \theta=0 . \nonumber$$

We next employ a technique that is useful for equations of the form

$$\ddot{\theta}+F(\theta)=0 \nonumber$$

when it is easy to integrate the function $F(\theta)$. Namely, we note that

$$\dfrac{d}{d t}\left[\dfrac{1}{2} \dot{\theta}^{2}+\int^{\theta(t)} F(\phi) d \phi\right]=[\ddot{\theta}+F(\theta)] \dot{\theta} \nonumber$$

For our problem, we multiply Equation (3.17) by $\dot{\theta}$,

$$\ddot{\theta} \dot{\theta}+\omega^{2} \sin \theta \dot{\theta}=0 \nonumber$$

and note that the left side of this equation is a perfect derivative. Thus,

$$\dfrac{d}{d t}\left[\dfrac{1}{2} \dot{\theta}^{2}-\omega^{2} \cos \theta\right]=0 \nonumber$$

Therefore, the quantity in the brackets is a constant. So, we can write

$$\dfrac{1}{2} \dot{\theta}^{2}-\omega^{2} \cos \theta=c . \nonumber$$

Solving for $\dot{\theta}$, we obtain

$$\dfrac{d \theta}{d t}=\sqrt{2\left(c+\omega^{2} \cos \theta\right)} \nonumber$$

This equation is a separable first order equation and we can rearrange and integrate the terms to find that

$$t=\int d t=\int \dfrac{d \theta}{\sqrt{2\left(c+\omega^{2} \cos \theta\right)}} . \nonumber$$

Of course, one needs to be able to do the integral. When one gets a solution in this implicit form, one says that the problem has been solved by quadratures. Namely, the solution is given in terms of some integral. In the appendix to this chapter we show that this solution can be written in terms of elliptic integrals and derive corrections to formula for the period of a pendulum.

The Stability of Fixed Points in Nonlinear Systems

We are now interested in studying the stability of the equilibrium solutions of the nonlinear pendulum. Along the way we will develop some basic methods for studying the stability of equilibria in nonlinear systems.

We begin with the linear differential equation for damped oscillations as given earlier in Equation (3.9). In this case, we have a second order equation of the form

$$x^{\prime \prime}+b x^{\prime}+\omega^{2} x . \nonumber$$

Using the methods of Chapter 2, this second order equation can be written as a system of two first order equations:

$$\begin{aligned} &x^{\prime}=y \\[4pt] &y^{\prime}=-b y-\omega^{2} x . \end{aligned} \nonumber$$

This system has only one equilibrium solution, $x=0, y=0$.

Turning to the damped nonlinear pendulum, we have the system

$$\begin{aligned} x^{\prime} &=y \\[4pt] y^{\prime} &=-b y-\omega^{2} \sin x . \end{aligned} \nonumber$$

This system also has the equilibrium solution, $x=0, y=0$. However, there are actually an infinite number of solutions. The equilibria are determined from $y=0$ and $-b y-\omega^{2} \sin x=0$. This implies that $\sin x=0$. There are an infinite number of solutions: $x=n \pi, n=0, \pm 1, \pm 2, \ldots$ So, we have an infinite number of equilibria, $(n \pi, 0), n=0, \pm 1, \pm 2, \ldots$

Next, we need to determine their stability. To do this we need a more general theory for nonlinear systems. We begin with the $n$-dimensional system

$$\mathbf{x}^{\prime}=\mathbf{f}(\mathbf{x}), \quad \mathrm{x} \in \mathrm{R}^{n} \nonumber$$

Here $\mathbf{f}: \mathrm{R}^{n} \rightarrow \mathrm{R}^{n}$. We define fixed points, or equilibrium solutions, of this system as points $\mathrm{x}^{*}$ satisfying $\mathbf{f}\left(\mathrm{x}^{*}\right)=\mathbf{0}$.

The stability in the neighborhood of fixed points can now be determined. We are interested in what happens to solutions of our system with initial conditions starting near a fixed point. We can represent a point near a fixed point in the form $\mathbf{x}=\mathbf{x}^{*}+\boldsymbol{\xi}$, where the length of $\boldsymbol{\xi}$ gives an indication of how close we are to the fixed point. So, we consider that initially, $|\boldsymbol{\xi}| \ll 1$.

As the system evolves, $\boldsymbol{\xi}$ will change. The change of $\boldsymbol{\xi}$ in time is in turn governed by a system of equations. We can approximate this evolution as follows. First, we note that

$$\mathbf{x}^{\prime}=\boldsymbol{\xi}^{\prime} \nonumber$$

Next, we have that

$$\mathbf{f}(\mathbf{x})=\mathbf{f}\left(\mathbf{x}^{*}+\boldsymbol{\xi}\right) \nonumber$$

We can expand the right side about the fixed point using a multidimensional version of Taylor’s Theorem. Thus, we have that

$$\mathbf{f}\left(\mathbf{x}^{*}+\boldsymbol{\xi}\right)=\mathbf{f}\left(\mathbf{x}^{*}\right)+D \mathbf{f}\left(\mathbf{x}^{*}\right) \boldsymbol{\xi}+O\left(|\boldsymbol{\xi}|^{2}\right) \nonumber$$

Here Df is the Jacobian matrix, defined as

$$D \mathbf{f}=\left(\begin{array}{cccc} \dfrac{\partial f_{1}}{\partial x_{1}} & \dfrac{\partial f_{1}}{\partial x_{2}} & \cdots & \dfrac{\partial f_{1}}{\partial x_{n}} \\[4pt] \dfrac{\partial f_{2}}{\partial x_{1}} & \ddots & \ddots & \vdots \\[4pt] \vdots & \ddots & \ddots & \vdots \\[4pt] \dfrac{\partial f_{n}}{\partial x_{1}} & \cdots & \cdots & \dfrac{\partial f_{n}}{\partial x_{n}} \end{array}\right) \nonumber$$

Noting that $\mathbf{f}\left(\mathbf{x}^{*}\right)=\mathbf{0}$, we then have that system (3.22) becomes

$$\xi^{\prime} \approx D \mathbf{f}\left(\mathbf{x}^{*}\right) \boldsymbol{\xi} \nonumber$$

It is this equation which describes the behavior of the system near the fixed point. We say that system (3.22) has been linearized or that Equation (3.23) is the linearization of system $(3.22)$. Example 3.4. As an example of the application of this linearization, we look at the system

$$\begin{aligned} &x^{\prime}=-2 x-3 x y \\[4pt] &y^{\prime}=3 y-y^{2} \end{aligned} \nonumber$$

We first determine the fixed points:

$$\begin{aligned} &0=-2 x-3 x y=-x(2+3 y) \\[4pt] &0=3 y-y^{2}=y(3-y) \end{aligned} \nonumber$$

From the second equation, we have that either $y=0$ or $y=3$. The first equation then gives $x=0$ in either case. So, there are two fixed points: $(0,0)$ and $(0,3)$

Next, we linearize about each fixed point separately. First, we write down the Jacobian matrix.

$$D \mathbf{f}(x, y)=\left(\begin{array}{cc} -2-3 y & -3 x \\[4pt] 0 & 3-2 y \end{array}\right) \nonumber$$

1. Case I $(0,0)$.

In this case we find that

$$D \mathbf{f}(0,0)=\left(\begin{array}{cc} -2 & 0 \\[4pt] 0 & 3 \end{array}\right) \nonumber$$

Therefore, the linearized equation becomes

$$\xi^{\prime}=\left(\begin{array}{cc} -2 & 0 \\[4pt] 0 & 3 \end{array}\right) \boldsymbol{\xi} \nonumber$$

This is equivalently written out as the system

$$\begin{aligned} &\xi_{1}^{\prime}=-2 \xi_{1} \\[4pt] &\xi_{2}^{\prime}=3 \xi_{2} \end{aligned} \nonumber$$

This is the linearized system about the origin. Note the similarity with the original system. We emphasize that the linearized equations are constant coefficient equations and we can use earlier matrix methods to determine the nature of the equilibrium point. The eigenvalues of the system are obviously $\lambda=-2,3$. Therefore, we have that the origin is a saddle point.

3. Case II $(0,3)$.

In this case we proceed as before. We write down the Jacobian matrix and look at its eigenvalues to determine the type of fixed point. So, we have that the Jacobian matrix is

$$D \mathbf{f}(0,3)=\left(\begin{array}{cc} -2 & 0 \\[4pt] 0 & -3 \end{array}\right) \nonumber$$

Here, we have the eigenvalues $\lambda=-2,-3$. So, this fixed point is a stable node. This analysis has given us a saddle and a stable node. We know what the behavior is like near each fixed point, but we have to resort to other means to say anything about the behavior far from these points. The phase portrait for this system is given in Figure 3.14. You should be able to find the saddle point and the node. Notice how solutions behave in regions far from these points.

We can expect to be able to perform a linearization under general conditions. These are given in the Hartman-Großman Theorem:

Theorem 3.5. A continuous map exists between the linear and nonlinear systems when $D \mathbf{f}\left(\mathbf{x}^{*}\right)$ does not have any eigenvalues with zero real part.

Generally, there are several types of behavior that one can see in nonlinear systems. One can see sinks or sources, hyperbolic (saddle) points, elliptic points (centers) or foci. We have defined some of these for planar systems. In general, if at least two eigenvalues have real parts with opposite signs, then the fixed point is a hyperbolic point. If the real part of a nonzero eigenvalue is zero, then we have a center, or elliptic point.

Nonlinear Population Models

We have already encountered several models of population dynamics. Of course, one could dream up several other examples. There are two standard types of models: Predator-prey and competing species. In the predator-prey model, one typically has one species, the predator, feeding on the other, the prey. We will look at the standard Lotka-Volterra model in this section. The competing species model looks similar, except there are a few sign changes, since one species is not feeding on the other. Also, we can build in logistic terms into our model. We will save this latter type of model for the homework.

The Lotka-Volterra model takes the form

$$\begin{aligned} &\dot{x}=a x-b x y, \\[4pt] &\dot{y}=-d y+c x y \end{aligned} \nonumber$$

In this case, we can think of $x$ as the population of rabbits (prey) and $y$ is the population of foxes (predators). Choosing all constants to be positive, we can describe the terms.

• ax: When left alone, the rabbit population will grow. Thus $a$ is the natural growth rate without predators.
• $-d y$ : When there are no rabbits, the fox population should decay. Thus, the coefficient needs to be negative.

• $-b x y$ : We add a nonlinear term corresponding to the depletion of the rabbits when the foxes are around.
• cxy: The more rabbits there are, the more food for the foxes. So, we add a nonlinear term giving rise to an increase in fox population.

The analysis of the Lotka-Volterra model begins with determining the fixed points. So, we have from Equation (3.34)

$$\begin{gathered} x(a-b y)=0, \\[4pt] y(-d+c x)=0 . \end{gathered} \nonumber$$

Therefore, the origin and $\left(\dfrac{d}{c} \dfrac{a}{b}\right)$ are the fixed points.

Next, we determine their stability, by linearization about the fixed points. We can use the Jacobian matrix, or we could just expand the right hand side of each equation in (3.34). The Jacobian matrix is $D f(x, y)=\left(\begin{array}{cc}a-b y & -b x \\[4pt] c y & -d+c x\end{array}\right)$. Evaluating at each fixed point, we have

$$\begin{gathered} D f(0,0)=\left(\begin{array}{cc} a & 0 \\[4pt] 0 & -d \end{array}\right), \\[4pt] D f\left(\dfrac{d}{c}, \dfrac{a}{b}\right)=\left(\begin{array}{cc} 0 & -\dfrac{b d}{c} \\[4pt] \dfrac{a c}{b} & 0 \end{array}\right) . \end{gathered} \nonumber$$

The eigenvalues of $(3.36)$ are $\lambda=a,-d$. So, the origin is a saddle point. The eigenvalues of (3.37) satisfy $\lambda^{2}+a d=0$. So, the other point is a center. In Figure $3.17$ we show a sample direction field for the Lotka-Volterra system.

Another way to linearize is to expand the equations about the fixed points. Even though this is equivalent to computing the Jacobian matrix, it sometimes might be faster.

Limit Cycles

So far we have just been concerned with equilibrium solutions and their behavior. However, asymptotically stable fixed points are not the only attractors. There are other types of solutions, known as limit cycles, towards which a solution may tend. In this section we will look at some examples of these periodic solutions.

Such solutions are common in nature. Rayleigh investigated the problem

$$x^{\prime \prime}+c\left(\dfrac{1}{3}\left(x^{\prime}\right)^{2}-1\right) x^{\prime}+x=0 \nonumber$$

in the study of the vibrations of a violin string. Van der Pol studied an electrical circuit, modelling this behavior. Others have looked into biological systems, such as neural systems, chemical reactions, such as Michaelis-Menton kinetics or systems leading to chemical oscillations. One of the most important models in the historical study of dynamical systems is that of planetary motion and investigating the stability of planetary orbits. As is well known, these orbits are periodic.

Limit cycles are isolated periodic solutions towards which neighboring states might tend when stable. A key example exhibiting a limit cycle is given by the system

$$\begin{aligned} &x^{\prime}=\mu x-y-x\left(x^{2}+y^{2}\right) \\[4pt] &y^{\prime}=x+\mu y-y\left(x^{2}+y^{2}\right) \end{aligned} \nonumber$$

It is clear that the origin is a fixed point. The Jacobian matrix is given as

$$\operatorname{Df}(0,0)=\left(\begin{array}{cc} \mu & -1 \\[4pt] 1 & \mu \end{array}\right) \nonumber$$

The eigenvalues are found to be $\lambda=\mu \pm i$. For $\mu=0$ we have a center. For $\mu<0$ we have a stable spiral and for $\mu>0$ we have an unstable spiral. However, this spiral does not wander off to infinity. We see in Figure $3.18$ that equilibrium point is a spiral. However, in Figure $3.19$ it is clear that the solution does not spiral out to infinity. It is bounded by a circle.

One can actually find the radius of this circle. This requires rewriting the system in polar form. Recall from Chapter 2 that this is done using

$$\begin{gathered} r r^{\prime}=x x^{\prime}+y y^{\prime}, \\[4pt] r^{2} \theta^{\prime}=x y^{\prime}-y x^{\prime} . \end{gathered} \nonumber$$

Inserting the system (3.39) into these expressions, we have

$$r r^{\prime}=\mu r^{2}-r^{4}, \quad r^{2} \theta^{\prime}=r^{2}, \nonumber$$

$$r^{\prime}=\mu r-r^{3}, \theta^{\prime}=1 . \nonumber$$

Of course, for a circle $r=$ const, therefore we need to look at the equilibrium solutions of Equation (3.43). This amounts to solving $\mu r-r^{3}=0$ for $r$. The solutions of this equation are $r=0, \pm \sqrt{\mu}$. We need only keep the one positive radius solution, $r=\sqrt{\mu}$. In Figures $3.18-3.19 \mu=0.4$, so we expect a circle with $r=\sqrt{0.4} \approx 0.63$. The $\theta$ equation just tells us that we follow the limit cycle in a counterclockwise direction.

Limit cycles are not always circles. In Figures 3.20-3.21 we show the behavior of the Rayleigh system (3.38) for $c=0.4$ and $c=2.0$. In this case we see that solutions tend towards a noncircular limit cycle.

The limit cycle for $c=2.0$ is shown in Figure $3.22$.

Can one determine ahead of time if a given nonlinear system will have a limit cycle? In order to answer this question, we will introduce some definitions.

We first describe different trajectories and families of trajectories. A flow on $R^{2}$ is a function $\phi$ that satisfies the following

1. $\phi(\mathbf{x}, t)$ is continuous in both arguments.
2. $\phi(\mathbf{x}, 0)=\mathbf{x}$ for all $\mathbf{x} \in R^{2}$
3. $\phi\left(\phi\left(\mathbf{x}, t_{1}\right), t_{2}\right)=\phi\left(\mathbf{x}, t_{1}+t_{2}\right)$.

The orbit, or trajectory, through $\mathbf{x}$ is defined as $\gamma=\{\phi(\mathbf{x}, t) \mid t \in I\}$. In Figure $3.23$ we demonstrate these properties. For $t=0, \phi(\mathbf{x}, 0)=\mathbf{x}$. Increasing $t$, one follows the trajectory until one reaches the point $\phi\left(\mathbf{x}, t_{1}\right)$. Continuing $t_{2}$ further, one is then at $\phi\left(\phi\left(\mathbf{x}, t_{1}\right), t_{2}\right)$. By the third property, this is the same as going from $\mathbf{x}$ to $\phi\left(\mathbf{x}, t_{1}+t_{2}\right)$ for $t=t_{1}+t_{2}$.

Having defined the orbits, we need to define the asymptotic behavior of the orbit for both positive and negative large times. We define the positive semiorbit through $\mathbf{x}$ as $\gamma^{+}=\{\phi(\mathbf{x}, t) \mid t>0\}$. The negative semiorbit through $\mathbf{x}$ is defined as $\gamma^{-}=\{\phi(\mathbf{x}, t) \mid t<0\}$. Thus, we have $\gamma=\gamma^{+} \cup^{-}$.

The positive limit set, or $\omega$-limit set, of point $\mathbf{x}$ is defined as

$$\Lambda^{+}=\left\{\mathbf{y} \mid \text { there exists a sequence of } t_{n} \rightarrow \infty \text { such that } \phi\left(\mathbf{x}, t_{n}\right) \rightarrow \mathbf{y}\right\} \nonumber$$

The $\mathbf{y}$ ’s are referred to as $\omega$-limit points. This is shown in Figure $3.24$.

Similarly, we define the negative limit set, or it alpha-limit sets, of point $\mathbf{x}$ is defined as $\Lambda^{-}=\left\{\mathbf{y} \mid\right.$ there exists a sequences of $t_{n} \rightarrow-\infty$ such that $\left.\phi\left(\mathbf{x}, t_{n}\right) \rightarrow \mathbf{y}\right\}$

and the corresponding $\mathbf{y}$ ’s are $\alpha$-limit points. This is shown in Figure $3.25$.

There are several types of orbits that a system might possess. A cycle or periodic orbit is any closed orbit which is not an equilibrium point. A periodic orbit is stable if for every neighborhood of the orbit such that all nearby orbits stay inside the neighborhood. Otherwise, it is unstable. The orbit is asymptotically stable if all nearby orbits converge to the periodic orbit.

A limit cycle is a cycle which is the $\alpha$ or $\omega$-limit set of some trajectory other than the limit cycle. A limit cycle $\Gamma$ is stable if $\Lambda^{+}=\Gamma$ for all $\mathbf{x}$ in some neighborhood of $\Gamma$. A limit cycle $\Gamma$ is unstable if $\Lambda^{-}=\Gamma$ for all $\mathbf{x}$ in some neighborhood of $\Gamma$. Finally, a limits cycle is semistable if it is attracting on one side and repelling on the other side. In the previous examples, we saw limit cycles that were stable. Figures $3.24$ and $3.25$ depict stable and unstable limit cycles, respectively.

We now state a theorem which describes the type of orbits we might find in our system.

Theorem 3.7. Poincaré-Bendixon Theorem Let $\gamma^{+}$be contained in a bounded region in which there are finitely many critical points. Then $\Lambda^{+}$is either

1. a single critical point;
2. a single closed orbit;
3. a set of critical points joined by heteroclinic orbits. [Compare Figures $3.27$ and ??.]

We are interested in determining when limit cycles may, or may not, exist. A consequence of the Poincaré-Bendixon Theorem is given by the following corollary.

Corollary Let $D$ be a bounded closed set containing no critical points and suppose that $\gamma^{+} \subset D$. Then there exists a limit cycle contained in $D$.

More specific criteria allow us to determine if there is a limit cycle in a given region. These are given by Dulac’s Criteria and Bendixon’s Criteria.

Dulac’s Criteria Consider the autonomous planar system

$$x^{\prime}=f(x, y), \quad y^{\prime}=g(x, y) \nonumber$$

and a continuously differentiable function $\psi$ defined on an annular region $D$ contained in some open set. If

$$\dfrac{\partial}{\partial x}(\psi f)+\dfrac{\partial}{\partial y}(\psi g) \nonumber$$

does not change sign in $D$, then there is at most one limit cycle contained entirely in $D$.

Bendixon’s Criteria Consider the autonomous planar system

$$x^{\prime}=f(x, y), \quad y^{\prime}=g(x, y) \nonumber$$

defined on a simply connected domain $D$ such that

$$\dfrac{\partial}{\partial x}(\psi f)+\dfrac{\partial}{\partial y}(\psi g) \neq 0 \nonumber$$

in $D$. Then there are no limit cycles entirely in $D$.

These are easily proved using Green’s Theorem in the plane. We prove Bendixon’s Criteria. Let $\mathbf{f}=(f, g)$. Assume that $\Gamma$ is a closed orbit lying in $D$. Let $S$ be the interior of $\Gamma$. Then

$$\begin{aligned} \int_{S} \nabla \cdot \mathbf{f} d x d y &=\oint_{\Gamma}(f d y-g d x) \\[4pt] &=\int_{0}^{T}(f \dot{y}-g \dot{x}) d t \\[4pt] &=\int_{0}^{T}(f g-g f) d t=0 \end{aligned} \nonumber$$

So, if $\nabla \cdot \mathbf{f}$ is not identically zero and does not change sign in $S$, then from the continuity of $\nabla \cdot \mathbf{f}$ in $S$ we have that the right side above is either positive or negative. Thus, we have a contradiction and there is no closed orbit lying in $D$

Example 3.8. Consider the earlier example in (3.39) with $\mu=1$.

$$\begin{aligned} &x^{\prime}=x-y-x\left(x^{2}+y^{2}\right) \\[4pt] &y^{\prime}=x+y-y\left(x^{2}+y^{2}\right) . \end{aligned} \nonumber$$

We already know that a limit cycle exists at $x^{2}+y^{2}=1$. A simple computation gives that

$$\nabla \cdot \mathbf{f}=2-4 x^{2}-4 y^{2} \nonumber$$

For an arbitrary annulus $a<x^{2}+y^{2}<b$, we have

$$2-4 b<\nabla \cdot \mathbf{f}<2-4 a . \nonumber$$

For $a=3 / 4$ and $b=5 / 4,-3<\nabla \cdot \mathbf{f}<-1$. Thus, $\nabla \cdot \mathbf{f}<0$ in the annulus $3 / 4<x^{2}+y^{2}<5 / 4$. Therefore, by Dulac’s Criteria there is at most one limit cycle in this annulus.

Example 3.9. Consider the system

$$\begin{aligned} &x^{\prime}=y \\[4pt] &y^{\prime}=-a x-b y+c x^{2}+d y^{2} . \end{aligned} \nonumber$$

Let $\psi(x, y)=e^{-2 d x}$. Then,

$$\dfrac{\partial}{\partial x}(\psi y)+\dfrac{\partial}{\partial y}\left(\psi\left(-a x-b y+c x^{2}+d y^{2}\right)\right)=-b e^{-2 d x} \neq 0 \nonumber$$

We conclude by Bendixon’s Criteria that there are no limit cycles for this system.

Nonautonomous Nonlinear Systems

In this section we discuss nonautonomous systems. Recall that an autonomous system is one in which there is no explicit time dependence. A simple example is the forced nonlinear pendulum given by the nonhomogeneous equation

$$\ddot{x}+\omega^{2} \sin x=f(t) . \nonumber$$

We can set this up as a system of two first order equations:

$$\begin{aligned} &\dot{x}=y \\[4pt] &\dot{y}=-\omega^{2} \sin x+f(t) \end{aligned} \nonumber$$

This system is not in a form for which we could use the earlier methods. Namely, it is a nonautonomous system. However, we introduce a new variable $z(t)=t$ and turn it into an autonomous system in one more dimension. The new system takes the form

$$\begin{aligned} &\dot{x}=y \\[4pt] &\dot{y}=-\omega^{2} \sin x+f(z) \\[4pt] &\dot{z}=1 \end{aligned} \nonumber$$

This system is a three dimensional autonomous, possibly nonlinear, system and can be explored using our earlier methods.

A more interesting model is provided by the Duffing Equation. This equation models hard spring and soft spring oscillations. It also models a periodically forced beam as shown in Figure 3.28. It is of interest because it is a simple system which exhibits chaotic dynamics and will motivate us towards using new visualization methods for nonautonomous systems.

The most general form of Duffing’s equation is given by

$$\ddot{x}+k \dot{x}+\left(\beta x^{3} \pm \omega_{0}^{2} x\right)=\Gamma \cos (\omega t+\phi) . \nonumber$$

This equation models hard spring $(\beta>0)$ and soft spring $(\beta<0)$ oscillations. However, we will use a simpler version of the Duffing equation:

$$\ddot{x}+k \dot{x}+x^{3}-x=\Gamma \cos \omega t . \nonumber$$

Let’s first look at the behavior of some of the orbits of the system as we vary the parameters. In Figures $3.29-3.31$ we show some typical solution plots superimposed on the direction field. We start with the the undamped $(k=0)$ and unforced $(\Gamma=0)$ Duffing equation,

$$\ddot{x}+x^{3}-x==0 . \nonumber$$

We can write this second order equation as the autonomous system

$$\begin{aligned} &\dot{x}=y \\[4pt] &\dot{y}=x\left(1-x^{2}\right) \end{aligned} \nonumber$$

We see there are three equilibrium points at $(0,0),(\pm 1,0)$. In Figure $3.29$ we plot several orbits for We see that the three equilibrium points consist of two centers and a saddle.

We now turn on the damping. The system becomes

$$\begin{aligned} &\dot{x}=y \\[4pt] &\dot{y}=-k y+x\left(1-x^{2}\right) . \end{aligned} \nonumber$$

In Figure $3.30$ we show what happens when $k=0.1$ These plots are reminiscent of the plots for the nonlinear pendulum; however, there are fewer equilibria. The centers become stable spirals.

Next we turn on the forcing to obtain a damped, forced Duffing equation. The system is now nonautonomous.

$$\begin{aligned} &\dot{x}=y \\[4pt] &\dot{y}=x\left(1-x^{2}\right)+\Gamma \cos \omega t \end{aligned} \nonumber$$

In Figure $3.31$ we only show one orbit with $k=0.1, \Gamma=0.5$, and $\omega=1.25$. The solution intersects itself and look a bit messy. We can imagine what we would get if we added any more orbits. For completeness, we show in Figure $3.32$ an example with four different orbits.

In cases for which one has periodic orbits such as the Duffing equation, Poincaré introduced the notion of surfaces of section. One embeds the orbit in a higher dimensional space so that there are no self intersections, like we saw in Figures $3.31$ and 3.32. In Figure $3.33$ we show an example where a simple orbit is shown as it periodically pierces a given surface.

In order to simplify the resulting pictures, one only plots the points at which the orbit pierces the surface as sketched in Figure 3.34. In practice, there is a natural frequency, such as $\omega$ in the forced Duffing equation. Then, one plots points at times that are multiples of the period, $T=\dfrac{2 \pi}{\omega}$. In Figure $3.35$ we show what the plot for one orbit would look like for the damped, unforced Duffing equation.

The more interesting case, is when there is forcing and damping. In this case the surface of section plot is given in Figure $3.36$. While this is not as busy as the solution plot in Figure 3.31, it still provides some interesting behavior. What one finds is what is called a strange attractor. Plotting many orbits, we find that after a long time, all of the orbits are attracted to a small region in the plane, much like a stable node attracts nearby orbits. However, this

set consists of more than one point. Also, the flow on the attractor is chaotic in nature. Thus, points wander in an irregular way throughout the attractor. This is one of the interesting topics in chaos theory and this whole theory of dynamical systems has only been touched in this text leaving the reader to wander of into further depth into this fascinating field.

Maple Code for Phase Plane Plots

For reference, the plots in Figures $3.29$ and $3.30$ were generated in Maple using the following commands:

This leads to what is known as a strange attractor.

The surface of section plots at the end of the last section were obtained using code from S. Lynch’s book Dynamical Systems with Applications Using Maple. The Maple code is given by

Appendix: Period of the Nonlinear Pendulum

In Section 3.5.1 we saw that the solution of the nonlinear pendulum problem can be found up to quadrature. In fact, the integral in Equation (3.19) can be transformed into what is know as an elliptic integral of the first kind. We will rewrite our result and then use it to obtain an approximation to the period of oscillation of our nonlinear pendulum, leading to corrections to the linear result found earlier.

We will first rewrite the constant found in (3.18). This requires a little physics. The swinging of a mass on a string, assuming no energy loss at the pivot point, is a conservative process. Namely, the total mechanical energy is conserved. Thus, the total of the kinetic and gravitational potential energies is a constant. Noting that $v=L \dot{\theta}$, the kinetic energy of the mass on the string is given as

$$T=\dfrac{1}{2} m v^{2}=\dfrac{1}{2} m L^{2} \dot{\theta}^{2} . \nonumber$$

The potential energy is the gravitational potential energy. If we set the potential energy to zero at the bottom of the swing, then the potential energy is $U=m g h$, where $h$ is the height that the mass is from the bottom of the swing. A little trigonometry gives that $h=L(1-\cos \theta)$. This gives the potential energy as

$$U=m g L(1-\cos \theta) . \nonumber$$

So, the total mechanical energy is

$$E=\dfrac{1}{2} m L^{2} \theta^{\prime 2}+m g L(1-\cos \theta) . \nonumber$$

We note that a little rearranging shows that we can relate this to Equation $(3.18)$

$$\dfrac{1}{2}\left(\theta^{\prime}\right)^{2}-\omega^{2} \cos \theta=\dfrac{1}{m L^{2}} E-\omega^{2}=c . \nonumber$$

We can use Equation (3.56) to get a value for the total energy. At the top of the swing the mass is not moving, if only for a moment. Thus, the kinetic energy is zero and the total energy is pure potential energy. Letting $\theta_{0}$ denote the angle at the highest position, we have that

$$E=m g L\left(1-\cos \theta_{0}\right)=m L^{2} \omega^{2}\left(1-\cos \theta_{0}\right) . \nonumber$$

Here we have used the relation $g=L \omega^{2}$.

Therefore, we have found that

$$\dfrac{1}{2} \dot{\theta}^{2}-\omega^{2} \cos \theta=\omega^{2}\left(1-\cos \theta_{0}\right) . \nonumber$$

Using the half angle formula,

$$\sin ^{2} \dfrac{\theta}{2}=\dfrac{1}{2}(1-\cos \theta) \nonumber$$

we can rewrite Equation $(3.57)$ as

$$\dfrac{1}{2} \dot{\theta}^{2}=2 \omega^{2}\left[\sin ^{2} \dfrac{\theta_{0}}{2}-\sin ^{2} \dfrac{\theta}{2}\right] \nonumber$$

Solving for $\theta^{\prime}$, we have

$$\dfrac{d \theta}{d t}=2 \omega\left[\sin ^{2} \dfrac{\theta_{0}}{2}-\sin ^{2} \dfrac{\theta}{2}\right]^{1 / 2} \nonumber$$

One can now apply separation of variables and obtain an integral similar to the solution we had obtained previously. Noting that a motion from $\theta=0$ to $\theta=\theta_{0}$ is a quarter of a cycle, then we have that

$$T=\dfrac{2}{\omega} \int_{0}^{\theta_{0}} \dfrac{d \phi}{\sqrt{\sin ^{2} \dfrac{\theta_{0}}{2}-\sin ^{2} \dfrac{\theta}{2}}} \nonumber$$

This result is not much different than our previous result, but we can now easily transform the integral into an elliptic integral. We define

$$z=\dfrac{\sin \dfrac{\theta}{2}}{\sin \dfrac{\theta_{0}}{2}} \nonumber$$

and

$$k=\sin \dfrac{\theta_{0}}{2} \nonumber$$

Then Equation (3.60) becomes

$$T=\dfrac{4}{\omega} \int_{0}^{1} \dfrac{d z}{\sqrt{\left(1-z^{2}\right)\left(1-k^{2} z^{2}\right)}} . \nonumber$$

This is done by noting that $d z=\dfrac{1}{2 k} \cos \dfrac{\theta}{2} d \theta=\dfrac{1}{2 k}\left(1-k^{2} z^{2}\right)^{1 / 2} d \theta$ and that $\sin ^{2} \dfrac{\theta_{0}}{2}-\sin ^{2} \dfrac{\theta}{2}=k^{2}\left(1-z^{2}\right)$. The integral in this result is an elliptic integral of the first kind. In particular, the elliptic integral of the first kind is defined

$$F(\phi, k) \equiv=\int_{0}^{\phi} \dfrac{d \theta}{\sqrt{1-k^{2} \sin ^{2} \theta}}=\int_{0}^{\sin \phi} \dfrac{d z}{\sqrt{\left(1-z^{2}\right)\left(1-k^{2} z^{2}\right)}} . \nonumber$$

In some contexts, this is known as the incomplete elliptic integral of the first kind and $K(k)=F\left(\dfrac{\pi}{2}, k\right)$ is called the complete integral of the first kind.

There are tables of values for elliptic integrals. Historically, that is how one found values of elliptic integrals. However, we now have access to computer algebra systems which can be used to compute values of such integrals. For small angles, we have that $k$ is small. So, we can develop a series expansion for the period, $T$, for small $k$. This is done by first expanding

$$\left(1-k^{2} z^{2}\right)^{-1 / 2}=1+\dfrac{1}{2} k^{2} z^{2}+\dfrac{3}{8} k^{2} z^{4}+O\left((k z)^{6}\right) \nonumber$$

Substituting this in the integrand and integrating term by term, one finds that

$$T=2 \pi \sqrt{\dfrac{L}{g}}\left[1+\dfrac{1}{4} k^{2}+\dfrac{9}{64} k^{4}+\ldots\right] \nonumber$$

This expression gives further corrections to the linear result, which only provides the first term. In Figure $3.37$ we show the relative errors incurred when keeping the $k^{2}$ and $k^{4}$ terms versus not keeping them. The reader is asked to explore this further in Problem 3.8.

Introduction

Until this point we have solved initial value problems. For an initial value problem one has to solve a differential equation subject to conditions on the unknown function and its derivatives at one value of the independent variable. For example, for $x=x(t)$ we could have the initial value problem

$$x^{\prime \prime}+x=2, \quad x(0)=1, \quad x^{\prime}(0)=0 \nonumber$$

In the next chapters we will study boundary value problems and various tools for solving such problems. In this chapter we will motivate our interest in boundary value problems by looking into solving the one-dimensional heat equation, which is a partial differential equation. for the rest of the section, we will use this solution to show that in the background of our solution of boundary value problems is a structure based upon linear algebra and analysis leading to the study of inner product spaces. Though technically, we should be lead to Hilbert spaces, which are complete inner product spaces.

For an initial value problem one has to solve a differential equation subject to conditions on the unknown function or its derivatives at more than one value of the independent variable. As an example, we have a slight modification of the above problem: Find the solution $x=x(t)$ for $0 \leq t \leq 1$ that satisfies the problem

$$x^{\prime \prime}+x=2, \quad x(0)=1, \quad x(1)=0 . \nonumber$$

Typically, initial value problems involve time dependent functions and boundary value problems are spatial. So, with an initial value problem one knows how a system evolves in terms of the differential equation and the state of the system at some fixed time. Then one seeks to determine the state of the system at a later time.

For boundary values problems, one knows how each point responds to its neighbors, but there are conditions that have to be satisfied at the endpoints. An example would be a horizontal beam supported at the ends, like a bridge. The shape of the beam under the influence of gravity, or other forces, would lead to a differential equation and the boundary conditions at the beam ends would affect the solution of the problem. There are also a variety of other types of boundary conditions. In the case of a beam, one end could be fixed and the other end could be free to move. We will explore the effects of different boundary value conditions in our discussions and exercises.

Let’s solve the above boundary value problem. As with initial value problems, we need to find the general solution and then apply any conditions that we may have. This is a nonhomogeneous differential equation, so we have that the solution is a sum of a solution of the homogeneous equation and a particular solution of the nonhomogeneous equation, $x(t)=x_{h}(t)+x_{p}(t)$. The solution of $x^{\prime \prime}+x=0$ is easily found as

$$x_{h}(t)=c_{1} \cos t+c_{2} \sin t \nonumber$$

The particular solution is easily found using the Method of Undetermined Coefficients,

$$x_{p}(t)=2 \nonumber$$

Thus, the general solution is

$$x(t)=2+c_{1} \cos t+c_{2} \sin t . \nonumber$$

We now apply the boundary conditions and see if there are values of $c_{1}$ and $c_{2}$ that yield a solution to our problem. The first condition, $x(0)=0$, gives

$$0=2+c_{1} \nonumber$$

Thus, $c_{1}=-2$. Using this value for $c_{1}$, the second condition, $x(1)=1$, gives

$$0=2-2 \cos 1+c_{2} \sin 1 \nonumber$$

This yields

$$c_{2}=\dfrac{2(\cos 1-1)}{\sin 1} . \nonumber$$

We have found that there is a solution to the boundary value problem and it is given by

$$x(t)=2\left(1-\cos t \dfrac{(\cos 1-1)}{\sin 1} \sin t\right) \nonumber$$

Boundary value problems arise in many physical systems, just as many of the initial values problems we have seen. We will see in the next section that boundary value problems for ordinary differential equations often appear in the solution of partial differential equations.

Partial Differential Equations

In this section we will introduce some generic partial differential equations and see how the discussion of such equations leads naturally to the study of boundary value problems for ordinary differential equations. However, we will not derive the particular equations, leaving that to courses in differential equations, mathematical physics, etc.

For ordinary differential equations, the unknown functions are functions of a single variable, e.g., $y=y(x)$. Partial differential equations are equations involving an unknown function of several variables, such as $u=u(x, y), u=$ $u(x, y), u=u(x, y, z, t)$, and its (partial) derivatives. Therefore, the derivatives are partial derivatives. We will use the standard notations $u_{x}=\dfrac{\partial u}{\partial x}, u_{x x}=\dfrac{\partial^{2} u}{\partial x^{2}}$, etc.

There are a few standard equations that one encounters. These can be studied in one to three dimensions and are all linear differential equations. A list is provided in Table 4.1. Here we have introduced the Laplacian operator, $\nabla^{2} u=u_{x x}+u_{y y}+u_{z z}$. Depending on the types of boundary conditions imposed and on the geometry of the system (rectangular, cylindrical, spherical, etc.), one encounters many interesting boundary value problems for ordinary differential equations.

Table 4.1. List of generic partial differential equations.

Let’s look at the heat equation in one dimension. This could describe the heat conduction in a thin insulated rod of length $L$. It could also describe the diffusion of pollutant in a long narrow stream, or the flow of traffic down a road. In problems involving diffusion processes, one instead calls this equation the diffusion equation.

A typical initial-boundary value problem for the heat equation would be that initially one has a temperature distribution $u(x, 0)=f(x)$. Placing the bar in an ice bath and assuming the heat flow is only through the ends of the bar, one has the boundary conditions $u(0, t)=0$ and $u(L, t)=0$. Of course, we are dealing with Celsius temperatures and we assume there is plenty of ice to keep that temperature fixed at each end for all time. So, the problem one would need to solve is given as

Another problem that will come up in later discussions is that of the vibrating string. A string of length $L$ is stretched out horizontally with both ends fixed. Think of a violin string or a guitar string. Then the string is plucked, giving the string an initial profile. Let $u(x, t)$ be the vertical displacement of the string at position $x$ and time $t$. The motion of the string is governed by the one dimensional wave equation. The initial-boundary value problem for this problem is given as

Connections to Linear Algebra

We have already seen in earlier chapters that ideas from linear algebra crop up in our studies of differential equations. Namely, we solved eigenvalue problems associated with our systems of differential equations in order to determine the local behavior of dynamical systems near fixed points. In our study of boundary value problems we will find more connections with the theory of vector spaces. However, we will find that our problems lie in the realm of infinite dimensional vector spaces. In this section we will begin to see these connections.

Introduction

In this chapter we will look at trigonometric series. Previously, we saw that such series expansion occurred naturally in the solution of the heat equation and other boundary value problems. In the last chapter we saw that such functions could be viewed as a basis in an infinite dimensional vector space of functions. Given a function in that space, when will it have a representation as a trigonometric series? For what values of $x$ will it converge? Finding such series is at the heart of Fourier, or spectral, analysis.

There are many applications using spectral analysis. At the root of these studies is the belief that many continuous waveforms are comprised of a number of harmonics. Such ideas stretch back to the Pythagorean study of the vibrations of strings, which lead to their view of a world of harmony. This idea was carried further by Johannes Kepler in his harmony of the spheres approach to planetary orbits. In the 1700 ’s others worked on the superposition theory for vibrating waves on a stretched spring, starting with the wave equation and leading to the superposition of right and left traveling waves. This work was carried out by people such as John Wallis, Brook Taylor and Jean le Rond d’Alembert.

In 1742 d’Alembert solved the wave equation

$$c^{2} \dfrac{\partial^{2} y}{\partial x^{2}}-\dfrac{\partial^{2} y}{\partial t^{2}}=0 \nonumber$$

where $y$ is the string height and $c$ is the wave speed. However, his solution led himself and others, like Leonhard Euler and Daniel Bernoulli, to investigate what "functions" could be the solutions of this equation. In fact, this lead to a more rigorous approach to the study of analysis by first coming to grips with the concept of a function. For example, in 1749 Euler sought the solution for a plucked string in which case the initial condition $y(x, 0)=h(x)$ has a discontinuous derivative! In 1753 Daniel Bernoulli viewed the solutions as a superposition of simple vibrations, or harmonics. Such superpositions amounted to looking at solutions of the form

$$y(x, t)=\sum_{k} a_{k} \sin \dfrac{k \pi x}{L} \cos \dfrac{k \pi c t}{L}, \nonumber$$

where the string extends over the interval $[0, L]$ with fixed ends at $x=0$ and $x=L$. However, the initial conditions for such superpositions are

$$y(x, 0)=\sum_{k} a_{k} \sin \dfrac{k \pi x}{L} . \nonumber$$

It was determined that many functions could not be represented by a finite number of harmonics, even for the simply plucked string given by an initial condition of the form

$$y(x, 0)=\left\{\begin{array}{cc} c x, & 0 \leq x \leq L / 2 \\[4pt] c(L-x), & L / 2 \leq x \leq L \end{array}\right. \nonumber$$

Thus, the solution consists generally of an infinite series of trigonometric functions.

Such series expansions were also of importance in Joseph Fourier’s solution of the heat equation. The use of such Fourier expansions became an important tool in the solution of linear partial differential equations, such as the wave equation and the heat equation. As seen in the last chapter, using the Method of Separation of Variables, allows higher dimensional problems to be reduced to several one dimensional boundary value problems. However, these studies lead to very important questions, which in turn opened the doors to whole fields of analysis. Some of the problems raised were

1. What functions can be represented as the sum of trigonometric functions?
2. How can a function with discontinuous derivatives be represented by a sum of smooth functions, such as the above sums?
3. Do such infinite sums of trigonometric functions a actually converge to the functions they represents?

Sums over sinusoidal functions naturally occur in music and in studying sound waves. A pure note can be represented as

$$y(t)=A \sin (2 \pi f t), \nonumber$$

where $A$ is the amplitude, $f$ is the frequency in hertz $(\mathrm{Hz})$, and $t$ is time in seconds. The amplitude is related to the volume, or intensity, of the sound. The larger the amplitude, the louder the sound. In Figure $5.1$ we show plots of two such tones with $f=2 \mathrm{~Hz}$ in the top plot and $f=5 \mathrm{~Hz}$ in the bottom one.

different amplitudes and frequencies. In Figure $5.2$ we see what happens when we add several sinusoids. Note that as one adds more and more tones with different characteristics, the resulting signal gets more complicated. However, we still have a function of time. In this chapter we will ask, "Given a function $f(t)$, can we find a set of sinusoidal functions whose sum converges to $f(t)$ ?"

Looking at the superpositions in Figure 5.2, we see that the sums yield functions that appear to be periodic. This is not to be unexpected. We recall that a periodic function is one in which the function values repeat over the domain of the function. The length of the smallest part of the domain which repeats is called the period. We can define this more precisely.

Definition 5.1. A function is said to be periodic with period $T$ if $f(t+T)=$ $f(t)$ for all $t$ and the smallest such positive number $T$ is called the period.

For example, we consider the functions used in Figure $5.2$. We began with $y(t)=2 \sin (4 \pi t)$. Recall from your first studies of trigonometric functions that one can determine the period by dividing the coefficient of $t$ into $2 \pi$ to get the period. In this case we have

$$T=\dfrac{2 \pi}{4 \pi}=\dfrac{1}{2} . \nonumber$$

Looking at the top plot in Figure $5.1$ we can verify this result. (You can count the full number of cycles in the graph and divide this into the total time to get a more accurate value of the period.)

Of course, this result makes sense, as the unit of frequency, the hertz, is also defined as $s^{-1}$, or cycles per second.

Returning to the superpositions in Figure 5.2, we have that $y(t)=$ $\sin (10 \pi t)$ has a period of $0.2 \mathrm{~Hz}$ and $y(t)=\sin (16 \pi t)$ has a period of $0.125 \mathrm{~Hz}$. The two superpositions retain the largest period of the signals added, which is $0.5 \mathrm{~Hz}$.

Our goal will be to start with a function and then determine the amplitudes of the simple sinusoids needed to sum to that function. First of all, we will see that this might involve an infinite number of such terms. Thus, we will be studying an infinite series of sinusoidal functions.

Secondly, we will find that using just sine functions will not be enough either. This is because we can add sinusoidal functions that do not necessarily peak at the same time. We will consider two signals that originate at different times. This is similar to when your music teacher would make sections of the class sing a song like "Row, Row, Row your Boat" starting at slightly different times.

We can easily add shifted sine functions. In Figure $5.3$ we show the functions $y(t)=2 \sin (4 \pi t)$ and $y(t)=2 \sin (4 \pi t+7 \pi / 8)$ and their sum. Note that this shifted sine function can be written as $y(t)=2 \sin (4 \pi(t+7 / 32))$. Thus, this corresponds to a time shift of $-7 / 8$.

We are now in a position to state our goal in this chapter.

Fourier Trigonometric Series

As we have seen in the last section, we are interested in finding representations of functions in terms of sines and cosines. Given a function $f(x)$ we seek a representation in the form

$$f(x) \sim \dfrac{a_{0}}{2}+\sum_{n=1}^{\infty}\left[a_{n} \cos n x+b_{n} \sin n x\right] \nonumber$$

Notice that we have opted to drop reference to the frequency form of the phase. This will lead to a simpler discussion for now and one can always make the transformation $n x=2 \pi f_{n} t$ when applying these ideas to applications.

The series representation in Equation (5.1) is called a Fourier trigonometric series. We will simply refer to this as a Fourier series for now. The set of constants $a_{0}, a_{n}, b_{n}, n=1,2, \ldots$ are called the Fourier coefficients. The constant term is chosen in this form to make later computations simpler, though some other authors choose to write the constant term as $a_{0}$. Our goal is to find the Fourier series representation given $f(x)$. Having found the Fourier series representation, we will be interested in determining when the Fourier series converges and to what function it converges.

From our discussion in the last section, we see that the infinite series is periodic. The largest period of the terms comes from the $n=1$ terms. The periods of $\cos x$ and $\sin x$ are $T=2 \pi$. Thus, the Fourier series has period $2 \pi$. This means that the series should be able to represent functions that are periodic of period $2 \pi$.

While this appears restrictive, we could also consider functions that are defined over one period. In Figure $5.4$ we show a function defined on $[0,2 \pi]$. In the same figure, we show its periodic extension. These are just copies of the original function shifted by the period and glued together. The extension can now be represented by a Fourier series and restricting the Fourier series to $[0,2 \pi]$ will give a representation of the original function. Therefore, we will first consider Fourier series representations of functions defined on this interval. Note that we could just as easily considered functions defined on $[-\pi, \pi]$ or any interval of length $2 \pi$.

Fourier Series Over Other Intervals

In many applications we are interested in determining Fourier series representations of functions defined on intervals other than $[0,2 \pi]$. In this section we will determine the form of the series expansion and the Fourier coefficients in these cases.

The most general type of interval is given as $[a, b]$. However, this often is too general. More common intervals are of the form $[-\pi, \pi],[0, L]$, or $[-L / 2, L / 2]$. The simplest generalization is to the interval $[0, L]$. Such intervals arise often in applications. For example, one can study vibrations of a one dimensional string of length $L$ and set up the axes with the left end at $x=0$ and the right end at $x=L$. Another problem would be to study the temperature distribution along a one dimensional rod of length $L$. Such problems lead to the original studies of Fourier series. As we will see later, symmetric intervals, $[-a, a]$, are also useful.

Given an interval $[0, L]$, we could apply a transformation to an interval of length $2 \pi$ by simply rescaling our interval. Then we could apply this transformation to our Fourier series representation to obtain an equivalent one useful for functions defined on $[0, L]$.

We define $x \in[0,2 \pi]$ and $t \in[0, L]$. A linear transformation relating these intervals is simply $x=\dfrac{2 \pi t}{L}$ as shown in Figure $5.5$. So, $t=0$ maps to $x=0$ and $t=L$ maps to $x=2 \pi$. Furthermore, this transformation maps $f(x)$ to a new function $g(t)=f(x(t))$, which is defined on $[0, L]$. We will determine the Fourier series representation of this function using the representation for $f(x)$

Recall the form of the Fourier representation for $f(x)$ in Equation (5.1):

$$f(x) \sim \dfrac{a_{0}}{2}+\sum_{n=1}^{\infty}\left[a_{n} \cos n x+b_{n} \sin n x\right] \nonumber$$

Inserting the transformation relating $x$ and $t$, we have

$$g(t) \sim \dfrac{a_{0}}{2}+\sum_{n=1}^{\infty}\left[a_{n} \cos \dfrac{2 n \pi t}{L}+b_{n} \sin \dfrac{2 n \pi t}{L}\right] . \nonumber$$

This gives the form of the series expansion for $g(t)$ with $t \in[0, L]$. But, we still need to determine the Fourier coefficients. Recall, that

$$a_{n}=\dfrac{1}{\pi} \int_{0}^{2 \pi} f(x) \cos n x d x \nonumber$$

We need to make a substitution in the integral of $x=\dfrac{2 \pi t}{L}$. We also will need to transform the differential, $d x=\dfrac{2 \pi}{L} d t$. Thus, the resulting form for our coefficient is

$$a_{n}=\dfrac{2}{L} \int_{0}^{L} g(t) \cos \dfrac{2 n \pi t}{L} d t \nonumber$$

Similarly, we find that

$$b_{n}=\dfrac{2}{L} \int_{0}^{L} g(t) \sin \dfrac{2 n \pi t}{L} d t \nonumber$$

We note first that when $L=2 \pi$ we get back the series representation that we first studied. Also, the period of $\cos \dfrac{2 n \pi t}{L}$ is $L / n$, which means that the representation for $g(t)$ has a period of $L$.

At the end of this section we present the derivation of the Fourier series representation for a general interval for the interested reader. In Table $5.1$ we summarize some commonly used Fourier series representations.

We will end our discussion for now with some special cases and an example for a function defined on $[-\pi, \pi]$.

Example 5.7. Let $f(x)=|x|$ on $[-\pi, \pi]$ We compute the coefficients, beginning as usual with $a_{0}$. We have

$$\begin{aligned} a_{0} &=\dfrac{1}{\pi} \int_{-\pi}^{\pi}|x| d x \\[4pt] &=\dfrac{2}{\pi} \int_{0}^{\pi}|x| d x=\pi \end{aligned} \nonumber$$

At this point we need to remind the reader about the integration of even and odd functions.

1. Even Functions: In this evaluation we made use of the fact that the integrand is an even function. Recall that $f(x)$ is an even function if $f(-x)=f(x)$ for all $x$. One can recognize even functions as they are symmetric with respect to the $y$-axis as shown in Figure 5.6(A). If one integrates an even function over a symmetric interval, then one has that

$$\int_{-a}^{a} f(x) d x=2 \int_{0}^{a} f(x) d x \nonumber$$

One can prove this by splitting off the integration over negative values of $x$, using the substitution $x=-y$, and employing the evenness of $f(x)$. Thus, Table 5.1. Special Fourier Series Representations on Different Intervals

Fourier Series on $[0, L]$

$$\begin{aligned} f(x) & \sim \dfrac{a_{0}}{2}+\sum_{n=1}^{\infty}\left[a_{n} \cos \dfrac{2 n \pi x}{L}+b_{n} \sin \dfrac{2 n \pi x}{L}\right] \\[4pt] a_{n} &=\dfrac{2}{L} \int_{0}^{L} f(x) \cos \dfrac{2 n \pi x}{L} d x . \quad n=0,1,2, \ldots, \\[4pt] b_{n} &=\dfrac{2}{L} \int_{0}^{L} f(x) \sin \dfrac{2 n \pi x}{L} d x . \quad n=1,2, \ldots \end{aligned} \nonumber$$

Fourier Series on $\left[-\dfrac{L}{2}, \dfrac{L}{2}\right]$

$$\begin{aligned} f(x) & \sim \dfrac{a_{0}}{2}+\sum_{n=1}^{\infty}\left[a_{n} \cos \dfrac{2 n \pi x}{L}+b_{n} \sin \dfrac{2 n \pi x}{L}\right] . \\[4pt] a_{n} &=\dfrac{2}{L} \int_{-\dfrac{L}{2}}^{\dfrac{L}{2}} f(x) \cos \dfrac{2 n \pi x}{L} d x . \quad n=0,1,2, \ldots, \\[4pt] b_{n} &=\dfrac{2}{L} \int_{-\dfrac{L}{2}}^{\dfrac{L}{2}} f(x) \sin \dfrac{2 n \pi x}{L} d x . \quad n=1,2, \ldots \end{aligned} \nonumber$$

Fourier Series on $[-\pi, \pi]$

$$\begin{gathered} f(x) \sim \dfrac{a_{0}}{2}+\sum_{n=1}^{\infty}\left[a_{n} \cos n x+b_{n} \sin n x\right] . \\[4pt] a_{n}=\dfrac{1}{\pi} \int_{-\pi}^{\pi} f(x) \cos n x d x . \quad n=0,1,2, \ldots, \\[4pt] b_{n}=\dfrac{1}{\pi} \int_{-\pi}^{\pi} f(x) \sin n x d x . \quad n=1,2, \ldots \end{gathered} \nonumber$$

$$\begin{aligned} \int_{-a}^{a} f(x) d x &=\int_{-a}^{0} f(x) d x+\int_{0}^{a} f(x) d x \\[4pt] &=-\int_{a}^{0} f(-y) d y+\int_{0}^{a} f(x) d x \\[4pt] &=\int_{0}^{a} f(y) d y+\int_{0}^{a} f(x) d x \\[4pt] &=2 \int_{0}^{a} f(x) d x \end{aligned} \nonumber$$

This can be visually verified by looking at Figure $5.6(\mathrm{~A})$.

3. Odd Functions: A similar computation could be done for odd functions. $f(x)$ is an odd function if $f(-x)=-f(x)$ for all $x$. The graphs of such functions are symmetric with respect to the origin as shown in Figure $5.6(\mathrm{~B})$. If one integrates an odd function over a symmetric interval, then one has that

We now continue with our computation of the Fourier coefficients for $f(x)=|x|$ on $[-\pi, \pi]$. We have

$$a_{n}=\dfrac{1}{\pi} \int_{-\pi}^{\pi}|x| \cos n x d x=\dfrac{2}{\pi} \int_{0}^{\pi} x \cos n x d x . \nonumber$$

Here we have made use of the fact that $|x| \cos n x$ is an even function. In order to compute the resulting integral, we need to use integration by parts,

$$\int_{a}^{b} u d v=\left.u v\right|_{a} ^{b}-\int_{a}^{b} v d u \nonumber$$

by letting $u=x$ and $d v=\cos n x d x$. Thus, $d u=d x$ and $v=\int d v=\dfrac{1}{n} \sin n x$. Continuing with the computation, we have

$$\begin{aligned} a_{n} &=\dfrac{2}{\pi} \int_{0}^{\pi} x \cos n x d x \\[4pt] &=\dfrac{2}{\pi}\left[\left.\dfrac{1}{n} x \sin n x\right|_{0} ^{\pi}-\dfrac{1}{n} \int_{0}^{\pi} \sin n x d x\right] \\[4pt] &=-\dfrac{2}{n \pi}\left[-\dfrac{1}{n} \cos n x\right]_{0}^{\pi} \\[4pt] &=-\dfrac{2}{\pi n^{2}}\left(1-(-1)^{n}\right) \end{aligned} \nonumber$$

Here we have used the fact that $\cos n \pi=(-1)^{n}$ for any integer $n$. This lead to a factor $\left(1-(-1)^{n}\right)$. This factor can be simplified as

$$1-(-1)^{n}=\left\{\begin{array}{l} 2, n \text { odd } \\[4pt] 0, n \text { even } \end{array}\right. \nonumber$$

So, $a_{n}=0$ for $n$ even and $a_{n}=-\dfrac{4}{\pi n^{2}}$ for $n$ odd.

Computing the $b_{n}$ ’s is simpler. We note that we have to integrate $|x| \sin n x$ from $x=-\pi$ to $\pi$. The integrand is an odd function and this is a symmetric interval. So, the result is that $b_{n}=0$ for all $n$.

Putting this all together, the Fourier series representation of $f(x)=|x|$ on $[-\pi, \pi]$ is given as

$$f(x) \sim \dfrac{\pi}{2}-\dfrac{4}{\pi} \sum_{n=1, \text { odd }}^{\infty} \dfrac{\cos n x}{n^{2}} \nonumber$$

While this is correct, we can rewrite the sum over only odd $n$ by reindexing. We let $n=2 k-1$ for $k=1,2,3, \ldots$ Then we only get the odd integers. The series can then be written as

$$f(x) \sim \dfrac{\pi}{2}-\dfrac{4}{\pi} \sum_{k=1}^{\infty} \dfrac{\cos (2 k-1) x}{(2 k-1)^{2}} . \nonumber$$

Throughout our discussion we have referred to such results as Fourier representations. We have not looked at the convergence of these series. Here is an example of an infinite series of functions. What does this series sum to? We show in Figure $5.7$ the first few partial sums. They appear to be converging to $f(x)=|x|$ fairly quickly.

Even though $f(x)$ was defined on $[-\pi, \pi]$ we can still evaluate the Fourier series at values of $x$ outside this interval. In Figure $5.8$, we see that the representation agrees with $f(x)$ on the interval $[-\pi, \pi]$. Outside this interval we have a periodic extension of $f(x)$ with period $2 \pi$.

Another example is the Fourier series representation of $f(x)=x$ on $[-\pi, \pi]$ as left for Problem 5.1. This is determined to be

$$f(x) \sim 2 \sum_{n=1}^{\infty} \dfrac{(-1)^{n+1}}{n} \sin n x . \nonumber$$

As seen in Figure $5.9$ we again obtain the periodic extension of our function. In this case we needed many more terms. Also, the vertical parts of the first plot are nonexistent. In the second plot we only plot the points and not the typical connected points that most software packages plot as the default style.

$$\pi=4\left[1-\dfrac{1}{3}+\dfrac{1}{5}-\dfrac{1}{7}+\ldots\right] \nonumber$$

Appendix: The Gibbs Phenomenon

We have seen that when there is a jump discontinuity in the periodic extension of our functions, whether the function originally had a discontinuity or developed one due to a mismatch in the values of the endpoints. This can be seen in Figures $5.9,5.11$ and 5.13. The Fourier series has a difficult time converging at the point of discontinuity and these graphs of the Fourier series show a distinct overshoot which does not go away. This is called the Gibbs phenomenon and the amount of overshoot can be computed.

In one of our first examples, Example $5.6$, we found the Fourier series representation of the piecewise defined function

$$f(x)=\left\{\begin{array}{c} 1, \quad 0<x<\pi \\[4pt] -1, \quad \pi<x<2 \pi \end{array}\right. \nonumber$$

to be

$$f(x) \sim \dfrac{4}{\pi} \sum_{k=1}^{\infty} \dfrac{\sin (2 k-1) x}{2 k-1} . \nonumber$$

In Figure $5.14$ we display the sum of the first ten terms. Note the wiggles, overshoots and under shoots near $x=0, \pm \pi$. These are seen more when we plot the representation for $x \in[-3 \pi, 3 \pi]$, as shown in Figure 5.15. We note that the overshoots and undershoots occur at discontinuities in the periodic extension of $f(x)$. These occur whenever $f(x)$ has a discontinuity or if the values of $f(x)$ at the endpoints of the domain do not agree.

One might expect that we only need to add more terms. In Figure $5.16$ we show the sum for twenty terms. Note the sum appears to converge better for points far from the discontinuities. But, the overshoots and undershoots are still present. In Figures $5.17$ and $5.18$ show magnified plots of the overshoot at $x=0$ for $N=100$ and $N=500$, respectively. We see that the overshoot persists. The peak is at about the same height, but its location seems to be getting closer to the origin. We will show how one can estimate the size of the overshoot.

We can study the Gibbs phenomenon by looking at the partial sums of general Fourier trigonometric series for functions $f(x)$ defined on the interval $[-L, L]$. Writing out the partial sums, inserting the Fourier coefficients and rearranging, we have

$$\begin{aligned} =& \dfrac{1}{2 L} \int_{-L}^{L} f(y) d y+\sum_{n=1}^{N}\left[\left(\dfrac{1}{L} \int_{-L}^{L} f(y) \cos \dfrac{n \pi y}{L} d y\right) \cos \dfrac{n \pi x}{L}\right.\\[4pt] &\left.+\left(\dfrac{1}{L} \int_{-L}^{L} f(y) \sin \dfrac{n \pi y}{L} d y\right) \sin \dfrac{n \pi x}{L}\right] \\[4pt] =&\left.\dfrac{1}{L} \int_{-L}^{L}\left\{\dfrac{1}{2}+\sum_{n=1}^{N}\left(\cos \dfrac{n \pi y}{L} \cos \dfrac{n \pi x}{L}+\sin \dfrac{n \pi y}{L} \sin \dfrac{n \pi x}{L}\right)\right\} f(y)\right\} f(y) d y \\[4pt] =& \dfrac{1}{L} \int_{-L}^{L}\left\{\dfrac{1}{2}+\sum_{n=1}^{N} \cos \dfrac{n \pi(y-x)}{L}\right\} \\[4pt] \equiv & \dfrac{1}{L} \int_{-L}^{L} D_{N}(y-x) f(y) d y . \end{aligned} \nonumber$$

We have defined

$$D_{N}(x)=\dfrac{1}{2}+\sum_{n=1}^{N} \cos \dfrac{n \pi x}{L}, \nonumber$$

which is called the $N$-th Dirichlet Kernel. We now prove

Introduction

In the last chapters we have explored the solution of boundary value problems that led to trigonometric eigenfunctions. Such functions can be used to represent functions in Fourier series expansions. We would like to generalize some of those techniques in order to solve other boundary value problems. A class of problems to which our previous examples belong and which have eigenfunctions with similar properties are the Sturm-Liouville Eigenvalue Problems. These problems involve self-adjoint (differential) operators which play an important role in the spectral theory of linear operators and the existence of the eigenfunctions we described in Section 4.3.2. These ideas will be introduced in this chapter.

In physics many problems arise in the form of boundary value problems involving second order ordinary differential equations. For example, we might want to solve the equation

$$a_{2}(x) y^{\prime \prime}+a_{1}(x) y^{\prime}+a_{0}(x) y=f(x) \nonumber$$

subject to boundary conditions. We can write such an equation in operator form by defining the differential operator

$$L=a_{2}(x) \dfrac{d^{2}}{d x^{2}}+a_{1}(x) \dfrac{d}{d x}+a_{0}(x) . \nonumber$$

Then, Equation (6.1) takes the form

$$L y=f \text {. } \nonumber$$

As we saw in the general boundary value problem (4.20) in Section 4.3.2, we can solve some equations using eigenvalue expansions. Namely, we seek solutions to the eigenvalue problem

$$L \phi=\lambda \phi \nonumber$$

with homogeneous boundary conditions and then seek a solution as an expansion of the eigenfunctions. Formally, we let

$$y=\sum_{n=1}^{\infty} c_{n} \phi_{n} \nonumber$$

However, we are not guaranteed a nice set of eigenfunctions. We need an appropriate set to form a basis in the function space. Also, it would be nice to have orthogonality so that we can easily solve for the expansion coefficients as was done in Section 4.3.2. [Otherwise, we would have to solve a infinite coupled system of algebraic equations instead of an uncoupled and diagonal system.]

It turns out that any linear second order operator can be turned into an operator that possesses just the right properties (self-adjointedness to carry out this procedure. The resulting operator is referred to as a Sturm-Liouville operator. We will highlight some of the properties of such operators and prove a few key theorems, though this will not be an extensive review of SturmLiouville theory. The interested reader can review the literature and more advanced texts for a more in depth analysis.

We define the Sturm-Liouville operator as

$$\mathcal{L}=\dfrac{d}{d x} p(x) \dfrac{d}{d x}+q(x) \nonumber$$

The Sturm-Liouville eigenvalue problem is given by the differential equation

$$\mathcal{L} u=-\lambda \sigma(x) u \nonumber$$

$$\dfrac{d}{d x}\left(p(x) \dfrac{d u}{d x}\right)+q(x) u+\lambda \sigma(x) u=0 \nonumber$$

for $x \in(a, b)$. The functions $p(x), p^{\prime}(x), q(x)$ and $\sigma(x)$ are assumed to be continuous on $(a, b)$ and $p(x)>0, \sigma(x)>0$ on $[a, b]$. If the interval is finite and these assumptions on the coefficients are true on $[a, b]$, then the problem is said to be regular. Otherwise, it is called singular.

We also need to impose the set of homogeneous boundary conditions

$$\begin{array}{r} \alpha_{1} u(a)+\beta_{1} u^{\prime}(a)=0 \\[4pt] \alpha_{2} u(b)+\beta_{2} u^{\prime}(b)=0 \end{array} \nonumber$$

The $\alpha$ ’s and $\beta$ ’s are constants. For different values, one has special types of boundary conditions. For $\beta_{i}=0$, we have what are called Dirichlet boundary conditions. Namely, $u(a)=0$ and $u(b)=0$. For $\alpha_{i}=0$, we have Neumann boundary conditions. In this case, $u^{\prime}(a)=0$ and $u^{\prime}(b)=0$. In terms of the heat equation example, Dirichlet conditions correspond to maintaining a fixed temperature at the ends of the rod. The Neumann boundary conditions would correspond to no heat flow across the ends, or insulating conditions, as there would be no temperature gradient at those points. The more general boundary conditions allow for partially insulated boundaries.

Another type of boundary condition that is often encountered is the periodic boundary condition. Consider the heated rod that has been bent to form a circle. Then the two end points are physically the same. So, we would expect that the temperature and the temperature gradient should agree at those points. For this case we write $u(a)=u(b)$ and $u^{\prime}(a)=u^{\prime}(b)$. Boundary value problems using these conditions have to be handled differently than the above homogeneous conditions. These conditions leads to different types of eigenfunctions and eigenvalues.

As previously mentioned, equations of the form (6.1) occur often. We now show that Equation (6.1) can be turned into a differential equation of SturmLiouville form:

$$\dfrac{d}{d x}\left(p(x) \dfrac{d y}{d x}\right)+q(x) y=F(x) \nonumber$$

Another way to phrase this is provided in the theorem:

Theorem 6.1. Any second order linear operator can be put into the form of the Sturm-Liouville operator (6.2).

The proof of this is straight forward, as we shall soon show. Consider the equation (6.1). If $a_{1}(x)=a_{2}^{\prime}(x)$, then we can write the equation in the form

$$\begin{aligned} f(x) &=a_{2}(x) y^{\prime \prime}+a_{1}(x) y^{\prime}+a_{0}(x) y \\[4pt] &=\left(a_{2}(x) y^{\prime}\right)^{\prime}+a_{0}(x) y \end{aligned} \nonumber$$

This is in the correct form. We just identify $p(x)=a_{2}(x)$ and $q(x)=a_{0}(x)$.

However, consider the differential equation

$$x^{2} y^{\prime \prime}+x y^{\prime}+2 y=0 . \nonumber$$

In this case $a_{2}(x)=x^{2}$ and $a_{2}^{\prime}(x)=2 x \neq a_{1}(x)$. The linear differential operator in this equation is not of Sturm-Liouville type. But, we can change it to a Sturm Liouville operator.

In the Sturm Liouville operator the derivative terms are gathered together into one perfect derivative. This is similar to what we saw in the first chapter when we solved linear first order equations. In that case we sought an integrating factor. We can do the same thing here. We seek a multiplicative function $\mu(x)$ that we can multiply through $(6.1)$ so that it can be written in Sturm-Liouville form. We first divide out the $a_{2}(x)$, giving

$$y^{\prime \prime}+\dfrac{a_{1}(x)}{a_{2}(x)} y^{\prime}+\dfrac{a_{0}(x)}{a_{2}(x)} y=\dfrac{f(x)}{a_{2}(x)} . \nonumber$$

Now, we multiply the differential equation by $\mu$ :

$$\mu(x) y^{\prime \prime}+\mu(x) \dfrac{a_{1}(x)}{a_{2}(x)} y^{\prime}+\mu(x) \dfrac{a_{0}(x)}{a_{2}(x)} y=\mu(x) \dfrac{f(x)}{a_{2}(x)} \nonumber$$

The first two terms can now be combined into an exact derivative $\left(\mu y^{\prime}\right)^{\prime}$ if $\mu(x)$ satisfies

$$\dfrac{d \mu}{d x}=\mu(x) \dfrac{a_{1}(x)}{a_{2}(x)} . \nonumber$$

This is formally solved to give

$$\mu(x)=e^{\int \dfrac{a_{1}(x)}{a_{2}(x)} d x} . \nonumber$$

Thus, the original equation can be multiplied by factor

$$\dfrac{\mu(x)}{a_{2}(x)}=\dfrac{1}{a_{2}(x)} e^{\int \dfrac{a_{1}(x)}{a_{2}(x)} d x} \nonumber$$

to turn it into Sturm-Liouville form.

In summary,

$$\begin{aligned} &\text { Equation (6.1), } \\[4pt] &\qquad a_{2}(x) y^{\prime \prime}+a_{1}(x) y^{\prime}+a_{0}(x) y=f(x) \end{aligned} \nonumber$$

can be put into the Sturm-Liouville form

$$\dfrac{d}{d x}\left(p(x) \dfrac{d y}{d x}\right)+q(x) y=F(x) \nonumber$$

where

$$\begin{aligned} p(x) &=e^{\int \dfrac{a_{1}(x)}{a_{2}(x)} d x} \\[4pt] q(x) &=p(x) \dfrac{a_{0}(x)}{a_{2}(x)} \\[4pt] F(x) &=p(x) \dfrac{f(x)}{a_{2}(x)} \end{aligned} \nonumber$$

Example 6.2. For the example above,

$$x^{2} y^{\prime \prime}+x y^{\prime}+2 y=0 . \nonumber$$

We need only multiply this equation by

$$\dfrac{1}{x^{2}} e^{\int \dfrac{d x}{x}}=\dfrac{1}{x}, \nonumber$$

to put the equation in Sturm-Liouville form:

$$\begin{aligned} 0 &=x y^{\prime \prime}+y^{\prime}+\dfrac{2}{x} y \\[4pt] &=\left(x y^{\prime}\right)^{\prime}+\dfrac{2}{x} y \end{aligned} \nonumber$$

Properties of Sturm-Liouville Eigenvalue Problems

There are several properties that can be proven for the (regular) SturmLiouville eigenvalue problem. However, we will not prove them all here. We will merely list some of the important facts and focus on a few of the properties.

1. The eigenvalues are real, countable, ordered and there is a smallest eigenvalue. Thus, we can write them as $\lambda_{1}<\lambda_{2}<\ldots$. However, there is no largest eigenvalue and $n \rightarrow \infty, \lambda_{n} \rightarrow \infty$.
2. For each eigenvalue $\lambda_{n}$ there exists an eigenfunction $\phi_{n}$ with $n-1$ zeros on $(a, b)$.
3. Eigenfunctions corresponding to different eigenvalues are orthogonal with respect to the weight function, $\sigma(x)$. Defining the inner product of $f(x)$ and $g(x)$ as

$$<f, g>=\int_{a}^{b} f(x) g(x) \sigma(x) d x \nonumber$$

then the orthogonality of the eigenfunctios can be written in the form

$$<\phi_{n}, \phi_{m}>=<\phi_{n}, \phi_{n}>\delta_{n m}, \quad n, m=1,2, \ldots \nonumber$$

5. The set of eigenfunctions is complete; i.e., any piecewise smooth function can be represented by a generalized Fourier series expansion of the eigenfunctions,

$$f(x) \sim \sum_{n=1}^{\infty} c_{n} \phi_{n}(x) \nonumber$$

where

$$c_{n}=\dfrac{<f, \phi_{n}>}{<\phi_{n}, \phi_{n}>} \nonumber$$

Actually, one needs $f(x) \in L_{\sigma}^{2}[a, b]$, the set of square integrable functions over $[a, b]$ with weight function $\sigma(x)$. By square integrable, we mean that $<f, f><\infty$. One can show that such a space is isomorphic to a Hilbert space, a complete inner product space.

6. Multiply the eigenvalue problem

$$\mathcal{L} \phi_{n}=-\lambda_{n} \sigma(x) \phi_{n} \nonumber$$

by $\phi_{n}$ and integrate. Solve this result for $\lambda_{n}$, to find the Rayleigh Quotient

$$\lambda_{n}=\dfrac{-\left.p \phi_{n} \dfrac{d \phi_{n}}{d x}\right|_{a} ^{b}-\int_{a}^{b}\left[p\left(\dfrac{d \phi_{n}}{d x}\right)^{2}-q \phi_{n}^{2}\right] d x}{<\phi_{n}, \phi_{n}>} \nonumber$$

The Rayleigh quotient is useful for getting estimates of eigenvalues and proving some of the other properties. Example 6.3. We seek the eigenfunctions of the operator found in Example 6.2. Namely, we want to solve the eigenvalue problem

$$\mathcal{L} y=\left(x y^{\prime}\right)^{\prime}+\dfrac{2}{x} y=-\lambda \sigma y \nonumber$$

subject to a set of boundary conditions. Let’s use the boundary conditions

$$y^{\prime}(1)=0, \quad y^{\prime}(2)=0 \nonumber$$

[Note that we do not know $\sigma(x)$ yet, but will choose an appropriate function to obtain solutions.]

Expanding the derivative, we have

$$x y^{\prime \prime}+y^{\prime}+\dfrac{2}{x} y=-\lambda \sigma y . \nonumber$$

Multiply through by $x$ to obtain

$$x^{2} y^{\prime \prime}+x y^{\prime}+(2+\lambda x \sigma) y=0 \nonumber$$

Notice that if we choose $\sigma(x)=x^{-1}$, then this equation can be made a Cauchy-Euler type equation. Thus, we have

$$x^{2} y^{\prime \prime}+x y^{\prime}+(\lambda+2) y=0 . \nonumber$$

The characteristic equation is

$$r^{2}+\lambda+2=0 . \nonumber$$

For oscillatory solutions, we need $\lambda+2>0$. Thus, the general solution is

$$y(x)=c_{1} \cos (\sqrt{\lambda+2} \ln |x|)+c_{2} \sin (\sqrt{\lambda+2} \ln |x|) . \nonumber$$

Next we apply the boundary conditions. $y^{\prime}(1)=0$ forces $c_{2}=0$. This leaves

$$y(x)=c_{1} \cos (\sqrt{\lambda+2} \ln x) . \nonumber$$

The second condition, $y^{\prime}(2)=0$, yields

$$\sin (\sqrt{\lambda+2} \ln 2)=0 \nonumber$$

This will give nontrivial solutions when

$$\sqrt{\lambda+2} \ln 2=n \pi, \quad n=0,1,2,3 \ldots \nonumber$$

In summary, the eigenfunctions for this eigenvalue problem are

$$y_{n}(x)=\cos \left(\dfrac{n \pi}{\ln 2} \ln x\right), \quad 1 \leq x \leq 2 \nonumber$$

and the eigenvalues are $\lambda_{n}=2+\left(\dfrac{n \pi}{\ln 2}\right)^{2}$ for $n=0,1,2, \ldots$

Note: We include the $n=0$ case because $y(x)=$ constant is a solution of the $\lambda=-2$ case. More specifically, in this case the characteristic equation reduces to $r^{2}=0$. Thus, the general solution of this Cauchy-Euler equation is

$$y(x)=c_{1}+c_{2} \ln |x| \nonumber$$

Setting $y^{\prime}(1)=0$, forces $c_{2}=0 . y^{\prime}(2)$ automatically vanishes, leaving the solution in this case as $y(x)=c_{1}$.

We note that some of the properties listed in the beginning of the section hold for this example. The eigenvalues are seen to be real, countable and ordered. There is a least one, $\lambda=2$. Next, one can find the zeros of each eigenfunction on $[1,2]$. Then the argument of the cosine, $\dfrac{n \pi}{\ln 2} \ln x$, takes values 0 to $n \pi$ for $x \in[1,2]$. The cosine function has $n-1$ roots on this interval.

Orthogonality can be checked as well. We set up the integral and use the substitution $y=\pi \ln x / \ln 2$. This gives

$$\begin{aligned} <y_{n}, y_{m}>&=\int_{1}^{2} \cos \left(\dfrac{n \pi}{\ln 2} \ln x\right) \cos \left(\dfrac{m \pi}{\ln 2} \ln x\right) \dfrac{d x}{x} \\[4pt] &=\dfrac{\ln 2}{\pi} \int_{0}^{\pi} \cos n y \cos m y d y \\[4pt] &=\dfrac{\ln 2}{2} \delta_{n, m} \end{aligned} \nonumber$$

The Eigenfunction Expansion Method

In section $4.3 .2$ we saw generally how one can use the eigenfunctions of a differential operator to solve a nonhomogeneous boundary value problem. In this chapter we have seen that Sturm-Liouville eigenvalue problems have the requisite set of orthogonal eigenfunctions. In this section we will apply the eigenfunction expansion method to solve a particular nonhomogenous boundary value problem.

Recall that one starts with a nonhomogeneous differential equation

$$\mathcal{L} y=f, \nonumber$$

where $y(x)$ is to satisfy given homogeneous boundary conditions. The method makes use of the eigenfunctions satisfying the eigenvalue problem

$$\mathcal{L} \phi_{n}=-\lambda_{n} \sigma \phi_{n} \nonumber$$

subject to the given boundary conditions. Then, one assumes that $y(x)$ can be written as an expansion in the eigenfunctions,

$$y(x)=\sum_{n=1}^{\infty} c_{n} \phi_{n}(x), \nonumber$$

and inserts the expansion into the nonhomogeneous equation. This gives

$$f(x)=\mathcal{L}\left(\sum_{n=1}^{\infty} c_{n} \phi_{n}(x)\right)=-\sum_{n=1}^{\infty} c_{n} \lambda_{n} \sigma(x) \phi_{n}(x) \nonumber$$

The expansion coefficients are then found by making use of the orthogonality of the eigenfunctions. Namely, we multiply the last equation by $\phi_{m}(x)$ and integrate. We obtain

$$\int_{a}^{b} f(x) \phi_{m}(x) d x=-\sum_{n=1}^{\infty} c_{n} \lambda_{n} \int_{a}^{b} \phi_{n}(x) \phi_{m}(x) \sigma(x) d x \nonumber$$

Orthogonality yields

$$\int_{a}^{b} f(x) \phi_{m}(x) d x=-c_{m} \lambda_{m} \int_{a}^{b} \phi_{m}^{2}(x) \sigma(x) d x \nonumber$$

Solving for $c_{m}$, we have

$$c_{m}=-\dfrac{\int_{a}^{b} f(x) \phi_{m}(x) d x}{\lambda_{m} \int_{a}^{b} \phi_{m}^{2}(x) \sigma(x) d x} . \nonumber$$

Example 6.11. As an example, we consider the solution of the boundary value problem

$$\begin{aligned} \left(x y^{\prime}\right)^{\prime}+\dfrac{y}{x} &=\dfrac{1}{x}, \quad x \in[1, e], \\[4pt] y(1) &=0=y(e) . \end{aligned} \nonumber$$

This equation is already in self-adjoint form. So, we know that the associated Sturm-Liouville eigenvalue problem has an orthogonal set of eigenfunctions. We first determine this set. Namely, we need to solve

$$\left(x \phi^{\prime}\right)^{\prime}+\dfrac{\phi}{x}=-\lambda \sigma \phi, \quad \phi(1)=0=\phi(e) . \nonumber$$

Rearranging the terms and multiplying by $x$, we have that

$$x^{2} \phi^{\prime \prime}+x \phi^{\prime}+(1+\lambda \sigma x) \phi=0 . \nonumber$$

This is almost an equation of Cauchy-Euler type. Picking the weight function $\sigma(x)=\dfrac{1}{x}$, we have

$$x^{2} \phi^{\prime \prime}+x \phi^{\prime}+(1+\lambda) \phi=0 . \nonumber$$

This is easily solved. The characteristic equation is

$$r^{2}+(1+\lambda)=0 . \nonumber$$

One obtains nontrivial solutions of the eigenvalue problem satisfying the boundary conditions when $\lambda>-1$. The solutions are

$$\phi_{n}(x)=A \sin (n \pi \ln x), \quad n=1,2, \ldots \nonumber$$

where $\lambda_{n}=n^{2} \pi^{2}-1$

It is often useful to normalize the eigenfunctions. This means that one chooses $A$ so that the norm of each eigenfunction is one. Thus, we have

$$\begin{aligned} 1 &=\int_{1}^{e} \phi_{n}(x)^{2} \sigma(x) d x \\[4pt] &=A^{2} \int_{1}^{e} \sin (n \pi \ln x) \dfrac{1}{x} d x \\[4pt] &=A^{2} \int_{0}^{1} \sin (n \pi y) d y=\dfrac{1}{2} A^{2} \end{aligned} \nonumber$$

Thus, $A=\sqrt{2}$

We now turn towards solving the nonhomogeneous problem, $\mathcal{L} y=\dfrac{1}{x}$. We first expand the unknown solution in terms of the eigenfunctions,

$$y(x)=\sum_{n=1}^{\infty} c_{n} \sqrt{2} \sin (n \pi \ln x) . \nonumber$$

Inserting this solution into the differential equation, we have

$$\dfrac{1}{x}=\mathcal{L} y=-\sum_{n=1}^{\infty} c_{n} \lambda_{n} \sqrt{2} \sin (n \pi \ln x) \dfrac{1}{x} \nonumber$$

Next, we make use of orthogonality. Multiplying both sides by $\phi_{m}(x)=$ $\sqrt{2} \sin (m \pi \ln x)$ and integrating, gives

$$\lambda_{m} c_{m}=\int_{1}^{e} \sqrt{2} \sin (m \pi \ln x) \dfrac{1}{x} d x=\dfrac{\sqrt{2}}{m \pi}\left[(-1)^{m}-1\right] . \nonumber$$

Solving for $c_{m}$, we have

$$c_{m}=\dfrac{\sqrt{2}}{m \pi} \dfrac{\left[(-1)^{m}-1\right]}{m^{2} \pi^{2}-1} . \nonumber$$

Finally, we insert our coefficients into the expansion for $y(x)$. The solution is then

$$y(x)=\sum_{n=1}^{\infty} \dfrac{2}{n \pi} \dfrac{\left[(-1)^{n}-1\right]}{n^{2} \pi^{2}-1} \sin (n \pi \ln (x)) . \nonumber$$

The Fredholm Alternative Theorem

Given that $L y=f$, when can one expect to find a solution? Is it unique? These questions are answered by the Fredholm Alternative Theorem. This theorem occurs in many forms from a statement about solutions to systems of algebraic equations to solutions of boundary value problems and integral equations. The theorem comes in two parts, thus the term "alternative". Either the equation has exactly one solution for all $f$, or the equation has many solutions for some $f^{\prime}$ ’s and none for the rest.

The reader is familiar with the statements of the Fredholm Alternative for the solution of systems of algebraic equations. One seeks solutions of the system $A x=b$ for $A$ an $n \times m$ matrix. Defining the matrix adjoint, $A^{*}$ through $<A x, y>=<x, A^{*} y>$ for all $x, y, \in \mathcal{C}^{n}$, then either

Classical Orthogonal Polynomials

We begin by noting that the sequence of functions $\left\{1, x, x^{2}, \ldots\right\}$ is a basis of linearly independent functions. In fact, by the Stone-Weierstrass Approximation Theorem this set is a basis of $L_{\sigma}^{2}(a, b)$, the space of square integrable functions over the interval $[a, b]$ relative to weight $\sigma(x)$. We are familiar with being able to expand functions over this basis, since the expansions are just power series representations of the functions,

$$f(x) \sim \sum_{n=0}^{\infty} c_{n} x^{n} . \nonumber$$

However, this basis is not an orthogonal set of basis functions. One can easily see this by integrating the product of two even, or two odd, basis functions with $\sigma(x)=1$ and $(a, b)=(-1,1)$. For example,

$$<1, x^{2}>=\int_{-1}^{1} x^{0} x^{2} d x=\dfrac{2}{3} \nonumber$$

Since we have found that orthogonal bases have been useful in determining the coefficients for expansions of given functions, we might ask if it is possible to obtain an orthogonal basis involving these powers of $x$. Of course, finite combinations of these basis element are just polynomials!

$\mathrm{OK}$, we will ask. "Given a set of linearly independent basis vectors, can one find an orthogonal basis of the given space?" The answer is yes. We recall from introductory linear algebra, which mostly covers finite dimensional vector spaces, that there is a method for carrying this out called the GramSchmidt Orthogonalization Process. We will recall this process for finite dimensional vectors and then generalize to function spaces.

Let’s assume that we have three vectors that span $\mathbf{R}^{3}$, given by $\mathbf{a}_{1}, \mathbf{a}_{2}$, and $\mathbf{a}_{3}$ and shown in Figure 7.1. We seek an orthogonal basis $\mathbf{e}_{1}, \mathbf{e}_{2}$, and $\mathbf{e}_{3}$, beginning one vector at a time.

First we take one of the original basis vectors, say $\mathbf{a}_{1}$, and define

$$\mathbf{e}_{1}=\mathbf{a}_{1} \nonumber$$

Of course, we might want to normalize our new basis vectors, so we would denote such a normalized vector with a "hat":

$$\hat{\mathbf{e}}_{1}=\dfrac{\mathbf{e}_{1}}{e_{1}}, \nonumber$$

where $e_{1}=\sqrt{\mathbf{e}_{1} \cdot \mathbf{e}_{1}}$.

Note that this is easily proven by writing the projection as a vector of length $a_{2} \cos \theta$ in direction $\hat{\mathbf{e}}_{1}$, where $\theta$ is the angle between $\mathbf{e}_{1}$ and $\mathbf{a}_{2}$. Using the definition of the dot product, $\mathbf{a} \cdot \mathbf{b}=a b \cos \theta$, the projection formula follows.

Combining Equations (7.1)-(7.2), we find that

$$\mathbf{e}_{2}=\mathbf{a}_{2}-\dfrac{\mathbf{a}_{2} \cdot \mathbf{e}_{1}}{e_{1}^{2}} \mathbf{e}_{1} \nonumber$$

It is a simple matter to verify that $\mathbf{e}_{2}$ is orthogonal to $\mathbf{e}_{1}$ :

$$\begin{aligned} \mathbf{e}_{2} \cdot \mathbf{e}_{1} &=\mathbf{a}_{2} \cdot \mathbf{e}_{1}-\dfrac{\mathbf{a}_{2} \cdot \mathbf{e}_{1}}{e_{1}^{2}} \mathbf{e}_{1} \cdot \mathbf{e}_{1} \\[4pt] &=\mathbf{a}_{2} \cdot \mathbf{e}_{1}-\mathbf{a}_{2} \cdot \mathbf{e}_{1}=0 \end{aligned} \nonumber$$

Now, we seek a third vector $\mathbf{e}_{3}$ that is orthogonal to both $\mathbf{e}_{1}$ and $\mathbf{e}_{2}$. Pictorially, we can write the given vector $\mathbf{a}_{3}$ as a combination of vector projections along $\mathbf{e}_{1}$ and $\mathbf{e}_{2}$ and the new vector. This is shown in Figure 7.3. Then we have,

$$\mathbf{e}_{3}=\mathbf{a}_{3}-\dfrac{\mathbf{a}_{3} \cdot \mathbf{e}_{1}}{e_{1}^{2}} \mathbf{e}_{1}-\dfrac{\mathbf{a}_{3} \cdot \mathbf{e}_{2}}{e_{2}^{2}} \mathbf{e}_{2} . \nonumber$$

Again, it is a simple matter to compute the scalar products with $\mathbf{e}_{1}$ and $\mathbf{e}_{2}$ to verify orthogonality.

We can easily generalize the procedure to the $N$-dimensional case.

The Rodrigues Formula

The first property that the Legendre polynomials have is the Rodrigues formula:

$$P_{n}(x)=\dfrac{1}{2^{n} n !} \dfrac{d^{n}}{d x^{n}}\left(x^{2}-1\right)^{n}, \quad n \in N_{0} . \nonumber$$

From the Rodrigues formula, one can show that $P_{n}(x)$ is an $n$th degree polynomial. Also, for $n$ odd, the polynomial is an odd function and for $n$ even, the polynomial is an even function.

As an example, we determine $P_{2}(x)$ from Rodrigues formula:

$$\begin{aligned} P_{2}(x) &=\dfrac{1}{2^{2} 2 !} \dfrac{d^{2}}{d x^{2}}\left(x^{2}-1\right)^{2} \\[4pt] &=\dfrac{1}{8} \dfrac{d^{2}}{d x^{2}}\left(x^{4}-2 x^{2}+1\right) \\[4pt] &=\dfrac{1}{8} \dfrac{d}{d x}\left(4 x^{3}-4 x\right) \\[4pt] &=\dfrac{1}{8}\left(12 x^{2}-4\right) \\[4pt] &=\dfrac{1}{2}\left(3 x^{2}-1\right) \end{aligned} \nonumber$$

Note that we get the same result as we found in the last section using orthogonalization.

One can systematically generate the Legendre polynomials in tabular form as shown in Table $7.2 .1$. In Figure $7.4$ we show a few Legendre polynomials.

Table 7.3. Tabular computation of the Legendre polynomials using the Rodrigues formula.

Three Term Recursion Formula

The classical orthogonal polynomials also satisfy three term recursion formulae. In the case of the Legendre polynomials, we have

$$(2 n+1) x P_{n}(x)=(n+1) P_{n+1}(x)+n P_{n-1}(x), \quad n=1,2, \ldots \nonumber$$

This can also be rewritten by replacing $n$ with $n-1$ as

$$(2 n-1) x P_{n-1}(x)=n P_{n}(x)+(n-1) P_{n-2}(x), \quad n=1,2, \ldots \nonumber$$

We will prove this recursion formula in two ways. First we use the orthogonality properties of Legendre polynomials and the following lemma.

Lemma 7.2. The leading coefficient of $x^{n}$ in $P_{n}(x)$ is $\dfrac{1}{2^{n} n !} \dfrac{(2 n) !}{n !}$.

Proof. We can prove this using Rodrigues formula. first, we focus on the leading coefficient of $\left(x^{2}-1\right)^{n}$, which is $x^{2 n}$. The first derivative of $x^{2 n}$ is $2 n x^{2 n-1}$. The second derivative is $2 n(2 n-1) x^{2 n-2}$. The $j$ th derivative is

$$\dfrac{d^{j} x^{2 n}}{d x^{j}}=[2 n(2 n-1) \ldots(2 n-j+1)] x^{2 n-j} \nonumber$$

Thus, the $n$th derivative is given by

$$\dfrac{d^{n} x^{2 n}}{d x^{n}}=[2 n(2 n-1) \ldots(n+1)] x^{n} \nonumber$$

This proves that $P_{n}(x)$ has degree $n$. The leading coefficient of $P_{n}(x)$ can now be written as

$$\begin{aligned} \dfrac{1}{2^{n} n !}[2 n(2 n-1) \ldots(n+1)] &=\dfrac{1}{2^{n} n !}[2 n(2 n-1) \ldots(n+1)] \dfrac{n(n-1) \ldots 1}{n(n-1) \ldots 1} \\[4pt] &=\dfrac{1}{2^{n} n !} \dfrac{(2 n) !}{n !} \end{aligned} \nonumber$$

In order to prove the three term recursion formula we consider the expression $(2 n-1) x P_{n-1}(x)-n P_{n}(x)$. While each term is a polynomial of degree $n$, the leading order terms cancel. We need only look at the coefficient of the leading order term first expression. It is

$$(2 n-1) \dfrac{1}{2^{n-1}(n-1) !} \dfrac{(2 n-2) !}{(n-1) !}=\dfrac{1}{2^{n-1}(n-1) !} \dfrac{(2 n-1) !}{(n-1) !}=\dfrac{(2 n-1) !}{2^{n-1}[(n-1) !]^{2}} . \nonumber$$

The coefficient of the leading term for $n P_{n}(x)$ can be written as

$$n \dfrac{1}{2^{n} n !} \dfrac{(2 n) !}{n !}=n\left(\dfrac{2 n}{2 n^{2}}\right)\left(\dfrac{1}{2^{n-1}(n-1) !}\right) \dfrac{(2 n-1) !}{(n-1) !} \dfrac{(2 n-1) !}{2^{n-1}[(n-1) !]^{2}} . \nonumber$$

It is easy to see that the leading order terms in $(2 n-1) x P_{n-1}(x)-n P_{n}(x)$ cancel.

The next terms will be of degree $n-2$. This is because the $P_{n}$ ’s are either even or odd functions, thus only containing even, or odd, powers of $x$. We conclude that

$$(2 n-1) x P_{n-1}(x)-n P_{n}(x)=\text { polynomial of degree } n-2 . \nonumber$$

Therefore, since the Legendre polynomials form a basis, we can write this polynomial as a linear combination of of Legendre polynomials:

$$(2 n-1) x P_{n-1}(x)-n P_{n}(x)=c_{0} P_{0}(x)+c_{1} P_{1}(x)+\ldots+c_{n-2} P_{n-2}(x) . \nonumber$$

Multiplying Equation $(7.17)$ by $P_{m}(x)$ for $m=0,1, \ldots, n-3$, integrating from $-1$ to 1 , and using orthogonality, we obtain

$$0=c_{m}\left\|P_{m}\right\|^{2}, \quad m=0,1, \ldots, n-3 . \nonumber$$

[Note: $\int_{-1}^{1} x^{k} P_{n}(x) d x=0$ for $k \leq n-1$. Thus, $\int_{-1}^{1} x P_{n-1}(x) P_{m}(x) d x=0$ for $m \leq n-3 .]$

Thus, all of these $c_{m}$ ’s are zero, leaving Equation (7.17) as

$$(2 n-1) x P_{n-1}(x)-n P_{n}(x)=c_{n-2} P_{n-2}(x) . \nonumber$$

The final coefficient can be found by using the normalization condition, $P_{n}(1)=1$. Thus, $c_{n-2}=(2 n-1)-n=n-1$.

The Generating Function

where $\theta$ is the angle between $\mathbf{r}_{1}$ and $\mathbf{r}_{2}$.

Typically, one of the position vectors is much larger than the other. Let’s assume that $r_{1} \ll r_{2}$. Then, one can write

$$\Phi \propto \dfrac{1}{\sqrt{r_{1}^{2}-2 r_{1} r_{2} \cos \theta+r_{2}^{2}}}=\dfrac{1}{r_{2}} \dfrac{1}{\sqrt{1-2 \dfrac{r_{1}}{r_{2}} \cos \theta+\left(\dfrac{r_{1}}{r_{2}}\right)^{2}}} \nonumber$$

Now, define $x=\cos \theta$ and $t=\dfrac{r_{1}}{r_{2}}$. We then have the tidal potential is proportional to the generating function for the Legendre polynomials! So, we can write the tidal potential as

$$\Phi \propto \dfrac{1}{r_{2}} \sum_{n=0}^{\infty} P_{n}(\cos \theta)\left(\dfrac{r_{1}}{r_{2}}\right)^{n} . \nonumber$$

The first term in the expansion is the gravitational potential that gives the usual force between the Earth and the moon. [Recall that the force is the gradient of the potential, $\mathbf{F}=\nabla\left(\dfrac{1}{r}\right)$.] The next terms will give expressions for the tidal effects

Now that we have some idea as to where this generating function might have originated, we can proceed to use it. First of all, the generating function can be used to obtain special values of the Legendre polynomials.

Example 7.3. Evaluate $P_{n}(0) . P_{n}(0)$ is found by considering $g(0, t)$. Setting $x=0$ in Equation (7.18), we have

$$g(0, t)=\dfrac{1}{\sqrt{1+t^{2}}}=\sum_{n=0}^{\infty} P_{n}(0) t^{n} \nonumber$$

We can use the binomial expansion to find our final answer. [See the last section of this chapter for a review.] Namely, we have

$$\dfrac{1}{\sqrt{1+t^{2}}}=1-\dfrac{1}{2} t^{2}+\dfrac{3}{8} t^{4}+\ldots \nonumber$$

Comparing these expansions, we have the $P_{n}(0)=0$ for $n$ odd and for even integers one can show (see Problem $7.10$ ) that

$$P_{2 n}(0)=(-1)^{n} \dfrac{(2 n-1) ! !}{(2 n) ! !} \nonumber$$

where $n ! !$ is the double factorial,

$$n ! !=\left\{\begin{array}{c} n(n-2) \ldots(3) 1, n>0, \text { odd } \\[4pt] n(n-2) \ldots(4) 2, n>0, \text { even } \\[4pt] 1 \quad n=0,-1 \end{array}\right. \nonumber$$

Example 7.4. Evaluate $P_{n}(-1)$. This is a simpler problem. In this case we have

$$g(-1, t)=\dfrac{1}{\sqrt{1+2 t+t^{2}}}=\dfrac{1}{1+t}=1-t+t^{2}-t^{3}+\ldots \nonumber$$

Therefore, $P_{n}(-1)=(-1)^{n}$.

We can also use the generating function to find recursion relations. To prove the three term recursion (7.14) that we introduced above, then we need only differentiate the generating function with respect to $t$ in Equation (7.18) and rearrange the result. First note that

$$\dfrac{\partial g}{\partial t}=\dfrac{x-t}{\left(1-2 x t+t^{2}\right)^{3 / 2}}=\dfrac{x-t}{1-2 x t+t^{2}} g(x, t) \nonumber$$

Combining this with

$$\dfrac{\partial g}{\partial t}=\sum_{n=0}^{\infty} n P_{n}(x) t^{n-1} \nonumber$$

we have

$$(x-t) g(x, t)=\left(1-2 x t+t^{2}\right) \sum_{n=0}^{\infty} n P_{n}(x) t^{n-1} \nonumber$$

Inserting the series expression for $g(x, t)$ and distributing the sum on the right side, we obtain

$$(x-t) \sum_{n=0}^{\infty} P_{n}(x) t^{n}=\sum_{n=0}^{\infty} n P_{n}(x) t^{n-1}-\sum_{n=0}^{\infty} 2 n x P_{n}(x) t^{n}+\sum_{n=0}^{\infty} n P_{n}(x) t^{n+1} \nonumber$$

Rearranging leads to three separate sums:

$$\sum_{n=0}^{\infty} n P_{n}(x) t^{n-1}-\sum_{n=0}^{\infty}(2 n+1) x P_{n}(x) t^{n}+\sum_{n=0}^{\infty}(n+1) P_{n}(x) t^{n+1}=0 \nonumber$$

Each term contains powers of $t$ that we would like to combine into a single sum. This is done by reindexing. For the first sum, we could use the new index $k=n-1$. Then, the first sum can be written

$$\sum_{n=0}^{\infty} n P_{n}(x) t^{n-1}=\sum_{k=-1}^{\infty}(k+1) P_{k+1}(x) t^{k} \nonumber$$

Using different indices is just another way of writing out the terms. Note that

$$\sum_{n=0}^{\infty} n P_{n}(x) t^{n-1}=0+P_{1}(x)+2 P_{2}(x) t+3 P_{3}(x) t^{2}+\ldots \nonumber$$

and

$$\sum_{k=-1}^{\infty}(k+1) P_{k+1}(x) t^{k}=0+P_{1}(x)+2 P_{2}(x) t+3 P_{3}(x) t^{2}+\ldots \nonumber$$

actually give the same sum. The indices are sometimes referred to as dummy indices because they do not show up in the expanded expression and can be replaced with another letter.

If we want to do so, we could now replace all of the $k$ ’s with $n$ ’s. However, we will leave the $k$ ’s in the first term and now reindex the next sums in Equation (7.21). The second sum just needs the replacement $n=k$ and the last sum we reindex using $k=n+1$. Therefore, Equation (7.21) becomes

$$\sum_{k=-1}^{\infty}(k+1) P_{k+1}(x) t^{k}-\sum_{k=0}^{\infty}(2 k+1) x P_{k}(x) t^{k}+\sum_{k=1}^{\infty} k P_{k-1}(x) t^{k}=0 . \nonumber$$

We can now combine all of the terms, noting the $k=-1$ term is automatically zero and the $k=0$ terms give

$$P_{1}(x)-x P_{0}(x)=0 . \nonumber$$

Of course, we know this already. So, that leaves the $k>0$ terms:

$$\sum_{k=1}^{\infty}\left[(k+1) P_{k+1}(x)-(2 k+1) x P_{k}(x)+k P_{k-1}(x)\right] t^{k}=0 \nonumber$$

Since this is true for all $t$, the coefficients of the $t^{k}$ ’s are zero, or

$$(k+1) P_{k+1}(x)-(2 k+1) x P_{k}(x)+k P_{k-1}(x)=0, \quad k=1,2, \ldots \nonumber$$

There are other recursion relations. For example,

$$P_{n+1}^{\prime}(x)-P_{n-1}^{\prime}(x)=(2 n+1) P_{n}(x) . \nonumber$$

This can be proven using the generating function by differentiating $g(x, t)$ with respect to $x$ and rearranging the resulting infinite series just as in this last manipulation. This will be left as Problem 7.4.

Another use of the generating function is to obtain the normalization constant. Namely, $\left\|P_{n}\right\|^{2}$. Squaring the generating function, we have

$$\dfrac{1}{1-2 x t+t^{2}}=\left[\sum_{n=0}^{\infty} P_{n}(x) t^{n}\right]^{2}=\sum_{n=0}^{\infty} \sum_{m=0}^{\infty} P_{n}(x) P_{m}(x) t^{n+m} \nonumber$$

Integrating from -1 to 1 and using the orthogonality of the Legendre polynomials, we have

$$\begin{aligned} \int_{-1}^{1} \dfrac{d x}{1-2 x t+t^{2}} &=\sum_{n=0}^{\infty} \sum_{m=0}^{\infty} t^{n+m} \int_{-1}^{1} P_{n}(x) P_{m}(x) d x \\[4pt] &=\sum_{n=0}^{\infty} t^{2 n} \int_{-1}^{1} P_{n}^{2}(x) d x \end{aligned} \nonumber$$

However, one can show that

$$\int_{-1}^{1} \dfrac{d x}{1-2 x t+t^{2}}=\dfrac{1}{t} \ln \left(\dfrac{1+t}{1-t}\right) \nonumber$$

Expanding this expression about $t=0$, we obtain

$$\dfrac{1}{t} \ln \left(\dfrac{1+t}{1-t}\right)=\sum_{n=0}^{\infty} \dfrac{2}{2 n+1} t^{2 n} \nonumber$$

Comparing this result with Equation $(7.27)$, we find that

$$\left\|P_{n}\right\|^{2}=\int_{-1}^{1} P_{n}(x) P_{m}(x) d x=\dfrac{2}{2 n+1} . \nonumber$$

Eigenfunction Expansions

Finally, we can expand other functions in this orthogonal basis. This is just a generalized Fourier series. A Fourier-Legendre series expansion for $f(x)$ on $[-1,1]$ takes the form

$$f(x) \sim \sum_{n=0}^{\infty} c_{n} P_{n}(x) . \nonumber$$

As before, we can determine the coefficients by multiplying both sides by $P_{m}(x)$ and integrating. Orthogonality gives the usual form for the generalized Fourier coefficients. In this case, we have

$$c_{n}=\dfrac{<f, P_{n}>}{\left\|P_{n}\right\|^{2}}, \nonumber$$

where

$$<f, P_{n}>=\int_{-1}^{1} f(x) P_{n}(x) d x \nonumber$$

We have just found $\left\|P_{n}\right\|^{2}=\dfrac{2}{2 n+1}$. Therefore, the Fourier-Legendre coefficients are

$$c_{n}=\dfrac{2 n+1}{2} \int_{-1}^{1} f(x) P_{n}(x) d x . \nonumber$$

Example 7.5. Expand $f(x)=x^{3}$ in a Fourier-Legendre series.

We simply need to compute

$$c_{n}=\dfrac{2 n+1}{2} \int_{-1}^{1} x^{3} P_{n}(x) d x . \nonumber$$

We first note that

$$\int_{-1}^{1} x^{m} P_{n}(x) d x=0 \quad \text { for } m<n \nonumber$$

This is simply proven using Rodrigues formula. Inserting Equation (7.12), we have

$$\int_{-1}^{1} x^{m} P_{n}(x) d x=\dfrac{1}{2^{n} n !} \int_{-1}^{1} x^{m} \dfrac{d^{n}}{d x^{n}}\left(x^{2}-1\right)^{n} d x \nonumber$$

Since $m<n$, we can integrate by parts $m$-times to show the result, using $P_{n}(1)=1$ and $P_{n}(-1)=(-1)^{n}$. As a result, we will have for this example that $c_{n}=0$ for $n>3$.

We could just compute $\int_{-1}^{1} x^{3} P_{m}(x) d x$ for $m=0,1,2, \ldots$ outright. But, noting that $x^{3}$ is an odd function, we easily confirm that $c_{0}=0$ and $c_{2}=0$. This leaves us with only two coefficients to compute. These are

$$c_{1}=\dfrac{3}{2} \int_{-1}^{1} x^{4} d x=\dfrac{3}{5} \nonumber$$

and

$$c_{3}=\dfrac{7}{2} \int_{-1}^{1} x^{3}\left[\dfrac{1}{2}\left(5 x^{3}-3 x\right)\right] d x=\dfrac{2}{5} \nonumber$$

Thus,

$$x^{3}=\dfrac{3}{5} P_{1}(x)+\dfrac{2}{5} P_{3}(x) . \nonumber$$

Of course, this is simple to check using Table $7.2 .1$ :

$$\dfrac{3}{5} P_{1}(x)+\dfrac{2}{5} P_{3}(x)=\dfrac{3}{5} x+\dfrac{2}{5}\left[\dfrac{1}{2}\left(5 x^{3}-3 x\right)\right]=x^{3} \nonumber$$

Well, maybe we could have guessed this without doing any integration. Let’s see,

$$\begin{aligned} x^{3} &=c_{1} x+\dfrac{1}{2} c_{2}\left(5 x^{3}-3 x\right) \\[4pt] &=\left(c_{1}-\dfrac{3}{2} c_{2}\right) x+\dfrac{5}{2} c_{2} x^{3} \end{aligned} \nonumber$$

Equating coefficients of like terms, we have that $c_{2}=\dfrac{2}{5}$ and $c_{1}=\dfrac{3}{2} c_{2}=\dfrac{3}{5}$. Example 7.6. Expand the Heaviside function in a Fourier-Legendre series.

The Heaviside function is defined as

$$H(x)=\left\{\begin{array}{l} 1, x>0 \\[4pt] 0, x<0 \end{array}\right. \nonumber$$

In this case, we cannot find the expansion coefficients without some integration. We have to compute

$$\begin{aligned} c_{n} &=\dfrac{2 n+1}{2} \int_{-1}^{1} f(x) P_{n}(x) d x \\[4pt] &=\dfrac{2 n+1}{2} \int_{0}^{1} P_{n}(x) d x, \quad n=0,1,2, \ldots \end{aligned} \nonumber$$

For $n=0$, we have

$$c_{0}=\dfrac{1}{2} \int_{0}^{1} d x=\dfrac{1}{2} . \nonumber$$

For $n>1$, we make use of the identity $(7.25)$ to find

$$c_{n}=\dfrac{1}{2} \int_{0}^{1}\left[P_{n+1}^{\prime}(x)-P_{n-1}^{\prime}(x)\right] d x=\dfrac{1}{2}\left[P_{n-1}(0)-P_{n+1}(0)\right] . \nonumber$$

Thus, the Fourier-Bessel series for the Heaviside function is

$$f(x) \sim \dfrac{1}{2}+\dfrac{1}{2} \sum_{n=1}^{\infty}\left[P_{n-1}(0)-P_{n+1}(0)\right] P_{n}(x) . \nonumber$$

We need to evaluate $P_{n-1}(0)-P_{n+1}(0)$. Since $P_{n}(0)=0$ for $n$ odd, the $c_{n}$ ’s vanish for $n$ even. Letting $n=2 k-1$, we have

$$f(x) \sim \dfrac{1}{2}+\dfrac{1}{2} \sum_{k=1}^{\infty}\left[P_{2 k-2}(0)-P_{2 k}(0)\right] P_{2 k-1}(x) . \nonumber$$

We can use Equation $(7.20)$

$$P_{2 k}(0)=(-1)^{k} \dfrac{(2 k-1) ! !}{(2 k) ! !}, \nonumber$$

to compute the coefficients:

$$\begin{aligned} f(x) & \sim \dfrac{1}{2}+\dfrac{1}{2} \sum_{k=1}^{\infty}\left[P_{2 k-2}(0)-P_{2 k}(0)\right] P_{2 k-1}(x) \\[4pt] &\left.=\dfrac{1}{2}+\dfrac{1}{2} \sum_{k=1}^{\infty}\left[(-1)^{k-1} \dfrac{(2 k-3) ! !}{(2 k-2) ! !}-(-1)^{k} \dfrac{(2 k-1) ! !}{(2 k) ! !}\right] P_{2 k-1}(x)\right) \\[4pt] &=\dfrac{1}{2}-\dfrac{1}{2} \sum_{k=1}^{\infty}(-1)^{k} \dfrac{(2 k-3) ! !}{(2 k-2) ! !}\left[1+\dfrac{2 k-1}{2 k}\right] P_{2 k-1}(x) \\[4pt] &=\dfrac{1}{2}-\dfrac{1}{2} \sum_{k=1}^{\infty}(-1)^{k} \dfrac{(2 k-3) ! !}{(2 k-2) ! !} \dfrac{4 k-1}{2 k} P_{2 k-1}(x) \end{aligned} \nonumber$$

The sum of the first 21 terms are shown in Figure 7.6. We note the slow convergence to the Heaviside function. Also, we see that the Gibbs phenomenon is present due to the jump discontinuity at $x=0$.

Gamma Function

Another function that often occurs in the study of special functions is the Gamma function. We will need the Gamma function in the next section on Bessel functions.

For $x>0$ we define the Gamma function as

$$\Gamma(x)=\int_{0}^{\infty} t^{x-1} e^{-t} d t, \quad x>0 . \nonumber$$

The Gamma function is a generalization of the factorial function. In fact, we have

$$\Gamma(1)=1 \nonumber$$

and

$$\Gamma(x+1)=x \Gamma(x) . \nonumber$$

The reader can prove this identity by simply performing an integration by parts. (See Problem 7.7.) In particular, for integers $n \in Z^{+}$, we then have

$$\Gamma(n+1)=n \Gamma(n)=n(n-1) \Gamma(n-2)=n(n-1) \cdots 2 \Gamma(1)=n ! \nonumber$$

We can also define the Gamma function for negative, non-integer values of $x$. We first note that by iteration on $n \in Z^{+}$, we have

$$\Gamma(x+n)=(x+n-1) \cdots(x+1) x \Gamma(x), \quad x<0, \quad x+n>0 \nonumber$$

Solving for $\Gamma(x)$, we then find

$$\Gamma(x)=\dfrac{\Gamma(x+n)}{(x+n-1) \cdots(x+1) x}, \quad-n<x<0 \nonumber$$

Note that the Gamma function is undefined at zero and the negative integers.

Example 7.7. We now prove that

$$\Gamma\left(\dfrac{1}{2}\right)=\sqrt{\pi} . \nonumber$$

This is done by direct computation of the integral:

$$\Gamma\left(\dfrac{1}{2}\right)=\int_{0}^{\infty} t^{-\dfrac{1}{2}} e^{-t} d t \nonumber$$

Letting $t=z^{2}$, we have

$$\Gamma\left(\dfrac{1}{2}\right)=2 \int_{0}^{\infty} e^{-z^{2}} d z \nonumber$$

Due to the symmetry of the integrand, we obtain the classic integral

$$\Gamma\left(\dfrac{1}{2}\right)=\int_{-\infty}^{\infty} e^{-z^{2}} d z \nonumber$$

which can be performed using a standard trick. Consider the integral

$$I=\int_{-\infty}^{\infty} e^{-x^{2}} d x \nonumber$$

Then,

$$I^{2}=\int_{-\infty}^{\infty} e^{-x^{2}} d x \int_{-\infty}^{\infty} e^{-y^{2}} d y . \nonumber$$

Note that we changed the integration variable. This will allow us to write this product of integrals as a double integral:

$$I^{2}=\int_{-\infty}^{\infty} \int_{-\infty}^{\infty} e^{-\left(x^{2}+y^{2}\right)} d x d y . \nonumber$$

This is an integral over the entire $x y$-plane. We can transform this Cartesian integration to an integration over polar coordinates. The integral becomes

$$I^{2}=\int_{0}^{2 \pi} \int_{0}^{\infty} e^{-r^{2}} r d r d \theta \nonumber$$

This is simple to integrate and we have $I^{2}=\pi$. So, the final result is found by taking the square root of both sides:

$$\Gamma\left(\dfrac{1}{2}\right)=I=\sqrt{\pi} . \nonumber$$

We have seen that the factorial function can be written in terms of Gamma functions. One can write the even and odd double factorials as

$$(2 n) ! !=2^{n} n !, \quad(2 n+1) ! !=\dfrac{(2 n+1) !}{2^{n} n !} \nonumber$$

In particular, one can write

$$\Gamma\left(n+\dfrac{1}{2}\right)=\dfrac{(2 n-1) ! !}{2^{n}} \sqrt{\pi} \nonumber$$

Another useful relation, which we only state, is

$$\Gamma(x) \Gamma(1-x)=\dfrac{\pi}{\sin \pi x} \nonumber$$

Hypergeometric Functions

Hypergeometric functions are probably the most useful, but least understood, class of functions. They typically do not make it into the undergraduate curriculum and seldom in graduate curriculum. Most functions that you know can be expressed using hypergeometric functions. There are many approaches to these functions and the literature can fill books. 1

In 1812 Gauss published a study of the hypergeometric series

$$\begin{aligned} y(x)=1 &+\dfrac{\alpha \beta}{\gamma} x+\dfrac{\alpha(1+\alpha)(1+\beta)}{2 ! \gamma(1+\gamma)} x^{2} \\[4pt] &+\dfrac{\alpha(1+\alpha)(2+\alpha) \beta(1+\beta)(2+\beta)}{3 ! \gamma(1+\gamma)(2+\gamma)} x^{3}+\ldots \end{aligned} \nonumber$$

Here $\alpha, \beta, \gamma$, and $x$ are real numbers. If one sets $\alpha=1$ and $\beta=\gamma$, this series reduces to the familiar geometric series

$$y(x)=1+x+x^{2}+x^{3}+\ldots . \nonumber$$

The hypergeometric series is actually a solution of the differential equation

$$x(1-x) y^{\prime \prime}+[\gamma-(\alpha+\beta+1) x] y^{\prime}-\alpha \beta y=0 \nonumber$$

This equation was first introduced by Euler and latter studied extensively by Gauss, Kummer and Riemann. It is sometimes called Gauss’ equation. Note that there is a symmetry in that $\alpha$ and $\beta$ may be interchanged without changing the equation. The points $x=0$ and $x=1$ are regular singular points. Series solutions may be sought using the Frobenius method. It can be confirmed that the above hypergeometric series results.

A more compact form for the hypergeometric series may be obtained by introducing new notation. One typically introduces the Pochhammer symbol, $(\alpha)_{n}$, satisfying (i) $(\alpha)_{0}=1$ if $\alpha \neq 0$. and (ii) $(\alpha)_{k}=\alpha(1+\alpha) \ldots(k-1+\alpha)$, for $k=1,2, \ldots$

Consider $(1)_{n}$. For $n=0,(1)_{0}=1$. For $n>0$,

$$(1)_{n}=1(1+1)(2+1) \ldots[(n-1)+1] \nonumber$$

This reduces to $(1)_{n}=n !$. In fact, one can show that

$$(k)_{n}=\dfrac{(n+k-1) !}{(k-1) !} \nonumber$$

for $k$ and $n$ positive integers. In fact, one can extend this result to noninteger values for $k$ by introducing the gamma function:

$$(\alpha)_{n}=\dfrac{\Gamma(\alpha+n)}{\Gamma(\alpha)} \nonumber$$

We can now write the hypergeometric series in standard notation as

${ }^{1}$ See for example Special Functions by G. E. Andrews, R. Askey, and R. Roy, 1999, Cambridge University Press.

$${ }_{2} F_{1}(\alpha, \beta ; \gamma ; x)=\sum_{n=0}^{\infty} \dfrac{(\alpha)_{n}(\beta)_{n}}{n !(\gamma)_{n}} x^{n} \nonumber$$

Using this one can show that the general solution of Gauss’ equation is

$$y(x)=A_{2} F_{1}(\alpha, \beta ; \gamma ; x)+B_{2} x_{2}^{1-\gamma} F_{1}(1-\gamma+\alpha, 1-\gamma+\beta ; 2-\gamma ; x) . \nonumber$$

By carefully letting $\beta$ approach $\infty$, one obtains what is called the confluent hypergeometric function. This in effect changes the nature of the differential equation. Gauss’ equation has three regular singular points at $x=0,1, \infty$. One can transform Gauss’ equation by letting $x=u / \beta$. This changes the regular singular points to $u=0, \beta, \infty$. Letting $\beta \rightarrow \infty$, two of the singular points merge.

The new confluent hypergeometric function is then given as

$${ }_{1} F_{1}(\alpha ; \gamma ; u)=\lim _{\beta \rightarrow \infty}{ }_{2} F_{1}\left(\alpha, \beta ; \gamma ; \dfrac{u}{\beta}\right) . \nonumber$$

This function satisfies the differential equation

$$x y^{\prime \prime}+(\gamma-x) y^{\prime}-\alpha y=0 . \nonumber$$

The purpose of this section is only to introduce the hypergeometric function. Many other special functions are related to the hypergeometric function after making some variable transformations. For example, the Legendre polynomials are given by

$$P_{n}(x)={ }_{2} F_{1}\left(-n, n+1 ; 1 ; \dfrac{1-x}{2}\right) . \nonumber$$

In fact, one can also show that

$$\sin ^{-1} x=x_{2} F_{1}\left(\dfrac{1}{2}, \dfrac{1}{2} ; \dfrac{3}{2} ; x^{2}\right) . \nonumber$$

The Bessel function $J_{p}(x)$ can be written in terms of confluent geometric functions as

$$J_{p}(x)=\dfrac{1}{\Gamma(p+1)}\left(\dfrac{z}{2}\right)^{p} e^{-i z}{ }_{1} F_{1}\left(\dfrac{1}{2}+p, 1+2 p ; 2 i z\right) . \nonumber$$

These are just a few connections of the powerful hypergeometric functions to some of the elementary functions that you know.

Appendix: The Binomial Expansion

In this section we had to recall the binomial expansion. This is simply the expansion of the expression $(a+b)^{p}$. We will investigate this expansion first for nonnegative integer powers $p$ and then derive the expansion for other values of $p$.

Lets list some of the common expansions for nonnegative integer powers.

$$\begin{aligned} &(a+b)^{0}=1 \\[4pt] &(a+b)^{1}=a+b \\[4pt] &(a+b)^{2}=a^{2}+2 a b+b^{2} \\[4pt] &(a+b)^{3}=a^{3}+3 a^{2} b+3 a b^{2}+b^{3} \\[4pt] &(a+b)^{4}=a^{4}+4 a^{3} b+6 a^{2} b^{2}+4 a b^{3}+b^{4} \end{aligned} \nonumber$$

We now look at the patterns of the terms in the expansions. First, we note that each term consists of a product of a power of $a$ and a power of $b$. The powers of $a$ are decreasing from $n$ to 0 in the expansion of $(a+b)^{n}$. Similarly, the powers of $b$ increase from 0 to $n$. The sums of the exponents in each term is $n$. So, we can write the $(k+1)$ st term in the expansion as $a^{n-k} b^{k}$. For example, in the expansion of $(a+b)^{51}$ the 6 th term is $a^{51-5} b^{5}=a^{46} b^{5}$. However, we do not know the numerical coefficient in the expansion.

We now list the coefficients for the above expansions.

This pattern is the famous Pascal’s triangle. There are many interesting features of this triangle. But we will first ask how each row can be generated.

We see that each row begins and ends with a one. Next the second term and next to last term has a coefficient of $n$. Next we note that consecutive pairs in each row can be added to obtain entries in the next row. For example, we have

With this in mind, we can generate the next several rows of our triangle.

Of course, it would take a while to compute each row up to the desired $n$. We need a simple expression for computing a specific coefficient. Consider

$$\begin{aligned} & n=0: \quad 1 \\[4pt] & n=1: \quad 1 \quad 1 \\[4pt] & n=2: \quad 1 \quad 2 \quad 1 \\[4pt] & n=3: \quad 1 \quad 3 \quad 3 \quad 1 \\[4pt] & n=3: \quad 1 \quad 3 \quad 3 \quad 1 \\[4pt] & n=4: \quad 1 \quad 4 \quad 6 \quad 4 \quad 1 \end{aligned} \nonumber$$

the $k$ th term in the expansion of $(a+b)^{n}$. Let $r=k-1$. Then this term is of the form $C_{r}^{n} a^{n-r} b^{r}$. We have seen the the coefficients satisfy

$$C_{r}^{n}=C_{r}^{n-1}+C_{r-1}^{n-1} \nonumber$$

Actually, the coefficients have been found to take a simple form.

$$C_{r}^{n}=\dfrac{n !}{(n-r) ! r !}=\left(\begin{array}{c} n \\[4pt] r \end{array}\right) \nonumber$$

This is nothing other than the combinatoric symbol for determining how to choose $n$ things $r$ at a time. In our case, this makes sense. We have to count the number of ways that we can arrange the products of $r$ b’s with $n-r a$ ’s. There are $n$ slots to place the $b$ ’s. For example, the $r=2$ case for $n=4$ involves the six products: $a a b b, a b a b, a b b a, b a a b, b a b a$, and bbaa. Thus, it is natural to use this notation. The original problem that concerned Pascal was in gambling.

So, we have found that

$$(a+b)^{n}=\sum_{r=0}^{n}\left(\begin{array}{c} n \\[4pt] r \end{array}\right) a^{n-r} b^{r} \nonumber$$

What if $a \gg b$ ? Can we use this to get an approximation to $(a+b)^{n}$ ? If we neglect $b$ then $(a+b)^{n} \simeq a^{n}$. How good of an approximation is this? This is where it would be nice to know the order of the next term in the expansion, which we could state using big $O$ notation. In order to do this we first divide out $a$ as

$$(a+b)^{n}=a^{n}\left(1+\dfrac{b}{a}\right)^{n} . \nonumber$$

Now we have a small parameter, $\dfrac{b}{a}$. According to what we have seen above, we can use the binomial expansion to write

$$\left(1+\dfrac{b}{a}\right)^{n}=\sum_{r=0}^{n}\left(\begin{array}{l} n \\[4pt] r \end{array}\right)\left(\dfrac{b}{a}\right)^{r} \nonumber$$

Thus, we have a finite sum of terms involving powers of $\dfrac{b}{a}$. Since $a \gg b$, most of these terms can be neglected. So, we can write

$$\left(1+\dfrac{b}{a}\right)^{n}=1+n \dfrac{b}{a}+O\left(\left(\dfrac{b}{a}\right)^{2}\right) \nonumber$$

note that we have used the observation that the second coefficient in the $n$th row of Pascal’s triangle is $n$.

Summarizing, this then gives

$$\begin{aligned} (a+b)^{n} &=a^{n}\left(1+\dfrac{b}{a}\right)^{n} \\[4pt] &=a^{n}\left(1+n \dfrac{b}{a}+O\left(\left(\dfrac{b}{a}\right)^{2}\right)\right) \\[4pt] &=a^{n}+n a^{n} \dfrac{b}{a}+a^{n} O\left(\left(\dfrac{b}{a}\right)^{2}\right) \end{aligned} \nonumber$$

Therefore, we can approximate $(a+b)^{n} \simeq a^{n}+n b a^{n-1}$, with an error on the order of $b a^{n-2}$. Note that the order of the error does not include the constant factor from the expansion. We could also use the approximation that $(a+b)^{n} \simeq a^{n}$, but it is not as good because the error in this case is of the order $b a^{n-1}$.

We have seen that

$$\dfrac{1}{1-x}=1+x+x^{2}+\ldots \nonumber$$

But, $\dfrac{1}{1-x}=(1-x)^{-1}$. This is again a binomial to a power, but the power is not a nonnegative integer. It turns out that the coefficients of such a binomial expansion can be written similar to the form in Equation (7.60).

This example suggests that our sum may no longer be finite. So, for $p$ a real number, we write

$$(1+x)^{p}=\sum_{r=0}^{\infty}\left(\begin{array}{l} p \\[4pt] r \end{array}\right) x^{r} . \nonumber$$

However, we quickly run into problems with this form. Consider the coefficient for $r=1$ in an expansion of $(1+x)^{-1}$. This is given by

$$\left(\begin{array}{c} -1 \\[4pt] 1 \end{array}\right)=\dfrac{(-1) !}{(-1-1) ! 1 !}=\dfrac{(-1) !}{(-2) ! 1 !} \text {. } \nonumber$$

But what is $(-1) ! ?$ By definition, it is

$$(-1) !=(-1)(-2)(-3) \cdots . \nonumber$$

This product does not seem to exist! But with a little care, we note that

$$\dfrac{(-1) !}{(-2) !}=\dfrac{(-1)(-2) !}{(-2) !}=-1 \text {. } \nonumber$$

So, we need to be careful not to interpret the combinatorial coefficient literally. There are better ways to write the general binomial expansion. We can write the general coefficient as

$$\begin{aligned} \left(\begin{array}{l} p \\[4pt] r \end{array}\right) &=\dfrac{p !}{(p-r) ! r !} \\[4pt] &=\dfrac{p(p-1) \cdots(p-r+1)(p-r) !}{(p-r) ! r !} \\[4pt] &=\dfrac{p(p-1) \cdots(p-r+1)}{r !} . \end{aligned} \nonumber$$

With this in mind we now state the theorem:

General Binomial Expansion The general binomial expansion for $(1+$ $x)^{p}$ is a simple generalization of Equation (7.60). For $p$ real, we have that

$$\begin{aligned} (1+x)^{p} &=\sum_{r=0}^{\infty} \dfrac{p(p-1) \cdots(p-r+1)}{r !} x^{r} \\[4pt] &=\sum_{r=0}^{\infty} \dfrac{\Gamma(p+1)}{r ! \Gamma(p-r+1)} x^{r} \end{aligned} \nonumber$$

Often we need the first few terms for the case that $x \ll 1$ :

$$(1+x)^{p}=1+p x+\dfrac{p(p-1)}{2} x^{2}+O\left(x^{3}\right) . \nonumber$$

The Method of Variation of Parameters

We are interested in solving nonhomogeneous second order linear differential equations of the form

$$a_{2}(x) y^{\prime \prime}(x)+a_{1}(x) y^{\prime}(x)+a_{0}(x) y(x)=f(x) . \nonumber$$

The general solution of this nonhomogeneous second order linear differential equation is found as a sum of the general solution of the homogeneous equation,

$$a_{2}(x) y^{\prime \prime}(x)+a_{1}(x) y^{\prime}(x)+a_{0}(x) y(x)=0, \nonumber$$

and a particular solution of the nonhomogeneous equation. Recall from Chapter 1 that there are several approaches to finding particular solutions of nonhomogeneous equations. Any guess would be sufficient. An intelligent guess, based upon the Method of Undetermined Coefficients, was reviewed previously in Chapter 1. However, a more methodical method, which is first seen in a first course in differential equations, is the Method of Variation of Parameters. Also, we explored the matrix version of this method in Section $2.8$. We will review this method in this section and extend it to the solution of boundary value problems.

While it is sufficient to derive the method for the general differential equation above, we will instead consider solving equations that are in SturmLiouville, or self-adjoint, form. Therefore, we will apply the Method of Variation of Parameters to the equation

$$\dfrac{d}{d x}\left(p(x) \dfrac{d y(x)}{d x}\right)+q(x) y(x)=f(x) \nonumber$$

Note that $f(x)$ in this equation is not the same function as in the general equation posed at the beginning of this section.

We begin by assuming that we have determined two linearly independent solutions of the homogeneous equation. The general solution is then given by

$$y(x)=c_{1} y_{1}(x)+c_{2} y_{2}(x) . \nonumber$$

In order to determine a particular solution of the nonhomogeneous equation, we vary the parameters $c_{1}$ and $c_{2}$ in the solution of the homogeneous problem by making them functions of the independent variable. Thus, we seek a particular solution of the nonhomogeneous equation in the form

$$y_{p}(x)=c_{1}(x) y_{1}(x)+c_{2}(x) y_{2}(x) . \nonumber$$

In order for this to be a solution, we need to show that it satisfies the differential equation. We first compute the derivatives of $y_{p}(x)$. The first derivative

$$y_{p}^{\prime}(x)=c_{1}(x) y_{1}^{\prime}(x)+c_{2}(x) y_{2}^{\prime}(x)+c_{1}^{\prime}(x) y_{1}(x)+c_{2}^{\prime}(x) y_{2}(x) . \nonumber$$

Without loss of generality, we will set the sum of the last two terms to zero. (One can show that the same results would be obtained if we did not. See Problem 8.2.) Then, we have

$$c_{1}^{\prime}(x) y_{1}(x)+c_{2}^{\prime}(x) y_{2}(x)=0 . \nonumber$$

Now, we take the second derivative of the remaining terms to obtain

$$y_{p}^{\prime \prime}(x)=c_{1}(x) y_{1}^{\prime \prime}(x)+c_{2}(x) y_{2}^{\prime \prime}(x)+c_{1}^{\prime}(x) y_{1}^{\prime}(x)+c_{2}^{\prime}(x) y_{2}^{\prime}(x) \nonumber$$

Expanding the derivative term in Equation (8.3),

$$p(x) y_{p}^{\prime \prime}(x)+p^{\prime}(x) y_{p}^{\prime}(x)+q(x) y_{p}(x)=f(x), \nonumber$$

and inserting the expressions for $y_{p}, y_{p}^{\prime}(x)$, and $y_{p}^{\prime \prime}(x)$, we have

$$\begin{aligned} f(x)=& p(x)\left[c_{1}(x) y_{1}^{\prime \prime}(x)+c_{2}(x) y_{2}^{\prime \prime}(x)+c_{1}^{\prime}(x) y_{1}^{\prime}(x)+c_{2}^{\prime}(x) y_{2}^{\prime}(x)\right] \\[4pt] &+p^{\prime}(x)\left[c_{1}(x) y_{1}^{\prime}(x)+c_{2}(x) y_{2}^{\prime}(x)\right]+q(x)\left[c_{1}(x) y_{1}(x)+c_{2}(x) y_{2}(x)\right] . \end{aligned} \nonumber$$

Rearranging terms, we find

$$\begin{aligned} f(x)=& c_{1}(x)\left[p(x) y_{1}^{\prime \prime}(x)+p^{\prime}(x) y_{1}^{\prime}(x)+q(x) y_{1}(x)\right] \\[4pt] &+c_{2}(x)\left[p(x) y_{2}^{\prime \prime}(x)+p^{\prime}(x) y_{2}^{\prime}(x)+q(x) y_{2}(x)\right] \\[4pt] &+p(x)\left[c_{1}^{\prime}(x) y_{1}^{\prime}(x)+c_{2}^{\prime}(x) y_{2}^{\prime}(x)\right] . \end{aligned} \nonumber$$

Since $y_{1}(x)$ and $y_{2}(x)$ are both solutions of the homogeneous equation. The first two bracketed expressions vanish. Dividing by $p(x)$, we have that

$$c_{1}^{\prime}(x) y_{1}^{\prime}(x)+c_{2}^{\prime}(x) y_{2}^{\prime}(x)=\dfrac{f(x)}{p(x)} . \nonumber$$

Our goal is to determine $c_{1}(x)$ and $c_{2}(x)$. In this analysis, we have found that the derivatives of these functions satisfy a linear system of equations (in the $c_{i}$ ’s):

This system is easily solved to give

$$\begin{aligned} c_{1}^{\prime}(x) &=-\dfrac{f(x) y_{2}(x)}{p(x)\left[y_{1}(x) y_{2}^{\prime}(x)-y_{1}^{\prime}(x) y_{2}(x)\right]} \\[4pt] c_{2}^{\prime}(x) &=\dfrac{f(x) y_{1}(x)}{p(x)\left[y_{1}(x) y_{2}^{\prime}(x)-y_{1}^{\prime}(x) y_{2}(x)\right]} \end{aligned} \nonumber$$

We note that the denominator in these expressions involves the Wronskian of the solutions to the homogeneous problem. Recall that

$$W\left(y_{1}, y_{2}\right)(x)=\left|\begin{array}{ll} y_{1}(x) & y_{2}(x) \\[4pt] y_{1}^{\prime}(x) & y_{2}^{\prime}(x) \end{array}\right| \nonumber$$

Furthermore, we can show that the denominator, $p(x) W(x)$, is constant. Differentiating this expression and using the homogeneous form of the differential equation proves this assertion.

$$\begin{aligned} \dfrac{d}{d x}(p(x) W(x))=& \dfrac{d}{d x}\left[p(x)\left(y_{1}(x) y_{2}^{\prime}(x)-y_{1}^{\prime}(x) y_{2}(x)\right)\right] \\[4pt] =&\left.y_{1}(x) \dfrac{d}{d x}\left(p(x) y_{2}^{\prime}(x)\right)\right)+p(x) y_{2}^{\prime}(x) y_{1}^{\prime}(x) \\[4pt] &\left.-y_{2}(x) \dfrac{d}{d x}\left(p(x) y_{1}^{\prime}(x)\right)\right)-p(x) y_{1}^{\prime}(x) y_{2}^{\prime}(x) \\[4pt] =&-y_{1}(x) q(x) y_{2}(x)+y_{2}(x) q(x) y_{1}(x)=0 \end{aligned} \nonumber$$

Therefore,

$$p(x) W(x)=\text { constant. } \nonumber$$

So, after an integration, we find the parameters as

$$\begin{aligned} &c_{1}(x)=-\int_{x_{0}}^{x} \dfrac{f(\xi) y_{2}(\xi)}{p(\xi) W(\xi)} d \xi \\[4pt] &c_{2}(x)=\int_{x_{1}}^{x} \dfrac{f(\xi) y_{1}(\xi)}{p(\xi) W(\xi)} d \xi \end{aligned} \nonumber$$

where $x_{0}$ and $x_{1}$ are arbitrary constants to be determined later.

Therefore, the particular solution of (8.3) can be written as

$$y_{p}(x)=y_{2}(x) \int_{x_{1}}^{x} \dfrac{f(\xi) y_{1}(\xi)}{p(\xi) W(\xi)} d \xi-y_{1}(x) \int_{x_{0}}^{x} \dfrac{f(\xi) y_{2}(\xi)}{p(\xi) W(\xi)} d \xi \nonumber$$

As a further note, we usually do not rewrite our initial value problems in self-adjoint form. Recall that for an equation of the form

$$a_{2}(x) y^{\prime \prime}(x)+a_{1}(x) y^{\prime}(x)+a_{0}(x) y(x)=g(x) . \nonumber$$

we obtained the self-adjoint form by multiplying the equation by

$$\dfrac{1}{a_{2}(x)} e^{\int \dfrac{a_{1}(x)}{a_{2}(x)} d x}=\dfrac{1}{a_{2}(x)} p(x) . \nonumber$$

This gives the standard form

$$\left(p(x) y^{\prime}(x)\right)^{\prime}+q(x) y(x)=f(x) \nonumber$$

where

$$f(x)=\dfrac{1}{a_{2}(x)} p(x) g(x) . \nonumber$$

With this in mind, Equation (8.13) becomes

$$y_{p}(x)=y_{2}(x) \int_{x_{1}}^{x} \dfrac{g(\xi) y_{1}(\xi)}{a_{2}(\xi) W(\xi)} d \xi-y_{1}(x) \int_{x_{0}}^{x} \dfrac{g(\xi) y_{2}(\xi)}{\left.a_{(} \xi\right) W(\xi)} d \xi . \nonumber$$

Example 8.1. Consider the nonhomogeneous differential equation

$$y^{\prime \prime}-y^{\prime}-6 y=20 e^{-2 x} . \nonumber$$

We seek a particular solution to this equation. First, we note two linearly independent solutions of this equation are

$$y_{1}(x)=e^{3 x}, \quad y_{2}(x)=e^{-2 x} . \nonumber$$

So, the particular solution takes the form

$$y_{p}(x)=c_{1}(x) e^{3 x}+c_{2}(x) e^{-2 x} \nonumber$$

We just need to determine the $c_{i}$ ’s. Since this problem is not in self-adjoint form, we will use

$$\dfrac{f(x)}{p(x)}=\dfrac{g(x)}{a_{2}(x)}=20 e^{-2 x} \nonumber$$

as seen above. Then the linear system we have to solve is

$$\begin{aligned} c_{1}^{\prime}(x) e^{3 x}+c_{2}^{\prime}(x) e^{-2 x} &=0 \\[4pt] 3 c_{1}^{\prime}(x) e^{3 x}-2 c_{2}^{\prime}(x) e^{-2 x} &=20 e^{-2 x} \end{aligned} \nonumber$$

Multiplying the first equation by 2 and adding the equations yields

$$5 c_{1}^{\prime}(x) e^{3 x}=20 e^{-2 x} \nonumber$$

or

$$c_{1}^{\prime}(x)=4 e^{-5 x} \nonumber$$

Inserting this back into the first equation in the system, we have

$$4 e^{-2 x}+c_{2}^{\prime}(x) e^{-2 x}=0, \nonumber$$

leading to

$$c_{2}^{\prime}(x)=-4 . \nonumber$$

These equations are easily integrated to give

$$c_{1}(x)=-\dfrac{4}{5} e^{-5 x}, \quad c_{2}(x)=-4 x . \nonumber$$

Therefore, the particular solution has been found as

$$\begin{aligned} y_{p}(x) &=c_{1}(x) e^{3 x}+c_{2}(x) e^{-2 x} \\[4pt] &=-\dfrac{4}{5} e^{-5 x} e^{3 x}-4 x e^{-2 x} \\[4pt] &=-\dfrac{4}{5} e^{-2 x}-4 x e^{-2 x} \end{aligned} \nonumber$$

Noting that the first term can be absorbed into the solution of the homogeneous problem. So, the particular solution can simply be written as

$$y_{p}(x)=-4 x e^{-2 x} \nonumber$$

This is the answer you would have found had you used the Modified Method of Undetermined Coefficients.

Example 8.2. Revisiting the last example, $y^{\prime \prime}-y^{\prime}-6 y=20 e^{-2 x}$.

The formal solution in Equation (8.13) was not used in the last example. Instead, we proceeded from the Linear System for Variation of Parameters earlier in this section. This is the more natural approach towards finding the particular solution of the nonhomogeneous equation. Since we will be using Equation (8.13) to obtain solutions to initial value and boundary value problems, it might be useful to use it to solve this problem.

From the last example we have

$$y_{1}(x)=e^{3 x}, \quad y_{2}(x)=e^{-2 x} \nonumber$$

We need to compute the Wronskian:

$$W(x)=W\left(y_{1}, y_{2}\right)(x)=\left|\begin{array}{cc} e^{3 x} & e^{-2 x} \\[4pt] 3 e^{3 x} & -2 e^{-2 x} \end{array}\right|=-5 e^{x} \nonumber$$

Also, we need $p(x)$, which is given by

$$p(x)=\exp \left(-\int d x\right)=e^{-x} . \nonumber$$

So, we see that $p(x) W(x)=-5$. It is indeed constant, just as we had proven earlier.

Finally, we need $f(x)$. Here is where one needs to be careful as the original problem was not in self-adjoint form. We have from the original equation that $g(x)=20 e^{-2 x}$ and $a_{2}(x)=1$. So,

$$f(x)=\dfrac{p(x)}{a_{2}(x)} g(x)=20 e^{-3 x} \nonumber$$

Now we are ready to construct the solution.

$$\begin{aligned} y_{p}(x) &=y_{2}(x) \int_{x_{1}}^{x} \dfrac{f(\xi) y_{1}(\xi)}{p(\xi) W(\xi)} d \xi-y_{1}(x) \int_{x_{0}}^{x} \dfrac{f(\xi) y_{2}(\xi)}{p(\xi) W(\xi)} d \xi \\[4pt] &=e^{-2 x} \int_{x_{1}}^{x} \dfrac{20 e^{-3 \xi} e^{3 \xi}}{-5} d \xi-e^{3 x} \int_{x_{0}}^{x} \dfrac{20 e^{-3 \xi} e^{-2 \xi}}{-5} d \xi \\[4pt] &=-4 e^{-2 x} \int_{x_{1}}^{x} d \xi+4 e^{3 x} \int_{x_{0}}^{x} e^{-5 x} d \xi \\[4pt] &=-\left.4 \xi e^{-2 x}\right|_{x_{1}} ^{x}-\left.\dfrac{4}{5} e^{3 x} e^{-5 \xi}\right|_{x_{0}} ^{x} \\[4pt] &=-4 x e^{-2 x}-\dfrac{4}{5} e^{-2 x}+4 x_{1} e^{-2 x}+\dfrac{4}{5} e^{-5 x_{0}} e^{3 x} \end{aligned} \nonumber$$

Note that the first two terms we had found in the last example. The remaining two terms are simply linear combinations of $y_{1}$ and $y_{2}$. Thus, we really have the solution to the homogeneous problem contained within the solution when we use the arbitrary constant limits in the integrals. In the next section we will make use of these constants when solving initial value and boundary value problems.

In the next section we will determine the unknown constants subject to either initial conditions or boundary conditions. This will allow us to combine the two integrals and then determine the appropriate Green’s functions.

Initial and Boundary Value Green’s Functions

We begin with the particular solution (8.13) of our nonhomogeneous differential equation (8.3). This can be combined with the general solution of the homogeneous problem to give the general solution of the nonhomogeneous differential equation:

$$y(x)=c_{1} y_{1}(x)+c_{2} y_{2}(x)+y_{2}(x) \int_{x_{1}}^{x} \dfrac{f(\xi) y_{1}(\xi)}{p(\xi) W(\xi)} d \xi-y_{1}(x) \int_{x_{0}}^{x} \dfrac{f(\xi) y_{2}(\xi)}{p(\xi) W(\xi)} d \xi \nonumber$$

As seen in the last section, an appropriate choice of $x_{0}$ and $x_{1}$ could be found so that we need not explicitly write out the solution to the homogeneous problem, $c_{1} y_{1}(x)+c_{2} y_{2}(x)$. However, setting up the solution in this form will allow us to use $x_{0}$ and $x_{1}$ to determine particular solutions which satisfies certain homogeneous conditions.

We will now consider initial value and boundary value problems. Each type of problem will lead to a solution of the form

$$y(x)=c_{1} y_{1}(x)+c_{2} y_{2}(x)+\int_{a}^{b} G(x, \xi) f(\xi) d \xi \nonumber$$

where the function $G(x, \xi)$ will be identified as the Green’s function and the integration limits will be found on the integral. Having identified the Green’s function, we will look at other methods in the last section for determining the Green’s function.

Boundary Value Green’s Function

We now turn to boundary value problems. We will focus on the problem

$$\begin{array}{r} \dfrac{d}{d x}\left(p(x) \dfrac{d y(x)}{d x}\right)+q(x) y(x)=f(x), \quad a<x<b \\[4pt] y(a)=0, \quad y(b)=0 \end{array} \nonumber$$

However, the general theory works for other forms of homogeneous boundary conditions.

Once again, we seek $x_{0}$ and $x_{1}$ in the form

$$y(x)=y_{2}(x) \int_{x_{1}}^{x} \dfrac{f(\xi) y_{1}(\xi)}{p(\xi) W(\xi)} d \xi-y_{1}(x) \int_{x_{0}}^{x} \dfrac{f(\xi) y_{2}(\xi)}{p(\xi) W(\xi)} d \xi \nonumber$$

so that the solution to the boundary value problem can be written as a single integral involving a Green’s function. Here we absorb $y_{h}(x)$ into the integrals with an appropriate choice of lower limits on the integrals.

We first pick solutions of the homogeneous differential equation such that $y_{1}(a)=0, y_{2}(b)=0$ and $y_{1}(b) \neq 0, y_{2}(a) \neq 0$. So, we have

$$\begin{aligned} y(a) &=y_{2}(a) \int_{x_{1}}^{a} \dfrac{f(\xi) y_{1}(\xi)}{p(\xi) W(\xi)} d \xi-y_{1}(a) \int_{x_{0}}^{a} \dfrac{f(\xi) y_{2}(\xi)}{p(\xi) W(\xi)} d \xi \\[4pt] &=y_{2}(a) \int_{x_{1}}^{a} \dfrac{f(\xi) y_{1}(\xi)}{p(\xi) W(\xi)} d \xi . \end{aligned} \nonumber$$

This expression is zero if $x_{1}=a$.

At $x=b$ we find that

$$\begin{aligned} y(b) &=y_{2}(b) \int_{x_{1}}^{b} \dfrac{f(\xi) y_{1}(\xi)}{p(\xi) W(\xi)} d \xi-y_{1}(b) \int_{x_{0}}^{b} \dfrac{f(\xi) y_{2}(\xi)}{p(\xi) W(\xi)} d \xi \\[4pt] &=-y_{1}(b) \int_{x_{0}}^{b} \dfrac{f(\xi) y_{2}(\xi)}{p(\xi) W(\xi)} d \xi \end{aligned} \nonumber$$

This vanishes for $x_{0}=b$.

So, we have found that

$$y(x)=y_{2}(x) \int_{a}^{x} \dfrac{f(\xi) y_{1}(\xi)}{p(\xi) W(\xi)} d \xi-y_{1}(x) \int_{b}^{x} \dfrac{f(\xi) y_{2}(\xi)}{p(\xi) W(\xi)} d \xi . \nonumber$$

We are seeking a Green’s function so that the solution can be written as one integral. We can move the functions of $x$ under the integral. Also, since $a<x<b$, we can flip the limits in the second integral. This gives

$$y(x)=\int_{a}^{x} \dfrac{f(\xi) y_{1}(\xi) y_{2}(x)}{p(\xi) W(\xi)} d \xi+\int_{x}^{b} \dfrac{f(\xi) y_{1}(x) y_{2}(\xi)}{p(\xi) W(\xi)} d \xi \nonumber$$

This result can be written in a compact form:

Properties of Green’s Functions

We have noted some properties of Green’s functions in the last section. In this section we will elaborate on some of these properties as a tool for quickly constructing Green’s functions for boundary value problems. Here is a list of the properties based upon our previous solution.

The Dirac Delta Function

We will develop a more general theory of Green’s functions for ordinary differential equations which encompasses some of the listed properties. The Green’s function satisfies a homogeneous differential equation for $x \neq \xi$,

$$\dfrac{\partial}{\partial x}\left(p(x) \dfrac{\partial G(x, \xi)}{\partial x}\right)+q(x) G(x, \xi)=0, \quad x \neq \xi \nonumber$$

When $x=\xi$, we saw that the derivative has a jump in its value. This is similar to the step, or Heaviside, function,

$$H(x)=\left\{\begin{array}{l} 1, x>0 \\[4pt] 0, x<0 \end{array}\right. \nonumber$$

In the case of the step function, the derivative is zero everywhere except at the jump. At the jump, there is an infinite slope, though technically, we have learned that there is no derivative at this point. We will try to remedy this by introducing the Dirac delta function,

$$\delta(x)=\dfrac{d}{d x} H(x) . \nonumber$$

We will then show that the Green’s function satisfies the differential equation

$$\dfrac{\partial}{\partial x}\left(p(x) \dfrac{\partial G(x, \xi)}{\partial x}\right)+q(x) G(x, \xi)=\delta(x-\xi) \nonumber$$

The Dirac delta function, $\delta(x)$, is one example of what is known as a generalized function, or a distribution. Dirac had introduced this function in the 1930’s in his study of quantum mechanics as a useful tool. It was later studied in a general theory of distributions and found to be more than a simple tool used by physicists. The Dirac delta function, as any distribution, only makes sense under an integral.

Before defining the Dirac delta function and introducing some of its properties, we will look at some representations that lead to the definition. We will consider the limits of two sequences of functions.

First we define the sequence of functions

$$f_{n}(x)=\left\{\begin{array}{l} 0,|x|>\dfrac{1}{n} \\[4pt] \dfrac{n}{2},|x|<\dfrac{1}{n} \end{array} .\right. \nonumber$$

This is a sequence of functions as shown in Figure 8.1. As $n \rightarrow \infty$, we find the limit is zero for $x \neq 0$ and is infinite for $x=0$. However, the area under each member of the sequences is one since each box has height $\dfrac{n}{2}$ and width $\dfrac{2}{n}$. Thus, the limiting function is zero at most points but has area one. (At this point the reader who is new to this should be doing some head scratching!)

The limit is not really a function. It is a generalized function. It is called the Dirac delta function, which is defined by

1. $\delta(x)=0$ for $x \neq 0$
2. $\int_{-\infty}^{\infty} \delta(x) d x=1$

Another example is the sequence defined by

$$D_{n}(x)=\dfrac{2 \sin n x}{x} \nonumber$$

We can graph this function. We first rewrite this function as

$$D_{n}(x)=2 n \dfrac{\sin n x}{n x} . \nonumber$$

Now it is easy to see that as $x \rightarrow 0, D_{n}(x) \rightarrow 2 n$. For large $x$, The function tends to zero. A plot of this function is in Figure 8.2. For large $n$ the peak grows and the values of $D_{n}(x)$ for $x \neq 0$ tend to zero as show in Figure 8.3.

We note that in the limit $n \rightarrow \infty, D_{n}(x)=0$ for $x \neq 0$ and it is infinite at $x=0$. However, using complex analysis one can show that the area is

$$\int_{-\infty}^{\infty} D_{n}(x) d x=2 \pi \nonumber$$

Thus, the area is constant for each $n$.

There are two main properties that define a Dirac delta function. First one has that the area under the delta function is one,

$$\int_{-\infty}^{\infty} \delta(x) d x=1 \nonumber$$

Integration over more general intervals gives

$$\int_{a}^{b} \delta(x) d x=1, \quad 0 \in[a, b] \nonumber$$

and

$$\int_{a}^{b} \delta(x) d x=0, \quad 0 \notin[a, b] . \nonumber$$

Another common property is what is sometimes called the sifting property. Namely, integrating the product of a function and the delta function "sifts" out a specific value of the function. It is given by

$$\int_{-\infty}^{\infty} \delta(x-a) f(x) d x=f(a) \nonumber$$

This can be seen by noting that the delta function is zero everywhere except at $x=a$. Therefore, the integrand is zero everywhere and the only contribution from $f(x)$ will be from $x=a$. So, we can replace $f(x)$ with $f(a)$ under the integral. Since $f(a)$ is a constant, we have that

$$\int_{-\infty}^{\infty} \delta(x-a) f(x) d x=\int_{-\infty}^{\infty} \delta(x-a) f(a) d x=f(a) \int_{-\infty}^{\infty} \delta(x-a) d x=f(a) \nonumber$$

Another property results from using a scaled argument, ax. In this case we show that

$$\delta(a x)=|a|^{-1} \delta(x) \nonumber$$

As usual, this only has meaning under an integral sign. So, we place $\delta(a x)$ inside an integral and make a substitution $y=a x$ :

$$\begin{aligned} \int_{-\infty}^{\infty} \delta(a x) d x &=\lim _{L \rightarrow \infty} \int_{-L}^{L} \delta(a x) d x \\[4pt] &=\lim _{L \rightarrow \infty} \dfrac{1}{a} \int_{-a L}^{a L} \delta(y) d y \end{aligned} \nonumber$$

If $a>0$ then

$$\int_{-\infty}^{\infty} \delta(a x) d x=\dfrac{1}{a} \int_{-\infty}^{\infty} \delta(y) d y \nonumber$$

However, if $a<0$ then

$$\int_{-\infty}^{\infty} \delta(a x) d x=\dfrac{1}{a} \int_{\infty}^{-\infty} \delta(y) d y=-\dfrac{1}{a} \int_{-\infty}^{\infty} \delta(y) d y \nonumber$$

The overall difference in a multiplicative minus sign can be absorbed into one expression by changing the factor $1 / a$ to $1 /|a|$. Thus,

$$\int_{-\infty}^{\infty} \delta(a x) d x=\dfrac{1}{|a|} \int_{-\infty}^{\infty} \delta(y) d y . \nonumber$$

Example 8.6. Evaluate $\int_{-\infty}^{\infty}(5 x+1) \delta(4(x-2)) d x$. This is a straight forward integration:

$$\int_{-\infty}^{\infty}(5 x+1) \delta(4(x-2)) d x=\dfrac{1}{4} \int_{-\infty}^{\infty}(5 x+1) \delta(x-2) d x=\dfrac{11}{4} \nonumber$$

A more general scaling of the argument takes the form $\delta(f(x))$. The integral of $\delta(f(x))$ can be evaluated depending upon the number of zeros of $f(x)$. If there is only one zero, $f\left(x_{1}\right)=0$, then one has that

$$\int_{-\infty}^{\infty} \delta(f(x)) d x=\int_{-\infty}^{\infty} \dfrac{1}{\left|f^{\prime}\left(x_{1}\right)\right|} \delta\left(x-x_{1}\right) d x \nonumber$$

This can be proven using the substitution $y=f(x)$ and is left as an exercise for the reader. This result is often written as

$$\delta(f(x))=\dfrac{1}{\left|f^{\prime}\left(x_{1}\right)\right|} \delta\left(x-x_{1}\right) . \nonumber$$

Example 8.7. Evaluate $\int_{-\infty}^{\infty} \delta(3 x-2) x^{2} d x$.

This is not a simple $\delta(x-a)$. So, we need to find the zeros of $f(x)=3 x-2$. There is only one, $x=\dfrac{2}{3}$. Also, $\left|f^{\prime}(x)\right|=3$. Therefore, we have

$$\int_{-\infty}^{\infty} \delta(3 x-2) x^{2} d x=\int_{-\infty}^{\infty} \dfrac{1}{3} \delta\left(x-\dfrac{2}{3}\right) x^{2} d x=\dfrac{1}{3}\left(\dfrac{2}{3}\right)^{2}=\dfrac{4}{27} . \nonumber$$

Note that this integral can be evaluated the long way by using the substitution $y=3 x-2$. Then, $d y=3 d x$ and $x=(y+2) / 3$. This gives

$$\int_{-\infty}^{\infty} \delta(3 x-2) x^{2} d x=\dfrac{1}{3} \int_{-\infty}^{\infty} \delta(y)\left(\dfrac{y+2}{3}\right)^{2} d y=\dfrac{1}{3}\left(\dfrac{4}{9}\right)=\dfrac{4}{27} \nonumber$$

More generally, one can show that when $f\left(x_{j}\right)=0$ and $f^{\prime}\left(x_{j}\right) \neq 0$ for $x_{j}$, $j=1,2, \ldots, n$, (i.e.; when one has $n$ simple zeros), then

$$\delta(f(x))=\sum_{j=1}^{n} \dfrac{1}{\left|f^{\prime}\left(x_{j}\right)\right|} \delta\left(x-x_{j}\right) . \nonumber$$

Example 8.8. Evaluate $\int_{0}^{2 \pi} \cos x \delta\left(x^{2}-\pi^{2}\right) d x$

In this case the argument of the delta function has two simple roots. Namely, $f(x)=x^{2}-\pi^{2}=0$ when $x=\pm \pi$. Furthermore, $f^{\prime}(x)=2 x$. Therefore, $\left|f^{\prime}(\pm \pi)\right|=2 \pi$. This gives

$$\delta\left(x^{2}-\pi^{2}\right)=\dfrac{1}{2 \pi}[\delta(x-\pi)+\delta(x+\pi)] . \nonumber$$

Inserting this expression into the integral and noting that $x=-\pi$ is not in the integration interval, we have

$$\begin{aligned} \int_{0}^{2 \pi} \cos x \delta\left(x^{2}-\pi^{2}\right) d x &=\dfrac{1}{2 \pi} \int_{0}^{2 \pi} \cos x[\delta(x-\pi)+\delta(x+\pi)] d x \\[4pt] &=\dfrac{1}{2 \pi} \cos \pi=-\dfrac{1}{2 \pi} \end{aligned} \nonumber$$

Finally, we previously noted there is a relationship between the Heaviside, or step, function and the Dirac delta function. We defined the Heaviside function as

$$H(x)=\left\{\begin{array}{l} 0, x<0 \\[4pt] 1, x>0 \end{array}\right. \nonumber$$

Then, it is easy to see that $H^{\prime}(x)=\delta(x)$.

Green’s Function Differential Equation

As noted, the Green’s function satisfies the differential equation

$$\dfrac{\partial}{\partial x}\left(p(x) \dfrac{\partial G(x, \xi)}{\partial x}\right)+q(x) G(x, \xi)=\delta(x-\xi) \nonumber$$

and satisfies homogeneous conditions. We have used the Green’s function to solve the nonhomogeneous equation

$$\dfrac{d}{d x}\left(p(x) \dfrac{d y(x)}{d x}\right)+q(x) y(x)=f(x) \nonumber$$

These equations can be written in the more compact forms

$$\begin{gathered} \mathcal{L}[y]=f(x) \\[4pt] \mathcal{L}[G]=\delta(x-\xi) \end{gathered} \nonumber$$

Multiplying the first equation by $G(x, \xi)$, the second equation by $y(x)$, and then subtracting, we have

$$G \mathcal{L}[y]-y \mathcal{L}[G]=f(x) G(x, \xi)-\delta(x-\xi) y(x) . \nonumber$$

Now, integrate both sides from $x=a$ to $x=b$. The left side becomes

$$\int_{a}^{b}[f(x) G(x, \xi)-\delta(x-\xi) y(x)] d x=\int_{a}^{b} f(x) G(x, \xi) d x-y(\xi) \nonumber$$

and, using Green’s Identity, the right side is

$$\int_{a}^{b}(G \mathcal{L}[y]-y \mathcal{L}[G]) d x=\left[p(x)\left(G(x, \xi) y^{\prime}(x)-y(x) \dfrac{\partial G}{\partial x}(x, \xi)\right)\right]_{x=a}^{x=b} \nonumber$$

Combining these results and rearranging, we obtain

$$y(\xi)=\int_{a}^{b} f(x) G(x, \xi) d x-\left[p(x)\left(y(x) \dfrac{\partial G}{\partial x}(x, \xi)-G(x, \xi) y^{\prime}(x)\right)\right]_{x=a}^{x=b} \nonumber$$

Next, one uses the boundary conditions in the problem in order to determine which conditions the Green’s function needs to satisfy. For example, if we have the boundary condition $y(a)=0$ and $y(b)=0$, then the boundary terms yield

$$\begin{aligned} y(\xi)=& \int_{a}^{b} f(x) G(x, \xi) d x-\left[p(b)\left(y(b) \dfrac{\partial G}{\partial x}(b, \xi)-G(b, \xi) y^{\prime}(b)\right)\right] \\[4pt] &+\left[p(a)\left(y(a) \dfrac{\partial G}{\partial x}(a, \xi)-G(a, \xi) y^{\prime}(a)\right)\right] \\[4pt] =& \int_{a}^{b} f(x) G(x, \xi) d x+p(b) G(b, \xi) y^{\prime}(b)-p(a) G(a, \xi) y^{\prime}(a) . \end{aligned} \nonumber$$

The right hand side will only vanish if $G(x, \xi)$ also satisfies these homogeneous boundary conditions. This then leaves us with the solution

$$y(\xi)=\int_{a}^{b} f(x) G(x, \xi) d x . \nonumber$$

We should rewrite this as a function of $x$. So, we replace $\xi$ with $x$ and $x$ with $\xi$. This gives

$$y(x)=\int_{a}^{b} f(\xi) G(\xi, x) d \xi . \nonumber$$

However, this is not yet in the desirable form. The arguments of the Green’s function are reversed. But, $G(x, \xi)$ is symmetric in its arguments. So, we can simply switch the arguments getting the desired result.

We can now see that the theory works for other boundary conditions. If we had $y^{\prime}(a)=0$, then the $y(a) \dfrac{\partial G}{\partial x}(a, \xi)$ term in the boundary terms could be made to vanish if we set $\dfrac{\partial G}{\partial x}(a, \xi)=0$. So, this confirms that other boundary value problems can be posed besides the one elaborated upon in the chapter so far.

We can even adapt this theory to nonhomogeneous boundary conditions. We first rewrite Equation (8.62) as

$$y(x)=\int_{a}^{b} G(x, \xi) f(\xi) d \xi-\left[p(\xi)\left(y(\xi) \dfrac{\partial G}{\partial \xi}(x, \xi)-G(x, \xi) y^{\prime}(\xi)\right)\right]_{\xi=a}^{\xi=b} \nonumber$$

Let’s consider the boundary conditions $y(a)=\alpha$ and $y^{\prime}(b)=$ beta. We also assume that $G(x, \xi)$ satisfies homogeneous boundary conditions,

$$G(a, \xi)=0, \quad \dfrac{\partial G}{\partial \xi}(b, \xi)=0 . \nonumber$$

in both $x$ and $\xi$ since the Green’s function is symmetric in its variables. Then, we need only focus on the boundary terms to examine the effect on the solution. We have

$$\begin{aligned} {\left[p(\xi)\left(y(\xi) \dfrac{\partial G}{\partial \xi}(x, \xi)-G(x, \xi) y^{\prime}(\xi)\right)\right]_{\xi=a}^{\xi=b} } &=\left[p(b)\left(y(b) \dfrac{\partial G}{\partial \xi}(x, b)-G(x, b) y^{\prime}(b)\right)\right] \\[4pt] &-\left[p(a)\left(y(a) \dfrac{\partial G}{\partial \xi}(x, a)-G(x, a) y^{\prime}(a)\right)\right] \\[4pt] &=-\beta p(b) G(x, b)-\alpha p(a) \dfrac{\partial G}{\partial \xi}(x, a) . \end{aligned} \nonumber$$

Therefore, we have the solution

$$y(x)=\int_{a}^{b} G(x, \xi) f(\xi) d \xi+\beta p(b) G(x, b)+\alpha p(a) \dfrac{\partial G}{\partial \xi}(x, a) . \nonumber$$

This solution satisfies the nonhomogeneous boundary conditions. Let’s see how it works. Example 8.9. Modify Example $8.4$ to solve the boundary value problem $y^{\prime \prime}=$ $x^{2}, \quad y(0)=1, y(1)=2$ using the boundary value Green’s function that we found:

$$G(x, \xi)=\left\{\begin{array}{l} -\xi(1-x), 0 \leq \xi \leq x \\[4pt] -x(1-\xi), x \leq \xi \leq 1 \end{array}\right. \nonumber$$

We insert the Green’s function into the solution and use the given conditions to obtain

$$\begin{aligned} y(x) &=\int_{0}^{1} G(x, \xi) \xi^{2} d \xi-\left[y(\xi) \dfrac{\partial G}{\partial \xi}(x, \xi)-G(x, \xi) y^{\prime}(\xi)\right]_{\xi=0}^{\xi=1} \\[4pt] &=\int_{0}^{x}(x-1) \xi^{3} d \xi+\int_{x}^{1} x(\xi-1) \xi^{2} d \xi+y(0) \dfrac{\partial G}{\partial \xi}(x, 0)-y(1) \dfrac{\partial G}{\partial \xi}(x) \\[4pt] &=\dfrac{(x-1) x^{4}}{4}+\dfrac{x\left(1-x^{4}\right)}{4}-\dfrac{x\left(1-x^{3}\right)}{3}+(x-1)-2 x \\[4pt] &=\dfrac{x^{4}}{12}+\dfrac{35}{12} x-1 \end{aligned} \nonumber$$

Of course, this problem can be solved more directly by direct integration. The general solution is

$$y(x)=\dfrac{x^{4}}{12}+c_{1} x+c_{2} . \nonumber$$

Inserting this solution into each boundary condition yields the same result.

We have seen how the introduction of the Dirac delta function in the differential equation satisfied by the Green’s function, Equation (8.59), can lead to the solution of boundary value problems. The Dirac delta function also aids in our interpretation of the Green’s function. We note that the Green’s function is a solution of an equation in which the nonhomogeneous function is $\delta(x-\xi)$. Note that if we multiply the delta function by $f(\xi)$ and integrate we obtain

$$\int_{-\infty}^{\infty} \delta(x-\xi) f(\xi) d \xi=f(x) \nonumber$$

We can view the delta function as a unit impulse at $x=\xi$ which can be used to build $f(x)$ as a sum of impulses of different strengths, $f(\xi)$. Thus, the Green’s function is the response to the impulse as governed by the differential equation and given boundary conditions.

In particular, the delta function forced equation can be used to derive the jump condition. We begin with the equation in the form

$$\dfrac{\partial}{\partial x}\left(p(x) \dfrac{\partial G(x, \xi)}{\partial x}\right)+q(x) G(x, \xi)=\delta(x-\xi) \nonumber$$

Now, integrate both sides from $\xi-\epsilon$ to $\xi+\epsilon$ and take the limit as $\epsilon \rightarrow 0$. Then,

$$\begin{aligned} \lim _{\epsilon \rightarrow 0} \int_{\xi-\epsilon}^{\xi+\epsilon}\left[\dfrac{\partial}{\partial x}\left(p(x) \dfrac{\partial G(x, \xi)}{\partial x}\right)+q(x) G(x, \xi)\right] d x &=\lim _{\epsilon \rightarrow 0} \int_{\xi-\epsilon}^{\xi+\epsilon} \delta(x-\xi) d x \\[4pt] &=1 \end{aligned} \nonumber$$

Since the $q(x)$ term is continuous, the limit of that term vanishes. Using the Fundamental Theorem of Calculus, we then have

$$\lim _{\epsilon \rightarrow 0}\left[p(x) \dfrac{\partial G(x, \xi)}{\partial x}\right]_{\xi-\epsilon}^{\xi+\epsilon}=1 . \nonumber$$

This is the jump condition that we have been using!

Series Representations of Green’s Functions

There are times that it might not be so simple to find the Green’s function in the simple closed form that we have seen so far. However, there is a method for determining the Green’s functions of Sturm-Liouville boundary value problems in the form of an eigenfunction expansion. We will finish our discussion of Green’s functions for ordinary differential equations by showing how one obtains such series representations. (Note that we are really just repeating the steps towards developing eigenfunction expansion which we had seen in Chapter 6.)

We will make use of the complete set of eigenfunctions of the differential operator, $\mathcal{L}$, satisfying the homogeneous boundary conditions:

$$\mathcal{L}\left[\phi_{n}\right]=-\lambda_{n} \sigma \phi_{n}, \quad n=1,2, \ldots \nonumber$$

We want to find the particular solution $y$ satisfying $\mathcal{L}[y]=f$ and homogeneous boundary conditions. We assume that

$$y(x)=\sum_{n=1}^{\infty} a_{n} \phi_{n}(x) . \nonumber$$

Inserting this into the differential equation, we obtain

$$\mathcal{L}[y]=\sum_{n=1}^{\infty} a_{n} \mathcal{L}\left[\phi_{n}\right]=-\sum_{n=1}^{\infty} \lambda_{n} a_{n} \sigma \phi_{n}=f \nonumber$$

This has resulted in the generalized Fourier expansion

$$f(x)=\sum_{n=1}^{\infty} c_{n} \sigma \phi_{n}(x) \nonumber$$

with coefficients

$$c_{n}=-\lambda_{n} a_{n} \nonumber$$

We have seen how to compute these coefficients earlier in the text. We multiply both sides by $\phi_{k}(x)$ and integrate. Using the orthogonality of the eigenfunctions,

$$\int_{a}^{b} \phi_{n}(x) \phi_{k}(x) \sigma(x) d x=N_{k} \delta_{n k} \nonumber$$

one obtains the expansion coefficients (if $\lambda_{k} \neq 0$ )

$$a_{k}=-\dfrac{\left(f, \phi_{k}\right)}{N_{k} \lambda_{k}}, \nonumber$$

where $\left(f, \phi_{k}\right) \equiv \int_{a}^{b} f(x) \phi_{k}(x) d x$.

As before, we can rearrange the solution to obtain the Green’s function. Namely, we have

$$y(x)=\sum_{n=1}^{\infty} \dfrac{\left(f, \phi_{n}\right)}{-N_{n} \lambda_{n}} \phi_{n}(x)=\int_{a}^{b} \underbrace{\sum_{n=1}^{\infty} \dfrac{\phi_{n}(x) \phi_{n}(\xi)}{-N_{n} \lambda_{n}}}_{G(x, \xi)} f(\xi) d \xi \nonumber$$

Therefore, we have found the Green’s function as an expansion in the eigenfunctions:

$$G(x, \xi)=\sum_{n=1}^{\infty} \dfrac{\phi_{n}(x) \phi_{n}(\xi)}{-\lambda_{n} N_{n}} . \nonumber$$