Download numerical analysis 1 1.1 Differential equations and more Exercises Differential Equations in PDF only on Docsity!
1.1 Differential equations
In this chapter we are going to study differential equations, with particular emphasis on how to solve them with computers. We assume that the reader has previously met differential equations, and so we’re going to review the most basic facts about them rather quickly.
A differential equation is an equation in an unknown function, say y(x), where the equation contains various derivatives of y and various known functions of x. The problem is to “find” the unknown function. The order of a differential equation is the order of the highest derivative that appears in it.
Here’s an easy equation of first order:
y′(x) = 0. (1. 1 .1)
The unknown function is y(x) = constant, and so we have solved the given equation (1.1.1).
The next one is a little harder: y′(x) = 2y(x). (1. 1 .2)
A solution will, now doubt, arrive after a bit of thought, namely y(x) = e^2 x. But, if y(x) is a solution of (1.1.2), then so is 10y(x), or 49. 6 y(x), or in fact cy(x) for any constant c. Hence y = ce^2 x^ is a solution of (1.1.2). Are there any other solutions? No there aren’t, because if y is any function that satisfies (1.1.2) then (ye−^2 x)′^ = e−^2 x(y′^ − 2 y) = 0,
and so ye−^2 x^ must be a constant, C.
In general, we can expect that if a differential equation is of the first order, then the most general solution will involve one arbitrary constant C. This is not always the case, since we can write down differential equations that have no solutions at all. We would have, for instance, a fairly hard time (why?) finding a real function y(x) for which
(y′)^2 = −y^2 − 2. (1. 1 .3)
There are certain special kinds of differential equations that can always be solved, and it’s often important to be able to recognize them. Among there are the “first-order linear” equations
y′(x) + a(x)y(x) = 0, (1. 1 .4)
where a(x) is a given function of x.
Before we describe the solution of these equations, let’s discuss the word linear. To say that an equation is linear is to say that if we have any two solutions y 1 (x) and y 2 (x) of the equation, then c 1 y 1 (x) + c 2 y 2 (x) is also a solution of the equation, where c 1 and c 2 are any two constants (in other words, the set of solutions forms a vector space).
Equation (1.1.1) is linear, in fact, y 1 (x) = 7 and y 2 (x) = 23 are both solutions, and so is 7 c 1 + 23c 2. Less trivially, the equation
y′′(x) + y(x) = 0 ( 1. 1 .5)
is linear. The linearity of (1.1.5) can be checked right from the equation itself, without knowing what the solutions are (do it!). For an example, though, we might note that y = sin x is a solution of (1.1.5), that y = cos x is another solution of (1.1.5), and finally, by linearity, that the function y = c 1 sin^ x^ +^ c 2 cos^ x^ is a solution, whatever the constants^ c 1 and^ c 2. Now let’s consider an instance of the first order linear equation (1.1.4):
y′(x) + xy(x) = 0. (1. 1 .6)
So we’re looking for a function whose derivative is −x times the function. Evidently y = e−x
(^2) / 2 will do, and the general solution is y(x) = ce−x
(^2) / 2 . If instead of (1.1.6) we had y′(x) + x^2 y(x) = 0, (1. 1 .7)
then we would have found the general solution ce−x
(^3) / 3 . As a last example, take y′(x) + ( cosx) y(x) = 0.
The right medicine is now y(x) = esin^ x. In the next paragraph we’ll give the general rule of which the above are three examples. The reader might like to put down the book at this point and try to formulate the rule for solving (1.1.4) before going on to read about it.
Ready? What we need is to choose some antiderivative A(x) of a(x), and then the solution is y(x) = ce−A(x).
Since that was so easy, next let’s put a more interesting right hand side into (1.1.4), by considering the equation y′(x) + a(x)y(x) = b(x) (. 1 .9)
where now b(x) is also a given function of x (Is (1.1.9) a linear equation? Are you sure?).
To solve (1.1.9), once again choose some antiderivative A(x) of a(x), and then note that we can rewrite (1.1.9) in the equivalent form
e−A(x)^
d dx
eA(x)y(x)
= b(x). (1. 1 .10)
Now if we multiply through by eA(x)^ we see that
d dx
eA(x)y(x)
= b(x)eA(x)
and so, if we integrate both sides,
eA(x)y(x) =
∫ (^) x b(t)eA(t)^ dt + const. ,
where on the right side, we mean any antiderivative of the function under the integral sign. Conse- quently
y(x) = e−A(x)
(∫ (^) x b(t)eA(t)^ dt + const.
As an example, consider the equation
y′^ + y x
= x + 1. (1. 1 .12)
We find that A(x) = log x, then from (1.1.11) we get
y(x) =
x
(∫ (^) x (t + 1)t dt + C
x^2 3
x 2
C
x
We may be doing a disservice to the reader by beginning with this discussion of certain types of differential equations that can be solved analytically, because it would be erroneous to assume that most, or even many, such equations can be dealt with by these techniques. Indeed, the reason for
Again, let’s see if there is a solution of the form y = eαx. This time, substitution into (1.2.4) and cancellation of the factor eαx^ leads to the quadratic equation
α^2 + 4α + 4 = 0 , (1. 2 .5)
whose two roots are identical, both being −2. Hence e−^2 x^ is a solution, and of course so is c 1 e−^2 x, but we don’t yet have the general solution because there is, so far, only one arbitrary constant. The difficulty, of course, is caused by the fact that the roots of (1.2.5) are not distinct.
In this case, it turns out that xe−^2 x^ is another solution of the differential equation (1.2.4) (verify this), and so the general solution is (c 1 + c 2 x)e−^2 x.
Suppose that we begin with an equation of third order, and that all three roots turn out to be the same. For instance, to solve the equation
y′′′^ + 3y′′^ + 3y′^ + y = 0 ( 1. 2 .6)
we would try y = eαx, and we would then be facing the cubic equation
α^3 + 3α^2 + 3α + 1 = 0 , (1. 2 .7)
whose “three” roots are all equal to −1. Now, not only is e−x^ a solution, but so are xe−x^ and x^2 e−x.
To see why this procedure works in general, suppose we have a linear differential equation with constant coeficcients, say
y(n)^ + a 1 y(n−1)^ + a 2 y(n−2)^ + · · · + any = 0 ( 1. 2 .8)
If we try to find a solution of the usual exponential form y = eαx, then after substitution into (1.2.8) and cancellation of the common factor eαx, we would find the polynomial equation
αn^ + a 1 αn−^1 + a 2 αn−^2 + · · · + an = 0. (1. 2 .9)
The polynomial on the left side is called the characteristic polynomial of the given differential equation. Suppose now that a certain number α = α∗^ is a root of (1.2.9) of multiplicity p. To say that α∗^ is a root of multiplicity p of the equation is to say that (α − α∗)p^ is a factor of the characteristic polynomial. Now look at the left side of the given differential equation (1.2.8). We can write it in the form
(Dn^ + a 1 Dn−^1 + a 2 Dn−^2 + · · · + an)y = 0 , (1. 2 .10)
in which D is the differential operator d/dx. In the parentheses in (1.2.10) we see the polynomial ϕ(D), where ϕ is exactly the characteristic polynomial in (1.2.9).
Since ϕ(α) has the factor (α − α∗)p, it follows that ϕ(D) has the factor (D − α∗)p, so the left side of (1.2.10) can be written in the form
g(D)(D − α∗)py = 0 , (1. 2 .11)
where g is a polynomial of degree n − p. Now it’s quite easy to see that y = xk^ eα
∗x satisfies (1.2.11) (and therefore (1.2.8) also) for each k = 0, 1 ,... , p − 1. Indeed, if we substitute this function y into (1.2.11), we see that it is enough to show that
(D − α∗)p(xkeα
∗x ) = 0 k = 0, 1 ,... , p − 1. (1. 2 .12)
However, (D−α∗)(xk^ e−α
∗x ) = kxk−^1 eα
∗x , and if we apply (D−α∗) again, (D−α∗)^2 (xke−α
∗x ) =
k(k − 1)xk−^2 eα
∗x , etc. Now since k < p it is clear that (D − α∗)p(xk^ e−α
∗x ) = 0, as claimed.
To summarize, then, if we encounter a root α∗^ of the characteristic equation, of multiplicity p, then corresponding to α∗^ we can find exactly p linearly independent solutions of the differential equation, namely
eα
∗x , xeα
∗x , x^2 eα
∗x ,... , xp−^1 eα
∗x .
Another way to state it is to say that the portion of the general solution of the given differential
equation that corresponds to a root α∗^ of the characteristic polynomial equation is Q(x)eα ∗x , where Q(x) is an arbitrary polynomial whose degree is one less than the multiplicity of the root α∗.
One last mild complication may arise from roots of the characteristic equation that are not real numbers. These don’t really require any special attention, but they do present a few options. For instance, to solve y′′^ + 4y = 0, we find the characteristic equation α^2 + 4 = 0, and the complex roots ± 2 i. Hence the general solution is obtained by the usual rule as
y(x) = c 1 e^2 ix^ + c 2 e−^2 ix.
This is a perfectly acceptable form of the solution, but we could make it look a bit prettier by using deMoivre’s theorem, which says that
e^2 ix^ = cos 2x + i sin 2x e−^2 ix^ = cos 2x − i sin 2x.
Then our general solution would look like
y(x) = (c 1 + c 2 ) cos 2x + (ic 1 − ic 2 ) sin 2x.
But c 1 and c 2 are just arbitrary constants, hence so are c 1 + c 2 and ic 1 − ic 2 , so we might as well rename them c 1 and c 2 , in which case solution would take the form
y(x) = c 1 cos 2x + c 2 sin 2x.
Here’s an example that shows the various possibilities:
y(8)^ − 5 y(7)^ + 17y(6)^ − 997 y(5)^ + 110y(4)^ − 531 y(3)^ + 765y(2)^ − 567 y′^ + 162y = 0. (1. 2 .13)
The equation was cooked up to have a characteristic polynomial that can be factored as
(α − 2)(α^2 + 9)^2 (α − 1)^3.
Hence the roots of the characteristic equation are 2 (simple), 3i (multiplicity 2), − 3 i (multiplicity 2), and 1 (multiplicity 3).
Corresponding to the root 2, the general solution will contain the term c 1 e^2 x. Corresponding to the double root at 3i we have terms (c 2 + c 3 x)e^3 ix^ in the solution. From the double root at − 3 i we get a contribution (c 4 + c 5 x)e−^3 ix, and finally from the triple root at 1 we get (c 6 + c 7 x + c 8 x^2 )ex. The general solution is the sum of these eight terms. Alternatively, we might have taken the four terms that come from 3i in the form
(c 2 + c 3 x) cos 3x + (c 4 + c 5 x) sin 3x.
Exercises 1.
- Obtain the general solutions of each of the following differential equations:
(a) y′′^ + 5y′^ + 6y = 0
When we take account of the given data y 0 = 1 and y 1 = 3, we get the two equations { 1 = c 1 + c 2 3 = (−2)c 1 + (−3)c 2
from which c 1 = 6 and c 2 = −5. Finally, we use these values of c 1 and c 2 in (1.3.3) to get
yn = 6(−2)n^ − 5(−3)n^ n = 0, 1 , 2 ,.... (1. 3 .5)
Equation (1.3.5) is the desired formula that represents the unique solution of the given difference equation together with the prescribed starting values.
Let’s step back a few paces to get a better view of the solution. Notice that the formula (1.3.5) expresses the solution as a linear combination of nth powers of the roots of the associated charac- teristic equation (1.3.2). When n is very large, is the number yn a large number or a small one? Evidently the powers of −3 overwhelm those of −2, and so the sequence will behave roughly like a constant times powers of −3. This means that we should expect the members of the sequence to alternate in sign and to grow rapidly in magnitude.
So much for the equation (1.3.1). Now let’s look at the general case, in the form of a linear difference equation of order p:
yn+p + a 1 yn+p− 1 + a 2 yn+p− 2 + · · · + apyn = 0. (1. 3 .6)
We try a solution of the form yn = αn, and after substituting and canceling, we get the characteristic equation αp^ + a 1 αp−^1 + a 2 αp−^2 + · · · + ap = 0. (1. 3 .7)
This is a polynomial equation of degree p, and so it has p roots, counting multiplicities, somewhere in the complex plane.
Let α∗^ be one of these p roots. If α∗^ is simple (i.e., has muiltiplicity 1) then the part of the general solution that corresponds to α∗^ is c(α∗)n. If, however, α∗^ is a root of multiplicity k > 1 then we must multiply the solution c(α∗)n^ by an arbitrary polynomial in n, of degree k − 1, just as in the corresponding case for differential equations we used an arbitrary polynomial in x of degree k − 1.
We illustrate this, as well as the case of complex roots, by considering the following difference equation of order five:
yn+5 − 5 yn+4 + 9yn+3 − 9 yn+2 + 8yn+1 − 4 yn = 0. (1. 3 .8)
This example is rigged so that the characteristic equation can be factored as
(α^2 + 1)(α − 2)^2 (α − 1) = 0 (1. 3 .9)
from which the roots are obviously i, −i, 2 (multiplicity 2), 1.
Corresponding to the roots i, −i, the portion of the general solution is c 1 in^ + c 2 (−i)n. Since
in^ = einπ/^2 = cos
( (^) nπ 2
( (^) nπ 2
and similarly for (−i)n, we can also take this part of the general solution in the form
c 1 cos
( (^) nπ 2
( (^) nπ 2
The double root α = 2 contributes (c 3 + c 4 n)2n, and the simple root α = 1 adds c 5 to the general solution, which in its full glory is
yn = c 1 cos
( (^) nπ 2
( (^) nπ 2
- (c 3 + c 4 n)2n^ + c 5. (1. 3 .11)
The five constants would be determined by prescribing five initial values, say y 0 , y 1 , y 2 , y 3 and y 4 , as we would expect for the equation (1.3.8).
Exercises 1.
- Obtain the general solution of each of the following difference equations:
(a) yn+1 = 3yn (b) yn+1 = 3yn + 2 (c) yn+2 − 2 yn+1 + yn = 0 (d) yn+2 − 8 yn+1 + 12yn = 0 (e) yn+2 − 6 yn+1 + 9yn = 1 (f) yn+2 + yn = 0
- Find the solution of the given difference equation that takes the prescribed initial values:
(a) yn+2 = 2yn+1 + yn; y 0 = 0; y 1 = 1 (b) yn+1 = αyn + β; y 0 = 1 (c) yn+4 + yn = 0; y 0 = 1; y 1 = −1; y 2 = 1; y 3 = − 1 (d) yn+2 − 5 yn+1 + 6yn = 0; y 0 = 1; y 1 = 2
- (a) For each of the difference equations in problems 1 and 2, evaluate
lim n→∞
yn+ yn
if it exists. (b) Formulate and prove a general theorem about the existence of, and value of the limit in part (a) for a linear difference equation with constant coefficients. (c) Reverse the process: given a polynomial equation, find its root of largest absolute value by computing from a certain difference equation and evaluating the ratios of consecutive terms. (d) Write a computer program to implement the method in part (c). Use it to calculate the largest root of the equation
x^8 = x^7 + x^6 + x^5 + · · · + 1.
1.4 Computing with difference equations
This is, after all, a book about computing, so let’s begin with computing from difference equations since they will give us a chance to discuss some important questions that concern the design of computer programs. For a sample difference equation we’ll use
yn+3 = yn+2 + 5yn+1 + 3yn (1. 4 .1)
together with the starting values y 0 = y 1 = y 2 = 1. The reader might want, just for practice, to find an explicit formula for this sequence by the methods of the previous section.
Suppose we want to compute a large number of these y’s in order to verify some property that they have, for instance to check that
lim n→∞
yn+ yn
One should not think that such programming methods are only for pocket calculators. As we progress through the numerical solution of differential equations we will see situations in which each of the quantities that appears in the difference equation will itself be an array (!), and that very large numbers, perhaps thousands, of these arrays will need to be computed. Even large computers might quake at the thought of using the first method above, rather than the second, or doing the calculation. Fortunately, it will almost never be necessary to save in memory all of the computed values simultaneously. Normally, they will be computed, and then printed or plotted, and never needed except in the calculation of their immediate successors.
Exercises 1.
- The Fibonacci numbers F 0 , F 1 , F 2 ,... are defined by the recurrence formula Fn+2 = Fn+1 + Fn for n = 0, 1 , 2 ,... together with the starting values F 0 = 1, F 1 = 1.
(a) Write out the first ten Fibonacci numbers. (b) Derive an explicit formula for the nth Fibonacci number Fn. (c) Evaluate your formula for n = 0, 1 , 2 , 3 , 4. (d) Prove directly from your formula that the Fibonacci numbers are integers (This is per- fectly obvious from their definition, but is not so obvious from the formula!). (e) Evaluate lim n→∞
Fn+ Fn (f) Suppose we change the starring value F 1 = 1 to something else, say F 1 = α. Show that there is exactly one value of α for which the sequence Fn approaches zero as n → ∞, and determine that value of α. (g) Write a computer program that will compute Fibonacci numbers and print out the limit in part (e) above, correct to six decimal places. (h) Write a computer program that will compute the first 40 members of the modified Fi- bonacci sequence, in which F 1 =^ α, the number that you found in (f) above. Do these computed numbers seem to be approaching zero? Explain carefully what you see and why it happens. (i) Modify the program of part (h) to run in higher (or double) precision arithmetic. Does it change any of your answers?
- Find the most general solution of each of the following difference equations:
(a) yn+1 − 2 yn + yn− 1 = 0 (b) yn+1 = 2yn (c) yn+2 + yn = 0 (d) yn+2 + 3yn+1 + 3yn + yn− 1 = 0
