Differential - Equations - Book by Daniel An
Differential - Equations - Book by Daniel An
Differential - Equations - Book by Daniel An
Engineers
Daniel An
November 6, 2017
This work is licensed under a Creative Commons “Attribution-
NonCommercial-ShareAlike 3.0 Unported” license.
1
Contents
Preface 4
2
8 Laplace Transform Method 139
8.1 Laplace Transform and Its Inverse . . . . . . . . . . . . . . . . . . . . . . . . . . . 139
8.2 Solving Differential Equations using Laplace Transforms . . . . . . . . . . . . . . 145
Appendices 153
3
Preface
This book is an outcome of my 9 years of teaching differential equations in SUNY Maritime Col-
lege. I have taught at least 2 courses of Differential equations per year and in some years I have
taught as many as 4 counting the summer course. Over many years of teaching I felt the need for
a text book that gives detailed step by step solutions, and this book aims to do just that. When
writing solutions to example problems, for most solutions I did not solve them on paper but solved
them as I was typing them on the computer. What resulted was solutions that may seem to take
too much space but including a lot of details that will be usually left out in a text book.
Another feature of this book is I tried not to use numbers for equations and reference previous
equations that way. This way, the readers do not have to go back and forth when trying to read a
solution of a problem. (However, this was not possible when solving partial differential equations
because it had so many parts to it.)
I do not plan to publish this book anywhere, but rather ride the Creative Commons band-
wagon. One of the reasons for this decision is because I am writing this book under the grant
received from New York State SUNY OER (Open Educational Resources) initiative. I thank the
SUNY Maritime Library department for applying and receiving this fund, which made this text-
book possible (and free).
Since this textbook is geared towards Engineers, most proofs are omitted except for a few
derivations that I feel is critical in order to understand the method of solutions.
Daniel An
Department of Science
SUNY Maritime College
4
First Order Differential Equations
1
1.1 Introduction to Differential Equations and Some Basic
Examples
Introduction to differential equations
A differential equation is an equation that contains derivatives. It’s that simple. For example,
y0 = x .
The above differential equation can be solved by integrating both sides. This yields the solution
y = 21 x 2 + C. There are two things you should notice about this solution. First is that the solution
is a function. In other words, when we say we are looking for a solution of a differential equation,
we are looking for a function that satisfies the differential equation. Second is that the solution
comes with a free parameter, which we denoted by C.
If the differential equation (sometimes abbreviated as DE) involves the second derivative, it
will have two free parameters because, loosely speaking, the solution will involve two integrals.
To state this property exactly, we need to define order of a differential equation. The order of a
differential equation is the highest order of derivative appearing in the differential equation. An
n-th order differential equation will have n free parameters.
5
P(t) be the population of the bacteria at time t, the rate of change P 0(t) will be proportional to
the population P(t). When one is proportional to another, we can make an equation by using a
proportionality constant, which in this context is:
P 0(t) = kP(t)
Here is another example. Newton’s famous equation F = ma can be thought as F = mv 0,
because acceleration is the derivative of velocity. Often the friction force on a moving object
is proportional to how fast it is moving, i.e. the faster it is moving the more friction it will
experience. In such a case, we can say F = −cv for some proportionality constant c (which
should be obtained via experiments). The reason for the negative sign is because the friction
force always acts against the velocity. Then the Newton’s equation becomes:
−cv = mv 0,
and we obtain a first order differential equation for the unknown function v(t). The method for
solving such differential equation will be explained later in this book.
There are many more examples. We will solve a few simple examples in the first chapter, and
then more complicated examples will be discussed as we learn.
Direct integration
As we have seen, the differential equation y 0 = x can be solved by integrating both sides. Some-
times differential equations come with an extra condition called initial condition, such as y(1) =
2. A differential equation with an initial condition is called initial value problem (often abbre-
viated as IVP). This means we are looking for a solution of y 0 = x (which we already know as
y = 21 x 2 + C) that satisfies y(1) = 2. Since y(1) = 21 (1)2 + C = 12 + C, we have the following
equation:
1
+C = 2
2
The solution is C = 32 and thus the solution of the initial value problem is:
1 3
y = x2 +
2 2
A differential equation that can be solved by simply integrating both sides is called direct
integration problem. Here is another example of the direct integration problem:
Example 1.1.1
3
y0 = , y(1) = 0
x3 +x
Solution) Integrate both sides to get:
3
∫
y= dx
x3 +x
To solve the integral on the right, we need partial fraction expansion. Note that x 3 +x = x(x 2 + 1).
So set:
3 A Bx + C
3
= + 2
x +x x x +1
6
Then, multiplying x 3 + x, we get:
3 3 3
∫ ∫
−3x
3
dx = + 2 dx = 3 ln x − ln x 2 + 1 + C
x +x x x +1 2
The last integral was done by using u-substitution with
u = x2 + 1
3 3
y = 3 ln x − ln x 2 + 1 + ln 2
2 2
Since integration is needed in order to solve direct integral problems, you will need to review
various methods of integrations such as integration by parts, integration by substitution, partial
fraction expansion, etc.
Example 1.1.2 A car weighs 2000kg and is going at 16m/s. If the brakes provide 8000N of friction
force, how long does it travel before stopping completely?
−8000 = 2000x 00
Since x 0 is the velocity and the initial velocity was 16 m/s, we have x 0(0) = 16. Also, let’s call
the initial position as 0 m(so that we would easily know the distance traveled once we know its
terminal position.) This means x(0) = 0. Notice that for second order differential equations, we
7
need two initial conditions. The differential equation −8000 = 2000x 00, solved for x 00 gives us
x 00 = −4. Integrating both sides, we get:
x 0 = −4t + C
x 0 = −4t + 16
Integrating this yet again, we have x = −2t 2 + 16t +C 2 . Then the initial condition x(0) = 0 implies
C 2 = 0 and we have:
x = −2t 2 + 16t
Now, we want to know when the car stops, i.e., when velocity=0. This means x 0 = −4t + 16 = 0.
The solution is t = 4, and at that time, the position is x(4) = −2(4)2 + 16 · 4 = 32. The answer is
32 m. (Note: One can use the well known physics shortcut v 2f − vi2 = 2a∆x to find the answer
easily. It will verify that our answer is correct.)
Natural growth
Think about a bacteria culture in a big tank. Denote P(t) as the population of bacteria at time
t. Take a moment to convince yourself that the bigger the population is, the faster it will grow
, because bacteria grow by dividing. Here, we are assuming here that there are large amount of
nutrients and space for bacteria to grow, because if there are limited resources, then too much
population may give rise to competition and thus could actually slow growth. (There is such a
model where that is also taken into consideration. That model is called the logistic growth model.)
The statement ”Rate of change of population is proportional to the population” can be translated
into: P 0(t) = kP(t), which means ”Rate of change of population is some constant k times the
population.” The above is a differential equation as you can see. However, to solve this one needs
to learn a new technique. Here is how it is solved. First, divide by P(t) both sides:
P 0(t)
=k
P(t)
Now, integrate both sides:
P 0(t)
∫ ∫
dt = kdt
P(t)
The integral on the right side is easy: kdt = kt + C However, on the left side, one has to do
∫
ln P(t) = kt + C
8
Exponentiating both sides, we get:
e ln P(t) = e kt+C
Using the fact that exponentiation cancels ln, we have:
P(t) = e kt eC
P(t) = C̃e kt
Now, since C is something we invented along the way of the solution, we can switch the names
and call C̃ as C instead (and vice versa). People who are waiting for our result wouldn’t care what
we call it, and in fact it is lot simpler just to say P(t) = Ce kt .
Here is an example problem that uses this result:
Example 1.1.3 At noon a bacteria culture of size 1000 was put into a tank. At 2 pm, there were
2000 bacteria in the tank. At what time will there be 1 million bacteria?
After finding such function P(t) we want to know what value of t will give us P(t) = 1000000.
We know from the above solution that P(t) = Ce kt Then P(0) = 1000 yields C = 1000. Also
P(2) = 2000 yields 2000 = 1000e 2k , which can be solved to reveal that k = 12 ln 2. Thus the
function P(t) is obtained:
1
P(t) = 1000e 2 (ln 2)t
Then we solve:
1
1000e 2 (ln 2)t = 1000000
1
e 2 (ln 2)t = 1000
1
(ln 2)t = ln 1000
2
2 ln 1000
t=
ln 2
The value is approximately 20 hours.
y 0 = f (x)д(y)
9
In other words, if y 0 is on the left side and you can make the right side to look like a product of a
function of x and a function of y, then it is called a separable equation. Such functions are easy
to solve because you can divide both sides by д(y)
1 0
y = f (x)
д(y)
(i.e. you can ”separate” x and y.), and then integrate both sides:
1 0
∫ ∫
y dx = f (x)dx
д(y)
Remember that y is really y(x), a function of x. The left side can be integrated by u-substitution
if one sets u = y(x), so that du = y 0(x)dx and the left side becomes:
1
∫
du
д(u)
Another way to explain the above method of solution is to use the Liebniz notation for y 0,
dy 1 0
namely, dx . Then д(y) y = f (x) can be written as:
1 dy
· = f (x)
д(y) dx
Then we integrate both sides:
1 dy
∫ ∫
· dx = f (x)dx
д(y) dx
After canceling the dx on the left side, we have:
1
∫ ∫
dy = f (x)dx
д(y)
This is the more common way of solving separable equations.
Example 1.2.1
y 0 = 3xy
dy
Step 1. Rewrite y 0, as dx .
dy
= 3xy
dx
Step 2. Divide both sides by y (or any function so that the right side only has a function of x)
1 dy
· = 3x
y dx
Step 3. Integrate both sides
1 dy
∫ ∫
· dx = 3xdx
y dx
10
Canceling dx, we have
1
∫ ∫
dy = 3xdx
y
The integral is:
3
ln |y| = x 2 + C
2
(When you integrate both sides, you should put the integral constant only on one side.)
Step 4. Solve for y We exponentiate both sides:
3 2
e ln |y| = e 2 x +C
11
More examples of Separable Differential Equations
Example 1.2.2
dy
= xy 2 + 4x + 2y 2 + 8
dx
Solution) Factor the right side by grouping:
Example 1.2.3
dy
= e x+2y , y(0) = 1
dx
Solution)
dy
= e x e 2y
dx
Divide by e 2y :
1
dy = e x dx
e 2y
12
Which is same as:
e −2ydy = e x dx
Integrate: ∫ ∫
e −2y
dy = e x dx
1
− e −2y = e x + C
2
Multiply by (-2) and take ln to get:
1
y = − ln (−2e x − 2C)
2
This is not yet the solution because we haven’t yet applied the initial condition y(0) = 1. Accord-
ing to the solution above,
1 1
y(0) = − ln (−2e 0 − 2C) = − ln (−2 − 2C)
2 2
Since we are given that y(0) = 1, we need this to be equal to 1.
1
− ln (−2 − 2C) = 1
2
ln (−2 − 2C) = −2
−2 − 2C = e −2
−2C = 2 + e −2
1
C = −1 − e −2
2
1 x
Plugging this back into y = − 2 ln (−2e − 2C) leads to the answer:
1
y = − ln (−2e x + 2 + e −2 )
2
Example 1.2.4
y 0 = xe 2x+y
Rewriting the right side as a product:
dy
= xe 2x ey
dx
Divide both sides by ey and multiply both sides by dx to get:
1
dy = xe 2x dx
ey
Then we integrate: ∫ ∫
e dy =
−y
xe 2x dx
The right side requires integration by parts, or tabular integration. (See appendix for explanation
of tabular integration.) The tabular integration is done as follows:
13
This yields 21 xe 2x − 14 e 2x + C and the result of the integration is:
1 1
−e −y = xe 2x − e 2x + C
2 4
Solving this for y gives us:
1 1
y = − ln (− xe 2x + e 2x − C)
2 4
Example 1.2.5
x(x + 1)(x + 2)y 0 = (3x − 3)y 2
Solution)
dy
x(x + 1)(x + 2) = (3x − 3)y 2
dx
Divide both sides by x(x + 1)(x + 2) and also by y 2 to get:
1 3x − 3
2
dy = dx
y x(x + 1)(x + 2)
Now we have to integrate both sides. The left side is simply:
1
∫ ∫
2
dy = y −2dy = −y −1
y
However, the right side requires something called partial fraction expansion:
3x − 3 A B C
= + +
x(x + 1)(x + 2) x x + 1 x + 2
We have to find A, B and C that satisfies this algebraic identity and then proceed to integrate. To
do this, we multiply x(x + 1)(x + 2) both sides to get:
What we have above is an algebraic identity, which means that the above equation has to hold
for all values of x. So whatever value we make x to equal to, the equation should hold, e.g. x = 0:
14
−3 = 2A + 0 + 0
3
A=−
2
Then for x = −1:
−6 = 0 − B + 0
B=6
Then for x = −2:
−9 = 0 + 0 + 2C
9
C=−
2
Using these values for A, B and C, we have an algebraic identity
3x − 3 −3 6 − 29
= 2 + +
x(x + 1)(x + 2) x x +1 x +2
Thus integrating the left side is same as integrating the right side:
∫ 3
3x − 3 −2 6 − 92
∫
dx = + + dx
x(x + 1)(x + 2) x x +1 x +2
3 9
= − ln |x | + 6 ln |x + 1| − ln |x + 2| + C
2 2
So we have:
3 9
−y −1 = − ln |x | + 6 ln |x + 1| − ln |x + 2| + C
2 2
Solving this for y, we get:
1
y= 3
2 ln |x | − 6 ln |x + 1| + 92 ln |x + 2| − C
1 2 1 2
y 0 − xy = Ce 2 x · x − x · Ce 2 x = 0
1 2
In the above, the solution y = Ce 2 x is called the general solution of the differential equation
1 2
y 0 − xy = 0. It means that any function that satisfies y 0 − xy = 0 can be written as y = Ce 2 x for
some value of C.
15
Suppose instead we had a initial value problem:
y 0 − xy = 0, y(0) = 3
1 2
Then the initial condition y(0) = 3 leads to y(0) = Ce 2 0 = C · 1 = C = 3. Hence the answer to
1 2
the initial value problem is y = 3e 2 x . Such an answer is called a specific solution, because it has
a specified value for the parameter C.
Solution) Since the ambient temperature is 400◦ F , the Newton’s equation becomes:
dT
= k(400 − T )
dt
Divide by (400 − T ) and integrate:
1
∫ ∫
dT = kdt
400 − T
− ln(400 − T ) = kt + C
Solving for T gives us:
T = 400 − e −kt−C
Now, using the trick of introducing C̃ and then erasing the tilde, we can write:
T = 400 − Ce −kt
Now, we apply the two conditions: a) Initial temperature was 32 : T (0) = 32 b) Temperature
became 50 in 2 minutes : T (2) = 50. T (0) = 32 implies that
32 = 400 − Ce 0
C = 400 − 32 = 368
Then T (2) = 50 implies that
50 = 400 − 368e −2k
Solving for k, we get:
1 350
k = − ln = 0.02507
2 368
16
Thus the temperature at time t is:
Example 1.2.8
(x 2 + 4)y 0 = 3xy
Solution)
dy
(x 2 + 4)= 3xy
dx
1 3x
dy = 2 dx
y x +4
1 3x
∫ ∫
dy = 2
dx
y x +4
3
ln y = ln (x 2 + 4) + C
2
3 2 +4)+C
y = e 2 ln (x
2 +4)3/2
y = C̃e ln (x
y = C̃(x 2 + 4)3/2
So the answer is
y = C(x 2 + 4)3/2
17
1.3 First Order Differential Equations
Introduction to first order linear differential equations
Suppose we have an equation, say,
x 2y = e 3x
Now, try differentiating the above using implicit differentiation (hopefully you still remember how
to do implicit differentiation).
x 2y = e 3x
0 0
The right side is just a regular chain rule, which yields 3e 3x . However, when you differentiate the
left side,
0 you have to know that y 2is actually
0 y(x), i.e. it is a function of x. So when you differentiate
2
x y , think of differentiating x f (x) , which means you have to differentiate using the product
rule:
x 2y = x 2 y + x 2 (y)0 = 2xy + x 2y 0
0 0
Definition 1.3.1 If a differential equation satifies the following properties, we call it a first order
linear differential equation.
• All appearances of y and y 0 have power of 1
• y and y 0 are not multiplied together
• y and y 0 are not nested inside a function such as cos, sin, exp, square roots … etc.
In general, a linear differential equation of any order is the above property holding for any
order derivative of y.
18
Definition 1.3.2 Any first order linear differential equation can be put into the form
y 0 + p(x)y = q(x)
Now, let’s trace back what we have done to solve the following problem:
2xy + x 2y 0 = 3e 3x
Now, notice that the left side is same as x 2y , if we apply the product rule to x 2y , we get
0 0
2xy + x 2y 0 (Again, we know this because we are reverse-engineering.) So we can write:
x 2y = 3e 3x
0
Notice that on the left side, the integral and the derivative will cancel each other. So we have:
∫
2
x y= 3e 3x dx
x 2y = e 3x + C
19
The Multiplier for the First Order Linear Differential Equation
Let’s say we are given a first order differential equation in the standard form:
y 0 + p(x)y = q(x)
The aim of this section is to find the multiplier M(x) that we can multiply on both sides of the
above equation so that the left side turns into a product rule form, i.e.
such that the left side M(x)y 0 +M(x)p(x)y is equal to (M(x)y)0. What does this requirement mean?
Well, first try to work out the product rule:
If we compare this with M(x)y 0 + M(x)p(x)y, we realize that what we want is M(x)p(x) = M 0(x).
If this is true, then we do have M(x)y 0 + M(x)p(x)y = M(x)y 0 + M 0(x)y = (M(x)y)0 and the left
side of M(x)y 0 + M(x)p(x)y = M(x)q(x) does become a product rule form. Let’s try to figure out
what M(x) should be:
M(x)p(x) = M 0(x)
dM
Mp(x) =
dx
1
p(x)dx = dM
M
1
∫ ∫
p(x)dx = dM
M
∫
p(x)dx = ln M + C
Now, the above suggests that there are many different possible M(x), because each time we choose
a different C we will have a different M(x). However, we are only looking for one such function.
So let’s set C = 0 and move on: ∫
p(x)dx = ln M
e ln M = e p(x)dx
∫
p(x)dx
∫
M(x) = e
So we have obtained an important result:
p(x)dx
∫
Definition 1.3.3 Multiplier. To solve y 0 + p(x)y = q(x), one has to compute M(x) = e and
multiply it to both sides. The result will be that the left side becomes same as (M(x)y)0
Example 1.3.2
y 0 + 2y = x
20
Solution) Since p(x) = 2,
2dx
= e 2x
∫
M(x) = e
(Notice that we did not put C when we integrated, because we only need one such function.)
Multiply this both sides:
e 2x y 0 + 2e 2x y = xe 2x
This turns the left side into a product rule: e 2x y = xe 2x So we integrate:
0
∫
2x
e y= xe 2x dx
(Notice that the derivative on the left side disappeared because it cancelled with the integral.) The
integral on the right can be computed using tabular integration as follows:
Solution) Since p(x) = −2x (Note: must include the negative sign if there is one!),
2
∫
M(x) = e (−2x)dx
= e −x
21
∫
−x 2 2
e y= 3xe −x dx
1
∫ ∫
−x 2 u
3xe dx = 3xe − du
2x
3 3 3 2
∫
= − eu du = − eu + C = − e −x + C
2 2 2
Using this result, we have:
2 3 2
e −x y = − e −x + C
2
Solving for y, we get the general solution:
3 2
y = − + Ce x
2
To get the specific solution satisfying y(0) = 1, we solve:
3 2
1 = − + Ce 0
2
3 5
C =1+ =
2 2
Answer:
3 5 2
y = − + ex
2 2
Example 1.3.4 Solve √
2xy 0 = y + 3 x 5
Solution) Since we see that y 0 has power of 1 and also y has power of 1, we know that this is a
first order linear differential equation. We solve it using the following steps:
Step 1 Change to standard form Divide by the coefficient of y 0, which is 2x:
1 3x 5/2
y0 = y+
2x 2x
Then make it into y 0 + p(x)y = q(x) form by moving terms:
1 3x 5/2
y − y=
0
2x 2x
Step 2 Identify the function p(x) and compute M(x) Since − 2x1 is in front of y in the standard
form, that is our p(x). Make sure you include the negative sign if there is one. Then
− 2x1 dx 1 1
x dx
∫ ∫
M(x) = e = e− 2
1
= e − 2 ln x = e ln x
−1/2
= x −1/2
22
In the above calculation study carefully how the expression was simplified. You must simplify
M(x) in order to use it. (If you don’t you’ll run into very complicated integrals.) Remember, e
and ln only cancel if they are directly next to each other. That’s why I had to move the − 12 away
before cancelling the two.
Step 3 Multiply the M(x) to the standard form and make the left side into a product rule form
I want to stress again that M(x) should be multiplied to the standard form, not the original given
differential equation. Thus:
1 3x 5/2
x −1/2y 0 − x −1/2 y = x −1/2
2x 2x
0 3x 5/2
x −1/2y = x −1/2
2x
Ideally, you would check that the two equations written above are indeed the same thing by
applying the product rule to the bottom equation, but as you become more confident in solving
differential equations this way, you can just write the above form by faith. You also have to
simplify the right side.
3 3
∫
x −1/2
y= xdx = x 2 + C
2 4
3
y = x 5/2 + Cx 1/2
4
Example 1.3.5
xy 0 − 5y = 2x, y(1) = 0
23
1
x −5y = − x −4 + C
2
1
y = − x + Cx 5
2
Now, we apply y(1) = 0
1 1
0 = y(1) = − · 1 + C(1)5 = − + C
2 2
1
C=
2
Answer:
1 1
y = − x + x5
2 2
When doing the above integral, a handy formula to know is the following:
f (x)
∫ 0
dx = ln f (x) + C
f (x)
24
∫
2 3/2
(x + 1) y= (x 2 + 1)1/2 · 6xdx
Using u = x 2 + 1, we get:
∫ ∫
2 1/2
(x + 1) · 6xdx = 3 u 1/2du = 2u 3/2 + C = 2(x 2 + 1)3/2 + C
So we have:
(x 2 + 1)3/2y = 2(x 2 + 1)3/2 + C
Answer is:
2(x 2 + 1)3/2 + C C
y= 2
= 2(x 2 + 1)1/2 + 2
x +1 x +1
− x1 1 u 1 2 u1
∫ ∫ ∫ ∫
e dx = e 3 x du = e du = −euudu
x3 x x
This needs integration by parts, or tabular integration.
1
∫
1
−euudu = −euu + eu + C = (1 − u)eu + C = (1 + )e − x + C
x
Thus we have:
1 −1
− x1
e y = 1+ e x +C
x
Now, the general solution is:
1 1
y =1+ + Ce x
x
25
To satisfy y(1) = 0,
1 1
0 = y(1) = 1 + + Ce 1 = 2 + eC
1
2
C=−
e
So the answer is:
1 2 1
y =1+ + − ex
x e
Example 1.4.1
x 2y + 5xy 2 + 3y 3 = 0
Find y 0 using implicit differentiation.
Solution)
x 2y + 5xy 2 + 3y 3 = 0
0 0 0
x 2y 0 + 2xy + 10xyy 0 + 5y 2 + 9y 2y 0 = 0
After this, one has to solve the above equation for y 0 to get the answer. Now, let’s try to look at
this a bit differently than what you were taught in Calculus I.
F (x, y) = 0
In our example above, it was x 2y + 5xy 2 + 3y 3 = 0, and the left side can be seen as F (x, y).
Now, although F (x, y) is a function of x and y, to differentiate implicitly is to treat y as implicitly
being a function of x. In other words, the function F , and variables y and x have the following
dependence:
26
A line drawn from F to x means F is a quantity that changes when x is changed. Other lines
can be interpreted similarly. Such dependence diagram leads to a chain rule of the form:
dF ∂F ∂F dy
= +
dx ∂x ∂y dx
You should read this chain rule as ”The rate of change of F when you change x is the sum of two
things. First is the rate of change of F when x is changed but y is fixed. Second is the rate of
change of F when change of x gives rise to change in y which in turn changes F .” If you see the
dependence diagram again, there are two ’routes’ from x to F : one going directly and another
∂F dy
going through y. The two terms in ∂F∂x + ∂y dx measures the change in F brought about by change
of x going through each of the two ’routes’. The magic of the multivariable chain rule is that if
you add these up, you can find the total rate of change of F brought about by change of x. What
does this have anything to do with the implicit differentiation? Well, differentiating both sides
F (x, y) = 0 by x is:
∂F ∂F dy
+ =0
∂x ∂y dx
where the left side is the chain rule. Let’s rewrite this equation using slightly different notation.
dy
First we use Fx for ∂F
∂x , and also likewise Fy . Then instead of dx we use y . The equation now reads:
0
Fx + Fyy 0 = 0
Fyy 0 = −Fx
Fx
y0 = −
Fy
F (x, y) = 0,
F (x, y) = 0 → Fx + Fyy 0 = 0
How does the above fact solve differential equation? Here is an example problem:
27
Solution) First, gather all the terms that have y 0 attached and factor the y 0 from them to get:
2xy + 5y 2 + (x 2 + 10xy + 9y 2 )y 0 = 0
Now what we aim to do is to find out what it was before the differentiation. In other words,
we are imagining the original equation as an equation of the form F (x, y) = 0 and then implicit
differentiation was applied to get the above equation. Since the derivative of F (x, y) = 0 is Fx +
Fyy 0 = 0, comparing this with the above equation reveals that:
Fx = 2xy + 5y 2
Fy = (x 2 + 10xy + 9y 2 )
Integrating Fx = 2xy + 5y 2 by x yields:
∫
F= 2xy + 5y 2dx = x 2y + 5y 2x + д(y)
Notice that instead of C, I wrote д(y). The reason for this is that since partial derivative Fx com-
puted by regarding y as a constant, the things that get deleted by derivatives could have been a
function of y. For similar reasons, if we integrate Fy = (x 2 + 10xy + 9y 2 ) we will need to put a
function of x: ∫
F = (x 2 + 10xy + 9y 2 )dy = x 2y + 5xy 2 + 3y 3 + h(x)
Now, this is similar to figuring out the details of an unknown object (in this case F (x, y)) from
two different reports (the two integral results we see above.) To get a clear picture, we have to
harmonize the two different accounts of F (x, y). If you compare them, you can see that indeed
the two descriptions do agree that F (x, y) must contain x 2y + 5xy 2 . Also, you see that what we
called д(y) in the first integral must be same as 3y 3 , while h(x) must be zero. This investigation
can conclude what F (x, y) is up to a constant, i.e.
F (x, y) = x 2y + 5xy 2 + 3y 3 + C
(If you studied hard in multivariable calculus, you will recognize this method as finding the po-
tential function of a conservative vector field problem.) Now that we found the F (x, y), we can
claim that the original equation before differentiation is
x 2y + 5xy 2 + 3y 3 + C = 0
To verify that indeed this equation is correct, just do implicit differentiation. WAIT! Didn’t we do
that in the beginning? Yes. Go back and check that this is indeed what we started with (except for
that C, which appears as the result of integration.) This somewhat convoluted approach is just a
way to say that solving exact differential equation is basically undoing the implicit differentiation.
I have not defined what exact differential equation is yet. It will be given in the next section.
28
Yet another introduction to Exact differential equations
(2xy)dx + (x 2 − 4y)dy = 0
This is an example of an exact differential equation. However, we haven’t defined what exact
differential equation is. That will be given after the solution of this problem. A little note about
the strange looking symbols dx and dy: If you divide both sides by dx, you get:
dy
(2xy) + (x 2 − 4y) =0
dx
dy
This looks more natural because dx has the meaning of the derivative of y, where as dx and dy
we had before is hard to understand. (Some calculus courses teach dx and dy as differentials.)
However, many differential equation textbooks present problems in the form of (2xy)dx + (x 2 −
4y)dy = 0, because in this form, it is easier to recognize whether the given differential equation
is exact.
To understand how to solve this equation, we need to review the multivariable chain rule for
the following dependence diagram:
In other words, we have a function F (x, y) and both x and y are functions of t. The chain rule
associated with the dependence diagram above is:
dF ∂F dx ∂F dy
= +
dt ∂x dt ∂y dt
If we multiply dt both sides, we get
dF = Fx dx + Fydy
29
In other words, if you take a function and differentiate it by x first and then y, you will get the
same answer as if you differentiated by y first and then x. Since we are claiming that Fx = 2xy
and Fy = (x 2 − 4y) in our question, we need to verify whether Fxy = Fyx does hold true or not.
Now, Fxy = ∂y Fx = ∂y (2xy) = 2x. Also, Fyx = ∂x Fy = ∂x (x 2 − 4y) = 2x. So we find that indeed
Fxy = Fyx does hold. What is surprising is that mathematician Henry Poincare proved that this is
the only condition one needs to verify in order to find out whether F (x, y) exists or not. (Actually,
there are some more delicate conditions in addition to this, but we will never have to deal with
other conditions.) Here is the statement:
Then if Py = Qx , we say that the differential equation is exact, or that the above is an exact
differential equation.
Once we verify that the differential equation is exact, then just remembering that Fx is the
one in front of dx and Fy is in front of dy will guide us to find F (x, y). In our example,
∫
Fx = 2xy → F = 2xydx = x 2y + д(y)
∫
2
Fy = (x − 4y) → F= (x 2 − 4y)dy = x 2y − 2y 2 + h(x)
F (x, y) = x 2y − 2y 2
x 2y − 2y 2 = C
3 2
2
(3x cos y + 5)dx + − x sin y + 3y dy = 0
2
Solution) Step 1. Check for exactness. The formula dF = Fx dx +Fydy tells us that Fx = (3x cos y+5)
and Fy = (− 32 x 2 sin y + 3y 2 ). We have to check whether Fxy = Fyx holds or not.
30
3 2
2
Fyx = ∂x Fy = ∂x − x sin y + 3y = −3x sin y
2
Thus we see that it is exact. Step 2. Integrate and figure out F (x, y)
3
∫
Fx = (3x cos y + 5) → F = (3x cos y + 5)dx = x 2 cos y + 5x + д(y)
2
3 2 3 2 3
∫
2
Fy = − x sin y + 3y → F= − x sin y + 3y dy = x 2 cos y + y 3 + h(x)
2
2 2 2
Comparing the two, we conclude that д(y) = y 3 and h(x) = 5x which means:
3
F = x 2 cos y + y 3 + 5x
2
Step 3. Set the function equal to C The answer is :
3 2
x cos y + y 3 + 5x = C
2
Example 1.4.4
5x 2 + 2y + 2xy 0 = 0
31
Then ∂d y(5x 2 + 2y) = 2 and ∂x 2x = 2. Therefore this differential equation is exact. This means
that there is a function F (x, y) which satisfies Fx = (5x 2 + 2y) and Fy = 2x.
5
∫
Fx = (5x + 2y) → F = (5x 2 + 2y)dx = x 3 + 2xy + д(y)
2
3
∫
Fy = (2x) → F = (2x)dy = 2xy + h(x)
Answer:
5
2xy + x 3 = C
3
If you solve this for y, you will see that both answers are the same.
∫
xy
Fy = (2y + xe ) → F= (2y + xe xy )dy = y 2 + e xy + h(x)
32
Higher Order Linear Differential Equations
2
2.1 Independence of Solutions and Wronskians
Introduction to Higher Order Linear Differential Equations
Order of a differential equations is the highest derivative that appears in the equation. Here are
some examples:
Higher order differential equations (i.e. equations whose order is 2 or above) are much harder
than first order differential equations. Therefore, we take a different approach than what we have
done for first order differential equations.
33
What justifies ”guessing”?
Thankfully, guessing the solutions of a linear differential equation is justified because mathemati-
cians proved two important theorems about it. They are:
• Existence of solutions
• Uniqueness of the solution given appropriate initial conditions
Existence allows you do go ahead and hunt for a solution. If we don’t know whether the solution
of a differential equation exists, then we may end up hunting for the solution forever. (Or wrongly
think that a false solution you found is a solution) Uniqueness allows you to say that the solution
you found is the only one, and thus you don’t have to find more.
y 00 = y
I know that y = e x is a function whose derivative is same as itself, so it would satisfy the above
equation. Then after trying many different functions, we find that y = e −x is another solution.
(Quick check: y 0 = −e −x . y 00 = e −x = y.) Whenever we find two different solutions like these,
take the linear combination of the two functions as:
y = C 1e x + C 2e −x
Then you can check that this is a solution to y 00 = y. This is a good thing, because we managed to
have two free parameters for our second order differential equation, which means we have found
the most general solution to this differential equation.
Note: 1. Later, we will learn what ’homogeneous’ means and that linear homogeneous differ-
ential equations always have the general solution in the form of linear combination of solutions.
2. Physicists often use the word superposition instead of linear combination.
With the success of the first question, let’s move on to a slightly harder question:
y 00 = −y
34
You may remember that sin x has second derivative that is negative of itself. So is cos x. So we
can quickly write down their linear combination and claim that the general solution is:
y = C 1 cos x + C 2 sin x
Did you notice that y = e ix is also a solution to y 00 = −y? This means that the linear combination
of three functions is also a solution:
y = C 1 cos x + C 2 sin x + C 3e ix
But wait, this has three parameters. So something is wrong here. Thinking about the above linear
combination for a while, you may come to question whether one of the functions is really other
functions written in disguise. And indeed, it is true. There is something called the Euler identity,
that says:
e ix = cos x + i sin x
So really what we thought was a new solution was in fact expressible as linear combination of the
other two functions. There is a language for such case. Since e ix was expressed as linear combina-
tion of sin x and cos x, we say ”e ix is NOT linearly independent from sin x and cos x”. Still another
way to say it is ”e ix is linearly independent on sin x and cos x”. The upshot of this discussion
is that whenever we have a set of functions that satisfy a differential equation, we should check
whether those functions are linearly independent before we make linear combination out of them
to make the general solution.
Definition 2.1.1 Functions f (x) and д(x) are said to be linearly dependent if there is any pair of
constants C 1 and C 2 (at least one of them non-zero) that satisfies
C 1 f (x) + C 2д(x) = 0.
Functions f (x), д(x) and h(x) are said to be linearly dependent if there is any triple of constants
C 1 ,C 2 and C 3 (at least one of them non-zero) that satisfies
In our previous example, the three functions e ix , sin x and cos x are linearly dependent because
35
Can this express all general solutions of y 00 = y? The answer is no because although above linear
combination looks like it has two free parameters, it really has just one. To see this, notice that:
y = C 1e x + C 2 (2e x ) = (C 1 + 2C 2 )e x
Since C 1 and C 2 are arbitrary, we can set C = (C 1 + 2C 2 ) and rewrite the above as:
y = Ce x
In other words, the linear combination only produces solutions that looks like multiples of e x ,
and thus it really has just one free parameter.
Actually, you should have noticed even from the beginning that f (x) = e x and д(x) = 2e x are
linearly dependent, because one function is a multiple of the other. Since
д(x) = 2f (x),
Example 2.1.1 Show that sin x and cos x are linearly independent.
Solution) We will prove by contradiction. Suppose they are linearly dependent. Then one is the
multiple of the other, so we can write:
sin x = C cos x
for some constant C. Now, if we divide both sides by cos x, then we have :
sin x
=C
cos x
tan x = C
But this last line is a contradiction because tan x is not a constant function. Therefore, it is im-
possible to find such C and therefore the two functions must be linearly independent.
This example shows that indeed the general solution of y 00 = −y is y = C 1 cos x + C 2 sin x
because we have established that cos x and sin x are linearly independent of each other.
Wronkian of functions
Wronskian is a determinant introduced by Jozef Hoene-Wronski in 1776. In the previous section,
we saw how cumbersome it is to prove that two functions are linearly independent. Wronskian
makes it easy to determine linear independence of multiple functions. To understand the idea
behind this, let us write the condition for linear dependence, using sin x and cos x:
C 1 cos x + C 2 sin x = 0
36
Now, the above is an equation of functions, which means that it is supposed to be true for any
value of x. Sometimes we call this (functional) identity, rather than an equation. There is a big
difference between algebraic equation and functional equality. If two sides are equal as functions,
their derivatives are equal to, e.g.
x 2 + 2x + 1 = (x + 1)2
is an identity because it holds for all values of x. Now, if we differentiate both sides, we get:
2x + 2 = 2(x + 1)
where the right side was computed using the chain rule. Notice that indeed the two sides are
same. On the other hand, think about the following equation:
x 2 = 2x
This equation is only true for x = 0 and x = 2. So it is not an identity. In this case, the derivatives
do not match:
2x , 2
Since C 1 cos x + C 2 sin x = 0 is an identity, we CAN differentiate it:
C 1 (− sin x) + C 2 cos x = 0
So saying that C 1 cos x + C 2 sin x = 0 holds leads to the following system of equations:
C 1 cos x + C 2 sin x = 0
C 1 (− sin x) + C 2 cos x = 0
You can view this as system of equations for C 1 and C 2 . Recall that you can write system of
equations as matrices. In our case
cos x sin x C1 0
=
− sin x cos x C2 0
Such equation can be solved by multiplying the inverse matrix both sides.
−1 −1
cos x sin x cos x sin x C1 cos x sin x 0
=
− sin x cos x − sin x cos x C2 − sin x cos x 0
In the left side, the inverse cancels the next matrix. On the right side, any matrix multiplied by a
zero matrix will produce zero matrix. Thus we have:
C1 0
=
C2 0
Now this means that the only pair of constants C 1 and C 2 that satisfies C 1 cos x + C 2 sin x = 0 is
C 1 = 0 and C 2 = 0. Which means that the two functions cos x and sin x are linearly independent!
37
(Go back to the previous section and read the definition of linear dependence if you don’t under-
stand why.) But wait, where did we use anything specific about cos x and sin x here? There is
one thing I forgot to mention. Inverse matrices don’t always exist. This calculation can only be
possible if the matrix
cos x sin x
− sin x cos x
has an inverse. Recall that the inverse of a matrix exists only when the determinant is nonzero. So
in order to justify our calculation, we have to compute the determinant as follows:
cos x sin x
2 2
− sin x cos x = cos x · cos x − sin x · (− sin x) = cos x + sin x = 1 , 0
Since the determinant is not zero, we know that the calculations were sound and the conclusion
that the two functions are linearly independent is correct. Now, for any two functions f (x) and
д(x), we can follow the exact same steps as above. You will find that at the end of the day, whether
or not the two functions are linearly independent depends on the following determinant:
f (x) д(x)
f (x) д0(x)
0
This matrix is called the Wronskian of two functions, and is denoted as W [f (x); д(x)]. If the
Wronskian is zero, then the two functions are linearly dependent, but if the Wronskian is non-
zero the two functions are linearly independent. (There is some requirement about the functions
being analytic, but we will avoid these technicalities.)
Notice that this time we had to involve the second derivative too. This is because determinants
are only defined for square matrices and by using the second derivative we can add one more
row to the matrix to make it a 3 × 3 matrix.
Let’s revisit some of the examples we investigated in the previous section:
Solution) x
e 2e x
x x = e x · 2e x − 2e x · e x = 0
W [e ; 2e ] = x
e 2e x
Thus the two functions are linearly dependent.
38
Solution)
cos x sin x e ix
= sin x(− sin xe ix + sin xe ix ) + cos x(cos xe ix − cos xe ix ) + e ix (cos x sin x − cos x sin x) = 0
Hence the three functions are linearly dependent.
Linear Operators
An operator is a function which takes a function as an input and gives you back another function
as an output. Derivative is an example of an operator. Partial derivatives are also operators. So
is a constant multiples. A linear operator is an operator that satisfies the following property:
In other words, the operator L has the property where it can be applied term by term and also
constant multiple can be brought in front.
Here is an example of a linear operator:
L = D2 + 5
where D is the operator that differentiates a function by x and 5 is the operator that multiplies
by 5. Here is how it operates on a function:
L = 5D 2 − 6xD + (e x − x)
is a linear operator. (I’ll leave you to check that it is indeed true, since it’s rather long.)
39
Why is linear operators important? It is important for us because linear differential equations
always contain a linear operator in it. For example, the differential equation:
y 00 − 3y 0 − 4y = 0
Since e rx can never be zero, the only way the above can be = 0 is when (r 2 − 3r − 4) = 0. Factoring
this as (r − 4)(x + 1) = 0, we find that r = −1, 4. This means that there are two solutions to this
differential equation: y = e −x and y = e 4x . You can also compute their Wronskian and find that
these two functions are linearly independent.
If we use the language of operators, y = e 4x being the solution to y 00 − 3y 0 − 4y = 0 is same
thing as saying L e 4x = 0, with L = D 2 − 3D − 4. Likewise L (e −x ) = 0. Here is why linearity is
important:
L c 1e 4x + c 2e −x = c 1L e 4x + c 2L (e −x ) = c 1 · 0 + c 2 · 0 = 0
Look at how the linearity of the operator was used to show that the linear combination of two
functions is again a zero if L is applied, i.e. c 1e 4x +c 2e −x is the general solution of y 00 − 3y 0 − 4y = 0
Conclusion
In other words, if your differential equation can be written in a form of Ly = 0 with a linear
operator L, then any linear combination of solutions is again a solution, due to the linear property
of the operator L.
ar 2 + br + c = 0
40
The above equation is called the characteristic equation of the differential equation ay 00 +by 0 +
cy = 0. Whenever we solve a differential equation of the form ay 00 + by 0 + cy = 0, we will imme-
diately write down the characteristic equation and solve it. The name ”characteristic equation”
comes from the fact that depending on whether the quadratic equation has two distinct roots, a
double root or complex roots, the solution of the differential equation takes different forms.
Let’s think about the solution of ar 2 + br + c = 0. Using the quadratic formula:
√
−b ± b 2 − 4ac
r=
2a
There are three different cases:
• When b 2 − 4ac > 0, the solutions are two real values r 1 and r 2 . This means e r 1x and e r 2x are
solutions. Then we take the linear combination of the two functions as the general solution:
y = C 1e r 1 x + C 2e r 2 x
(One should check that the two solutions are linearly independent by computing the Wron-
skian.)
• When b 2 − 4ac < 0, the quadratic formula produces r = α ± βi. In this case the general
solution is
y = e αx (C 1 cos βx + C 2 sin βx)
• When b 2 − 4ac = 0, then the quadratic formula produces just one solution r = r 1 . Hence we
only have one solution e r 1x . This is bad because we need two linearly independent solutions
in order to write the general solution. Thankfully, someone figured out that in this case,
xe r 1x also becomes one solution. So the general solution is:
y = C 1e r 1x + C 2xe r 1x
r 2 − 6r + 8 = 0
(r − 4)(r − 2) = 0
r = 4, 2
Hence the general solution is:
y = C 1e 4x + C 2e 2x
41
Now we want the specific solution that satisfies the initial condition y(0) = 1, y 0(0) = −6. Since
y 0 = C 1 · 4e 4x + C 2 · 2e 2x , we have:
−6 = y 0(0) = C 1 · 4e 4·0 + C 2 · 2e 2·0 = 4C 1 + 2C 2
Also:
1 = y(0) = C 1e 4·0 + C 2e 2·0 = C 1 + C 2
So we have a system of equations:
4C 1 + 2C 2 = −6
C1 + C2 = 1
Solving C 1 + C 2 = 1 for C 2 as C 2 = 1 − C 1 and then plugging it into 4C 1 + 2C 2 = −6, we get:
4C 1 + 2(1 − C 1 ) = −6
2C 1 + 2 = −6
C 1 = −4
C 2 = 1 − (−4) = 5
After finding the constants, put it back into the general solution:
y = (−4)e 4x + 5e 2x
42
Second order constant coefficient homogeneous linear differential equa-
tions: two complex roots
Example 2.2.3 Solve
y 00 − 6y 0 + 10y = 0, y(0) = 3, y 0(0) = −2
r 2 − 6r + 10 = 0
You can think of 3 ± i = (3) ± (1)i, so that α = 3 and β = 1 in the above formula. Thus the
general solution is:
y = e 3x (C 1 cos x + C 2 sin x)
The derivative is:
using the product rule. Using y(0) = 3, y 0(0) = −2 and the fact that sin 0 = 0 and cos 0 = 1 we
have:
3 = y(0) = C 1
−2 = y 0(0) = 3C 1 + C 2
Plugging in C 1 = 3 into the second yields C 2 = −11.
So the answer is:
y = e 3x (3 cos x − 11 sin x)
43
Second order constant coefficient homogeneous linear differential equa-
tions: Factor by grouping
Example 2.2.4 Solve
3y 00 − 50 − 2y = 0, y(0) = 1, y 0(0) = −2
Solution) The characteristic equation is:
3r 2 − 5r − 2 = 0
This can be factored using ’factor by grouping’. First, multiply the first number and the last
number (which are 3 and −2) to get −6. Then think of two numbers that multiply to −6 and also
adds to the middle number −5. We find that -6 and 1 satisfies this property. Now split the middle
number as the sum of the two like this:
3r 2 + (−6 + 1)r − 2 = 0
3r 2 − 6r + r − 2 = 0
Then you factor by grouping the first two terms and the rest:
(3r 2 − 6r ) + (r − 2) = 0
3r (r − 2) + (r − 2) = 0
(3r + 1)(r − 2) = 0
Therefore the solution of the characteristic equation is
1
r =− , 2
3
Hence the general solution is:
1
y = C 1e − 3 x + C 2e 2x
Now we want the specific solution that satisfies the initial condition y(0) = 1, y 0(0) = −2. Since
1
y 0 = − 31 C 1e − 3 x + 2C 2e 2x , we have:
1 1 1
−2 = y 0(0) = − C 1e − 3 ·0 + 2C 2e 2·0 = − C 1 + 2C 2
3 3
Also, 1 = y(0) = C 1 + C 2 . So we have a system of equations:
1
− C 1 + 2C 2 = −2
3
C1 + C2 = 1
Multiply 3 to − 13 C 1 + 2C 2 = −2 to get −C 1 + 6C 2 = −6 then add this to C 1 + C 2 = 1. Then in the
left side, −C 1 and C 1 will cancel leaving you with 7C 2 = −5. Thus C 2 = − 75 . Plugging this into
C 1 + C 2 = 1 yields:
5
C1 − = 1
7
5 12
C1 = 1 + =
7 7
Thus the answer is:
12 1 5
y = e − 3 x − e 2x
7 7
44
2.3 Method of Undetermined Coefficients
Introduction to the Method of Undetermined Coefficients
We will slowly learn this method by looking into various examples.
Example 2.3.1
y 00 − 2y 0 − 3y = 5e 7x
Solution) Since the right side has 5e 7x , we can predict that if we plug in Ae 7x with some constant
A in the left side, it will produce the 5e 7x on the right. Method of Undetermined Coefficients is
basically the method to plug in Ae 7x into y to figure out what A is. The word ”Undetermined
Coefficients” is about the undetermined A which we figure out by plugging it into y. Now,
yp = Ae 7x
yp0 = A · 7e 7x
yp00 = A · 49e 7x
and then:
yp00 − 2yp0 − 3yp = 49Ae 7x − 14Ae 7x − 3Ae 7x
= 32Ae 7x
5
Thus in order for this to equal to 5e 7x , we need 32A = 5. Solving it gives us A = 32 , and we have
5 7x
yp = 32 e .
After finding the particular solution, we need the complementary solution to go with it. Write
the characteristic equation:
r 2 − 2r − 3 = 0
(r − 3)(r + 1) = 0 → r = 3, −1
yc = C 1e 3x + C 2e −x
So the final answer is:
5 7x
y = yc + yp = C 1e 3x + C 2e −x + e
32
Solution) If you try just yp = A sin 2x, then you’ll notice that the middle term on the left, −2y 0 will
produce a cos 2x with some multiple. Since there is no cos 2x on the right, we see this is going to
cause trouble. To remedy this, you want to include cos 2x into the yp like the following:
yp = A sin 2x + B cos 2x
45
Hopefully, there is a combination of A and B that will make any cos term disappear and produce
exactly the 4 sin 2x when we plug it into the differential equation. Differentiate this twice so we
can plug it into the left side of the differential equation.
yp00 − 2yp0 − 3y
−7A + 4B = 4
−4A − 7B = 0
Solving the second for B, we have B = − 74 A. Plugging this into the first, you get −7A − 16
7A = 4.
49 + 16
− A=4
7
28
A=−
65
Then B = (− 47 ) · (− 28
65 ) =
16
65 . So
28 16
yp = − cos 2x + sin 2x
65 65
Thus
28 16
y = yc + yp = C 1e 3x + C 2e −x − cos 2x + sin 2x
65 65
r 2 − 6r + 9 = 0
(r − 3)(r − 3) = 0
r = 3, 3
yc = C 1e 3x + C 2xe 3x
46
Now let’s think about yp . At first, you might think plugging in yp = Ae 3x + Be 4x will produce
the right side: 3e 3x − e 4x . However, if you plug in Ae 3x into y, you’ll get zero on the left side.
Why? That’s because this term is actually included in the complementary solution. When this
happens, we say that there is a duplication. Being a complementary solution means it’s a solution
to the homogeneous equation, which means the right side is zero. So what do you do? When the
particular solution has any overlapping terms with the complementary, multiply x to those terms
until there is no overlap. In the above case, yp = Axe 3x + Be 4x doesn’t work because of duplication.
(Notice that xe 3x is included in the complementary solution yc ). So you have to multiply another
x to it, which gives you yp = Ax 2e 3x + Be 4x . Now, we differentiate this twice in order to plug it
into the differential equation.
= A(2x + 3x 2 )e 3x + 4Be 4x
yp00 = A(2 + 6x)e 3x + A(6x + 9x 2 )e 3x + 16Be 4x
= A(2 + 12x + 9x 2 )e 3x + 16Be 4x
Plugging these into the differential equation, we have:
= 2Ae 3x + Be 4x
3
For this to equal to 3e 3x − e 4x , we need A = 2 and B = −1. Therefore,
3
yp = x 2e 3x − e 4x
2
Thus the final answer is:
3
y = yc + yp = C 1e 3x + C 2xe 3x + x 2e 3x − e 4x
2
step 1. Find yc . (You do this first, so that you can recognize any duplication if it happens.)
step 2. Set yp = A( ) + B( ) +C( ) + · · · · · · where the blanks are filled with functions related
to the functions f (x), д(x), h(x), · · · appearing on the right side of the differential equation. In this
list, include all the f (x), д(x), h(x), · · · and their higher derivatives (ignoring coefficients). Here
are some examples:
47
• If the right side is e 3x , use yp = Ae 3x .
• If the right side is sin 3x, use yp = A sin 3x + B cos 3x. (You need cos 3x because it is the
derivative of sin 3x, ignoring the coefficients.)
• If the right side is x 3 , use yp = Ax 3 + Bx 2 + Cx + D.
• If the right side is xe 2x cos 4x, use y = Axe 2x cos 4x +Be 2x cos 4x +Cxe 2x sin 4x +De 2x sin 4x
• If the right side is 3 sin x + cos x, use yp = A sin x + B cos x.
step 3. Check for duplication, and multiply by x to any term that is duplicated. Keep mul-
tiplying by x until you have just avoided duplication. (Warning: Don’t multiply more than you
have to.) For example, let’s say you have xe 2x cos 4x in the complementary solution. and let’s say
the right side is xe 2x cos 4x + x 3 . Then first, you write
yp = Axe 2x cos 4x + Be 2x cos 4x + Cxe 2x sin 4x + De 2x sin 4x + Ex 3 + Fx 2 + Gx + H
However, you see that xe 2x cos 4x is duplicated. Then all the functions that are related to this
function by derivatives must be multiplied by x. So the proper form of yp is:
yp =x Axe 2x cos 4x + Be 2x cos 4x + Cxe 2x sin 4x + De 2x sin 4x
+ Ex 3 + Fx 2 + Gx + H
step 4. Find all coefficients A, B, · · ·
step 5. Then you’ll have your yp , which means y = yc + yp is your final answer.
48
Find the appropriate form of a particular solution when using the method
of undetermined coefficients. Second example.
Example 2.3.5
y 00 + 4y 0 + 4y = 7x 2e −2x + 5e 2x − sin x
a) Find yc
b) What is the form of yp appropriate for this differential equation if using the method of undeter-
mined coefficients?
r 2 + 4r + 4 = 0
(r + 2)(r + 2) = 0
r = −2, −2
Thus yc = C 1e −2x + C 2xe −2x
Now, for 7x 2e −2x we need (Ax 2 + Bx + C)e −2x , for 5e 2x we need De 2x and for − sin x we need
E sin x + F cos x. Hence we have:
However, if you expand x(Ax 2 + Bx + C)e −2x , it’s Ax 3e −2x + Bx 2e −2x + Cxe −2x and that last term
is still duplicated in the complementary solution. Therefore, another x should be multiplied:
Solution) This can be done with integration by parts twice. But it’s hard! Here is a completely
new way to look at it.
Think of the direct integration problem:
y 0 = e 2x cos 3x
This question is equivalent to the integration problem above. However, we now have a different
method to solve this because we can use the method of undetermined coefficients.
So put yp = Ae 2x cos 3x + Be 2x sin 3x Then
49
So in order for this to equal to e 2x cos 3x, we need: 2A + 3B = 1 −3A + 2B = 0 Now solving this
2
system gives us A = −313 and B = 13 . Plugging this back into yp , we see that
−3 2x 2
yp = e cos 3x + e 2x sin 3x
13 13
Just putting the C (which by the way is the complementary solution to the above direct integral
differential equation) gives us the answer:
−3 2x 2
e cos 3x + e 2x sin 3x + C
13 13
Solution) This is a question that cannot be solved using method of undetermined coefficients,
because if you think about all the derivatives of tan x, it never ends: (tan x)0 = sec2 x, (sec2 x)0 =
2 sec2 x tan x, · · · . This means you would need infinitely many terms for yp when you prepare for
the method of undetermined coefficients. Hence that approach will fail. So we have to do another
approach, which is called the Variation of parameters method.
Variation of parameters method also requires one to start with the complementary solution.
Write the characteristic equation: r 2 + 1 = 0 r = ±i Hence yc = C 1 cos x + C 2 sin x.
The variation of parameters method is a method where we change the parameters C 1 and
C 2 into functions of x, like the following:
Now, because cos x and sin x are complementary solutions, when you plug in this form of yp into
y, many terms on the left will cancel. That’s the basic idea of variation of parameters.
Let’s differentiate:
Now, this is already complicated. So we put another additional assumption that part of the above
is zero, namely u 10 (x) cos x + u 20 (x) sin x = 0. Then yp0 becomes much simpler:
Why is the assumption u 10 (x) cos x + u 20 (x) sin x = 0 okay to have? Well, if you plug in yp =
u 1 (x) cos x +u 2 (x) sin x into the given differential equation, you get a single equation involving two
unknown functions u 1 (x) and u 2 (x). This means another equation must be provided in order to
find out exactly what u 1 (x) and u 2 (x) are. So we make the assumption u 10 (x) cos x +u 20 (x) sin x = 0
to provide us this extra equation. (Remember, solving differential equations is basically an art of
guessing. So if this approach works, it works.) Now let’s differentiate the simplified yp0 to get:
50
Plugging in, we get:
yp00 + yp = −u 10 (x) sin x + u 20 (x) cos x,
which should equal to tan x because that’s what’s on the right side of our differential equation.
Hence we have the following system of equations:
(Since we are looking for a particular solution, we don’t need C in the integration.) Then we have:
51
Formula for Variation of Parameters
Suppose we are given a differential equation:
ay 00 + by 0 + cy = f (x)
Also let y1 and y2 be linearly independent solutions of the homogeneous equation ay 00+by 0+cy = 0
Then ’Variation of Parameters Method’ leads to:
y2 f (x) y1 f (x)
∫ ∫
yp = −y1 dx + y2 dx,
W W
where
y1 y2
W = 0 0
y1 y2
which we know as the Wronskian of y1 and y2 . The derivation of the above formula is straight-
forward but long, so we will skip it. Here is an example showing how to use the formula.
y 00 − y = x
r 2 − 1 = 0 → r = 1, −1
y = yc + yp = C 1e x + C 2e −x − x
52
2.5 Mechanical Vibrations
Spring-Mass-Dashpot System
Suppose a spring with spring constant k is attached to the wall and a block of mass m, where the
block is also attached to the dashpot with friction coefficient c. (See picture below.)
We will assume that the floor is frictionless and the only friction force is caused by the dashpot.
Using Newton’s law F = ma we aim to derive an equation that governs the motion of the
Í
block. In the picture, there are two kinds of forces: force of the spring and the friction force of
the dashpot. Hence we have:
Fs + Fd = ma
For the spring, there is the Hooke’s law, which states that
Fs = −kx
where x is the displacement of the block from the equilibrium position. Notice that the − sign
is due to the fact that the spring always inserts force opposite of the displacement direction. For
the dashpot, the friction force tries to counter movement of the block. For example, if v = 0 then
there is no movement and therefore the dashpot doesn’t exert any force. When v is large, i.e.
if the block tries move fast, the dashpot will exert a huge friction force in order to counter the
velocity. This fact can be written as the following equation
Fd = −cv
where the negative sign means that the friction always tries to counter the direction of motion.
If we replace Fs and Fd with the equations we found, we get:
−kx − cv = ma
mx 00 + cx 0 + kx = 0
53
This is the differential equation for the motion of the block. When describing vibration of an
object, often one can find an analogue of a spring-mass-dashpot system. For example, the me-
chanical vibration of a car as it runs on the road has the mass as the mass of the car and the
struts as springs (and various mechanical friction as dashpot). So the above equation has a quite
universal application for describing mechanical vibration.
FE + Fs + Fd = ma
FE − kx − cv = ma
FE = mx 00 + cx 0 + kx
This last equation is the differential equation to use when there is an external force acting on the
block. Usually FE is not constant but a force that varies with time. So one writes FE (t) and we
have:
FE (t) = mx 00 + cx 0 + kx
54
Undamped mechanical vibration
The master equation for mechanical vibration is:
F (t) = mx 00 + cx 0 + kx,
where F (t) is the external force, m is the mass, c is the friction coefficient and k is the spring
constant. Undamped motion is when there is no friction in the system, i.e. c = 0. Let’s investigate
when there is no external force, so that we have:
mx 00 + kx = 0
To solve this, first write the characteristic equation
mr 2 + k = 0
k
r2 = −
m
r
k
r =± −
m
r
k
r =± i
m
If we set r
k
ω= ,
m
then r = ±ωi, and the solution is:
x(t) = C 1 cos ωt + C 2 sin ωt
Note that both cos ωt and sin ωt are periodic functions with period
2π
T = .
ω
(One way to understand the above formula is to think about the special case when ω = 1. Then
cos ωt = cos t which we know to have the period of 2π .) Now, since frequency f = T1 , we can
further write:
1 2π
T = = .
f ω
Another useful formula you can obtain from above is
ω = 2π f
This formula tells us what ω to choose when the q frequency is given. ω is called the angular
frequency. Solving this for f and also using ω = mk , we have
r
1 k
f =
2π m
Above formula states that the frequency of the mechanical vibration depends on√the mass and
spring constant. √More specifically, frequency is inversely proportional to the m while it is
proportional to k. This can be easily understood if you think about a heavy object (one with
big value for m), for it will move slowly, causing the frequency to go down.
55
Putting C 1 cos ωt + C 2 sin ωt into a single cosine function
Graphs of A sin(at +b) and A cos(at +b) are called ’sinusoidal’ because they have the same shape
as the graph of sin t. When two sinusoidal functions of same angular frequency are added, the
resulting function is again a sinusoidal function with the same angular frequency.
Say we have a sum of a sine and a cosine function of the same angular frequency, i.e.
A cos ωt + B sin ωt
This way of writing makes it easy to see that the result of the sum of sinusoidal functions of
the same period is again a sinusoidal function. It also tells you the amplitude C of the sum and
also the phase angle α. So when writing your answer to mechanical vibrations, always make it
into a single cosine function using the above formula.
Example 2.5.1 Suppose a block of 4 kg is attached to a spring with spring constant 1N/m. There is
no external force or friction in the system. If the spring was compressed by 2m and was pushed with
1m/s initial velocity to the right, find period, amplitude and phase angle of the motion.
4x 00 + 1 · x = 0
4r 2 + 1 = 0
1
r2 = −
4
1
r =± i
2
1 1
x(t) = C 1 cos t + C 2 sin t
2 2
−2 = x(0) = C 1 · 1 + C 2 · 0
56
C 1 = −2
1 1 1 1
x 0(t) = − C 1 sin t + C 2 cos t
2 2 2 2
1
1 = x 0(0) = C 2
2
C2 = 2
p √ √
C = (−2)2 + 22 = 8 = 2 2
−1 2 π 3π
α = tan +π = − +π =
−2 4 4
√ 1 3π
x = 2 2 cos t −
2 4
√ 3π
Answer: Amplitude is 2 2. Phase angle is 4 . Period is:
2π
T = = 4π
1/2
Solution)
F (t) = mx 00 + cx 0 + kx,
with m = 1/2 and c = 1. Also, by the Hooke’s law F = kx, we can write 100 = k · 2 and we get
k = 50. There is no external force indicated, so we have:
1 00
x + x 0 + 50x = 0
2
with x(0) = 1, v(0) = x 0(0) = −5
1 2
r + r + 50 = 0
2
r 2 + 2r + 100 = 0
(r + 1)2 + 99 = 0
√
r = −1 ± 99i
Hence √ √
x(t) = e −t (C 1 cos 99t + C 2 sin 99t)
1 = x(0) = C 1
√ √
x 0(t) = −e −t (C 1 cos 99t + C 2 sin 99t)
57
√ √ √ √
+e −t (− 99C 1 sin 99t + 99C 2 cos 99t)
√
−5 = x 0(0) = −C 1 + C 2 99
Hence C 2 = −4
√ Thus:
99
√ 4 √
x(t) = e −t (cos
99t − √ C 2 sin 99t)
99
Now we have to put the thing inside the parenthesis as a single cosine function:
s 2 s
115
−4
C = 12 + √ = √
99 99
−4
α = tan −1
√ = −0.3822
99
Hence, s
115 √
x(t) = √ e −t cos( 99t − (−0.3822))
99
Note that this is not a periodic function because√e −t multiplied to above means that the amplitude
is decreasing. However, we still call the ω = 99 as the pseudo-angular frequency and use
that to compute the pseudoperiod:
2π
T = √
99
Also, what the question calls the ’time lag’ is what is usally called the phase shift:
α
δ=
ω
(in other words, phase angle divided by the angular frequency gives you the phase shift) So the
time lag calculated is =-0.038 (seconds). Finally, in order to know when the block crosses the
equilibrium, one uses the fact that cosine is zero at π /2, 3π /2, 5π /2, · · · . Therefore, the first four
times block crosses the equilibrium can be calculated when t satisfies:
√
99t + 0.3822 = π /2, 3π /2, 5π /2, 7π /2
Solution)
1 · x 00 + 4x 0 + 8x = 3 cos(5t)
r 2 + 4r + 8 = 0
58
(r + 2)2 + 4 = 0
r = −2 ± 2i
xc = e −2t (C 1 cos 2t + C 2 sin 2t)
We use the method of undetermined coefficients to find xp :
xp = A cos 5t + B sin 5t
mx 00 + kx = F 0 cos ωt .
Now there is the natural angular frequency of the spring-mass system, given by
r
k
ω0 = .
m
(Notice that heavier the block is, slower it will move, and thus lower frequency it will have.)
On the other hand, F 0 cos ωt is the external force with angular frequency ω. When ω = ω 0 , i.e.
when the external force’s angular frequency matches the natural frequency of the system, we say
59
that there is resonance in the system. Let’s look at two examples comparing with and without
resonance:
Case of no resonance: x + 4x = sin t
Since m = 1 and k = 4, the natural angular frequency is ω0 = 4/1 = 2. Since the external
p
force sin t = sin 1 · t has angular frequency 1, it’s angular frequency is different from the natural
angular frequency. So we know that this is a case of no resonance.
Let us solve this using method of undetermined coefficients. Plugging in x + p = A sin t on
the left side, we find that xp = 31 sin t is a particular solution. Since the characteristic equation
r 2 + 4 = 0 has the solution r = ±2i, the complementary solution is C 1 cos 2t + C 2 sin 2t. Thus the
general solution is x = xc + xp = C 1 cos 2t + C 2 sin 2t + 13 sin t. Important aspect of the solution is
that all the functions above are bounded functions.
Case of resonance: x 00 + 4x = sin 2t
This time, because there is duplication (when doing method of undetermined coefficients),
xp = t(A cos 2t + B sin 2t) must be used. Plugging this into the left side gives:
In other words, if there is resonance, the resulting motion grows without bound. In the real
world, this leads to bridges falling down, glasses breaking etc.
mx 00 + cx 0 + kx = F 0 cos ωt
Then we have learned that xc has a dying amplitude and thus is called the transient solution
(written as xt r ) while xp = C cos(ωt − α) is a periodic solution whose amplitude doesn’t get
60
smaller and thus is called the steady periodic solution (written as xs p). Using the method of
undetermined coefficients, you can find that
Using C = A2 + B 2 we get:
F0
C= p
(k − mω 2 )2 + (cω)2
If the value of ω is taken to make this amplitude to be maximized, we say that the external
frequency has been tuned to practical resonance. It is the frequency when xs p has the largest
amplitude, assuming F 0 is fixed. To maximize C, we need the denominator to be minimized,
which in turn is minimized when the radicand is minimized. Since minimum must occur at the
critical point, we differentiate with respect to ω:
(k − mω 2 )2 + (cω)2 = 0
0
2(k − mω 2 )(−2mω) + 2c 2ω = 0
Find the angular frequency of the external force which creates practical resonance for the
system with m=1, c=4 and k=10 solution)
F0 F0
C= p = p
(k − mω 2 )2 + (cω)2 (10 − ω 2 )2 + (4ω)2
Then
(10 − ω 2 )2 + (4ω)2 = 0
0
61
Now the voltage drop across an inductor is VL = −L dI dt because inductors try to counter the
changes in current. Voltage drop across a resistor is VR = −IR which is called the Ohm’s law.
Lastly, the voltage drop across a capacitor is VC = − C1 Q where Q is the charge. If the circuit is
provided with the voltage of E(t), then using the Kirchoff’s loop law (which states that the sum
of voltages around a loop is zero), we have the following:
dI 1
E(t) − L − IR − Q = 0
dt C
or
dI 1
E(t) = L + IR + Q
dt C
In the above, there are two unknown functions of time t. One is the current I (t) and the other is
the charge Q(t). There is a relation between the two, because the current is the rate of the change
of the charge, i.e. I (t) = Q 0(t). So if we differentiate the differential equation, we get:
1
E 0(t) = LI 00 + RI 0 + I
C
This equation is the master equation for LCR circuits.
Notice that this is a second order linear differential equation. You can compare this with the
mechanical vibrations equation and find an analogy between x(t) and I (t). In fact, what L is do-
ing is exactly what m is doing in mechanical vibration. R corresponds to c and C1 corresponds to
k. The E 0(t) acts like the external force. Therefore, everything you learned about the mechan-
ical vibration applies to the LCR circuits. We can talk about under-damped, over-damped, and
critically damped. We can talk about frequency, period and pseudo period. We can talk about
transient solutions and steady periodic solutions, etc. In LCR circuits, the phase angle is an im-
portant quantity because it tells you how much the current is lagging or leading the input voltage.
Note about the lagging/leading: If the current peaks a little before the voltage, then we say the
current ’leads’ the voltage. The graph of the current then would be shifted left, which means the
phase angle is negative. There could be a little confusion here because many books just write
C cos(ωt + ϕ) so that what we are calling as −α is instead called ϕ. In other words, if ϕ is positive,
so that α is negative, then the current ’leads’ the voltage.
Example 2.5.4 Suppose R= 30 Ω, L= 10 H, C = 0.02 F and E(t) = 50 sin 2t V. Find amplitude of the
current and also how much the current is lagging or leading the voltage.
62
and
60/37
α = tan−1 = 1.406 rad
10/37
So the current is lagging the voltage by 1.406 radians. The peak current is 1.644 amperes.
63
Series Solutions to Differential Equations
3
3.1 Power Series Solutions
Introduction to Series Solutions
Consider the differential equation:
y 00 − 3xy 0 + 2y = 0
This differential equation is quite different than what we have done so far, because there is a
coefficient −3x in front of y 0 which makes it non-constant. In the constant coefficient case, we
made a ’guess’ that the answer is in some format of e rx . However, when the coefficient is not
constant but rather a function of x, the solution quickly becomes very complicated and usually
ends up being a function which cannot be expressed using polynomials, trigs, exponentials, logs
or any other functions that we are used to. Thus we need to use a completely different approach
to this differential equation. The way we are going to solve it is by representing the solution as a
power series.
64
First, the equal sign holds for all values of x as long as the right side converges. This means,
for x = 1, we have:
1 1 1
e = e 1 = 1 + 1 + 12 + 13 + 14 + · · ·
2! 3! 4!
1 2 1 3 1 4
≈ 1 + 1 + 1 + 1 + 1 ≈ 2.70833
2! 3! 4!
In other words, if we took the first 5 terms of the series we get the above approximation. Con-
sidering that e = 2.7182818..., the above is a pretty good approximation. In fact, if you take more
terms, more accurate the approximation will be.
Second, you have to know how this power series is obtained. Suppose all that you know is
that e x can be written as a power series. Then you may write:
e x = a 0 + a 1x + a 2x 2 + a 3x 3 + · · ·
However, we don’t know what the coefficients a 0 , a 1 , a 2 , · · · are. So we have to figure out what
they are. The idea is to understand that the above equation is an equation for functions, which
means (a) it should hold for any value of x and (b) even if you differentiate both sides the two
sides are still equal. To get a 0 , we plug in x = 0 to get:
e 0 = a 0 + a 1 · 0 + a 2 · 02 + a 3 · 03 + · · ·
1 = a0
Then we differentiate both sides:
e x = a 1 + 2a 2x + 3a 3x 2 + 4a 4x 3 + · · ·
e 0 = a 1 + 2a 2 · 0 + 3a 3 · 02 + 4a 4 · 03 + · · ·
1 = a1
Differentiate more:
e x = 2a 2 + 3 · 2a 3x + 4 · 3a 4x 2 + 5 · 4a 5x 3 + · · ·
Again, plug in x = 0 to get:
e 0 = 2a 2 + 3 · 2a 3 · 0 + 4 · 3a 4 · 02 + 5 · 4a 5 · 03 + · · ·
1 = 2a 2
1
a2 =
2
We can continue this until there is a pattern. Taylor series (or Maclaurin series) is basically this
packaged as an easy-to-use formula.
65
Back to the differential equation
Once again we owe to the mathematicians who proved the existence of the solutions of differential
equations such as:
y 00 − 3xy 0 + 2y = 0
Because we know that the general solution exists, we can think about its Taylor series. In other
words, we take y to be:
∞
an x n
Õ
y=
n=0
n=1
∞ ∞ ∞
n−2 n−1
2an x n
Õ Õ Õ
= n(n − 1)an x + (−3x)nan x +
n=2 n=1 n=0
∞ ∞ ∞
n−2 n
2an x n
Õ Õ Õ
= n(n − 1)an x + −3nan x +
n=2 n=1 n=0
Shifting
Notice that in the term n=2 n(n−1)an x n−2 , the exponent of x n−2 is n−2, while for the other to ’s
Í∞ Í
they are n. For the reason that will be clear later, you want to match the exponents over x, and
that technique is called shifting of indices. Here is what it is: Take the n=2 n(n − 1)an x n−2 . Since
Í∞
we want the exponent to be n, we replace the n in x n−2 by n + 2, so that x n−2 becomes x n+2−2 = x n .
But when you do that, you have to change every n in the summation n=2 n(n − 1)an x n−2 , even
Í∞
including the n = 2 in the bottom.
n = 2 is transformed like this: Since we are replacing n by n + 2, it becomes n + 2 = 2. But then
if you solve for n, you have n = 0. After replacing n by n + 2, the summation n=2 n(n − 1)an x n−2
Í∞
66
becomes + 2)(n + 1)an+2x n . Now this shifting is complete. After the shifting, this is what
Í∞
n=0 (n
we have:
∞ ∞ ∞
(n + 2)(n + 1)an+2x n + −3nan x n + 2an x n
Õ Õ Õ
y 00 − 3xy 0 + 2y =
n=0 n=1 n=0
There are three summations you see here. Summations are things that produce terms, e.g. n=0
Í∞
(n+
2)(n+1)an+2x n produces 2·1·a 2x 0 when n = 0. Then it produces 3·2·a 3x 1 when n = 1. So basically
you have three term-producing factories working. There is a difference in when they start their
work though. When n = 0, only the first and the last factories (i.e. summations) produce terms.
But then, starting from n = 1 and afterwards, all factories are working to produce terms. Let’s
make use of this observation by separating out the terms produced when n = 0 as follows:
∞ ∞ ∞
= 2 · 1 · a 2x 0 + 2a 0x 0 + (n + 2)(n + 1)an+2x n + −3nan x n + 2an x n
Õ Õ Õ
Then the three sums can be put together as a single sum as:
∞
= 2 · 1 · a 2x 0 + 2a 0x 0 + ((n + 2)(n + 1)an+2x n + −3nan x n + 2an x n )
Õ
n=1
Now, this has to equal to zero, because y 00 − 3xy 0 + 2y = 0 that we are trying to solve has zero
on the right. When you have a power series organized in terms of power of x which has to equal
to zero, it is a mathematical theorem that each coefficient has to equal to zero (the reason being
that the Taylor series of zero is 0 + 0 · x + 0 · x 2 · · · .) and hence we have the following identities:
2 · 1 · a 2 + 2a 0 = 0
67
This is called the recurrence relation because it could be used over and over again. Let’s think
about the consequences of the recurrence relation. When n = 1,
1
a3 = a1
6
Then when n = 2,
4
a4 = a2
12
But since a 2 = −a 0 we have:
1
a4 = − a0
3
Then when n = 3,
7 7
a5 = a3 = a1
20 120
You can see that all the an can be represented by using either a 0 or a 1 (or sometimes both). This
is good because we know that at the end of the day there should be two free parameters. a 0 and
a 1 will end up being the free parameters.
Already we have some idea of the power series solution for y:
y = a 0 + a 1x + a 2x 2 + a 3x 3 + a 4x 4 + a 5x 5 + · · ·
1 1 7
= a 0 + a 1x − a 0x 2 + a 1x 3 − a 0x 4 + a 1x 5 + · · ·
6 3 120
The above is a significant result because a 0 and a 1 can be determined if the initial condition
is given. When they are determined, the above will tell us the few terms of the power series
expansion (a.k.a. Taylor series) of y, which can be used to find approximate value of y when x is
near 0.
x 2y 00 + 3y 00 + 4xy 0 + y = 0
68
(When doing power series solutions, the above is always used.) Plugging into the left side of
(x 2 + 3)y 00 + 4xy 0 + y = 0 we get:
∞ ∞ ∞ ∞
x2 n(n − 1)an x n−2 + 3 n(n − 1)an x n−2 + 4x nan x n−1 + an x n
Õ Õ Õ Õ
n=2 n=0
Hence we have:
∞ ∞ ∞ ∞
n(n − 1)an x n + 3(n + 2)(n + 1)an+2x n + 4nan x n + an x n
Õ Õ Õ Õ
We look at each summation and see that the second and fourth summations begin at n = 0, and
that starting from n = 2 and afterwards all summations produce terms. So we write:
∞
0 1
(n(n − 1)an + 3(n + 2)(n + 1)an+2 + 4nan + an ) x n
Õ
= (6a 2 + a 0 )x + (18a 3 + 4a 1 + a 1 )x +
n=0
Since the above has to equal to zero on the right side of (x 2 + 3)y 00 + 4xy 0 + y = 0, we obtain the
following equations:
6a 2 + a 0 = 0
18a 3 + 4a 1 + a 1 = 0
n(n − 1)an + 3(n + 2)(n + 1)an+2 + 4nan + an = 0
The first two equations yield a 2 = − 61 a 0 and a 3 = − 18
5
a 1 . The third equation, which holds for
n ≥ 2, solved for an+2 gives us:
n 2 + 3n + 1
an+2 = − an , n ≥ 2
3(n + 2)(n + 1)
This last equation is called the recurrence relation. We use the recurrence relation to get the
following results:
11 11
n = 2 → a4 = − a2 = a0
36 216
19 19
n = 3 → a5 = − a3 = a1
60 216
and so on. So the answer, the first six non-zero terms of y is:
1 5 11 19
y = a 0 + a 1x − a 0x 2 − a 1x 3 + a 0x 4 + a 1x 5 + · · ·
6 18 216 216
69
Power series solution when initial condition is given
Example 3.1.2
y 00 + xy 0 + y = 0, y(0) = 1, y 0(0) = −2
Solution) Let’s think about what the initial conditions mean. Since y = a 0 +a 1x +a 2x 2 +a 3x 3 +· · · ,
the initial condition y(0) = 1 means 1 = y(0) = a 0 + a 1 · 0 + a 2 · 02 + a 3 · 03 + · · · = a 0 . In other
words, y(0) gives you the value of a 0 . Also, y 0 = a 1 + 2a 2x + 3a 3x 2 + 4a 4x 3 + · · · and y 0(0) = −2
means −2 = y 0(0) = a 1 + 2a 2 · 0 + 3a 3 · 02 + 4a 4 · 03 + · · · = a 1 . Therefore, we find that a 0 = 1 and
a 1 = −2. Note: If the value of y 00(0) is given, then a 2 will be equal to half of that. (Why?)
Other than this part, the rest is solved exactly like before.
∞
an x n
Õ
y=
n=0
∞
nan x n−1
Õ
y =
0
n=1
∞
n(n − 1)an x n−2
Õ
y 00 =
n=2
Plugging these into the left side of y 00 + xy 0 + y = 0 yields:
∞ ∞ ∞
n−2 n−1
an x n
Õ Õ Õ
n(n − 1)an x +x nan x +
n=2 n=1 n=0
∞ ∞ ∞
n−2 n
an x n
Õ Õ Õ
= n(n − 1)an x + nan x +
n=2 n=1 n=0
∞ ∞ ∞
(n + 2)(n + 1)an+2x n + nan x n + an x n
Õ Õ Õ
=
n=0 n=1 n=0
∞
= (2a 2 + a 0 )x 0 + ((n + 2)(n + 1)an+2 + nan + an ) x n
Õ
n=1
Therefore,
1 1
2a 2 + a 0 = 0 → a 2 = − a 0 = −
2 2
(n + 2)(n + 1)an+2 + nan + an = 0
1
→ an+2 = − an , (n ≥ 1)
n+2
Using the recurrence relation above:
1 2
n = 1 → a3 = − a1 =
3 3
1 1
n = 2 → a4 = − a2 = −
4 8
70
1 2
n = 3 → a5 = − a3 = −
5 15
Hence the first 6 non-zero terms of the solution are:
1 2 1 2
y = 1 − 2x − x 2 + x 3 − x 4 − − x 5 · · ·
2 3 8 15
Here, if r = 0, it’s the usual power series. By a Frobenius solution, we are talking about
solutions where r is negative or a fraction and a 0 , 0. Those solutions cannot be expressed using
just a power series. When do a differential equation have power series solutions and when do a
differential equation have a Frobenius solution? The answer to this requires the understanding
of singular and ordinary points of a differential equation, and we will discuss it later. Before that,
let’s look at an example where Frobenius solution is required:
Example 3.2.1
2x 2y 00 + 3xy 0 − (x 2 + 1)y = 0
∞
(n + r )an x n+r −1
Õ
y0 =
n=0
71
∞
(n + r )(n + r − 1)an x n+r −2
Õ
y 00 =
n=0
Important: In the y 0,
no longer you start the summation from n = 1. That’s because in a Frobe-
nius solution, the series does not start from the constant term, and therefore the first term does
not become zero when differentiated. Same goes for y 00.
Now, let’s plug in to the left side:
∞ ∞ ∞ ∞
2 00 2 2 n+r −2 n+r −1 2
an x − an x n+r
n+r
Õ Õ Õ Õ
2x y +3xy −x y−y = 2x
0
(n+r )(n+r −1)an x +3x (n+r )an x −x
n=0 n=0 n=0 n=0
∞ ∞ ∞ ∞
2(n + r )(n + r − 1)an x n+r + 3(n + r )an x n+r − an x n+r +2 − an x n+r
Õ Õ Õ Õ
=
n=0 n=0 n=0 n=0
Let us shift the third summation by replacing n by n − 2.
∞ ∞ ∞ ∞
2(n + r )(n + r − 1)an x n+r + 3(n + r )an x n+r − an−2x n+r − an x n+r
Õ Õ Õ Õ
=
n=0 n=0 n=2 n=0
2r 2 + r − 1 = 0
(2r − 1)(r + 1) = 0
1
r = , −1
2
Each value of r will lead to different (i.e. linearly independent) Frobenius solution provided that
they do not differ by an integer. If they do differ by an integer, the higher value of r always
72
provides a Frobenius solution, while the lower value may or may not give us a Frobenius solution.
Since solving linear second order differential equation requires finding two linearly independent
solutions, if the second solution is not of the Frobenius form, it has to be found in some other
way. Luckily, it turns out that in such a case, the second solution can be found in the following
form:
∞
bn x k+r 2 ,
Õ
y2 = Cy1 ln x +
n=0
where y1 is the Frobenius solution corresponding to the larger root r 1 , and r 2 is the smaller root.
Such a solution is called an extended Frobenius solution .
2a 1 = 0
1 1 3
2(n + )(n − ) + 3n + − 1 an = an−2 , n ≥ 2
2 2 2
So we have a 1 = 0 and
2n2 + 3n an = an−2 , n ≥ 2
1
an = an−2 , n ≥ 2
n(2n + 3)
This is the recurrence relation.
Let’s try to find the first four non-zero terms using the recurrence relation:
1 1
a2 = a0 = a0
2·7 14
1
a3 = a1 = 0
3·9
(since a 1 = 0. In fact, all the odd items will be zero.)
1 1
a4 = a2 = a0
4 · 11 616
1 1
a6 = a4 = a0
6 · 15 55440
73
Hence:
1 1 1 1 2 1 4 1
1/2 2 4 6 1/2 6
y =x a 0 + a 0x + a 0x + a 0x . . . = a 0x 1+ x + x + x ...
14 616 55440 14 616 55440
is the solution obtained. Now, since a 0 is arbitrary we can choose a 0 = 1 and find that one solution
obtained for r = 12 is:
1 2 1 4 1
1/2 6
y1 = x 1+ x + x + x ...
14 616 55440
But since we are solving second order linear differential equation, we need two linearly indepen-
dent solutions. So we have to find the second solution by putting r = −1. We go back to
1
an = an−2 , n ≥ 2
2n2 − 3n
This is the recurrence relation.
Let’s try to find the first four non-zero terms using the recurrence relation:
1
a2 = a0
2
1
a3 = a1 = 0
9
(since a 1 = 0. In fact, all the odd items will be zero.)
1 1
a4 = a2 = a0
20 40
1 1
a6 = a4 = a0
54 2160
Hence:
1 1 1 1 2 1 4 1 6
2 4 6
y =x −1
a 0 + a 0x + a 0x + a 0x . . . = a 0x −1
1+ x + x + x ...
2 40 2160 2 40 2160
Again, set a 0 = 1 to get the second solution:
1 2 1 4 1 6
y2 = x −1
1+ x + x + x ...
2 40 2160
74
Now take the linear combination of the two solutions y1 and y2 to get the full solution:
1 2 1 4 1 1 2 1 4 1 6
1/2 6
y = C 1x 1+ x + x + x . . . + C 2x −1
1+ x + x + x ...
14 616 55440 2 40 2160
Note: The fact that y1 and y2 are linearly independent can be most easily proved by looking at
the asymptotic behavior of the function near zero, but the proof will be omitted since asymptotic
behavior lies beyond the scope of this book.
Why power series solutions not good for singular differential equation
We will learn that the equation 2x 2y 00 + 3xy 0 − (x 2 + 1)y = 0 has x = 0 as a singular point. But
before we learn about that, we will try to see what goes wrong if one tried to find a power series
solution to this equation.
Let us try to look for the power series form of y as a solution to 2x 2y 00 + 3xy 0 − (x 2 + 1)y = 0.
This means:
∞
an x n
Õ
y=
n=0
∞
nan x n−1
Õ
y0 =
n=1
∞
n(n − 1)an x n−2
Õ
y 00 =
n=2
Plugging this in, we get:
∞ ∞ ∞ ∞
2 00 2 2 n−2 n−1 2 n
an x n
Õ Õ Õ Õ
2x y + 3xy − x y − y = 2x
0
n(n − 1)an x + 3x nan x −x an x −
n=2 n=1 n=0 n=0
∞ ∞ ∞ ∞
n n n+2
an x n
Õ Õ Õ Õ
= 2n(n − 1)an x + 3nan x − an x −
n=2 n=1 n=0 n=0
∞ ∞ ∞ ∞
2n(n − 1)an x n + 3nan x n − an−2x n − an x n
Õ Õ Õ Õ
=
n=2 n=1 n=2 n=0
∞
= a 0x 0 + (3a 1 + a 1 )x 1 + (2n(n − 1)an + 3nan − an−2 − an ) x n
Õ
n=2
Then we have:
a0 = 0
3a 1 + a 1 = 0 → a 1 = 0
2n(n − 1)an + 3nan − an−2 − an = 0 → 2n 2 + n − 1 an = an−2 , n ≥ 2
75
Let’s plug in values of n:
1
n = 2 → a2 = a0 = 0
9
1
n = 3 → a3 = a1 = 0
20
1
n = 4 → a4 = a2 = 0
35
As you can see, all an will turn out to be zero, yielding the solution y = 0. Such answer is called
the trivial solution because any homogeneous differential equation will have zero as one answer
(which is a useless answer in most cases.)
2x 2y 00 + 3xy 0 − (x 2 + 1)y = 0
y 00 + P(x)y 0 + Q(x)y = 0
if you divide by what’s in front of y 00 (which in the above case we had 2x 2 ). Comparing, we see
3x 3 x 2 +1
that P(x) = 2x 2 = 2x and Q(x) = − 2x 2 . You must include the sign when you figure out P(x) and
Q(x).
There are some important concepts you have to know before we go for the indicial equations:
• When, after simplifying P(x) and Q(x), both P(0) and Q(0) are well defined , we say x = 0
is a ordinary point of the differential equation.
• When, after simplifying P(x) and Q(x), if any or both of P(0) and Q(0) are undefined (i.e.
x = 0 causes the denominator to become zero) then we say x = 0 is a singular point of
the differential equation.
If p(0) and q(0) are well defined, we say the differential equation has a regular singular point
at x = 0. If not, x = 0 is said to be an irregular singular point. In other words, singular points
come in two flavors - regular and irregular.
Here is what’s known: If x = 0 is ordinary, then the differential equation can be solved by
power series. If x = 0 is regular singular, then the differential equation can be solved by Frobenius
method. If it is irregular singular, then Frobenius nor the power series method will work.
76
Indicial equation
Once p(x) and q(x) are calculated, the indicial equation is:
r (r − 1) + p(0)r + q(0) = 0
3 1
r (r − 1) + r + = 0
2 2
Multiplying by 2 and simplifying gives us 2r 2 + r − 1 = 0, exactly same as what we had before.
Example 3.2.2
2x 2y 00 + (7x − x 2 )y 0 + 2y = 0
7x − x 2 0 1
y 00 + y + 2y = 0
2x 2 x
x(7−x) 1
We see that P(x) = 2x 2
and Q(x) = x2
. Hence
7−x 7
p(x) = xP(x) = , q(x) = x 2Q(x) = 1 → p0 = p(0) = , q 0 = q(0) = 1
2 2
This makes the indicial equation
7 5
r (r − 1) + r + 1 = 0 → r 2 + r + 1 = 0
2 2
Multiplying the equation by 2, we have 2r 2 + 5r + 2 = 0. Using the quadratic formula, we have
√
−5 ± 25 − 4 · 2 · 2 −5 ± 3 1
r= = → r = −2 , −
2·2 4 2
Since the two numbers −2 and − 12 do not differ by a constant, both of them yield a Frobenius
solution. Let us find the Frobenius solution for each case.
But before we do so, we will rewrite the differential equation 2x 2y 00 + (7x − x 2 )y 0 + 2y = 0 as
2x 2y 00 + 7xy 0 − x 2y 0 + 2y = 0
77
∞ ∞
r n
an x n+(−2)
Õ Õ
y =x an x =
n=0 n=0
We differentiate twice:
∞
(n − 2)an x n−3
Õ
y0 =
n=0
∞
(n − 2)(n − 3)an x n−4
Õ
y 00 =
n=0
(Note that contrary to the case when deriving power series solutions, during Frobenius method
we don’t change n = 0 to n = 1 or n = 2 when taking derivatives.)
Now plug in these formulas into 2x 2y 00 + 7xy 0 − x 2y 0 + 2y = 0:
∞ ∞ ∞ ∞
2 n−4 n−3 2 n−3
an x n−2 = 0
Õ Õ Õ Õ
2x (n − 2)(n − 3)an x + 7x (n − 2)an x −x (n − 2)an x +2
n=0 n=0 n=0 n=0
∞ ∞ ∞ ∞
2(n − 2)(n − 3)an x n−2 + 7(n − 2)an x n−2 − (n − 2)an x n−1 + 2 an x n−2 = 0
Õ Õ Õ Õ
Now, when n = 0, the shifted summation does not produce any terms because that one begins
with n = 1, while the other summations produce terms. We can separate out the terms produced
when n = 0 and rewrite it as:
2(−2)(−3)a 0x −2 + 7(−2)a 0x −2 + 2a 0x −2
∞ ∞ ∞ ∞
2(n−2)(n−3)an x n−2 + 7(n−2)an x n−2 − (n−3)an−1x n−2 +2 an x n−2 = 0
Õ Õ Õ Õ
+
n=1 n=1 n=1 n=1
But 2(−2)(−3)a 0 x −2
+ 7(−2)a 0 x −2
+ 2a 0 x −2
= (12a 0 − 14a 0 + 2a 0 = 0. Hence it cancels away.
)x −2
Also, since all the summations now begin at n = 1, we can just use a single summation to express
all the sums at once:
∞
2(n − 2)(n − 3)an x n−2 + 7(n − 2)an x n−2 − (n − 3)an−1x n−2 + 2an x n−2 = 0
Õ
n=1
78
∞
(2(n − 2)(n − 3)an + 7(n − 2)an − (n − 3)an−1 + 2an ) x n−2 = 0
Õ
n=1
Thus we must have:
2n 2 − 3n an = (n − 3)an−1
n−3
an = an−1 , (n ≥ 1)
n(2n − 3)
This leads to:
a 1 = −21 · (−1)a 0 = 2a 0
−1 1 1
a2 = a 1 = − a 1 = − · 2a 0 = −a 0
2·1 2 2
0
a3 = a2 = 0
3·3
But then a 4 = (somethinд)a 3 = 0 and a 5 is zero for the same reason, and we end up with all an = 0
after n ≥ 3.
Hence the solution is:
y = x −2 a 0 + a 1x + a 2x 2 + a 3x 3 + · · · = x −2 a 0 + 2a 0x − a 0x 2 + 0 + 0 + · · ·
1 2
2
y = a 0x 1 + 2x − x = a 0 2 + − 1
−2
x x
Case r = − 21
∞ ∞
r n
an x n+(−1/2)
Õ Õ
y =x an x =
n=0 n=0
We differentiate twice:
∞
1
n − an x n−3/2
Õ
y =0
n=0
2
∞
1 3
n − an x n−5/2
Õ
y =
00
n−
n=0
2 2
Now plug in these formulas into 2x 2y 00 + 7xy 0 − x 2y 0 + 2y = 0:
∞ ∞ ∞ ∞
1 3 1 1
2 n−5/2 n−3/2 2 n−3/2
an x n−1/2 = 0
Õ Õ Õ Õ
2x n− n − an x +7x n − an x −x n − an x +2
n=0
2 2 n=0
2 n=0
2 n=0
79
∞ ∞ ∞ ∞
1 3 1 1
n−1/2 n−1/2
n − an x n+1/2 + 2an x n−1/2 = 0
Õ Õ Õ Õ
2 n− n − an x + 7 n − an x −
n=0
2 2 n=0
2 n=0
2 n=0
We see that the second from the last summation needs shifting of indices. For that summation
only, we replace n by n − 1 so that:
∞ ∞
1 3
n+1/2
n − an−1x n−1/2
Õ Õ
n − an x →
n=0
2 n=1
2
Now, when n = 0, the shifted summation does not produce any terms because that one begins
with n = 1, while the other summations produce terms. We can separate out the terms produced
when n = 0 and rewrite it as:
1 3 1
2(− )(− )a 0x −1/2 + 7(− )a 0x −1/2 + 2a 0x −1/2
2 2 2
∞ ∞ ∞ ∞
1 3 1 3
n−1/2 n−1/2 n−1/2
+ 2an x n−1/2 =
Õ Õ Õ Õ
+ 2 n− n − an x + 7 n − an x − n − an−1x
n=1
2 2 n=1
2 n=1
2 n=1
But 2(− 12 )(− 32 )a 0x −1/2 + 7(− 12 )a 0x −1/2 + 2a 0x −1/2 = ( 23 a 0 − 72 a 0 + 2a 0 )x −1/2 = 0. Hence it cancels
away. Also, since all the summations now begin at n = 1, we can just use a single summation to
express all the sums at once:
∞
1 3 1 3
n−1/2 n−1/2
− n − an−1x n−1/2 + 2an x n−1/2 = 0
Õ
2 n− n − an x + 7 n − an x
n=1
2 2 2 2
∞
1 3 1 3
n − an + 7 n − an − n − an−1 + 2an x n−1/2 = 0
Õ
2 n−
n=1
2 2 2 2
Thus we must have:
1 3 1 3
2 n− n − an + 7 n − an − n − an−1 + 2an = 0
2 2 2 2
80
4n2 + 6n an = (2n − 3)an−1
2n − 3
an = an−1 , (n ≥ 1)
2n(2n + 3)
This leads to:
1
a 1 = −12 · 5a 0 = − a 0
10
1 1 1 1 1
a2 = a1 = a1 = · − a0 = − a0
4·7 28 28 10 280
3 1 1 1 1
a3 = a2 = a2 = · − a0 = − a0
6·9 18 18 280 5040
5 1
a4 = a2 = − a0
8 · 11 88704
Hence the solution is:
1 1 1 1
2 3 2 3 4
y =x −1/2
a 0 + a 1x + a 2x + a 3x + · · · = x a 0 − a 0x −
−1/2
a 0x − a 0x − a 0x + · · ·
10 280 5040 88704
1 1 2 1 3 1
4
y = a 0x −1/2
1− x − x − x − x +···
10 280 5040 88704
Now taking the linear combinations of the two solutions obtained for r = −2 and r = − 12 , we
get:
1 2 1 1 2 1 3 1
4
y = C 1 2 + − 1 + C 2x −1/2
1− x − x − x − x +···
x x 10 280 5040 88704
Case when the two roots of the indicial equation differ by an integer
Example 3.2.3
x 2y 00 + 5xy 0 + (3 − x)y = 0
r (r − 1) + 5r + 3 = 0 → r 2 + 4r + 3 = 0
81
∞ ∞
r n
an x n+(−1)
Õ Õ
y =x an x =
n=0 n=0
We differentiate twice:
∞
(n − 1)an x n−2
Õ
y0 =
n=0
∞
(n − 1)(n − 2)an x n−3
Õ
y 00 =
n=0
(Note that contrary to the case when deriving power series solutions, during Frobenius method
we don’t change n = 0 to n = 1 or n = 2 when taking derivatives.)
Now rewrite the given equation x 2y 00 + 5xy 0 + (3 − x)y = 0 as x 2y 00 + 5xy 0 + 3y − xy = 0 and
plug in these formulas:
∞ ∞ ∞ ∞
x2 (n − 1)(n − 2)an x n−3 + 5x (n − 1)an x n−2 + 3 an x n−1 − x an x n−1 = 0
Õ Õ Õ Õ
Now, when n = 0, only the first three summations produce terms and the last summation is
at rest. (The last one only begins producing terms when n = 1.) We can separate out the terms
produced when n = 0 and rewrite it as:
(−1)(−2)a 0x −1 + 5 · (−1)a 0x −1 + 3a 0x −1
∞ ∞ ∞ ∞
(n − 1)(n − 2)an x n−1 + 5 (n − 1)an x n−1 + 3 an x n−1 − an−1x n−1 = 0
Õ Õ Õ Õ
+
n=1 n=1 n=1 n=1
But (−1)(−2)a 0 x −1
+ 5 · (−1)a 0 x −1
+ 3a 0 x −1
= (2a 0 − 5a 0 + 3a 0 = 0. Hence it cancels away.
)x −1
Also, since all the summations now begin at n = 1, we can just use a single summation to express
all the sums at once:
∞
(n − 1)(n − 2)an x n−1 + 5(n − 1)an x n−1 + 3an x n−1 − an−1x n−1 = 0
Õ
n=1
82
∞
((n − 1)(n − 2)an + 5(n − 1)an + 3an − an−1 ) x n−1 = 0
Õ
n=1
Thus we must have:
2
a3 = a0
4!5!
In general, the pattern is:
2
an = a 0 , (n ≥ 1)
(n + 1)!(n + 2)!
2
Notice that the above is true even when n = 0, because (0+1)!(0+2)! a 0 = 22 a 0 = a 0 . So we can take
the above as a solution for all an . Thus the solution is:
∞ ∞
2 n−1 2
x n−1
Õ Õ
y= a 0x = a0
n=0
(n + 1)!(n + 2)! n=0
(n + 1)!(n + 2)!
Case r = −3
∞ ∞
y = xr an x n = an x n+(−3)
Õ Õ
n=0 n=0
We differentiate twice:
∞
(n − 3)an x n−4
Õ
y0 =
n=0
∞
(n − 3)(n − 4)an x n−5
Õ
y 00 =
n=0
83
Now rewrite the given equation x 2y 00 + 5xy 0 + (3 − x)y = 0 as x 2y 00 + 5xy 0 + 3y − xy = 0 and
plug in these formulas:
∞ ∞ ∞ ∞
2 n−5 n−4 n−3
an x n−3 = 0
Õ Õ Õ Õ
x (n − 3)(n − 4)an x + 5x (n − 3)an x +3 an x −x
n=0 n=0 n=0 n=0
∞ ∞ ∞ ∞
n−3 n−3 n−3
an x n−2 = 0
Õ Õ Õ Õ
(n − 3)(n − 4)an x + 5 (n − 3)an x +3 an x −
n=0 n=0 n=0 n=0
We see that the last summation needs shifting of indices. For that last summation only, we
replace n by n − 1 so that:
∞ ∞
an x n−2 → an−1x n−3
Õ Õ
n=0 n=1
Note that the n = 1 in the bottom is due to changing n = 0 to n − 1 = 0 which is equivalent to
n = 1. The resulting summation is
∞ ∞ ∞ ∞
(n − 3)(n − 4)an x n−3 + 5 (n − 3)an x n−3 + 3 an x n−3 − an−1x n−3 = 0
Õ Õ Õ Õ
(−3)(−4)a 0x −3 + 5 · (−3)a 0x −3 + 3a 0x −3
∞ ∞ ∞ ∞
n−3 n−3 n−3
an−1x n−3 = 0
Õ Õ Õ Õ
+ (n − 3)(n − 4)an x +5 (n − 3)an x +3 an x −
n=1 n=1 n=1 n=1
n=1
Thus we must have:
84
1
an = an−1 , (n ≥ 1)
n(n − 2)
Until now, this solution may look exactly like the solution we wrote for r = −1. How-
1
ever, notice that an = n(n−2) an−1 is not defined for n = 2, because we get a 0 in the denominator.
That is why there is no Frobenius solution corresponding to r = −3.
85
System of Differential Equations
4
4.1 Introduction to System of Differential Equations
Introduction to system of differential equations (part I)
Suppose we have a system of differential equations:
(
x + 2y = x 0
4x + 3y = y 0
with initial condition x(0) = 1 , y(0) = 3. In the above, x and y are functions of t and the
derivatives are with respect to t. This means we are looking for two functions x(t) and y(t) that
satisfy the above differential equations and initial conditions.
We can put the system into a matrix notation as follows:
0
x 1 2 x
=
y 4 3 y
After the system is changed into the matrix form, it can be solved using a related problem in
matrices called the eigenvector problem.
a
86
For a given matrix A, an eigenvector of A is a non-zero column vector v satisfying Av = λv
for some constant λ. For example, if we consider the eigenvector of the matrix appearing in the
above system of differential equations, it means we are looking to find numbers a, b and λ that
satisfy:
1 2 a a
=λ
4 3 b b
In the above, vector v is the column vector made of a and b, and the condition that v should be
non-zero means that at least one of a and b should be non-zero. The value of λ that satisfies the
above equation is called the eigenvalue associated to the eigenvector v.
Let’s try to find what a, b and λ can satisfy the above matrix equation. First, multiply the left
side:
a + 2b λa
=
4a + 3b λb
Now, the above is equivalent to the following system:
(
a + 2b = λa
4a + 3b = λb
We can treat this equation as an equation of a and b. Move λa to the left and group it with a. Also
move λb to the left and group it with 3b. We end up with:
(
(1 − λ)a + 2b = 0
4a + (3 − λ)b = 0
1−λ 2 a 0
=
4 3−λ b 0
This has an easy solution. Just multiply the inverse matrices to both sides:
−1 −1
1−λ 2 1−λ 2 a 1−λ 2 0
=
4 3−λ 4 3−λ b 4 3−λ 0
Then there is a cancellation on the left side, while on the right, anything times zero matrix is zero:
a 0
=
b 0
BUT WAIT! We are not supposed to have both a and b equal to zero!
So what went wrong? The only way this solution could go wrong is when the inverse matrix
does not exist. (Since we want at least one of a and b to be non-zero, we want it to go wrong -
then there is a chance that we could have at least one of a and b to be non-zero. )
So when does a matrix NOT have an inverse? That’s when the determinant is zero, i.e.
1 − λ 2
4 =0
3 − λ
87
Computing the determinant, we have:
(1 − λ)(3 − λ) − 2 · 4 = 0
λ2 − 4λ − 5 = 0
(λ − 5)(λ + 1) = 0
So far our conclusion means this: The only way we could have a non-zero a and b is when
λ = 5 or λ = −1
Now what we need is to figure out what a and b could be for each value of λ.
• Case λ = 5
• Case λ = −1
88
unknowns a and b. It means we have infinitely many such pairs. (This will always be the case,
when looking for eigenvectors. ) By trial and error, you can quickly see that a = 1 and b = −1
is one such pair satisfying 2a + 2b = 0. Hence we found that for eigenvalue λ = −1, we can take
the eigenvector to be:
1
−1
1
1
−1
1 5t
u1 (t) = e
2
and
1 −t
u2 (t) = e ,
−1
you can check that they separately satisfy the matrix form of the system of differential equations.
Then we take the linear combinations of the two solutions to write the general solution:
x 1 5t 1 −t
= C 1 u1 (t) + C 2 u2 (t) = C 1 e + C2 e
y 2 −1
89
• Step 2. Find the eigenvalues λ 1 and λ 2 , and also their associated eigenvectors, so that they
satisfy:
k l a1 a
= λ1 1
m n b1 b1
and likewise for λ 2 , a 2 and b2 .
x a 1 λ 1t a
= C1 e + C 2 2 e λ 2t
y b1 b2
Let’s see what we have achieved. How many derivatives do you see? Only one! So we have in
front of us a first order differential equation. However, there is a problem with just having the
above. Since we have two unknown functions y1 and y2 , we need two differential equations. What
is another relationship between y1 and y2 ? Well, here is the obvious relationship: y10 = y 0 = y2 .
Therefore we have the following system of differential equations:
(
y20 − 2y2 + 3y1 = x 3 + 5x
y10 = y2
Now, it is customary to rewrite it so that y10 and y20 appear left. So the final answer is:
(
y10 = y2
y20 = −3y1 + 2y2 + x 3 + 5x
90
Example 4.1.1 Rewrite the differential equation as a system of first order differential equations:
y 000 − 3y 00 + 2y 0 − y = 0
Solution) We define y1 = y, y2 = y 0 and y3 = y 00. (Rule of the thumb is that a 3rd order differ-
ential equation will need three functions y1 , y2 and y3 , while second order differential equation
would need just two.) Now, y 000 = (y 00)0 = y30 . Therefore:
y30 − 3y3 + 2y2 − y1 = 0
Then we add in two more relations (two more because we need a total of three differential equa-
tions in all) which are:
y10 = y2
y20 = y3
Writing everything as a system of differential equations, we have:
y 0 = y2
1
y20 = y3
y 0 = y1 − 2y2 + 3y3
3
4 3
We have learned that finding eigenvalues and eigenvectors is the key to solving system of
differential equations. Here are the steps in finding eigenvalues and eigenvectors.
• Step 1. Subtract λ from the diagonal and set the determinant equal to zero, i.e.
1 − λ 2
4 =0
3 − λ
We have already solved this in the first section. The solution is λ = 5, −1.
• Step 2. For each solution, plug it back into λ and look for a null vector.
91
System of differential equations example
Example 4.1.2 Solve:
7 −4
X =
0
X
5 −2
Solution) In the above, X is a column vector of functions of t, i.e.
x(t)
X=
y(t)
We have to first find the eigenvalues and eigenvectors. So we follow the steps outlined in the
previous section:
7 − λ
−4 =0
5 −2 − λ
(7 − λ)(−2 − λ) − (−4) · 5 = λ2 − 5λ + 6 = 0
(λ − 3)(λ − 2) = 0
λ = 2, 3
• Case λ = 2
5 −4 a 0
=
5 −4 b 0
Multiplying the top row of the first matrix by the column vector, we get:
5a − 4b = 0
(You never have to worry about the second line of the matrix, because it will always give you
another equation equivalent to what you get from the first line.) We just need one solution for
5a − 4b = 0, which we choose to be a = 4, b = 5. (TIP: Switch the coefficients around and flip
one of the signs. In this example, the coefficients are 5 and -4. Switching them gives you -4 and
5. Now change the sign of -4 to make 4. Then we have 4 and 5. That’s how we got a = 4, b = 5.)
• Case λ = 3
4 −4 a 0
=
5 −5 b 0
The first line multiplied to the column vector yields 4a −4b = 0. Dividing by 4 we get a −b = 0.
It’s easy to see that a = 1, b = 1 is a solution.
Since we have obtained the eigenvectors for the eigenvalues, we put them together to get the
solution of the system of differential equations:
x(t) 4 2t 1 3t
= C1 e + C2 e
y(t) 5 1
This solution is equivalent to saying:
(
x(t) = 4C 1e 2t + C 2e 3t
y(t) = 5C 2e 2t + C 2e 3t
92
System of differential equations example with initial conditions
Example 4.1.3 Solve: (
x 10 = 4x 1 − 3x 2
x 20 = 6x 1 − 7x 2
with initial conditions x 1 (0) = 1, x 2 (0) = −1.
(4 − λ)(−7 − λ) − (−3) · 6 = λ2 + 3λ − 10 = 0
(λ + 5)(λ − 2) = 0
λ = −5, 2
• Case λ = −5
9 −3 a 0
=
6 −2 b 0
The first line multiplied to the column vector yields 9a − 3b = 0. Dividing by 3 we get
3a − b = 0. It’s easy to see that a = 1, b = 3 is a solution.
• Case λ = 2
2 −3 a 0
=
6 −9 b 0
The first line multiplied to the column vector yields 2a − 3b = 0. A possible solution is
a = 3, b = 2 is a solution.
x(t) 1 −5t 3 2t
= C1 e + C2 e
y(t) 3 2
1 x(0) 1 0 3 0
= = C1 e + C2 e
−1 y(0) 3 2
This is equivalent to (
1 = C 1 + 3C 2
−1 = 3C 1 + 2C 2
93
Solving this system gives us: C 1 = − 75 , C 2 = 47 . Thus after plugging the constants into the general
solution, the solution to the initial problem is:
x 1 (t) 5 1 −5t 4 3 2t
=− e + e
x 2 (t) 7 3 7 2
94
End Point Value Problems
5
5.1 End Point Value Problems
Introduction to endpoint value problems
Example 5.1.1
y 00 − y = 0, y(0) = 2, y(4) = 1
The above problem seems very similar to what we have been solving before. However, instead of
giving the value of y 0(0), y(4) = 1 is given. In other words, the value of the function y(x) is given
at two different points. Since this often happens when solving a physics problem with conditions
given at two endpoints, such a problem is called an endpoint value problem as opposed to the
case when the value of y(0) and y 0(0) are given in which case we call initial value problem. (Some
initial value problems provide the value of y(a) and y 0(a) at some point other than a = 0.)
Let’s try to solve this to see what difference it makes. For this question, it’s not much different.
We use the characteristic equation
r 2 − 1 = 0 → r = 1, −1
Which means
y = C 1e x + C 2e −x
Now, applying endpoint conditions y(0) = 2, y(4) = 1 gives us:
(
C1 + C2 = 2
C 1e 4 + C 2e −4 = 1
Multiplying e 4 to the second equation turns this into:
(
C1 + C2 = 2
C 1e 8 + C 2 = e 4
95
Then we subtract both sides to get:
C 1 (1 − e 8 ) = 2 − e 4
2 − e4
C1 =
1 − e8
2 − e4
C2 = 2 − C1 = 2 −
1 − e8
Thus the solution is:
2 − e4 x 2 − e 4 −x
y= e + 2− e
1 − e8 1 − e8
So far you may not see anything different from an initial value problem. However, think about
the next problem:
Example 5.1.2
y 00 + y = 0, y(0) = 0, y(π ) = 1
If you compare this problem with the first one, there is not much of a difference other than
a negative sign changed to positive and change of the end point conditions. However, you’ll see
this equation behaves quite differently. Let’s solve it:
r 2 + 1 = 0 → r = ±i → y = C 1 cos x + C 2 sin x
Example 5.1.3
y 00 + y = 0, y(0) = 0, y(π ) = 0
Here is another question very similar to the previous example, but with a very different conclu-
sion. The general solution is the same, so let’s just look at the endpoint conditions:
96
Now, we see that just C 1 = 0 alone satisfies both endpoint conditions. Therefore C 2 can be
anything, and the solution satisfying this problem is:
y = C 2 sin x
What does this mean? It means that y = 10 sin x, y = 100 sin x or y = −12125 sin x are all
solutions. There are infinite choices for C 2 and thus there are infinite solutions for this question.
Lesson we learn here is that for endpoint problems there may be infinite number of solutions
even though two conditions are given. In other words, while there is the uniqueness theorem for
initial value problems, there is no uniqueness theorem for endpoint value problems.
Q: For what value of λ does y 00 + λy = 0, y(0) = 0, y(L) = 0 have infinitely many solutions?
This is a special case of what people call the Sturm-Liouville problem. Since y = 0 always
satisfies the above, we know that either there is one solution or infinitely many solutions. Having
a non-zero solution (also called a non-trivial solution) y = f (x) immediately gives us infinitely
many solutions because any constant multiple of it, i.e. y = C f (x) will be again a solution. In
other words, the above question is equivalent to the following one:
Q: For what value of λ does y 00 + λy = 0, y(0) = 0, y(L) = 0 have solution other than the trivial
solution x(0) = 0 ?
It turns out that only for a few values of λ such non-trivial solution will exist. Sturm-
Liouville problem is about finding that value of λ that allows a non-trivial solution y (and
also find those non-trivial functions y).
Often y 00 + λy = 0 is written as −y 00 = λy. This is because the actual Sturm-Liouville problem
is finding the value of λ where the following equation has a non-trivial solution:
d d
− p(x) + r (x) y = λy,
dx dx
subject to various endpoint conditions such as y(0) = 0, y(L) = 0. The equation −y 00 = λy is
obtained when p(x) = 1 and r (x) = 0 in the above problem. That’s why −y 00 = λy is a special
case of the Sturm-Liouville problem. We will limit our attention to this special case only, with
two kinds of endpoint conditions: either {y(0) = 0, y(L) = 0} or {y 0(0) = 0, y 0(L) = 0}.
In any case, Sturm-Liouville problem always takes the form of Ly = λy where L is a second
order linear differential operator of the form stated above (which is called the Sturm-Liouville
97
operator). When a value λ yields a non-trivial solution y, we say λ is an eigenvalue and a
nontrivial solution y an eigenfunction associated with λ.
There are several important theorems regarding Sturm-Liouville problems, most of which are
beyond the scope of this book. However, there are two that we will make use of, without giving
proofs:
• For each eigenvalue, the eigenfunctions are unique up to a constant multiple, i.e. if f (x) and
д(x) are both eigenfunctions of the same eigenvalue, then f (x) = cд(x) for some constant
c. An eigenvalue having such property is called simple. In other words, all eigenvalues of
the Sturm-Liouville operator is simple.
Solution) Note that λ can be negative, zero or positive. However, the mathematical theorem
on Sturm-Liouville operator guarantees us that λ cannot be negative, so we will only consider
the case when λ is zero or positive.
a) case λ = 0
In this case the equation reads:
x 00 = 0.
Integrating both sides x 00dt∫ = 0dt, we get x 0 = C 1 , where C 1 is an integration constant. Then
∫ ∫
x = C 1t + C 2
x(5) = C 1 · 5 + 0 = 5C 1 → 5C 1 = 0
Thus we have both C 1 and C 2 zero and we get x(t) = 0. This is a trivial solution, and therefore
λ = 0 is not an eigenvalue.
b) case λ > 0
√
r 2 + λ = 0 → r 2 = −λ → r = ± λi
Thus √ √
x(t) = C 1 cos λt + C 2 sin λt
Now we apply the endpoint conditions:
98
√ √ √
0 = x(5) = 0 · cos( λ5) + C 2 sin( λ5) → C 2 sin( λ5) = 0
√ √
Now, suppose
√ sin( λ5) , 0. Then we are forced to have C 2 = 0 but
√ then x(t) = C 1 cos λt +
C 2 sin λt leads to x(t) = 0, a trivial solution. So√ we must have sin( λ5) = 0.√Since sin t is zero
for t = nπ with n an integer, the only way sin( λ5) = 0 can happen is when λ5 = nπ . Solving
for λ we get
n2π 2
λ=
25
with n = 1, 2, 3, · · · .
With these values of λ, we have
nπ
x(t) = C 2 sin , n = 1, 2, 3, · · ·
5
So the answer is that the eigenvalues and eigenfunctions are
n2π 2 nπ
λ= , x(t) = sin for n = 1, 2, 3, · · · .
25 5
Some comments
Note that we dropped C 2 in front of sin nπ 5 when writing eigenfunctions. Since any non-zero
multiple of an eigenfunction is an eigenfunction, it is customary to write just one eigenfunction
to represent a class of eigenfunctions who are multiples of each other.
If you are curious to know what happens if λ < 0, you can check directly to confirm that
indeed there are only trivial solutions in this case. If λ < 0, then set λ = −k 2 for some positive
number k. Then the equation reads:
x 00 − k 2x = 0
The characteristic equation of the above is r 2 − k 2 = 0 whose solutions are r = k and r = −k.
Then we have the solution:
x(t) = C 1e kt + C 2e −kt .
Applying x(0) = 0 yields:
x(0) = C 1 · 1 + C 2 · 1 → C 1 + C 2 = 0,
e 1 0C 1 − C 1 = 0 → (e 10 − 1)C 1 = 0
So C 1 = e 100−1 = 0. We also have C 2 = −C 1 = 0. Thus we have both C 1 and C 2 zero and we get
x(t) = 0. This is a trivial solution, and therefore λ = 0 is not an eigenvalue.
When solving eigenvalue-eigenfunction problems for Sturm-Liouville operators, we will al-
ways think about two cases — first λ = 0 and second λ > 0. In the following example, we will
find that λ = 0 does give us a non-trivial solution and thus λ = 0 is an eigenvalue.
99
An example where the boundary conditions are given by derivatives.
Example 5.2.2 x 00 + λx = 0, x 0(0) = x 0(4) = 0 Find eigenvalues and eigenfunctions.
Solution)
a) case λ = 0
In this case the equation reads:
x 00 = 0.
Integrating both sides x 00dt = 0dt, we get x 0 = C 1 , where C 1 is an integration constant. Then
∫ ∫
x = C 1t + C 2
x(0) = C 1 → C 1 = 0
But then x(t) = C 2 and x 0(t) = 0. This means x 0(4) is satisfied automatically.
Thus we see that x(t) = C 2 satisfies the endpoint conditions, and indeed this is a non-trivial
solution (provided that we don’t choose C 2 = 0). Therefore λ = 0 is an eigenvalue with the
eigenfunction x(t) = 1.
b) case λ > 0
√
r 2 + λ = 0 → r 2 = −λ → r = ± λi
Thus √ √
x(t) = C 1 cos λt + C 2 sin λt
√ √ √ √
We also get x 0(t) = −C 1 λ sin λt + C 2 λ cos λt . Now we apply the endpoint conditions:
√ √ √
0 = x 0(0) = −C 1 λ sin 0 + C 2 λ cos 0 = C 2 λ → C 2 = 0
100
Curiously, if you think about the case when n = 0, the above function produces
x(t) = C 1 cos 0 = C 1 ,
which is the eigenfunction for the case when λ = 0. Thus, we can put n = 0, 1, 2, 3, · · · instead.
So the answer is that the eigenvalues and eigenfunctions are
n2π 2 nπ
λ= , x(t) = cos for n = 0, 1, 2, 3, · · · .
16 4
101
Fourier Series
6
6.1 Fourier Series
Even and odd functions:
An even function is a function whose graph is symmetric with respect to the y-axis. Another
way to put it is that if you draw the graph and fold it along the y-axis, the left and right will
match evenly. An odd function is a function whose graph is symmetric with respect to the
origin. Another way to put it is that if you draw the graph and rotate it 180 degrees with respect
to the origin, the graph is still the same. For example, the graph of y = cos x is even, while the
graph of y = sin x is odd. See the picture below:
Other examples of even functions are x 2 , x 4 , or any x n with n an even number. Any constant
multiple of an even function is even, e.g. 17x 4 . You can take sum or difference of even functions
and you get an even function, e.g. x 2 + x 4 , x 8 − 17x 4 .
102
Other examples of odd functions are x 3 , x 5 , or any x n with n an odd number. Any constant
multiple of an odd function is even, e.g. 16x 3 . You can take sum or difference of odd functions
and you get an odd function, e.g. x 3 + x, x − 16x 3 .
Most functions are neither even nor odd. y = ln x, y = e x , and y = x 2 + 6x + 7 are all examples
of functions which are neither symmetric with respect to the y-axis nor the origin.
If you multiply an odd function with an even function, you get an odd function. This can
be easily seen if you consider the example x 3 · x 2 = x 5 . In similar fashion, you can see that
multiplying two odd functions or multiplying two even functions will give you an even function.
Why are even or odd important? It is because of the following observations:
• Integrating an odd function over a symmetric domain gives you zero. In other words, if f (x)
is an odd function, then ∫ a
f (x)dx = 0
−a
In the above integral, the interval [−a, a] is symmetric over the origin, so we say it’s a
symmetric domain. (On the other hand [−2, 10] would be an example of a domain that is
not symmetric.)
• If you integrate an even function over a symmetric domain, you can simplify the calculation
by doubling the integral over the half of the domain, i.e.
∫ a ∫ a
f (x)dx = 2 f (x)dx
−a 0
Orthogonality relations
Since sin nx is an odd function and cos mx is an even function for n and m positive integers, their
product is an odd function. Which means:
∫ π
sin nx cos mxdx = 0
−π
π
(
0, ifn , m
∫
cos nx cos mxdx =
−π π , ifn = m
You just have to use the identity cos A · cos B = 12 (cos(A + B) + cos(A − B)) to prove the above.
These integral results are called orthogonality relations.
103
One more formula that we will use frequently when computing Fourier Series is the following:
sin nπ = 0
cos nπ = (−1)n
You should draw the graph of sine and cosine to convince yourself that indeed sine is zero at
multiples of π while cosine alternates between 1 and -1 at multiples of π .
Here is another example. This time, sum of some sine functions are trying to mimick a step
function:
Not impressed? Well, if we add more sine functions, the result is more convincing:
104
Now all this begs the question: ”How can we figure out the coefficients of these sine or cosine
functions used above?”
Let’s take a look at the function y = x 2 and try to figure out its Fourier series from scratch.
First thing to do is to try to represent x 2 as a sum of sine or cosine functions, i.e.:
Now, x 2 is an even function, so it can only be sum of even functions. So all those sine functions
shouldn’t be there, or another way to say it will be that bn will be zero for all n. Therefore we
have:
x 2 = a 0 + a 1 cos x + a 2 cos 2x + a 3 cos 3x + · · ·
Then let’s integrate both sides over [−π , π ].
∫ π ∫ π
2
x dx = a 0 + a 1 cos x + a 2 cos 2x + a 3 cos 3x + · · · dx
−π −π
Using
π
1 1
∫
cos nxdx = (sin nπ − sin −nπ ) = (0 − 0) = 0,
−π n n
∫π
we see that most of the integral of the right side is zero except for −π a 0dx = a 0 · 2π
On the other hand, the left side is:
∫ π ∫ π
2 1 2π 3
x dx = 2 x 2dx = 2 · π 3 − 03 =
−π 0 3 3
So setting both sides equal gives us:
2π 3
= a 0 · 2π
3
2π 3 1 π2
a0 = · =
3 2π 3
105
Look at what we have! We found what a 0 is!
Now comes the ingenuous part. We multiply cos x both sides:
and integrate it over [−π , π ]. Good thing that we already have these integrals calculated, in the
form the orthogonality relations:
∫ π (
0, ifn , m
cos nx cos mxdx =
−π π , ifn = m
On the other hand, the integral on the left side can be done via integral by parts (or tabular
integration), to get:
∫ π ∫ π
2 π
x cos xdx = 2 x 2 cos xdx = 2 x 2 sin x + 2x cos x − 2 sin x 0 = −4π
−π 0
−4π = a 1 π
a 1 = −4
We’ve only done this for a 0 and a 1 , but the rest is same.
∞
a0 Õ
f (x) = + (an cos nx + bn sin nx)
2 n=0
106
(Notice the denominator 2 under a 0 . This change will be explained in a moment.) Just like
the Taylor series approximation, if you truncate the right side, you get an approximation of the
function f (x), and the more terms you include on the right side, the better the approximation
is. The only exception to this are points where the function f (x) is not continuous. There is a
mathematical theorem which proves that if the right side is taken to be the infinite series, the
two sides will be equal everywhere except at points of jump discontinuity and also possibly at
the endpoints of the interval (−π and π in our case). So this is the way the equal sign of the above
formula should be understood. (The mathematically elaborate terminology for this is equal almost
everywhere or equal in L 2 sense.)
As we have seen before, multiplying cos nx both sides and integrating both sides over the
interval [−π , π ] enables one to find an . Going through this process, one can find the formula:
1 π
∫
an = f (x) cos nxdx
π −π
Now, a 0 is slightly different because we don’t have to multiply anything. However, since cos 0x =
cos 0 = 1, we may still use the above formula for a 0 . The reason there is the denominator 2 under
a 0 is to make the above formula work for n = 0. (If you want to know why, look back on the
previous calculation done on x 2 and see that when a 0 was calculated, we had to divide by 2π .)
Here is the formula for bn :
1 π
∫
bn = f (x) sin nxdx
π −π
The formula is obtained in the same manner, with only difference being multiplying sin nx
instead of cos nx.
π
(
0, ifn , m
∫
sin nx sin mxdx =
−π π , ifn = m
∫ π (
0, ifn , m
cos nx cos mxdx =
−π π , ifn = m
The above integral can be split as follows:
∫ π
= 4 cos 3x · 3 sin x + 4 cos 3x · 4x + 4 cos 3x · 5x 3
−π
107
−4 cos 3x · 10 cos 3x + 4 cos 3x · 2 cos 8xdx
∫ π ∫ π ∫ π
= 12 cos 3x sin xdx + 16x cos 3xdx + 20x 3 cos 3xdx
−π −π −π
∫ π ∫ π
+ −40 cos 3x cos 3xdx + 8 cos 3x cos 8xdx
−π −π
∫π ∫π
From orthogonality relations, we know that −π 12 cos 3x sin xdx = 12 −π cos 3x sin 1xdx = 0
∫π ∫π ∫π
and −π 8 cos 3x cos 8xdx = 0, while −π −40 cos 3x cos 3xdx = −40 −π cos 3x cos 3xdx = −40π .
∫π
For −π 16x cos 3xdx, we know that 16x is an odd function while cos 3x is an even function.
Thus the∫ product is odd and since it is integrated over the symmetric domain [−π , π ], the integral
π
is zero. −π 20x 3 cos 3xdx is zero for the same reason.
Therefore, we found that all the integrals are zero except the one where cos 3x is multiplied
by itself. In that case the answer is:
∫ π ∫ π
−40 cos 3x cos 3xdx = −40 cos 3x cos 3xdx = −40π
−π −π
108
π
1
∫
= (x 2 − 3x) cos nxdx
π −π
Now, using the fact that even function times even function is an even function and also that odd
function times odd function is odd, we see that the above is equal to:
2 π 2
∫
= x cos nxdx
π 0
Then using the integration by parts, we have:
π
2 x2 2x 2 4
= sin nx + 2 cos nx − 3 sin nx = 2 (−1)n
π n n n 0 n
step 3 Find bn
1 π
∫
bn = f (x) sin nxdx
π −π
2 π
∫
= −3x sin nxdx
π 0
(Note: x 2 sin nx is an odd function, and so its integral will be zero.) By integration by parts, we
have: π
2 −3x 6
−3
=− cos nx + 2 sin nx = (−1)n
π n n 0 n
step 4 Put everything into the series:
∞
π2 Õ 4 6
2 n n
x − 3x = + (−1) cos nx + (−1) sin nx
3 n=1 n2 n
109
interval [−π , π ], some integer multiples of periods of these trigonometric functions fit exactly.
Adding such functions will give us a function that has period of 2π . To explain why, think about
a combination of events that repeat every week and events that repeat every day and events
that repeat every hour. Then every week will have the same schedule and the combination of all
events will have a period of one week (i.e. every week will be a repeat of the previous week.)
Now, notice that the formulas for an and bn are
1 π
∫
an = f (x) cos nxdx
π −π
1 π
∫
bn = f (x) sin nxdx
π −π
which only makes use of the value of the function f (x) on the interval [−π , π ]. Thus when you
have a function f (x), not necessarily periodic with period 2π , then the Fourier Series
∞
a0 Õ
f (x) = + (an cos nx + bn sin nx)
2 n=0
are meant to hold only true on the interval [−π , π ]. Here is another way to think about this. Since
the right side is periodic with period 2π , the Fourier series takes a function’s value on the interval
[−π , π ] and makes a period extension of this restricted function.
All this observation can now be used to come up with a formula for the Fourier series on the
interval [−L, L]. First, we introduce L into the formula cos nπx nπx
L and sin L so that these functions
will have integer multiple periods that fit exactly in the interval [−L, L].
The result is
∞
a0 Õ nπx nπx
f (x) = + an cos + bn sin
2 n=0 L L
and the integrals for finding the coefficients are also changed accordingly:
1 L nπx
∫
an = f (x) cos dx
L −L L
1 L nπx
∫
bn = f (x) sin dx
L −L L
1
The L factor is introduced because the orthogonality relations are
∫ L (
nπx mπx 0, ifn , m
cos cos =
−L L L L, ifn = m
and
L
(
nπx mπx 0, ifn , m
∫
sin sin =
−L L L L, ifn = m
which I leave the reader to check. You should also check that the Fourier series is exactly same
as the previous one if you set L = π , just as it should.
110
Calculation of Fourier Series
Example 6.1.3 Find the Fourier series of
f (x) = |x | + 2x on [−3, 3]
Solution) Before we calculate, note that the interval given is [−3, 3], so that L = 3. Also
note that |x | is an even function and 2x is an odd function. Thus when integrating with the odd
function sin mπ 3 x, the even part |x | will yield zero while when integrating with the even function
mπ
cos 3 x, the odd part 2x will yield zero. We compute the integrals as follows.
3 3
1 0 · πx 1
∫ ∫
a0 = (|x | + 2x) cos dx = (|x | + 2x) · 1dx
3 −3 3 3 −3
3 3 3
1 2 2 1 2
∫ ∫
= |x |dx = xdx = x =3
3 −3 3 0 3 2 0
Then for an and bn we need to perform integration by parts (or tabular integration) to get:
where the left and right side agree (except for a few points) on the interval [−3, 3].
111
Fourier series of a piecewise function
Example 6.1.4 Find the Fourier series of
(
−2, if − 5 < x < 0
f (x) =
2, if0 ≤ x ≤ 5
Solution) The graph is almost same as
except for the fact that the two open circles of the bottom line should be closed. However, Fourier
series only gives equality away from jump discontinuity and also possibly away from the end
points, it doesn’t matter whether we are looking for Fourier series with the closed circle or the
open circle.
Note that this function is an odd function. Therefore we know that an ’s vanish. So we will
only compute bn .
Note that we could have simplified the calculation by just using the fact that f (x)sin nπx
5 is an
even function and therefore we can just double the integration on the half interval [0, 5]. However,
the above integration was performed to demonstrate how to deal with piecewise functions.
112
6.2 Fourier Sine Series and Cosine Series
Even and Odd extension of a function
Example 6.2.1 Find the even extension of the function f (x) = x which is defined on [0, 4].
An even function is symmetric with respect to the y-axis, so if the function is extended as an
even function, its graph would be:
As you can see, the newly added left part is a straight line whose slope is (-1) and whose
y-intercept is zero. So it is easy to see that the function on the left is −x.
Thus we conclude that the extended function is a piecewise function whose formula is given
as: (
x, 0≤x ≤4
f (x) =
−x , −4 ≤ x < 0
satisfies the requirement f (x) = f (−x). Check that this formula indeed gives you the same
answer for the example in the beginning. Let’s apply this to a second example:
113
Example 6.2.2 f (x) = e x on [0, 3]. What is the even extension?
Answer: (
ex , 0 ≤ x ≤ 3
f (x) = −x
e , −3 ≤ x < 0
Odd extension
A function is odd if it satisfies f (x) = −f (−x). If a function f (x) is defined on the interval [0, L],
you can extend it to an odd function by formula:
(
f (x) , 0 ≤ x ≤ L
f (x) =
−f (−x) , −L ≤ x < 0
−f (−x) , −L ≤ x < 0
agrees with the original function f (x) on (0, L] but may not agree at x = 0.
114
where the coefficients are obtained by:
L
2 nπ
∫
an = f (x) cos xdx
L 0 L
Note that in the above formula for an , we are using the fact that f (x) cos nπ
L x is even, to take the
double value of integration to be on [0, L].
If f (x) is an odd function, then the cos terms and also the a20 term (being even functions)
won’t be there and you’ll only get the sin terms, i.e.:
∞
Õ nπ
f (x) = bn sin x
n=1
L
Note that in the above formula for bn , we are using the fact that f (x) sin nπ
L x is even, to take the
double value of integration to be on [0, L].
115
nπx 5 2 nπx 5
!
2 5
= (x + 2) · sin + cos
5 nπ 5 nπ 5
0
10
= ((−1)n − 1)
n2π 2
As done above, it is important to find a 0 separately. Also once you get the coefficients, you
should write down the sum for the final answer:
∞
9 Õ 10 nπ
f (x) = + 2 2
((−1)n − 1) cos x
2 n=1 n π 5
116
f (x) = 2x + 1 on [0, 4]
Solution) Since we are taking an even extension, we are to use:
∞
a0 Õ nπ
f (x) = an cos x
2 n=1 L
16
= 2 2 ((−1)n − 1)
n π
Thus the answer is:
∞
16 nπ
((−1)n − 1) cos x
Õ
2x + 1 = 5 + 2 2
n=1
n π 4
which holds true on the interval [0, 4].
Solution) Since the function is defined on the interval (0, 4], we have L = 4. When computing the
integral of a piecewise function, it is necessary to split the integral.
2 4
∫ 2 ∫ 4
1
∫
a0 = f (x)dx = f (x)dx + f (x)dx
4 0 2 0 2
1 2 1 4
∫ ∫
= 3dx + x + 1dx
2 0 2 2
2 4
3 1 2 1
= x + x + x = 7
2 4 2
0 2
117
4 2 4
2 nπx 1 nπx 1 nπx
∫ ∫ ∫
an = f (x) cos dx = f (x) cos dx + f (x) cos dx
4 0 4 2 0 4 2 2 4
The first integral is:
nπx 2
2
1 nπx 1 4
∫
3 cos dx = · 3 · sin
2 4 2 nπ 4
0
0
6 nπ
= sin
nπ 2
The second integral is:
4
1 nπx
∫
(x + 1) cos dx
2 2 4
nπx 4
1 4 nπx 16
= (x + 1) · sin + 2 2 cos
2 nπ 4 n π 4
2
6 nπ 8 8 nπ
=− sin + 2 2 (−1)n − 2 2 cos
nπ 2 n π n π 2
Adding the two we get:
6 nπ 6 nπ 8 8 nπ
an = sin − sin + 2 2 (−1)n − 2 2 cos
nπ 2 nπ 2 n π n π 2
8 n 8 nπ
(−1) − = cos
n2π 2 n2π 2 2
nπ
The cos 2 takes the value of 1, 0 or -1. So we should investigate it in two cases. If n is even,
then n = 2m for some m, which makes cos nπ m nπ
2 = cos(mπ ) = (−1) . If n is odd, then cos 2 = 0
(You should draw the graph of cos in order to convince yourself.)
For each case we have:
• Case n = 2m
8 8 nπ
an = (−1)n − cos
n2π 2 n2π 2 2
8 8
= 2 2
(−1)2m − (−1)m
(2m) π (2m)2 π 2
2 2
= 2 2 − 2 2 (−1)m
m π m π
2
= 2 2 (1 − (−1)m )
m π
• Case n = 2m − 1 (i.e. n is odd)
118
8 n 8 nπ
an = (−1) − cos
n2π 2 n2π 2 2
8
= 2 2 (−1)2m−1 − 0
n π
8
=−
(2m − 1)2 π 2
Answer:
∞
7 Õh 8 (2m − 1)πx 2 m 2mπx i
f (x) = + − cos + (1 − (−1) ) cos
2 m=1 (2m − 1)2 π 2 4 m2π 2 4
2 3 nπ 2 1 nπ 2 3 nπ
∫ ∫ ∫
bn = f (x) sin xdx = f (x) sin xdx + f (x) sin xdx
3 0 3 3 0 3 3 1 3
2 1 nπ 2 3 nπ
∫ ∫
= x sin xdx + 1 · sin xdx
3 0 3 3 1 3
1 3
2 3 nπ 9 nπ 2 3 nπ
= − x cos x + 2 2 sin x + − cos x
3 nπ 3 n π 3 3 nπ 3
0 1
1 3
2 nπ 6 nπ 2 nπ
= − x cos x + 2 2 sin x + − cos x
nπ 3 n π 3 nπ 3
0 1
2 nπ 6 nπ 2 2 nπ
=− cos + 2 2 sin −0− (−1)n + cos
nπ 3 n π 3 nπ nπ 3
119
6 nπ 2
= sin − (−1)n
n2π 2 3 nπ
Since we are told to find 3 non-zero terms, we compute b1 , b2 and b3 using the above formula:
√
6 π 2 1 3 3 2
b1 = 2 sin − (−1) = 2 +
π 3 π π π
√
6 2π 2 2 3 3 1
b2 = 2 sin − (−1) = −
4π 3 2π 4π 2 π
2 2
b3 = 0 − (−1)3 =
3π 3π
Thus the first three non-zero terms of the sine series is:
√ √
3 3 2 π 3 3 1 2π 2 3π
f (x) = 2
+ sin x + 2
− sin x + sin x + · · ·
π π 3 4π π 3 3π 3
Note: The conditions y(0) = 0 and y(2) = 0 are endpoint value conditions, and the differential
equation y 00 + 4y = t is thought to hold for 0 < t < 2 (i.e. in between the endpoints 0 and 2). The
goal is to find a continuous function that satisfies both the differential equation and the endpoint
conditions. Although this can be solved using the usual method of undetermined coefficients, let
us learn how to find the solution using Fourier series. (The advantage of Fourier series solution
is that the right side, t, could be replaced by any piecewise continuous function and we can still
solve it.)
Solution) Notice that the trigonometric functions that satisfy y(0) = 0 and y(2) = 0 are sin
functions of certain period. More specifically, sin 2nπ
2 t satisfies the requirement. The idea of the
solution is to put y as sum of these functions, i.e.
∞
Õ 2nπ
y= cn sin t
n=1
2
Then it will automatically satisfy the endpoint value conditions. Let us differentiate in order to
plug this into the differential equation.
∞
Õ n2π 2 2nπ
y =
00
− cn sin t
n=1
4 2
120
∞
n2π 2 2nπ
Õ
= − + 4 cn sin t
n=1
4 2
Now this has to be matched with the right side, which is t. But when you set them equal, the
conclusion is not clear:
∞
n2π 2 2nπ
Õ
− + 4 cn sin t =t
n=1
4 2
In other words, we can’t figure out what cn should be. Thus we take the Fourier sine series of t.
∞
Õ 2nπ
t= bn sin t
n=1
2
where
2 2
2 2nπ 2nπ 4 nπ
∫
−2
bn = t sin t= t cos t + 2 2 sin t
2 0 2 nπ 2 n π 2
0
4
=− (−1)n
nπ
So we have:
∞
4 2nπ
(−1)n sin
Õ
t= − t
n=1
nπ 2
Now, since both sides look similar, we can compare the coefficients to deduce the following:
∞ ∞
n2π 2 2nπ 4 2nπ
− (−1)n sin
Õ Õ
− + 4 cn sin t= t
n=1
4 2 n=1
nπ 2
n2π 2 4
− + 4 cn = − (−1)n
4 nπ
Solving this for cn gives us:
4
− nπ (−1)n
cn = 2 2
− n 4π + 4
Plugging this back into y gives us the answer:
∞ 4 n
Õ − nπ (−1) 2nπ
y= sin t
n2 π 2
n=1 − 4 + 4
2
121
Partial Differential Equations
7
7.1 Heat Equation
Introduction to heat equation
Consider the following situation. A thin, long rod is inside some insulation with exposed ends,
so that the only way heat can enter or escape is at the end points. Our goal is to figure out the
temperature distribution u(x, t) which means ”temperature at point x at time t”. There is some
initial temperature distribution at time t = 0. Let’s call it u 0 (x) so that u(x, 0) = u 0 (x). (This
condition is called the initial condition .) Also let the two endpoints touch ice so that the end
points are always kept to 0 degree Celsius. If we denote x = 0 for the position of the left end
point and x = L for the position of the right end point (which means that the length of the rod is
L), then having temperature of zero at those points at all time t means u(0, t) = 0 and u(L, t) = 0.
(This condition is called the boundary conditions .) Here is a picture to help you understand
the situation:
In order to figure out the temperature distribution u(x, t), we need an equation that governs the
change of temperature, which is called the heat equation . Usually, differential equation text
books start with the Newton’s cooling law to derive the one dimensional heat equation. However,
we will approach it more heuristically by thinking of an example. Let’s say the entire rod was
122
at the room temperature 20 degrees Celsius (=68 degrees Fahrenheit). If the rod is 2 meters long,
the graph of the temperature distribution will be as follows:
Then, attach the ice at the end points so that the temperature at end points quickly drop to 0
degree Celsius. A brief moment later, we will have the following temperature distribution:
Then as time goes by, the temperature will overall drop down to zero:
What kind of behavior can we deduce from the above discussion? Well, we see that whether
the temperature rises or falls at a point is related to the shape of the temperature graph near by.
123
For example, at the center of the last picture above, the temperature is bent like an upside down
U, and the temperature at the center decreases. In mathematics, when a graph is bent like a U,
we say the graph is concave up. So we know that the increase/decrease of the temperature at a
point is related to concavity. Recall that concavity is measured by the second derivative. A way
to remember the relationship between concave up/down and the sign of the second derivative is
shown in the table below:
The smiley face ’smiles’ in a ’concave up’ way and the frowning face ’frowns’ in a ’concave down’
way. Not only that, you should convince yourself that the more concave the function is at a point,
the faster the rate of change of the function will be. What is ’rate of change’ in math? It is the
derivative. But we have to be careful because there is two kinds of derivative at work. Since u(x, t)
is a function of x and t, there is the derivative with respect to x and a derivative with respect to
t. ’Rate of change’ is the change of temperature with respect to time, so it should be thought as
the partial derivative
∂u
(x, t).
∂t
On the other hand, ’concavity’ is the shape of the temperature when drawn with the x-axis (which
are the three temperature graphs we have seen before). So concavity should be measured by
derivative of x as:
∂ 2u
(x, t)
∂x 2
Now, the previous discussion says that the two partial derivatives are related in the following
∂ 2u
way: If ∂x 2 (x, t) is positive, then the temperature graph is shaped like a U, and the temperature at
the point will increase, i.e. ∂u ∂t (x, t) is positive. Same kind of statement can be said for the negative
case. Not only that, since more concavity means faster rate of change, we can write:
∂u ∂ 2u
(x, t) ∝ 2 (x, t)
∂t ∂x
which is equivalent to saying:
∂u ∂u
(x, t) = k (x, t)
∂t ∂t
for some positive constant k. The above is what we call the heat equation. Since this differential
equation involves partial differentiation, it is an example of a partial differential equation. It is one
of the simplest partial differential equations, and it can be solved by hand in many easy cases.
124
Putting everything together, the heat equation along with the initial condition and boundary
conditions, we have the full set of equations describing the behavior of the temperature distribu-
tion of a 1-dimensional rod:
∂u ∂ 2u
(x, t) = k 2 (x, t) , 0 < x < L, t > 0
∂t ∂x
u(x, 0) = u 0 (x) , 0 < x < L
u(0, t) = u(L, t) = 0 , t > 0
There is another type of boundary condition that is also used often. That is the case when the
end points are insulated. In that case we say that there is no change in the temperature across
the end points, i.e. the derivative in the x direction is zero at the end points. This is written as:
∂u ∂ 2u
(x, t) = k 2 (x, t) , 0 < x < L, t > 0
∂t ∂x
u(x, 0) = u 0 (x) , 0 < x < L
∂u ∂u
(0, t) = (L, t) = 0 , t > 0
∂x ∂x
As you can see, the last line has the boundary condition that describes insulated end points. Such
boundary conditions are called Neumann boundary conditions.
When instead of insulation, the end points are given some fixed value (such as zero as we
have in the previous one) for all time, we call it Dirichlet boundary conditions.
125
can’t walk. Ones that look like cars would have handles that don’t do anything.) The reason for
this discussion will be apparent later in the solution of the heat equation.
Now, let’s actually solve one of these equations. Say we have:
Example 7.1.1
∂u ∂ 2u
(x, t) = 3 2 (x, t) , 0 < x < L, t > 0
∂t ∂x
u(0, t) = u(2, t) = 0 , t > 0
u(x, 0) = 4 , 0 < x < L
The above means that the end points are set to 0 degrees Celsius and the initial temperature
was 4 degrees everywhere. k value is given as 3, and it controls how rapidly the heat spreads.
Step 1 We solve this by using a method called separation of variables. The idea is to
(A) only consider the heat equation and the boundary conditions at first (i.e. ignore the u(x, 0) = 4
for a moment.)
(B) assume that the solution takes the form of:
Which means that the function u(x, t) is a product of a function of x and a function of t. This
is where the name ’separation of variables’ come from. Let’s compute the partial derivatives in
order to plug into the heat equation:
∂u
(x, t) = f (x) · д0(t)
∂t
∂ 2u
(x, t) = f 00(x) · д(t)
∂x 2
Then the heat equation becomes:
∂u ∂ 2u
(x, t) = 3 2 (x, t) → f (x) · д0(t) = 3f 00(x) · д(t)
∂t ∂x
Here is the IMPORTANT STEP. You should divide both sides by k f (x)д(t). (Many books, and some
of my later sections use X (x) and T (t) instead of f (x) and д(t). In that case you would divide both
sides by kXT .) The result is:
f (x) · д0(t) 3f 00(x) · д(t) д0(t) f 00(x)
= → =
3f (x)д(t) 3f (x)д(t) 3д(t) f (x)
Now comes an interesting argument. First, you see that the above equation involves two inde-
pendent variables x and t. So let’s assign it as a function H (x, t):
д0(t) f 00(x)
H (x, t) := =
3д(t) f (x)
д 0 (t)
But then, looking at H (x, t) = 3д(t) , we see that H (x, t) does not depend on x. On the other hand,
f 00 (x)
looking at H (x, t) = f (x) , we see that H (x, t) does not depend on t. What does that say about
126
H (x, t)? It means that it’s a constant function, as we have discussed in the beginning of this
section. For a reason that will be apparent later, we will call that constant −λ. So we have the
д 0 (t) f 00 (x) д 0 (t) f 00 (x)
following equality: −λ = 3д(t) = f (x) This leads to two equations: −λ = 3д(t) and −λ = f (x) This
last equation turns into
f 00(x) + λ f (x) = 0
This looks familiar. When we solved the eigenfunction-eigenvalue problem for Endpoint Value
Problem, we had x 00(t) + λx(t) = 0. The above question is same, if we take f to be x and x to be
t. In fact, the whole point of eigenfunction-eigenvalue problem was that it was needed in order
to solve PDEs. We will come back to this in a moment.
Step 2 Next we consider the boundary conditions. Recall that u(x, t) = f (x)д(t). The bound-
ary condition u(0, t) = u(2, t) = 0 then turns into f (0)д(t) = f (2)д(t) = 0. We see from f (0)д(t) =
0 that either f (0) = 0 or д(t) = 0. However, д(t) = 0 leads to u(x, t) = f (x)д(t) = f (x) · 0 = 0.
So that part is not that interesting. (It is certainly a solution satisfying the heat equation and the
boundary conditions, albeit something anyone can easily see without having to do all this work.
This solution is called the trivial solution.) Therefore, in order to have useful solutions, we must
have f (0) = 0. For the same reason, we use f (2)д(t) = 0 to get f (2) = 0.
Step 3 Now take the previous equation f 00(x)+λ f (x) = 0 along with the boundary conditions
f (0) = 0 and f (2) = 0. This is indeed the eigenfunction-eigenvalue problem with end point
conditions that we have seen before. We know that when the end points are set to zero, the only
2 2
kind of non-trivial solutions are f (x) = sin nπ 2 x with the eigenvalue λ = n 4π for n = 1, 2, 3, · · · .
д 0 (t) 2 2
For д(t), we had −λ = 3д(t) . Plugging in λ = n 4π and multiplying 3 both sides, we have:
n 2 π 2 д0(t)
−3 =
4 д(t)
Integrating both sides, we have:
n2π 2
−3 t + C̃ = ln д(t)
4
Then exponentiating, we have:
n2 π 2
д(t) = Ce −3 4 t
Since we have found out what f (x) and д(t) could be, we know what u(x, t) is:
nπ n2 π 2
u(x, t) = C f (x)д(t) = C sin x e −3 4 t
2
What you should notice in the above solution is that this is not just one solution but infinitely
many solutions for n = 1, 2, 3, · · · .
Step 4 Think about how in homogeneous differential equations when we could take linear
combinations of solutions to form a solution. The same is true for our case too. (You can check this
by assuming that u 1 (x, t) and u 2 (x, t) both satisfy the heat equation and the boundary conditions
and then trying to plug in u(x, t) = c 1u 1 (x, t)+c 2u 2 (x, t). ) The only difference in this case would be
that since we have infinitely many solutions, adding them all gives us the infinite series, written
as:
∞ nπ 3n 2 π 2
x e− 4 t
Õ
u(x, t) = cn sin
n=1
2
127
You can think of the above as the general solution for our heat equation. The only thing we
haven’t included in our consideration is the initial condition u(x, 0) = 4. We will do so in the next
section.
Step 5 We now have:
∞ nπ 3n 2 π 2
x e− 4 t
Õ
u(x, t) = cn sin
n=1
2
To match this with u(x, 0) = 4, we will plug in t = 0 into the above equation. Since e 0 = 0, the
entire exponential factor disappears when t = 0 is plugged in. Setting it equal to 4:
∞
Õ nπ
4= cn sin x
n=1
2
This looks familiar. This is the Fourier Series, or Fourier sine series to be exact. Since we need
u(x, 0) = 4 to hold when 0 < x < 2, we need to find the Fourier sine series of 4 in order to make
sense out of the above equation. We compute:
2 2 nπ
∫
bn = 4 · sin x dx
2 0 2
nπ 2
2
=4· − cos x
nπ 2
0
8
=− ((−1)n − 1)
nπ
8
= (1 − (−1)n )
nπ
Step 6 Okay, so we now know that:
∞
8 nπ
(1 − (−1)n ) sin
Õ
4= x
n=1
nπ 2
from the Fourier sine series. If you look closely, you see that in place of cn we had before, we now
8
have nπ (1 − (−1)n ). Thus we found that:
8
cn = (1 − (−1)n )
nπ
Plugging this back into the u(x, t) formula in the beginning of this section, we get the final answer:
∞ nπ 3n 2 π 2
x e− 4 t
Õ
u(x, t) = cn sin
n=1
2
128
Plot of the heat equation
Using a computer one can plot the graph of
∞ nπ n2 π 2
x e −3 4 t
Õ
u(x, t) = cn sin
n=1
2
For example, I used a programming language called python, with the pylab package, to draw the
graph of u(x, t) for various values of t :
As you can see, the temperature will eventually become zero everywhere.
For those who are curious to know, here is the python (with pylab) code used to generate the
above graph. (People who know matlab will find that it is very similar to matlab)
from p y l a b i m p o r t ∗
x= l i n s p a c e ( 0 , 2 , 1 0 0 1 )
t=linspace (0 ,1 ,1001)
fs =zeros ((1001 ,1001))
f o r n in range ( 1 , 1 0 0 0 , 2 ) :
fs [n , : ] = 1 6 ∗ sin ( n∗ pi ∗ x / 2 ) / ( n∗ pi )
T= z e r o s ( ( 1 0 0 1 , 1 0 0 1 ) )
f o r n in range ( 0 , 1 0 0 1 ) :
T [ n , : ] = exp ( − 3 ∗ n ∗ n ∗ p i ∗ p i ∗ t / 4 )
u= d o t ( t r a n s p o s e ( T ) , f s )
plot (x , u [ 1 , : ] )
plot (x , u [100 ,:])
plot (x , u [200 ,:])
plot (x , u [300 ,:])
plot (x , u [400 ,:])
plot (x , u [500 ,:])
129
Heat Equation with Neumann boundary conditions
Example 7.1.2 Solve the heat equation with Neumann boundary conditions:
Since T (t) = 0 would give a trivial solution, we must have T (t) , 0 and the above forces X 0(0) =
X 0(8) = 0. Then we have the following eigenfunction-eigenvalue problem:
X 00(x)
−λ = → X 00(x) + λX (x) = 0
X (x)
with the end point condition X 0(0) = X 0(8) = 0. We have solved such problems before. When
the end point conditions say the derivative is zero, then we have the cosine functions as the
eigenfunctions. We have
nπ n2π 2
X (x) = cos x , λ=
8 64
T 0 (t)
for n = 0, 1, 2, 3, · · · . Also, since −λ = 7T (t) , solving it gives us: T (t) = Ce −7λt . Using the eigenval-
ues we found, we obtain:
7n 2 π 2
T (t) = Ce − 64 t
130
We have the solutions: nπ 7n 2 π 2
u(x, t) = C cos x e − 64 t
8
for n = 0, 1, 2, 3, · · · . We move on to take the linear combination of these solutions make this into
a series:
∞ nπ 7n 2 π 2
x e − 64 t
Õ
u(x, t) = cn cos
n=0
8
nπ 7n 2 π 2
Now, when n = 0, we have cos 64 t = cos 0 · e 0 = 1. Therefore we just have c 0 by itself
8 x e−
if n = 0. So we have:
∞ nπ 7n 2 π 2
x e − 64 t
Õ
u(x, t) = c 0 + cn cos
n=1
8
Let’s apply the initial condition
and
8
2 nπ
∫
an = (−x + 6) cos x dx
8 0 8
nπ 8
2 8 nπ 64
= (−x + 6) · sin x − 2 2 cos x
8 nπ 8 n π 8
0
16 16
2
=−
2
((−1)n − 1) = 2 2 (1 − (−1)n )
n π n π
Thus the Fourier Cosine law tells us:
∞
4 Õ 16 n
nπ
−x + 6 = + (1 − (−1) ) cos x
2 n=1 n2 π 2 8
4 16
Comparing, we see that c 0 = 2 = 2 while cn = n2 π 2
(1 − (−1)n ). Thus the general solution is:
∞
16 n
nπ 7n 2 π 2
x e − 64 t
Õ
u(x, t) = 2 + 2 2
(1 − (−1) ) cos
n=1
n π 8
131
Heat equation solution that do not require Fourier Series
Think about the following example:
Example 7.1.3
∂u ∂ 2u
(x, t) = 0.2 2 (x, t) , 0 < x < L, t > 0
∂t ∂x
u(0, t) = u(3, t) = 0 , t > 0
u(x, 0) = cos(6πx) − cos(2πx) + 7 , 0 < x < L
Since T (t) = 0 would give a trivial solution, we must have T (t) , 0 and the above forces X 0(0) =
X 0(3) = 0. Then we have the following eigenfunction-eigenvalue problem:
X 00(x)
−λ = → X 00(x) + λX (x) = 0
X (x)
with the end point condition X 0(0) = X 0(3) = 0. We have solved such problems before. When
the end point conditions say the derivative is zero, then we have the cosine functions as the
eigenfunctions. We have
nπ n2π 2
X (x) = cos x , λ=
3 9
T 0 (t)
for n = 0, 1, 2, 3, · · · . Also, since −λ = 0.2T (t) , solving it gives us: T (t) = Ce −0.2λt . Using the
eigenvalues we found, we obtain:
0.2n 2 π 2
t
T (t) = Ce − 9
132
We have the solutions: nπ 0.2n 2 π 2
u(x, t) = C cos x e− 9 t
3
for n = 0, 1, 2, 3, · · · . Just as before we take the linear combination:
∞ nπ 0.2n 2 π 2
x e− 9 t
Õ
u(x, t) = c 0 + cn cos
n=1
3
However when we look at the initial condition u(x, 0) = cos(6πx) − cos(2πx) + 7 and set it
equal with t = 0:
Õ∞ nπ
cos(6πx) − cos(2πx) + 7 = c 0 + cn cos x
n=1
3
We see that already we can make c 0 = 7 and then the two other terms cos(6πx) − cos(2πx) are
already sums of cosines. Let’s try to find out when nπ nπ
3 = 6π or when 3 = 2π . We have n = 18
and n = 6. In other words, if we take c 0 = 7, c 6 = −1, c 18 = 1 and set all other cn as zero then
the above equality holds. This is a case where we did not have to find the Fourier cosine series
in order to find the coefficients, because the initial condition is already given to us as a Fourier
cosine series.
Plugging in c 0 = 7, c 6 = −1, c 18 = 1 and set all other cn as zero, we have the answer:
2 2t
u(x, t) = 7 − cos (2πx) e −0.8π t + cos (6πx) e −7.2π
133
What does this mean? It means that the force is propotional to the bending, i.e. concavity:
∂ 2u
F∝ (x, t)
∂x 2
Here, u(x, t) has the meaning of ’vertical displacement of the string at position x at time t ’. Since
F = ma, we know that above statement means acceleration is proportional to concavity too:
∂ 2u
a∝ (x, t)
∂x 2
But then, acceleration is the second derivative in time, so we write:
∂ 2u ∂ 2u
(x, t) ∝ (x, t)
∂t 2 ∂x 2
When we write this as an equation, we put the proportionality constant as c 2 to write:
∂ 2u 2
2∂ u
(x, t) = c (x, t)
∂t 2 ∂x 2
The above is a partial differential equation called the wave equation . The reason for using c 2
for the proportionality constant is because the solution of the equation becomes easier to write
if we put it that way. (It also stresses the fact that c 2 is positive.)
Since it is a second order differential equation in time, we need to give two initial conditions:
u(x, 0) = f (x)
and
∂u
(x, 0) = д(x)
∂t
The first gives the initial configuration of the string and the second gives the initial velocity of
the string at each point x.
Also, the fixed end points of the string means that the vertical displacement is always zero at
the end points. So we have:
u(0, t) = u(L, t) = 0
where L is the length of the string.
You can also have Neumann boundary conditions at the end points. Such conditions describe
wind pipe instruments like flute or pipe organ. For wind instruments, u(x, t) has the meaning of
pressure (compared to the atmospheric pressure) at point x at time t. If the wind instrument has
open ends, then the pressure is zero (i.e. same as the atmospheric pressure) at the end points, so
we have Dirichlet boundary conditions (set to zero). However, if one end is closed (as instruments
like flute is, where the end close to the mouth piece is closed) then the derivative of u(x, t) should
be set to zero at that end point, giving us the Neumann boundary condition at that end.
To summarize, a vibrating string or a wind instrument with both ends open follows the partial
differential equation and boundary conditions:
∂ 2u ∂ 2u
(x, t) ∝ (x, t) , 0 < x < L, t > 0
∂t 2 ∂x 2
∂u
u(x, 0) = f (x) , (x, 0) = д(x) , 0 < x < L
∂t
u(0, t) = u(L, t) = 0 , t > 0
134
Solving Wave Equation using Fourier Series
Suppose a string 4 in. long is tied at both ends. We pulle the string at the middle to make the
following shape:
Then we release this string, meaning that the initial velocity of every part of the string will be
zero. This can be written as yt (x, 0) = 0 for all x in the interval [0, 4]. We still have to provide a
value for c 2 as explained in the previous section. Let’s assume that we are given the value c 2 = 9.
Solution) We use the separation of variables method again. If you do not understand this
method, please review how the heat equation was solved by separation of variables.
step 1. We will find set of solutions that satisfy the wave equation and the boundary condi-
tions.
Put
y(x, t) = X (x)T (t)
Plugging this into the wave equation, we get
135
The equation XX (x)
00(x)
= −λ gives us X 00(x) + λX (x) = 0. Now, observe that the boundary conditions
y(0, t) = y(4, t) = 0 in separation of variables is X (0)T (t) = X (4)T (t) = 0. Since T (t) = 0 leads to
trivial solution, we are forced to say X (0) = X (4) = 0. Thus we have the eigenvalue-eigenfunction
problem:
X 00(x) + λX (x) = 0 , X (0) = X (4) = 0.
The solution of this problem is (see the chapter Endpoint Value problems for details of how eigen-
functions and eigenvalues are found)
nπx
X (x) = C sin , n = 1, 2, 3, 4, · · ·
4
with
n2π 2
λ= , n = 1, 2, 3, 4, · · ·
16
Then T9T (t)
00(t)
= −λ can now be written as (after multiplying 9T (t),
9n2 π 2
T 00(t) = − T (t).
16
2 2
After moving everything to the left as T 00(t) + 9n16π T (t) = 0, we can write the characteristic
equation as
9n 2 π 2
r2 + =0
16
3nπ
r =± i
4
3nπ 3nπ
T (t) = C 1 cos t + C 2 sin t
4 4
Now, y(x, t) = X (x)T (t) becomes
nπx 3nπ 3nπ nπx 3nπ nπx 3nπ
y(x, t) = C sin C 1 cos t + C 2 sin t = C·C 1 sin cos t+C·C 2 sin sin t
4 4 4 4 4 4 4
Since C · C 1 and C · C 2 are arbitrary constants, we can reassign the name of them as An and
Bn respectively. The reason for the subscript n is because the above is actually an infinite set of
solutions for n = 1, 2, 3, 4, · · · . If we can take linear combinations of all of them, then we end up
with the following series:
∞ ∞
Õ nπx 3nπ Õ nπx 3nπ
y(x, t) = An sin cos t+ Bn sin sin t
n=1
4 4 n=1
4 4
The linear combination satisfies the wave equation ytt = 9yxx because it is linear. You should
also check that it satisfies the boundary condition y(0, t) = y(4, t) = 0, , (t > 0).
step 2. We now put the initial value problem into consideration. First we consider y(x, 0) =
4 − |2 − x | , (0 < x < 4). Setting t = 0 in the solution obtained in step 1, we have:
∞ ∞
nπx 3nπ nπx 3nπ
Õ Õ
An sin cos ·0 + Bn sin sin · 0 = 4 − |2 − x |
n=1
4 4 n=1
4 4
136
The second summation becomes zero because sin 0 = 0. So using cos 0 = 1 we get:
∞
Õ nπx
An sin = 4 − |2 − x |
n=1
4
This means that in order to find out the values of An ’s, we need to find the Fourier sine
expansion of the right side. We use the Fourier sine formula with L = 4 to get
2 4 nπ
∫
bn = (4 − |2 − x |) sin xdx
4 0 4
∫ 2 ∫ 4
1 nπ 1 nπ
= (4 − |2 − x |) sin xdx + (4 − |2 − x |) sin xdx
2 0 4 2 2 4
∫ 2 ∫ 4
1 nπ 1 nπ
= (4 − (2 − x)) sin xdx + (4 − (x − 2)) sin xdx
2 0 4 2 2 4
∫ 2 ∫ 4
1 nπ 1 nπ
= (2 + x) sin xdx + (6 − x) sin xdx
2 0 4 2 2 4
Perform the integration by parts (or tabular integration) for each integral to get:
2 4
1 4 nπ 16 nπ 1 4 nπ 16 nπ
= − (2 + x) cos x + 2 2 sin x + − (6 − x) cos x − 2 2 sin x
2 nπ 4 n π 4 0 2 nπ 4 n π 4 2
1 16 nπ 16 nπ 8 1 8 16 nπ 16 nπ
= − cos + 2 2 sin − − + − cos nπ − − cos − 2 2 sin
2 nπ 2 n π 2 nπ 2 nπ nπ 2 n π 2
8 nπ 8 nπ 4 4 8 nπ 8 nπ
=− cos + 2 2 sin + − cos nπ + cos + 2 2 sin
nπ 2 n π 2 nπ nπ nπ 2 n π 2
16 nπ 4 4
= 2 2 sin + − cos nπ
n π 2 nπ nπ
Now, depending on whether n is even (which can be expressed as n = 2m for some integer
m) or odd (which can be expressed as n = 2m + 1 for some integer m), the above can be further
simplified using sin 2mπ
2 = sin mπ = 0 and sin
(2m+1)π
2 = (−1)m . If n = 2m, we get:
16 nπ 4 4 16 2mπ 4 4 4 4
2 2
sin + − cos nπ = 2 2 sin + − cos 2mπ = 0 + − =0
n π 2 nπ nπ n π 2 nπ nπ nπ nπ
If n = 2m + 1, we get:
16 nπ 4 4 16 (2m + 1)π 4 4 16 8
2 2
sin + − cos nπ = 2 2 sin + − (−1) = 2 2 (−1)m +
n π 2 nπ nπ n π 2 nπ nπ n π nπ
If we replace all the remaining n by 2m + 1 we have
16 8
2 2
(−1)m +
(2m + 1) π (2m + 1)π
137
The Fourier sine series we obtain is therefore
∞ ∞
nπ 16 8 (2m + 1)π
m
Õ Õ
4 − |2 − x | = bn sin x = 2 2
(−1) + sin x
n=1
4 m=0
(2m + 1) π (2m + 1)π 4
Note that the summation starts from m = 0 because then n = 2m + 1 = 2 · 0 + 1 = 1. Now this
means for An ’s :
∞ ∞
nπx Õ 16 8 (2m + 1)π
m
Õ
An sin = 2 2
(−1) + sin x
n=1
4 m=0
(2m + 1) π (2m + 1)π 4
This yields that if n is even, An = 0. For n odd, we have:
16 8
A2m+1 = 2 2
(−1)m +
(2m + 1) π (2m + 1)π
We still have the other initial condition yt (x, 0) = 0. If we differentiate (by t)
∞ ∞
Õ nπx 3nπ Õ nπx 3nπ
y(x, t) = An sin cos t+ Bn sin sin t
n=1
4 4 n=1
4 4
to get
∞ ∞
Õ 3nπ nπx 3nπ Õ 3nπ nπx 3nπ
y(x, t) = −An sin sin t+ Bn sin cos t.
n=1
4 4 4 n=1
4 4 4
Then
∞
Õ 3nπ nπx
y(x, 0) = Bn sin
n=1
4 4
Setting this equal to yt (x, 0) = 0 we have
∞
Õ 3nπ nπx
Bn sin =0
n=1
4 4
Which suggests that we must find the Fourier sine series of the right side. However, the function
0 of the right side has the Fourier sine series 0. Hence
3nπ
Bn =0
4
and we end up with Bn = 0 for all n. So the second summation in y(x, t) completely disappears.
We finally arrive at the solution
∞
Õ nπx 3nπt
y(x, t) = An sin cos
n=1
4 4
∞
16 8 (2m + 1)πx 3(2m + 1)πt
m
Õ
= 2 2
(−1) + sin cos
m=0
(2m + 1) π (2m + 1)π 4 4
138
Laplace Transform Method
8
8.1 Laplace Transform and Its Inverse
Introduction to Laplace Transforms
A transform is an operation where it takes in a function and returns a function of a different
variable. Among such transforms, one transform called Laplace transform is used for solving
differential equations. It is defined as:
∫ ∞
L { f (t)} = e −st f (t)dt
0
where s is assumed to be bigger than some positive number (just a technicality assumed in order
to make the integral converge.) The best way to understand what the above formula is doing is
to try out an example: ∫ ∞
2
L t = e −st t 2dt
0
We use integration by parts (or tabular integration) to do the integration. The result is:
2 ∞
t −st 2t −st 2 −st
= − e − 2 e − 3e
s s s
0
2
s3
139
So the result is that
2
L t2 = 3
s
As you can see, a function of t was plugged in and the result is a function of s.
You can calculate similarly to get the following results:
1
L {t } =
s2
1
L {1} =
s
In general, you have:
n!
L {t n } =
s n+1
where n! means n! = 1 · 2 · 3 · · · n
There are results of the Laplace transform of other functions. These are provided to you as a
table shown below:
140
Take dv = f 0(t)dt
∫ so that v =
∫ f (t). Also u = e so that du = −se dt. Then the integration by
−st −st
If you look closely at the integral, you see that it is same as L { f (t)}. Therefore we have what we
wanted to prove:
L { f 0(t)} = −f (0) + sL { f (t)} = sF (s) − f (0)
Now, part of the above formula, L { f 0(t)} = −f (0)+sL { f (t)}, can be used to obtain a formula
for the second derivative:
Which becomes:
L { f 00(t)} = s 2 F (s) − f 0(0) − s f (0)
If you keep going, a general pattern emerges for the n-th derivative of f (x):
n o
L f (t) = s n F (s) − f (n−1) (0) − s f (n−2) (0) − s 2 f (n−3) (0) · · · − s n−1 f (0)
(n)
What the above formula shows us is that Laplace transform takes derivatives and changes it
into multiplication by s (and then some more…) .
This means that if you have a differential equation, taking the Laplace transform will make it
into an algebraic equation which does not involve any derivatives.
Here is an example:
Example 8.1.1
y 00 − 5y 0 + 4y = cos 3t , y(0) = 1, y 0(0) = 0
Take the Laplace transform of the left:
= L {y 00 } + −5L {y 0 } + 4L {y}
= s 2Y (s) − y 0(0) − sy(0) − 5 (sY (s) − y(0)) + 4Y (s)
= s 2Y (s) − 0 − s · 1 − 5sY (s) + 5 · 1 + 4Y (s) = (s 2 − 5s + 4)Y (s) − s + 5
141
For the right side, you can look at the Table of Laplace Transforms and find that:
s
L {cos at } =
s2 + a2
s
Thus the right side becomes s 2 +32
and the entire differential equation becomes:
s
(s 2 − 5s + 4)Y (s) − s + 5 =
s2 + 32
As you can see, there is no derivative involved in the above equation. It has been changed into
an algebraic equation. So solving it is really easy:
s
(s 2 − 5s + 4)Y (s) = s − 5 +
s2 + 32
s −5 s
Y (s) = + 2
s2 − 5s + 4 (s − 5s + 4)(s 2 + 9)
However, this is NOT what we want. We need y(t), and the above is the Laplace transform of
y(t). Thus we have to know how to undo the Laplace transform (also known as ”taking the inverse
Laplace transform”) That is major thing we will need to learn how to do.
Using tabular integration (see appendix if you do not know what tabular integration is.):
we find that:
1 −st a −st a2
∫ ∫
e −st
sin atdt = − e sin at − 2 e cos at − 2 e −st sin atdt
s s s
142
a2 1 −st a −st
∫ ∫
e sin atdt + 2
−st
e sin atdt = − e sin at − 2 e cos at
−st
s s s
a2 1 −st a −st
∫
1+ 2 e sin atdt = − e sin at − 2 e cos at
−st
s s s
2
s + a2 1 −st a −st
∫
e sin atdt = − e sin at − 2 e cos at
−st
s2 s s
s2 1 −st a 2 −st s a
∫
e −st
sin atdt = 2 2
− e sin at − 2 e cos at = − 2 2
e −st sin at − 2 e −st cos at
s +a s s s +a s + a2
Thus we have:
∞
∞
s a
∫
L {sin at } = e −st sin atdt = − 2 e −st sin at − 2 e −st cos at
s +a 2 s +a 2
0
0
Because of e −st , all terms become 0 when t = ∞ is plugged in. Also, using the fact that cos 0 = 1
and sin 0 = 0, we get:
∞
s a a
L {sin at } = − 2 e −st sin at − 2 e −st cos at = 2
s +a 2 s +a 2 s + a2
0
L e at f (t) = F (s − a)
Proof
∫ ∞
L e at f (t) = e −st e at f (t)dt
0
∫ ∞
= e −(s−a)t f (t)dt
0
∫∞
On the other hand, since F (s) = 0
e −st f
(t)dt, replacing s by s − a becomes:
∫ ∞
F (s − a) = e −(s−a)t f (t)dt
0
143
Then
s
L e bt cos at =
(s − b)2 + a 2
Note that some tables have a and b changed, which is not a big deal because they are just names:
s
L e at cos bt =
(s − a)2 + b 2
Similarly, we have:
a b
L {sin at } = → L e at sin bt =
s2 +a 2 (s − a)2 + b 2
and:
n! n!
L {t n } = → L e at t n =
s n+1 (s − a)n+1
s 2 + 4s + 13 = (s + 2)2 + 9 = (s + 2)2 + 32
144
So we must have:
5s + 3 = A(s + 2) + 3B
Since 5s + 3 = As + 2A + 3B, we must have 5 = A and 3 = 2A + 3B. The second one, after plugging
in A = 5, becomes −7 = 3B. Thus we have B = − 73
This means:
5s + 3 s +2 7 3
2 2
=5 2 2
−
(s + 2) + 3 (s + 2) + 3 3 (s + 2)2 + 32
Then taking the inverse transform:
s +2 7 −1 3
2 2
L −1
f rac5s + 3(s + 2) + 3 = 5L −1
− L
(s + 2)2 + 32 3 (s + 2)2 + 32
7
= 5e −2t cos 3t − e −2t sin 3t
3
y 00 − 3y 0 + 2y = 0 , y(0) = 1, y 0(0) = −1
s 2 − 3s + 2 Y (s) − s + 4 = 0
145
Multiply the denominator (s − 1)(s − 2) both sides:
s − 4 = A(s − 2) + B(s − 1)
One important thing to know is that in the above equation, the two sides are equal as functions
of s. It means the the following holds:
• If you expand the right side and organize it according to order of s, i.e. (A + B)s − 2A − B,
then the coefficients of the polynomial on the left s − 1 should match with that of the right,
i.e. (A + B) = 1 and −2A − B = −1.
• Even if you differentiate both sides, you still have equality as functions of s The above three
facts can be used in different ways to find the unknown constants A and B. For now, we
will just use the first property. However, there are some cases where the other ones are
useful.
−2 = 0 + B · 1
−3 = A · (−1) + 0
So we have A = 3. (Why did we choose the values 2 and 1 to plug in? That’s because it makes
some terms become 0 when plugged in. This will be the key reason for choosing values to plug
in.)
Now, we had:
s −4 A B
Y (s) = = +
(s − 1)(s − 2) s − 1 s − 2
So plugging in A = 3 and B = −2 gives us:
3 2
Y (s) = −
s −1 s −2
Which means
3 2
y(t) = L {Y (s)} = L
−1 −1
−
s −1 s −2
1 1
= 3L −1
− 2L −1
s −1 s −2
1
Then looking up the Laplace transform table, we find that L e at = s−a . Thus, applied to above,
we get:
y(t) = 3e t − 2e 2t
146
Why do we need Laplace Transforms?
You may not be impressed about Laplace transforms so far, because we already know how to solve
differential equations in the previous example. There are two reasons why Laplace transform is
used. First, it solves the initial value problem directly (you don’t have to go through the general
solution first). Second, in electrical circuits, often a switch that can turn on and off is included.
The differential equation for such circuits involve a function which quickly jumps from 0 V to
9 V (or whatever volts), and those functions are hard to work with just using the method of
undetermined coefficients. Laplace transforms are perfect for that situation. Indeed, if you see
differential equation being used in Electrical Circuits, you’ll mainly learn how to solve it using
Laplace transforms.
147
• When you have the denominator of the form (s − a)2 + b 2 , you should write the numerator
as B(s − a) + bC. This is different from what you are used to doing when using partial de-
composition for integration. You will find that this way makes the inverse Laplace transform
much easier.
Plug in s = −2 to get:
13 = A · 13
So A = 1. Now plug in s = −4 to get:
−15 = 9A − 6C
−24 = −6C
So we have C = 4.
Now, we see that there aren’t any more values to plug into s to make a few terms disappear.
In that case, you have two options. You can plug in some other value that’s easy to compute, say
s = 0, or you can expand both sides and organize both according to order of s and compare the
coefficients. The easiest method often turns out to be lazily doing the second. How? Well, just
pretend that you’ve done it. Think about what will be the highest order term on the left once
organized. (In our case it’s −s 2 ) Then try to think about what would be the corresponding term
on the right. (In our case, the s 2 will appear when (s + 4)2 is computed and also when (s + 4) gets
multiplied to (s + 2). This will result in As 2 + Bs 2 ) Using this, we obtain that:
(−1)s 2 = (A + B)s 2
Thus the coefficients of s 2 tell us that −1 = A + B, and we get B = −2. The result is the following
decomposition:
1 −2(s + 4) + 3 · 4
Y (s) = +
s +2 (s + 4)2 + 32
Which means
1 (s + 4) 3
y(t) = L {Y (s)} = L
−1 −1
−2 +4
s +2 (s + 4)2 + 32 (s + 4)2 + 32
1 (s + 4) 3
=L −1
− 2L −1
+ 4L −1
= e −2t − 2e −4t cos 3t + 4e −4t sin 3t
s +2 (s + 4)2 + 32 (s + 4)2 + 32
148
Solution) Applying the Laplace transform to the differential equation, we get:
2
s 2Y (s) − y 0(0) − sy(0) − 3 (sY (s) − y(0)) + 2Y (s) = 6
s2 + 22
12
s 2 − 3s + 2 Y (s) − s + 3 =
s2 + 4
Solving for Y (s) gives us:
s −3 12
Y (s) = + 2
s2 − 3s + 2 (s + 4)(s 2 − 3s + 2)
We see that the denominator s 2 − 3s + 2 can be factored as (s − 2)(s − 1), so we rewrite it as:
s −3 12
Y (s) = + 2
(s − 2)(s − 1) (s + 4)(s − 2)(s − 1)
Now take partial fraction decomposition.
s −3 12 As + 2B C D
+ 2 = 2 + +
(s − 2)(s − 1) (s + 4)(s − 2)(s − 1) s +4 s −2 s −1
Plug in s = 1 to get:
2
2 = −5D → D = −
5
Plug in s = 2 to get:
1
4 = 8C → C =
2
Plug in s = 0 to get:
3
0 = 4B − 4C − 8D → B = C + 2D = −
10
Now we’ve run out of good numbers to plug in. (One possibility is to plug in s = 2i so that
s 2 + 4 = 0. But that’s too confusing.) We then compare the highest order terms, i.e. the s 3 term
of left and right:
s 3 = As 3 + Cs 3 + Ds 3
9
This means 1 = A + C + D. Thus A = 1 − C − D = 10 and plugging all we have found back to the
decomposition gives us:
9 s 3 2 1 1 2 1
Y (s) = · 2 − · 2 + · − ·
10 s + 4 10 s + 4 2 s − 2 5 s − 1
Then looking at the Laplace Transform table tells us that:
9 3 1 2
y(t) = cos 2t − sin 2t + e 2t − e t
10 10 2 5
149
Example 8.2.4
y 00 + 4y 0 + 13y = 0 , y(0) = 0, y 0(0) = 2
(s 2 + 4s + 13)Y (s) = 2
2
Y (s) = 2
s + 4s + 13
Since s 2 + 4s + 13 = (s + 2)2 + 32 , we set up the partial fraction expansion as:
2 3A B(s + 2)
= +
s2 2
+ 4s + 13 (s + 2) + 32 (s + 2)2 + 32
Note that the 3 in the 3A is due to the term 32 in the denominator. Putting the partial fraction
expansion form this way makes it much easier to perform the inverse Laplace Transforms.
Multiplying (s + 2)2 + 32 both sides, we get:
2 = 3A + B(s + 2)
3 · 32
( )
2 −1 3 2
y(t) = L −1
= L = e −2t sin 3t
(s + 2)2 + 32 3 2
(s + 2) + 32 3
Example 8.2.5
y 00 + 4y 0 + 13y = e 2t , y(0) = 0, y 0(0) = 2
150
1
(s 2 + 4s + 13)Y (s) = 2 +
s −2
2 1
Y (s) = 2 +
s + 4s + 13 (s − 2)(s 2 + 4s + 13)
Since s 2 + 4s + 13 = (s + 2)2 + 32 , we set up the partial fraction expansion as:
2 1 3A B(s + 2) C
+ = + +
s 2 + 4s + 13 (s − 2)(s 2 + 4s + 13) (s + 2)2 + 32 (s + 2)2 + 32 (s − 2)
Note that the 3 in the 3A is due to the term 32 in the denominator. Putting the partial fraction
expansion form this way makes it much easier to perform the inverse Laplace Transforms.
Multiplying (s − 2)(s 2 + 4s + 13) both sides, we get:
0 = Bs 2 + Cs 2
1
Thus B + C = 0 and we have B = −C = − 25 .
46 1 1
Having A = 75 , B = − 25 and C = 25 means
46 1 1
2 1 3 · 75 − 25 (s + 2) 25
+ = + +
s 2 + 4s + 13 (s − 2)(s 2 + 4s + 13) (s + 2)2 + 32 (s + 2)2 + 32 (s − 2)
Therefore the answer is:
46 1 1
( )
3 · 75 − 25 (s + 2) 25 46 1 1
y(t) = L −1
2 2
+ 2 2
+ = e −2t sin 3t − e −2t cos 3t + e 2t
(s + 2) + 3 (s + 2) + 3 (s − 2) 75 25 25
Example 8.2.6
y 00 + 5y 0 + 4y = 17 cos t , y(0) = −1, y 0(0) = 2
151
17s
(s 2 + 5s + 4)Y (s) = −3 − s +
+1 s2
17s
(s + 4)(s + 1)Y (s) = −3 − s + 2
s +1
−3 − s 17s
Y (s) = + 2
(s + 4)(s + 1) (s + 1)(s + 4)(s + 1)
We set up the partial fraction expansion as:
−3 − s 17s As + B C D
+ 2 = 2 + +
(s + 4)(s + 1) (s + 1)(s + 4)(s + 1) s +1 s +4 s +1
Multiplying (s 2 + 1)(s + 4)(s + 1) both sides, we get:
(−3 − s)(s 2 + 1) + 17s = (As + B)(s + 4)(s + 1) + C(s 2 + 1)(s + 1) + D(s 2 + 1)(s + 4)
−s 3 = As 3 + Cs 3 + Ds 3
−3 − s 17s − 21 s + 52 1 − 72
+ = 2 + +
(s + 4)(s + 1) (s 2 + 1)(s + 4)(s + 1) s +1 s +4 s +1
Therefore the answer is:
1 5 7
( )
− 2 s + 2 1 − 1 5 7
y(t) = L −1 2
+ + 2 = − cos t + sin t + e −4t − e −t
s +1 s +4 s +1 2 2 2
152
Appendices
153
Integration by Parts using Tabular Integration
A
Integration by parts
Suppose you are trying to integrate (x 2 +3x −1)e 4x dx. You should recognize that this is a product
∫
of two functions
∫ which is not∫in any form of a chain rule. Thus one needs to use the integration by
parts formula udv = uv − vdu. Now, one chooses which one to be u and which one to be dv
by using the prescription LIATE (meaning Log, Inverse Trig, Algebraic, Trig, and Exponential).
In our example, (x 2 + 3x − 1) is Algebraic while e 4x is Exponential. Since A appears before E in
LIATE, we choose (x 2 + 3x − 1) as u and e 4x dx as dv (Notice that dx is included when choosing
it as dv.) Then we do the following computation:
u = (x 2 + 3x − 1) → du = (2x + 3)dx
1
dv = e 4x dx → v = e 4x
4
Thus the integration by parts reads:
1 4x 1 4x
∫ ∫
2 4x 2
(x + 3x − 1)e dx = (x + 3x − 1) e − e (2x + 3)dx
4 4
But now, you see there is another integral you have to compute which requires another integra-
tion by parts. After the second integration by parts, you get: (x 2 + 3x − 1) 41 e 4x − (2x + 3) 16
1 4x
e +
1 4x
32 e + C as the answer.
Tabular integration
Here is another way to do the integration by parts. This is called the tabular integration. Anything
you can solve using integration by parts can be solved using tabular integration and vice-versa.
Tabular integration is usually much quicker, and is easier to do once you learn it. First make a
table as follows:
154
Then choose one function to be differentiated and another to be integrated, using LIATE as before,
like this:
(What do we do if it never will become zero? Then you have to know when to stop - which is the
only tricky business with tabular.) Then you integrate the integration column:
155
The first column is for alternating signs. You always start with ’+’ then alternate:
Now you are ready to get the result. Basically, you draw curves from left to right which goes
straight to the right and then bends down one row, as these red curves shown below. Then the
last line has to be drawn from right to left, as the blue line is shown below. The red curves are
multiplied, and that’s all. However, the blue line must be multiplied and then integrated:
156
In the above, since the blue lines multiply to zero, the integral just becomes C (the integration
constant) and then if we add up everything we obtained, we get the same exact answer:
1 1 1
(x 2 + 3x − 1) e 4x − (2x + 3) e 4x + e 4x + C
4 16 32
Here is another example where the last line does not become zero. We aim to integrate
x ln xdx. Using LIATE, we determine that ln x has to be differentiated while x has to be in-
∫
Notice that ln x differentiates to x1 and none of its higher derivatives will become zero. So we have
to know when to stop, and the answer is ’we stop at the second line’. This is the only delicate part
about the tabular integration: knowing when to stop. For integrals involving ln x, we stop at the
second line. For integrals involving exponential and trigonometric function multiplied together, we
stop at the third line. Just memorize this and you will have no problem in applying the tabular
integration. Now let’s draw the curves and lines:
1 1
∫
x ln xdx = x 2 ln x − x 2 + C
2 4
157
Fourier series over [0, A]
B
In this book, all the formulas for Fourier series was written for a function defined on the interval
[−L, L], which means we get a periodic function of period 2L. This was done so that the introduc-
tion of Fourier sine series and Fourier cosine series will be easy and then use them to solve partial
differential equations. However, when Fourier series is used for signal processing and filtering,
often it is used for functions defined on [0, A]. Since this will have period A, the A can be thought
as 2L. The formula for Fourier Series looks quite different in this case:
2 A 2πnx
∫
an = cos f (x)dx,
A 0 A
2 A 2πnx
∫
bn = sin f (x)dx,
A 0 A
∞
a0 Õ 2πnx 2πnx
f (x) = + an cos + bn sin .
2 n=1 A A
Often students of engineering will learn the formula with [−L, L] and then see this other
formula and find it confusing. To add even more to this confusion, sometimes instead of a 0 , we
use a constant term c which is computed as:
1 A
∫
c= f (x)dx
A 0
and then we have:
∞
2πnx 2πnx
Õ
f (x) = c + an cos + bn sin .
n=1
A A
In other words, the a20 is replaced by c (so that there is no longer 2 in the denominator.)
Although these formulas may look different, the theory and derivation of these rules are done
in the exactly same manner as the Fourier series on [−L, L]
158
Index
amplitude, 56 order, 5
angular frequency, 55 ordinary point, 76
over-damped, 54
boundary condition, 118
phase angle, 56
characteristic equation, 41
phase shift, 58
critically damped, 54
power series, 64
differential equation, 5 power series solution, 68
pseudo-angular frequency, 58
eigenfunction, 98 pseudoperiod, 58
eigenvalue, 87, 98
eigenvector, 87 regular singular point, 76
endpoint value problem, 95 resonance, 60
even function, 102
separable equation, 9
exact, 30
separation of variables, 122
exact differential equation, 30
shifting of indices, 66
existence of solutions, 34
singular point, 76
first order linear differential equation, 18 spring-mass-dashpot system, 54
Frobenius Method, 71 standard form, of first order linear DE , 19
Frobenius solution, extended, 73 steady periodic solution, 59
Sturm-Liouville problem, 97
heat equation, 118 superposition, 34
indicial equation, 77 Table of Laplace Transforms, 136
initial condition of ODE, 6, 34 transient solution, 59
initial condition of PDE, 118 trivial solution, 123
initial value problem, 6
irregular singular point, 76 under-damped, 54
uniqueness of solutions, 34
linear combination, 34
linear operator, 39 variation of parameters, 50
linearly dependent, 35
linearly independent, 35 wave equation, 130
Wronskian, 38
mechanical vibration, 54
multiplier for first order linear DE, 20
natural angular frequency, 59
odd function, 102
159