The first semester of this course is about ordinary differential equations. The ......
where dj = y(j)(tk) for j = 1,2,··· ,N. In this case the LTE is O(hN+1) and the GTE is
... one-step Runge-Kutta methods which are related to Taylor's method but do not
...
1
Method for Ordinary Differential Equations
This chapter will introduce the reader to the terminology and notation of differential equations. Students will also be reminded of some of the elementary solution methods they are assumed to have encountered in an undergraduate course on the subject. At the conclusion of this review one should have an idea of what it means to ‘solve’ a differential equation and some confidence that they could construct a solution to some simple and special types of differential equations.
1.1
Definitions and Notation
Differential equations are divided into two classes, ordinary and partial. An ordinary differential equation (ODE) involves one independent variable and derivatives with respect to that variable. A partial differential equation (PDE) involves more than one independent variable and corresponding partial derivatives. The first semester of this course is about ordinary differential equations. The phrase ‘differential equation’ will usually mean ordinary differential equation unless the context makes it clear that a partial differential equation is intended. The order of the highest derivative in the differential equation is the order of the equation. Example (1.1.1)
(b)
d3 y dy − 6x + 10y = 0 3rd order ODE 3 dx dx x (t) + sin x(t) = 0 2nd order ODE
(c)
y + λy = 0
2nd order ODE
(d)
x + x = 0
1st order ODE
(a)
(e) (f )
x2
∂u ∂u − =0 ∂x ∂y ∂2u ∂2u ∂ 2u = + ∂t2 ∂x2 ∂y 2
1st order PDE 2nd order PDE
From the equation it is usually clear which is the independent variable and which is dependent. For example, in (a) it is implied that y = y(x) while in (b) the dependence of x upon t is explicitly noted. In (c) and (d) the choice of the independent variable is arbitrary but as a matter of choice we will usually regard y = y(x) and x = x(t). The latter notation is usually reserved for the discussion of dynamical systems in which case one should think of x(t) as denoting the state of a system at time t. In contrast the notation y = y(x) might be used in the context of a boundary value problem when x has the interpretation of a spatial coordinate. Of course this is much ado 1
about nothing and the student should be prepared to see differential equations written in a variety of ways. For the partial differential equation in (e) u = u(x, y) while in (f), u = u(x, y, t). We will often denote partial derivatives by subscripts in which case (e) could be written ux − uy = 0 and similarly for (f) utt = uxx + uyy . Differential equations are divided into two other classes, linear and nonlinear. An nth order linear differential equation can be put into the form a0 y (n) + a1 y (n−1) + · · · + an−1 y + an y = b where b and the coefficients ak depend on at most the independent variable. If b = 0 the equation is said to be homogeneous, otherwise it is called nohomogeneous. Example (1.1.2)
y = y + 1
linear, nonhomogeneous
x sin t = x
linear, homogeneous
(y )2 = x + y
nonlinear
General differential equations of order 1 and 2 may be written as F (x, y, y ) = 0
(1.1.1)
F (x, y, y , y ) = 0 with the obvious notation for higher order equations. Such equations are often assumed to be solvable for the highest derivative and then written as y = f (x, y)
(1.1.2)
y = f (x, y, y ). A differentiable function y = φ(x) is a solution of (1.1.1) on an interval J if F (x, φ(x), φ (x)) = 0, x ∈ J. Usually a first order equation has a family of solutions y = φ(x, c) depending on a single parameter c. A second order equation usually has a two-parameter family of solutions y = ϕ(x, c1 , c2 ) and so on. These parameters are like constants of integration and may be determined by specifying initial conditions corresponding to some ‘time’ when a process starts. Example (1.1.3)
A two-parameter family of solutions of x = 6 is x = 3t2 + c1 t + c2
The solution that satisfies the initial conditions x(0) = 1, x (0) = −1 is x = 3t2 − t + 1. 2
(1.1.3)
Because solutions may involve one or more integrations, the graph of a solution is called an integral curve. The integral curves in the previous example are parabolas. A family of functions is a complete solution of a differential equation if every member of the family satisfies the differential equation and every solution of the differential equation is a member of the family. The family of functions given by (1.1.3) is the complete solution of y = 6. For the student who is familiar with the terminology of a general solution we should state at this point that the notion of a general solution will only be used in the context of linear equations. The appropriate definitions will be introduced in Chapter 3 at which time the distinction between a general and complete solution will be made clear. Example (1.1.4)
The family y = ϕ(x, c) = √
1 2x + c
solves y + y 3 = 0. The family is not complete since the trivial solution, y(x) = 0 (also denoted y ≡ 0), cannot be obtained for any value of c. In equation (1.1.2) the left hand member dy/dx denotes the slope of a solution curve whereas the right hand member f (x, y) gives the value of the slope at (x, y). From this point of view the solution of an ordinary differential equations is a curve that flows along directions specified by f (x, y). The use of dy/dx breaks down when the tangent is vertical but the geometric interpretation remains meaningful. To deal with this matter note that a smooth curve in the (x, y) plane may be described locally by y = ϕ(x) or x = ψ(y). If dy/dx is meaningless, dx/dy may be permissible. Consequently, it may be advantageous to treat x or y as independent or dependent variables in formulating a differential equation. These remarks motivate the following discussion. A differential equation is said to be in differential form if it is written as M (x, y)dx + N (x, y)dy = 0.
(1.1.4)
A differentiable function y = φ(x) is a solution of (1.1.4) on an interval if the substitutions y = φ(x), dy = φ (x)dx make (1.1.4) true. A differentiable function x = ψ(y) is a solution if the substitutions x = ψ(y), dx = ψ (y)dy make (1.1.4) true. An equation G(x, y) = c, c constant
(1.1.5)
is said to furnish an implicit solution of (1.1.2) if (1.1.5) can be solved for y = φ(x) or x = ψ(y) in a neighborhood of a point (x, y) satisfying (1.1.5) and φ or ψ is a solution in the sense defined above. When (1.1.5) furnishes an implicit solution the graph of this equation is also called an integral curve. 3
There is some ambiguity as to what one is expected to produce when asked to solve a differential equation. Often the differential equation is regarded as solved when one has derived an equation such as (1.1.5). However, the context of the problem may make it clear that one needs to produce an explicit formula for the solution, say y = φ(x) or x = ψ(y). The next two examples will expound upon this issue, but in general the specific problem at hand will usually dictate as to what is meant by the phrase, ‘solution’. Example (1.1.5) Consider the equation in differential form xdx + ydy = 0. This equation is easily solved by writing it as a ‘pure derivative’, i.e., 1 d(x2 + y 2 ) = 0. 2 When written in this form it is apparent that possible implicit solutions are x2 + y 2 = C.
(1.1.6)
If C < 0, the locus is empty. If C = 0 the locus consists of the point (0,0). This case does not correspond to a solution since it does not describe a differentiable curve of the form y = ϕ(x) or x = ψ(y). If C > 0, the integral curves are circles centered at the origin. Example (1.1.6)
Consider the three equations (a)
dy x dx y =− (b) =− (c) xdx + ydy = 0. dx y dy x
(a) This form of the equation implicitly requires a solution y = φ(x) and no integral curve can contain a point (x, y) where y = 0. (b) This requires a solution y = ϕ(x) and no integral curve can contain a point (x, y) where x = 0. (c) This equation allows for solutions of the form x = ψ(y), or y = ϕ(x) and leads to the family of circles (1.1.6) obtained in the previous example. In particular, in the neighborhood of any point (x, y) on the circle we can express y as a differentiable function of x or vice versa. One could obtain from (1.1.6) additional functions y = ϕ(x) by taking y > 0 on part of the circle and y < 0 on the remainder (as in (c) below). Such a function is not differentiable and hence not a solution.
(a)
y=φ(x)
(b)
x=ψ(y)
4
(c)
y=χ(x)
Fig. 1.1.1. Solutions defined by an integral curve are displayed in (a) and (b). Though the graph of χ(x) lies on an integral curve, it is not a solution of the differential equation. Having introduced the notion of a solution and integral curve, we now will review some of the most elementary methods for solving ordinary differential equations.
1.2
Examples of Explicit Solution Techniques
(a) Separable Equations. A differential equation is separable if it can be wrtitten in the form F (x, y, y ) =
dy − f (x)/g(y) = 0. dx
The differential equation is solved by ‘separating’ the variables and performing the integrations g(y)dy = f (x)dx.
Example (1.2.1)
Solve dy/dx = y 2 /x, x = 0.
We separate the variables and integrate to obtain 1 1 dy = dx y2 x or y=−
1 ln |x| + C
(1.2.1)
In the above calculation division by y 2 assumed y = 0. The function y = 0 is also a solution, referred to as a singular solution. Therefore, the family given in (1.2.1) is not complete. (b) Homogeneous Equations. A differential equation is homogeneous if it can put in the form F (x, y, y ) = y − f (y/x) = 0. In this case let z = y/x so dy dz =z+ = f (z) dx dx or
dz f (z) − z = , dx x
which is separable.
5
(c) Exact Equations. The differential equation F (x, y, y ) = 0 is exact if it can be written in differential form M (x, y)dx + N (x, y)dy = 0 where
∂N ∂M = . ∂y ∂x
Recall that if M, N are C 1 on a simply connected domain then there exists a function P such that ∂P ∂P = M, = N. ∂x ∂y Thus the exact differential equation may be written as a pure derivative dP = M dx + N dy = 0 and hence implicit solutions may be obtained from P (x, y) = C.
Example (1.2.2)
Solve (3x2 + y 3 ey )dx + (3xy 2 ey + xy 3 ey + 3y 2 )dy = 0
Here My = 3y 2 ey + y 3 ey = Nx . Thus there exists P such that
∂P ∂P = M, = N. ∂x ∂y
It is simplest to obtain P from P =
Px dx = x3 + xy 3 ey + h(y).
Differentiate with respect to y to obtain ∂P = 3xy 2 ey + xy 3 ey + h (y) = N ∂y and so h (y) = 3y 2 or h(y) = y 3 . Thus P (x, y) = x3 + y 3 ey x + y 3 and so x3 + y 3 ey x + y 3 = C provides an implicit solution. 6
(d) Integrating Factors. Given M (x, y)dx + N (x, y)dy = 0, µ = µ(x, y) is called an integrating factor if ∂ ∂ (µM ) = (µN ) ∂y ∂x The point of this definition is that if µ is an integrating factor, µM dx + µN dy = 0 is an exact differential equation. Here are three suggestions for finding integrating factors: 1) Try to determine m, n so that µ = xm y n is an integrating factor. 2) If
M y − Nx = Q(x), N µ(x) = exp Q(x)dx
then
is an integrating factor. 3) If
Nx − M y = R(y), M µ(y) = exp R(y)dy
then
is an integrating factor. (e) First Order Linear Equations. The general form of a first order linear equation is a0 (x)y + a1 (x) = b(x) We assume p =
a1 b ,f = are continuous on some interval and rewrite the equation as a0 a0 y (x) + p(x)y(x) = f (x).
By inspection we see that
µ=e
p(x)dx
is an integrating factor for the homogeneous equation. From (1.2.2) we obtain d p(x)dx y) = f (x)e p(x)dx (e dx
7
(1.2.2)
and so −
y(x) = e
Example (1.2.3)
p(x)dx
f (x)e
p(x)dx
dx + c .
Find all solutions of xy − y = x2 .
If x = 0, rewrite the equation as 1 y − y = x x and so µ = e− If x > 0,
(1/x)dx
=
(1.2.3) 1 . |x|
d y =1 dx x
and so y = x2 + Cx. If x < 0, d dx
−1 y x
= −1
and y = x2 + Cx. Hence y = x2 + Cx gives a differentiable solution valid for all x. Now consider the initial value problem xy − y = x2 y(0) = y0 . If y0 = 0, we see that there are infinitely many solutions of this problem whereas if y0 = 0, there are no solutions. The problem in this latter case arises from the fact that when the equation is put in the form of (1.2.3) the coefficient p(x) is not continuous at x = 0. The significance of this observation will be elaborated upon in the next chapter. (e) Reduction in Order. If in the differential equation the independent variable is missing, that is, the equation is of the form F (y, y , · · · , y (n) ) = 0, set
dy , dx d dv dv dy dv (y ) = = = v, y = dx dx dy dx dy d dv dv dv d dv dy v = +v y = dx dy dy dx dy dy dx v = y =
8
(1.2.4)
=v
dv dy
2 +
2 2d v v , dy 2
etc. These calculations show that with this change of variables we obtain a differential equation of one less order with y regarded as the independent variable, i.e., the differential equation (1.3.1) becomes G(y, v, v , · · · v (n−1) ) = 0. If the dependent variable is missing, i.e., F (x, y , y , · · · y (n) = 0, again set v =
dy = y to obtain dx G(x, v, v , · · · v (n−1) ) = 0.
(f) Linear Constant Coefficient Equations. A homogeneous linear differential equation with constant real coefficients of order n has the form y (n) + a1 y (n−1) + · · · + an y = 0. d and write the above equation as dx
P (D)y ≡ Dn + a1 D(n−1) + · · · + an y = 0.
We introduce the notation D =
By the fundamental theorem of algebra we can write P (D) = (D − r1 )m1 · · · (D − rk )mk (D2 − 2α1 D + α12 + β12 )p1 · · · (D2 − 2α D + α 2 + β 2 )p , where
k j=1
mj + 2
pj = n.
j=1
Lemma 1.1. The general solution of (D − r)k y = 0 is
y = c1 + c2 x + · · · + ck x(k−1) erx and the general solution of (D2 − 2αD + α2 + β 2 )k y = 0 is
y = c1 + c2 x + · · · + ck x(k−1) eαx cos(βx) + d1 + d2 x + · · · + dk x(k−1) eαx sin(βx). Proof. Note first that (D − r)erx = D(erx ) − rerx = rerx − rerx = 0 and for k > j
(D − r) xj erx = D xj erx − r xj erx = jxj−1 erx .
9
Thus we have (D − r)
k
j rx
xe
j rx = (D − r) (D − r) x e
= j(D − r)k−1 xj−1 erx = · · · = k−1
= j! (D − r)k−j (erx ) . Therefore, each function xj erx , for j = 0, 1, · · · , (k − 1), is a solution of the equation and by the fundamental theory of algebra these functions are linearly independent, i.e., 0=
k
cj xj−1 erx = erx
j=1
k
cj xj−1 , for all x
j=1
implies c1 = c2 = · · · = ck = 0. Note that each factor (D2 − 2αD + α2 + β 2 ) corresponds to a pair of complex conjugate roots r = α ± iβ. In the above calculations we did not assume that r is real so that for a pair of complex roots we must have solutions e(α±iβ)x = eiβx eαx = eαx (cos(βx) + i sin(βx)) , and any linear combination of these functions will also be a solution. In particular the real and imaginary parts must be solutions since 1 1 αx [e (cos(βx) + i sin(βx))] + [eαx (cos(βx) − i sin(βx))] = eαx cos(βx) 2 2 1 1 αx [e (cos(βx) + i sin(βx))] − [eαx (cos(βx) − i sin(βx))] = eαx sin(βx) 2i 2i Combining the above results we find that the functions
y = c1 + c2 x + · · · + cn x(n−1) eαx cos(βx) and
y = d1 + d2 x + · · · + dn x(n−1) eαx sin(βx).
are solutions and, as is shown in Section 1.6, these solutions are linearly independent.
The general solution of P (D)y = 0 is given as a linear combination of the solutions for each real root and each pair of complex roots. Let us consider an example which is already written in factored form
(D + 1)3 (D2 + 4D + 13) y = 0 The term (D + 1)3 gives a part of the solution as (c1 + c + 2x + c3 x2 )e−x 10
and the term (D2 + 4D + 13) corresponds to complex roots with α = −2 and β = 3 giving the part of the solution c4 e−2x cos(3x) + c5 e−2x sin(3x). The general solution is y = (c1 + c + 2x + c3 x2 )e−x + c4 e−2x cos(3x) + c5 e−2x sin(3x). (g) Equations of Euler (and Cauchy) Type. A differential equation is of Euler type if F (x, y , y , · · · , y (n) ) = F (y, xy , · · · xn y (n) ) = 0. Set x = et so that
dy dy dx = = y x, dt dx dt d2 y dy dy dx d − x) − y x = = (y x + y x − y x = y x2 , dt2 dt dt dx dt etc. y˙ =
In this way the differential equation can be put into the form G(y, y, ˙ · · · , y (n) ) = 0 and now the method of reduction may be applied. An important special case of this type of equation is the so-called Cauchy’s equation , xn y (n) + · · · + a1 xy + a2 y = 0, or
P (D)y ≡ xn Dn + a1 xn−1 D(n−1) + · · · + an y = 0.
With x = et and D ≡
dy we have dt
dy dt 1 dy dy = = ⇒ xDy = Dy, dx dt dx x dt d 1 dy 1 d2 y dy 2 D y= − = 2 ⇒ x2 D2 y = D(D − 1)y, dx x dt x dt2 dt .. . r r x D y = D(D − 1)(D − 2) · · · (D − r + 1)y. Dy =
Thus we can write
P (D)y = P (D)y = D(D − 1) · · · (D − n + 1) + a1 D(D − 1) · · · (D − n + 2) + · · · + an−1 D + an y = 0.
11
The second order case ax2 y + bxy + cy = 0.
(1.2.5)
will arise numerous times throughout the course. We assume the coefficients are constant and x > 0 and as an alternative to the above approach we seek a solution in the form y = xr . In this way one obtains the following quadratic for the exponent r, ar2 + (b − a)r + c = 0.
(1.2.6)
There are three cases to be considered: 1) If (1.2.6) has distinct real roots r1 , r2 then we obtain solutions xr1 , xr2 . 2) If r1 = α + iβ, r2 = α − iβ then solutions may be written as xα+iβ = xα eiβ ln x = xα (cos(β ln x) + i sin(β ln x)) and similarly xα−iβ = xα (cos(β ln x) − i sin(β ln x)). Observe that a linear combination of solutions of (1.3.2) is again a solution and so we obtain the solutions 1 α+iβ (x + xα−iβ ) = xα cos(β ln x) 2 and 1 α+iβ − xα−iβ ) = xα sin(β ln x). (x 2i 3) If (1.2.6) has repeated roots then (b − a)2 = 4ac and r1 = (a − b)/2a. We seek a second solution as y = v(x)xr1 and observe that v must satisfy Set w = v to get
axv + av = 0. xw + w = 0
and so w = c1 /x, and v = c1 ln x + c2 . Thus in the case of repeated roots we obtain the solutions xr1 ln x, xr1 . One might try to verify that c1 xr1 +c2 xr2 , xα (c1 cos(β ln x)+c2 sin(β ln x)), and c1 xr1 ln x+ c2 xr1 are complete families of solutions of Euler’s equation in each of the three cases respectively. That this is indeed the case will be a trivial consequence of a more general result for linear equations that we will prove in Chapter 3. 12
(h) Equations of Legendre Type. A slightly more general class of equations than the Cauchy equations is the Legendre equations which have the form
P (D)y ≡ (ax + b)n Dn + an−1 (ax + b)n−1 D(n−1) + · · · + a1 (ax + b)D + a0 y = 0. These equations, just as Cauchy equations, can be reduced to linear constant coefficient equations. For these equations we use the substitution (ax + b) = et which, using the chain rule just as above, gives: dy dy dt a dy = = ⇒ (ax + b)Dy = aDy, dx dt dx (ax + b) dt (ax + b)2 D2 y = a2 D(D − 1)y, .. . r r r (ax + b) D y = a D(D − 1)(D − 2) · · · (D − r + 1)y. Dy =
Thus we can write
P (D)y = P (D)y = an D(D − 1) · · · (D − n + 1) 1
+ a an−1
D(D − 1) · · · (D − n + 2) + · · · + an−1 D + an y = 0.
As an example consider the equation
(x + 2)2 D2 − (x + 2)D + 1 y = 0 Set (x + 2) = et ; then the equation can be written as D(D − 1) − D + 1 y = 0 or
D2 − 2D + 1 y = (D − 1)2 y = 0.
The general solution to this problem is y(t) = C1 et + C2 tet and we can readily change back to the x variable using t = ln(x + 2) to obtain
y = (x + 2) C1 + C2 ln(x + 2) . (i) Equations of Bernoulli Type. A Bernoulli equation is a differential equation of the form y + p(x)y = q(x)y n . 13
I can be shown that the substitution v = y 1−n changes the Bernoulli equation into the linear differential equation v (x) + (1 − n)p(x)v = (1 − n)q(x). The special cases n = 0, n = 1 should be considered separately. As an example consider the differential equation y + y = y a+1 where a is a nonzero constant. This equation is separable so we can separate variables to obtain an implicit solution or we can use Bernoulli’s procedure to derive the explicit solution y = (1 + ceax )−1/a . The student should check that both results give this solution. (j) Equations of Clairaut Type. An equation in the form
y = xy + f (y )
for any function f is called a Clairaut equation. While, at first glance, it would appear that such equation could be rather complicated, it turns out that this is not the case. As can be readily verified, the general solution is given by y = Cx + f (C) where C is a constant. As an example the equation has solution y = Cx +
√
y = xy +
4 + (y )2
4 + C 2.
Sometimes one can transform an equation into Clairaut form. For example, consider the equation y = 2xy + 6y 2 (y )2 . If we multiply the equation by y 2 we get y 3 = 3xy 2 y + 6y 4 (y )2 . Now use the transformation v = y 3 , which implies v = 3y 2 y , to write the equation as 2 v = xv + (v )2 3 whose solution is v = Cx + 23 C 2 which gives y 3 = Cx + 23 C 2 .
14
(k) Other First Order and Higher Degree Equations. A first order differential equation may have higher degree. Often an equation is given in polynomial form in the variable y and in this case we refer to the equation as having order n if n is the highest order power of y that occurs in the equation. If we write p = y for notational convienence, then such an equation can be written as pn + gn−1 (x, y)pn−1 + · · · + g1 (x, y)p + g0 (x, y) = 0.
(1.2.7)
It may be possible to solve such equations using one of the following methods (see [?]). Equation Solvable for p: It may happen that (1.2.7) can be factored into (p − F1 )(p − F2 ) · · · (p − Fn ) = 0 in which case we can solve the resulting first order and first degree equations y = F1 (x, y), y = F2 (x, y), · · · , y = Fn (x, y). This will lead to solutions f1 (x, y, c) = 0, f2 (x, y, c) = 0, · · · , fn (x, y, c) = 0 and the general solution is the product of the solutions since the factored equation can be rewritten in any form (i.e., the ordering of the terms does not matter). Thus we have f1 (x, y, c)f2 (x, y, c) · · · fn (x, y, c) = 0. For example, p4 − (x + 2y + 1)p3 + (x + 2y + 2xy)p2 − 2xyp = 0 can be factored into p(p − 1)(p − x)(p − 2y) = 0 resulting in the equations y = 0, y = 1, y = x, y = 2y. These equations yield the solutions y − c = 0, y − x − c = 0, 2y − x2 − c = 0, y − ce2x = 0 giving the solution (y − c)(y − x − c)(2y − x2 − c)(y − ce2x ) = 0.
15
Equation Solvable for y: In this case we can write the equation G(x, y, y ) = 0 in the form y = f (x, p). In this case differentiate this equation with respect to x to obtain, p=
dy dp dp = fx + fp = F (x, p, ). dx dx dx
This is an equation for p which is first order and first degree. Solving this equation to obtain φ(x, p, c) = 0 we then use the original equation y = f (x, p) to try to eliminate the p dependence and obtain Φ(x, y, c) = 0. Consider the example y = 2px + p4 x2 . We differentiate with respect to x to obtain p = 2x which can be rewritten as
dp dp + 2p + 2p4 x + 4p3 x2 dx dx
dp p + 2x dx
(1 − 2p3 x) = 0.
An analysis of singular solutions shows that we can discard the factor (1−2p3 x) and we have dp (p+2x ) = 0 which implies xp2 = c. If we write the original equation as (y−p4 x2 ) = 2px dx and square both sides we have (y − p4 x2 )2 = 4p2 x2 . With this and xp2 = c we can eliminate p to obtain (y − c2 )2 = 4cx. Equation Solvable for x: In this case an equation G(x, y, y ) = 0 can be written as x = f (y, p). We proceed by differentiating with respect to y to obtain 1 dp dx dp = = fy + fp = F (y, p, ) p dy dy dy dp . Solving this equation to obtain φ(x, y, c) = 0 we dy then use the original equation y = f (x, p) to try to eliminate the p dependence and obtain Φ(x, y, c) = 0. which is first order and first degree in
As an example consider p3 − 2xyp + 4y 2 = 0 which we can write as 2x = Differentiating with respect to y gives 2 2p dp p2 +4 = − p y dy y 2 or
dp p − 2y dy
1 y dp − 2 p p dy
(2y 2 − p3 ) = 0.
The term (2y 2 − p3 ) gives rise to singular solutions and we consider only dp =0 p − 2y dy 16
p2 4y + . y p
which has solution p2 = cy. We now use this relation and the original equation to eliminate p. First we have 4y 2x − c = p which implies (2x − c)2 = 16y
y 16y = p2 c
and finally 16y = c(2x − c)2 .
1.3
Linear Nonhomogeneous Problems
Variation of parameters Consider a nonhomogeneous second order linear equation L(y) = y + a1 y + a2 y = β and suppose {y1 , y2 } is a fundamental set. Then c1 y1 (t)+c2 y2 (t) is a general solution of L(y) = 0. A method due to Lagrange, called the method of variation of parameters, for solving L(y) = β is based on the idea of seeking a solution as yp (t) = c1 (t)y1 (t) + c2 (t)y2 (t). Then
yp = c1 y1 + c2 y2 + c1 y1 + c2 y2 .
To simplify the algebra, we impose the auxiliary condition c1 y1 + c2 y2 = 0. Then
yp = c1 y1 + c2 y2 + c1 y1 + c2 y2 .
If we substitute into L(y) = β, we want c1 (t)(y1 + a1 y1 + a2 y1 ) + c2 (t)(y2 + a1 y2 + a0 y2 ) + c1 y1 + c2 y2 = β(t). Note that the two parenthesized expressions are zero because y1 and y2 are solutions of the homogeneous equation. Thus we need to solve c1 y1 + c2 y2 = 0 c1 y1 + c2 y2 = β. By Cramer’s Rule c1 (t) =
−y2 (t)β1 (t) , W (y1 , y2 )(t)
c2 (t) = 17
y1 (t)β(t) . W (y1 , y2 )(t)
Thus a particular solution is given as t t y2 (s)β(s) y1 (s)β(s) ds + y2 (t) ds yp (t) = −y1 (t) W (s) W (s) t0 t0 t y1 (s)y2 (s) − y1 (s)y2 (s) β(s) ds = W (y1 , y2 )(s) t0 t g(x, s)β(t) ds. = t0
g(t, s) is called a Fundamental Solution. The same method works, if not as smoothly, in the general case. Consider the nth order case L(y) = y (n) + a1 y (n−1) + · · · + an−1 y (1) + an y = β and let {y1 , . . . , yn } be a fundamental set of solutions of the homogeneous problem L(y) = 0. As in the last section, given a basis of solutions {yj }nj=1 of the homogeneous problem Ly = 0 we seek a solution of L(y) = β in the form yp (t) = u1 (t)y1 (t) + · · · + un (t)yn (t). We seek a system of equations that can be solved to find u1 , . . . , un . To this end we note that by applying the product rule to yp and, collecting terms carefully, we can conclude that u1 y1 + u2 y2 + · · · + un yn = 0
⇒ yp = u1 y1 + u2 y2 + · · · + un yn
u1 y1 + u2 y2 + · · · + un yn = 0
⇒ yp = u1 y1 + u2 y2 + · · · + un yn
u1 y1 + u2 y2 + · · · + un yn = 0
⇒ yp = u1 y1 + u2 y2 + · · · + un yn
.. .
.. .
⇒
u1 y1
+ u2 y2
+ · · · + un yn(n−2) = 0 ⇒ yp(n−1) = u1 y1
u1 y1
+ u2 y2
+ · · · + un yn(n−1) = β ⇒ yp(n) = u1 y1 + u2 y2 + · · · + un yn(n) + β
(n−2)
(n−1)
(n−2)
(n−1)
(n−1)
(n)
(n−1)
+ u2 y 2 (n)
which implies L(yp ) = u1 L(y1 ) + u2 L(y2 ) + · · · + un L(yn ) + β = β. Now we note that the system of equations becomes y1 y2 ··· yn y y2 ··· yn 1 . .. .. . . . ··· . (n−1)
y1
(n−1)
y2
0 0 u2 = . . .. .. . un (n−1) β · · · yn 18
u1
+ · · · + un yn(n−1)
The determinant of the coefficient matrix is nonvanishing since it is the Wronskian W (t) of a set of linearly independent solutions to an n-order linear differential equation (see equation (1.5.4) in Section 1.5). Applying Kramer’s rule we can write the solutions as uk (t) =
Wk (t) , k = 1, · · · , n W (t)
where Wk (t) is the determinant of the matrix obtained from the coefficient matrix by replacing the T
T kth column i.e., yk yk · · · yk(n−1) by the vector 0 0 · · · β If we define n yk (t)Wk (s) g(t, s) = W (s) k=1 then a particular solution of L(y) = β is
t
yp =
g(t, s) β(s) ds. t0
Method of Undetermined Coefficients As we have already learned, the method of variation of parameters provides a method of representing a particular solution to a nonhomogeneous linear problem Ly = y (n) + a1 y (n−1) + · · · + a(n−1) y (1) + an y = f in terms of a basis of solutions {yj }nj=1 of the linear homogeneous problem. In the special case in which the operator L has constant coefficients, we have just seen that it is possible to construct such a basis of solutions for the homogeneous problem. Thus given any f we can write out a formula for a particular solution in integral form t yp (t) = g(t, s)f (s) ds. t0
Unfortunately, the method of variation of parameters often requires much more work than is needed. As an example consider the problem Ly = y + y + y + y = 1 y(0) = 0, y (0) = 1, y (0) = 0. Example 1.1. For the homogeneous problem we have (D3 + D2 + D + 1)y = (D + 1)(D2 + 1)y = 0 so we can take
y1 = cos t, y2 = sin t, y3 = e−t . 19
cos t sin t e−t − sin t cos t −e−t W (t) = −t − cos t − sin t e
Thus the wronskian is
and we can apply Abel’s theorem to obtain
W (t) = W (0)e−
t 0
1 ds
1 0 1 0 1 −1 −t e = 2e−t . = −1 0 1
Thus by the variation of parameters formula yp = u1 y1 + u2 y2 where 0 sin t e−t 1 1 t 0 cos t −e−t = − u1 = e (cos t + sin t), 2 2 −t 1 − sin t e which implies 1 u1 (t) = (cos(t) − sin(t)). 2 Similarly, we obtain 1 1 u2 = (cos t − sin t), u3 (t) = et , 2 2 which imply 1 1 u2 = (cos t + sin t), u3 (t) = et , 2 2 So we get y p = u1 y 1 + u2 y 2 1 1 1 = (cos t − sin t) cos t + (sin t + cos t) sin t + et e−t 2 2 2 =1 Wow! In retrospect it would have been easy to see that yp = 1 is a particular solution. Now the general solution is y = 1 + c1 cos t + c2 sin t + c3 e−t and we can apply the initial conditions to determine the constants which yields 1 y = 1 + (sin t − cos t − e−t ). 2 20
We note that if we were to apply the method for finding a particlar solution with the properties yp (0) = 0, yp (0) = 0, yp (0) = 0, (as given in the proof of the n order case), then we would get 1 yp (t) = 1 − (cos t + sin t + e−t ). 2 We note that the second term is part of the homogeneous solution so it can be excluded. In any case this is a lot of work to find such a simple particular solution. In case the function f is a linear combination of: 1. polynomials, 2. polynomials times exponentials or, 3. polynomials times exponentials times sine or cosine, i.e., if f is a solution of a linear constant coefficient homogeneous differential equation, one can apply the method of undetermined coefficients. The method goes as follows: n
1. Let L = P (D) = D +
n
aj Dn−j
j=1
2. Let Ly = 0 have general solution yh =
n
cj yj .
j=1
3. Assume that M = Q(D) = Dm +
m
bj Dm−j is a constant coefficient linear differential
j=1
operator such that M f = 0. = M L is a polynomial constant coefficient differential operator. 4. Then L 5. Notice that if yp is any solution of Ly = f and yh is the general solution of the homogeneous problem, then we have h + yp ) = 0. L(y
(1.3.1)
6. On the other hand we can write the general solution of this problem by simply factoring L and applying the results of the previous section. 7. Note that the solution yh =
n
cj yj is part of this general solution.
j=1
21
8. So let us denote the general solution by y=
n
cj yj +
j=1
m
dj wj .
(1.3.2)
j=1
9. Now we also know, by the variation of parameters formula, that there exists a particular solution of Lyp = f . 10. But since we know the general solution of (1.3.1) is (1.3.2), this means that any particular solution must also be a part of the full general solution of the large homogeneous problem, n m cj yj + dj wj . i.e., yp = j=1
j=1
11. We know that Lyh = 0 so we can choose yp =
m
dj wj .
j=1
12. We now substitute this expression into the original equation and solve for the constants {dj }. Example 1.2. Example: Ly = (D2 − 2D + 2)y = t2 et sin(t): The general solution of the homogeneous equation Ly = 0 is yh = c1 et cos(t) + c2 et sin(t) According to the above disscussion we seek a differential operator M so that M (t2 et sin(t)) = 0. We immediately choose M = (D2 − 2D + 2)3 and we need to compute the general solution to the homogeneous problem M Ly = (D2 − 2D + 2)3 (D2 − 2D + 2)y = (D2 − 2D + 2)4 y = 0, which implies y = (c1 + c2 t + c3 t2 + c4 t3 )et cos(t) + (d1 + d2 t + d3 t2 + d4 t3 )et sin(t). If we now remove the part of this function corresponding to the solution of the homogeneous problem Ly = 0, we have yp = (at3 + bt2 + ct)et cos(t) + (dt3 + f t2 + gt)et sin(t) After a lengthy calculation the first derivative
y = (d + a) t3 + (3 a + b + f ) t2 + (2 b + c + g) t + c et cos(t)
+ (−a + d) t3 + (−b + f + 3 d) t2 + (−c + g + 2 f ) t + g et sin(t) 22
and the second derivative
y = 2 dt3 + (6 d + 6 a + 2 f ) t2
+ (6 a + 4 b + 4 f + 2 g) t + 2 b + 2 c + 2 g et cos(t)
+ − 2 at3 + (−6 a + 6 d − 2 b) t2 + (−4 b + 4 f + 6 d − 2 c) t + 2 f + 2 g − 2 c et sin(t)
Plugging all of this into the equation yields
y − 2y + 2y = 6 dt2 + (6 a + 4 f ) t + 2 g + 2 b et cos(t) +
−6 at2 + (−4 b + 6 d) t − 2 c + 2 f et sin(t) Equating coefficients with the right hand side leads to the equations 6d = 0 6a + 4f = 0 2g + 2b = 0 −6 a = 1 −4 b + 6 d = 0 −2 c + 2 f = 0 which have the solutions a = −1/6, b = 0, c = 1/4, d = 0, f = 1/4 and g = 0. Thus, a particular solution of the equation is 3 t t2 t y= − + et cos(t) + et sin(t) 6 4 4 The following table contains a guide for generating a particular solution when one applies the method of undetermined coefficients. In particular, consider Ly = f. If f (t) =
seek yp (t).
Pm (t) = c0 tm + · · · + cm
xs (a0 tm + · · · + am )
Pm (t)eαt
ts (a0 tm + · · · + am )eαt s αt t e (a0 tm + · · · + am ) cos βt m +(b0 t + · · · + bm ) sin βt
αt
Pm (t)e
sin βt cos βt
Example 1.3. Returning to Example 1.1 we see that the operator M = D annihilates f = 1 so we seek a particular solution yp = 1. 23
1.4
Some Numerical Methods for ODEs
Method: Collocation Method (a Weighted Residual Method) Given a differential equation L[u] = 0 for u(ξ) with ξ ∈ Q (Q some domain) and boundary conditions B[u] = 0, we seek an approximate solution u(ξ) = w(ξ; α) where α = {αj }N j=1 is a set of parameters and B[w] = 0 for all choices of α. We try to determine the α by requiring that the equation be satisfied by w at a set of points S = {ξj }N j=1 ⊂ Q. Thus we arrive at a system of N equations in the N unknowns α: L[w(ξj , α)] = 0, j = 1, · · · , N.
Example 1.4. Consider L[u] = u + u + x = 0, u(0) = 0, u(1) = 0
(1.4.1)
Suppose we choose to approximate the solution by w(x) = α1 x(1 − x) + α2 x(1 − x2 ). Note that w satisfies the boundary conditions. A straightforward calculation gives: L[w(x)] = −α1 (2 − x + x2 ) − α2 (5x + x3 ) + x. Applying our collocation method at x1 = 1/3 and x2 = 2/3, we obtain the 2 × 2 system (48/27)α1 + (46/27)α2 = 1/3 (48/27)α1 + (98/27)α2 = 2/3 which has solutions α1 = 9/416 and α2 = 9/52. Thus we obtain w(x) =
9 9 x(1 − x) + x(1 − x2 ). 416 52
The exact solution to this problem is u(x) =
sin(x) −x sin(1)
and the maximum difference between y and w on [0, 1] occurs at x = .7916 with a value of .00081. Method: Taylor Series This method yields an approximate solution to a differential equation near a single point a. The method is applicable for differential equations in which all expressions are analytic functions of the variables. The idea is to use the differential equation to obtain the coefficients in a Taylor series about the point a. We note that this is really an application of the famous Cauchy-Kovalevskaya theorem.
24
Example 1.5. We seek an approximate solution to y (x) = F (x, y), y(a) = y0 in the form y(x) =
N y (j) (a) j=0
j!
(x − a)j .
The idea is to use the differential equation and initial condition to find the Taylor coefficients. y(a) = y0 y =F (x, y) implies y (a) = F (a, y0 ) y (x) =Fx (x, y) + Fy (x, y)yx implies y (a) = Fx (a, y0 ) + Fy (a, y0 )F (a, y0 ), etc.
Consider y = x2 − y 2 , y(0) = 1 A straightforward calculation yields: y = x2 − y 2 , y = 2x − 2yy y = 2 − 2(y )2 − 2yy y = −6y y − 2yy from which we immediately obtain for a = 0 y (0) = −1, y (0) = 2 y (0) = −4 y (0) = 20. Thus we obtain the approximate solution 2 2 4 20 x − x3 + x4 2! 3! 4! 2 5 = 1 − x + x2 − x3 + x4 3 6
y =1−x+
25
Euler’s Method: This method, which can be derived many different ways, is not particularly practical but is by far the simplest method and the starting point for many more complicated methods. Unlike the methods discussed in the earlier section, this method does not produce a function that can be evaluated at any point but rather is gives approximate values to the solution at points in +1 the form {(tj , yj )}M j=1 where y(tj ) ≈ yj , a = t1 < t2 < · · · < tM < tM +1 = b. While any mesh can be used, we will select uniform mesh tk = a + h(k − 1), k = 1, · · · , (M + 1), h =
(b − a) . M
Assuming that y, y and y are continuous we can apply Taylor’s theorem at t1 to obtain y(t) = y(t1 ) + y (t1 )(t − t1 ) + y (c1 )
(t − t1 )2 , c1 ∈ [t1 , t]. 2
Using y (t1 ) = f (t1 , y(t1 )) and h = (t2 − t1 ) we get an approximation for y at t2 y(t2 ) = y(t1 ) + hf (t1 , y(t1 )) + y (c1 )
h2 . 2
For h sufficiently small we write y2 = y1 + hf (t1 , y1 ) to obtain an approximate value at t2 . Repeating, we have tk+1 = tk + h, yk+1 = yk + hf (tk , yk ), k = 1, · · · (M + 1) Euler’s Method. Assuming that at each step we begin at exactly y(tk ) (which won’t be quite true), we define the local truncation error (LTE) to be the error obtained at a given step. In this case the LTE is y (ck )
h2 . 2
Summing the LTE’s all the way to tM we get, roughly anyway, the so-called Global Truncation Error (GTE). In this case we have M
h2 y (2) (c) (b − a) (2) 2 y (ck ) ≈ y (c)M h = M h = O(h). 2 2 M k=1 (2)
Thus we see that the GTE, E(h), for Euler’s method is E(h) ≈ Ch. From this we have
h 1 1 E ≈ C(h/2) = Ch ≈ E(h). 2 2 2 26
So cutting the step size in half roughly reduces the GTE by a factor of 1/2. Note that if the equation is y = f (t) with y(a) = 0 then this method give y(b) ≈
M
f (tk )h
k=1
b
f (t) dt.
which is a Riemann sum approximation to the integral a
Modified Euler’s Method: Once again we consider (??) and apply the fundamental theorem of calculus to get t2 t2 f (t, y(t)) dt = y (t) dt = y(t2 ) − y(t1 ). t1
t1
This implies that
t2
y(t2 ) = y(t1 ) +
f (t, y(t)) dt. t1
Now use the trapezoid rule to evaluate the integral with step size h = t2 − t1 to obtain h y(t2 ) ≈ y(t1 ) + [f (t1 , y(t1 )) + f (t2 , y(t2 ))]. 2 Unfortunately, we do not know y(t2 ) on the right side, so we use, for example, Euler’s method to approximate it: y(t2 ) ≈ y(t1 ) + hf (t1 , y(t1 )). We thus obtain a numerical procedure by denoting y(tk ) by yk for each k we set pk = yk + hf (tk , yk ), tk+1 = tk + h, h yk + 1 = yk + [f (tk , yk ) + f (tk+1 , pk )] . 2 Note that if the equation were y = f (t) with y(0) = 0 then we would have h y(b) ≈ [f (tk ) + f (tk+1 )] 2 k=1 M
which is exactly the trapezoid rule for numerical quadrature. In this case the LTE for the trap rule is −y (ck ) and the GTE is −
M k=1
y (2) (ck )
h3 12
h3 y (2) (c)(b − a) 2 ≈− h = O(h2 ). 12 12
27
Thus we see that the GTE E(h) for Modified Euler Method is E(h) ≈ Ch2 . 1 1 h ≈ C(h/2)2 = Ch2 ≈ E(h). E 2 4 2
From this we have
So cutting the step size in half roughly reduces the GTE by a factor of 1/4. This is sometimes called Huen’s method. See the codes a9_2.m and huen.m. Taylor’s Method: Assuming that y has derivatives of all order we could, again by Taylor’s theorem, for any N write y(t + h) = y(t) + y (t)h + with a LTE of
y (2) (t)h2 y (N ) (t)hN + ··· + 2 N!
(1.4.2)
y (N +1) (c)hN +1 . (N + 1)!
We can adapt this result to each interval [tk , tk+1 ] to obtain a numerical procedure dN hN d2 h2 d3 h3 + + ··· + 2 6 N!
yk+1 = yk + d1 h +
where dj = y (j) (tk ) for j = 1, 2, · · · , N . In this case the LTE is O(hN +1 ) and the GTE is EN (h) ≈ ChN = O(hN ). Thus we have h 1 EN = N EN (h). 2 2 Computing the values dk for a specific example is not to bad. Consider, for example, (see the exercises) t−y 2 2−t+y = 4 −2 + t − y = 8 2−t+y = 16
y = y (2) y (3) y (4)
from which we readily compute the values dj . There is a Taylor’s matlab code named taylor.m taken from John Mathews book Chapter 9. The file a9_3.m runs an example. Runge-Kutta Methods: 28
The dissadvantage of Taylor’s method is obvious. We have to do some analytic calculations on every problem. To eliminate this difficulty and still have a higher order method we consider the one-step Runge-Kutta methods which are related to Taylor’s method but do not require separate calculations. We give a brief discussion of this method in the second order case. The fourth order case is more common but the algebra is very messy. In the second order case the idea is to look for a formula yk+1 = yk + hF (tk , yk , h, f ) with F in the form F (t, y, h, f ) = γ1 f (t, y) + γ2 f (t + αh, y + βhf (t, y). where γ1 , γ2 , α, β are to be determined. We apply Taylor’s theorem in two variable (t, y) applied to the second term (the term containing γ2 ) to obtain an approximation through the second derivative terms:
2 1 2 1 2 2 F (t, y, h, f ) = γ1 f (t, y)+γ2 f (t, y)+h αft +βf fy +h α ftt +αβfty f + β f fyy +O(h3 ) 2 2 where ft =
∂f ∂f , fy = , etc. ∂t ∂y
Similarly, we have y = f (t, y) y = ft + fy y = ft + fy f y = ftt + 2fty f + fyy f 2 + ft fy + fy2 f. Now if we expand y(t) in a Taylor polynomial expansion about tk we have y(tn + h) = y(tn ) +
y (tn ) y (tn ) 2 y (tn ) 3 h+ h + h + O(h4 ). 1 2 6
(k)
Let us denote y (k) (tn ) = yn , then the LTE is LT E =y(tn+1 ) − y(tn ) − hF (tn , y(tn ), h; f ) h2 h3 =hyn(1) + yn(2) + yn(3) + O(h4 ) − hF (tn , y(tn ), h; f ) 2 6 =h[1 − γ1 − γ2 ]f + h2 [(1/2 − γ2 α)ft + (1/2 − γ2 β)fy f ] + h3 [(1/6 − 1/2γ2 α2 )ft t + (1/3 − γ2 αβ)fty f + (1/6 − 1/2γ2 β 2 )fyy f 2 + 1/6fy ft + fy2 f ] + O(h4 ) For general f we cannot eliminate the third term but by choosing γ1 + γ2 = 1 γ2 α = 1/2 γ2 β = 1/2 29
we can eliminate the first two terms. This implies that the LTE is O(h3 ). The system of equations is underdetermined but we can write γ1 = 1 − γ2 , α = β =
1 2γ2
for arbitrary γ2 . Special Cases: For γ2 = 1/2 we get the modified Euler’s Method h yn+1 = yn + [f (tn , yn ) + f (tn+1 , yn + hf (tn , yn ))]. 2 For γ2 = 1 we get the so-called Midpoint Method h yn+1 = yn + [f (tn + h/2, yn ) + f (tn+1 , yn + (h/2)f (tn , yn ))]. 2 We now consider choosing γ2 in order to minimize the LTE at the O(h3 ) level. Note that LT E = c(f, γ2 )h3 + O(h4 ) where
1
1
1 1 1 1 1 1 2 2 c(f, γ2 ) = − − − ftt + fty f + fyy f + fy ft + fy f . 6 8γ2 3 4γ2 6 8γ2 6 6
Now applying the Cauchy-Schwartz inequality we get |c(f, γ2 )| ≤ c1 (f )c2 (γ2 ),
2 4 f + fy2 ft2 + fy4 f 2 c1 (f ) = ftt2 + fty2 f 2 + fyy
1/2 ,
1/2
1 1 2 1 1 2 1 1 2 1 − − − + + + . c2 (γ2 ) = 6 8γ2 3 4γ2 6 8γ2 18 We can compute the minimum of c2 (·) to obtain √ γ2 = 3/4, with c2 (3/4) = 1/ 18. The resulting method is usually refered to as Huen’s Method which takes the form yk+1 = yk +
h 2 2 f (tk , yk ) + 3f (tk + h, yn + hf (tk , yk )) . 4 3 3
30
There is a Huen’s matlab code named huen.m in Chapter 9 taken from John Mathews book. This same analysis can be carried out for “higher order” methods but the algebra becomes very messy. One method that was very popular until about 1970 was the 4th order Runge-Kutta (RK4) method with a LTE of O(h5 ). If f does not depend on y this method reduces to simpson’ rule for quadrature. The RK4 method goes like this f1 f2 f3 f4
= f (tk , yk ) = f (tk + h/2, yk + h/2f1 ) = f (tk + h/2, yk + h/2f2 ) = f (tk + h, yk + hf3 ) h yk+1 = yk + (f1 + 2f2 + 2f3 + f4 ). 6
Assignment 3: 1. Program Example 1.4 using Maple, plot y and w and check the error given in the example. 2. Try using two different functions for example to approximate the function w in Example 1.4. 3. Try using different points x1 and x2 , e.g. x1 = 1/4 and x2 = 3/4 in Example 1.4. 4. Use Maple to carry out the Taylor method up to order n for y = F (x, y) with y(a) = y0 . Apply your program with n = 4, 6, 8 and plot the results. Use a) F (x, y) = y 2 , y(1) = 1. How do your answers compare with the exact solution on [1, 2]? b) Same question with F (x, y) = 3x2 y, y(0) = 1 on [0, 1]. 5. Apply Euler’s method to the following: a) y = t2 − y on [0, 2] with y(0) = 1. The exact answer is y = −e−t + t2 − 2t + 2. Use h = .2, .1, .05. 1 b) y = on [0, 1] with y(0) = 1. Note that the exact answer doesn’t exist on this 1 − t2 interval. Try the methods anyway with h = .1, .05, .025. c) y = e−2t − 2y on [0, 2] with y(0) = 1/10 and h = .2, .1, .05. Here y = e−2t /10 + te . −2t
6. Apply the midpoint method, the modified euler method and Huen’s method for the following problem with the same h in each case. y = −y + t + 1, 0 ≤ t ≤ 1, y(0) = 1. What do you observe, i.e., how do they compare? 7. Solve the equation y = 1 + y 2 using any method above and Taylor’s method of order 2, 4 and 6 on [0, 5] with y(0) = 1 and compare your answers to y = tan(t − π/4). 31
1.5
Wronskian and Linear Independence
Consider an nth order linear nonhomogeneous equation y (n) + a1 (t)y (n−1) + a2 (t)y (n−2) + · · · + an (t)y = β(t)
(1.5.1)
assuming ai (t), bi (t) ∈ C[a, b]. Let
dn (·) dn−1 + a (·) + · · · + an (·). 1 dtn dtn−1 Then (1.5.1) may be written as the linear nonhomogeneous system L(·) =
L(y) = β.
(1.5.2)
The operator L is said to be linear since L(cy1 (t) + c2 y2 (t)) = c1 L(y1 (t)) + c2 L(y2 (t)). It easily follows that the set of solutions of the homogeneous linear system L(y) = 0
(1.5.3)
is a vector space. The Wronskian of {yj (t)}nj=1 is defined as y1 (t) ··· yn (t)
.. .. .. W (t) ≡ W {yj (t)}nj=1 (t) = . . . (n−1) (n−1) y1 (t) · · · yn (t)
.
(1.5.4)
Theorem 1.1. [Abel’s Formula] If y1 , . . . yn are solutions of (1.5.3) , and t0 ∈ (a, b), then t a1 (s)ds . W (t) = W (t0 ) exp − t0
Thus the Wronskian of {y1 , · · · yn } is never 0 or identically 0. Proof. We compute y1 ··· yn ··· yn y1 W (t) = .. .. .. . . . (n−1) (n−1) y · · · y1 1
y1 y 1 + y1 .. . y (n−1) 1
···
yn
···
yn
··· yn .. .. . . (n−1) · · · yn
32
y1 ··· yn y11 .. .. .. + ··· + . . . (n−2) (n−2) y · · · yn 1 (n) (n) y · · · yn 1
y1 y1 .. . = (n−2) y1 n (n−j) αj y1 − j=1
···
yn
··· .. . ···
yn .. .
···
··· y1 y1 ··· = .. .. . . −a (t)y (n−1) · · · 1 1
(n−2)
yn n − αj yn(n−j) j=1
yn yn .. . (n−1) −a1 (t)yn
= −a1 (t)W (t).
t a1 (ts)ds W (t) = K exp −
Hence
t0
or, more explicitly,
t a1 (s)ds . W (t) = W (t0 ) exp − t0
Definition 1.1. A collection of functions {yi (t)}ki=1 is linearly independent on (a, b) if k
ci yi (t) = 0 for all x ∈ (a, b) ⇒ cj = 0 for j = 1, · · · , n.
i=1
Otherwise we say the set {yi (t)} is linearly dependent. Theorem 1.2. Suppose y1 , . . . , yn are solutions of (1.5.2). If the functions are linearly dependent on (a, b) then W (t) = 0 for all t ∈ (a, b). Conversely, if there is an t0 ∈ (a, b) so that W (t0 ) = 0, then W (t) = 0 for all t ∈ (a, b) and the yi (t) are linearly dependent on (a, b). Proof. If the {yi (t)} are linearly dependent, then there are ci (t) not all zero such that ci yi (t) = 0 for all t ∈ (a, b) i
⇒ Hence, defining
(k)
ci yi (t) = 0 for all t and any k.
y1 (t) ··· yn (t) c1 .. . . . M (t) = .. .. .. , C = . , (n−1) (n−1) cn y1 (t) · · · yn (t) 33
the system can be written as M (t)C = 0 and since C = 0 we see that M is singular and therefore W (t) = det(M (t))0 for all x ∈ (a, b). Conversely, if det(M (t)) = W (t0 ) = 0 then M (t)C = 0 has a nontrivial solution. For this choice of ci ’s, let y(t) = ci yi (t). Then
y(t0 ) = 0, y (t0 ) = 0, · · · , y (n−1) (t0 ) = 0.
and since y is a solution of Ly = 0, from the uniqueness part of the fundamental existence uniqeness theorem, we must have y(t) = 0 for all x ∈ (a, b). Example 1.6. Consider y1 = t2 , y2 = t|t|. Then 2 t t|t| = 2t3 sgn(t) − 2t2 |t| ≡ 0. W (t) = 2x 2x sgn(t) However, y1 (t), y2 (t) are not linearly dependent. For suppose c1 y1 (t) + c2 y2 (t) = 0, for all t. Then c1 t + c2 |t| = 0 for all t. If t > 0, c1 t + c2 t = 0 ⇒ c1 = −c2 while for t < 0, c1 t − c2 t = 0 ⇒ c 1 = c2 . Hence c1 = c2 = 0 and so y1 , y2 are linearly independent on any interval (a, b) containing 0. Thus y1 , y2 are not solutions of a linear homogeneous 2nd order equation on (a, b).
1.6
Completeness of Solution for Constant Coefficient ODE
We have already learned, in Section 1.2, how to find a set of n solutions to any homogeneous equation of the form Ly = 0 with L = D(n) + a1 D(n−1) + · · · + an−1 D + an . Namely, we factor the operator into a product of factors (D − r)k and (D2 − 2αD + α2 + β 2 )m . Having done this we
34
simply observe that the general solution of the associated homogeneous problem for each of these types of operators is easy to write out. Namely, we have (D − r) y = 0
⇒ y=
k
k
cj tj−1 ert
(1.6.1)
j=1
(D − 2αD + α + β ) y = 0 2
2
⇒ y=
2 m
k
cj tj−1 eαt cos(βt)
j=1
+
k
cj tj−1 eαt sin(βt)
(1.6.2)
j=1
In the case that the coefficients ai are constant, it is possible to describe the solutions explicitly by simply solving the homogeneous equation for each factor and adding these terms together. What we have not proved is that all such solutions give a basis for the null space of L, i.e., we have not shown that the soutions are linearly independent. To show that these solutions are linearly independent is not really difficult but to do it completely rigorously and carefully is a bit lengthy. First we note that Lemma 1.2. If λ = α + βi is a real (β = 0 ) or complex number, then k j−1 y= cj t eλt j=1
is the complete solution of (D − λ)k y = 0. Proof. Showing that the solutions of (D − λ)k y = 0 are linearly independent amounts to showing that k cj tj−1 eλt = 0 for all t ∈ R ⇒ cj = 0, j = 1, 2, · · · , k. j=1
But, on noting that eλt = 0 and dividing, this result is obvious from the fundamental theorem of algebra which says that a polynomial of degree k has exactly k zeros. Lemma 1.3. If λ1 = λ2 are two complex numbers and p(t) =
k
j−1
cj t
and q(t) =
j=1
dj tj−1 ,
j=1
are two polynomials, then p(t)eλ1 t = q(t)eλ2 t for all t ∈ R ⇒ p(t) = 0,
35
q(t) = 0.
Proof. To see that this is true we first multiply both sides of the equation by e−λ1 t so that p(t) = q(t)e(λ2 −λ1 )t for all t ∈ R. Now consider the cases in which α < 0, α > 0 and α = 0 where (λ2 − λ1 ) ≡ α + βi. If α < 0 then (using L’Hospital’s rule in the first term) lim q(t)e(λ2 −λ1 )t = 0 while
t→+∞
lim p(t) = ±∞ (as ck is pos. or neg.).
t→+∞
So that we must have p(t) ≡ 0 and then q(t) ≡ 0. If α > 0 we repeat the same argument with the first limit replace by t → −∞. Finally, in the case α = 0 we divide both sides of the equation by q(t) and collect real and imaginary parts to obtain r1 (t) + ir2 (t) =
p(t) = eβit = cos(βt) + i sin(βt) q(t)
where r1 (t) and r1 (t) are rational functions with real coefficients. Equating real and imaginary parts we see that this would imply that r1 (t) = cos(βt), r2 (t) = sin(βt) which is impossible unless r1 (t) = 0 and r2 (t) = 0 since the right side has infinitely many zeros while the left can have only a finite number. This in turn implies that p(t) = 0 and also q(t) = 0. Lemma 1.4. If < > 0, λ1 = λ2 are real or complex numbers and
(D − λ2 ) p(t)eλ1 t = 0 where p(t) is a polynomial, then p(t) ≡ 0.
Proof. We know that every solution of (D − λ2 ) y = 0 can be written as y = q(t)eλ2 t for some polynomial q(t) of degree at most (< − 1). So the equestion is whether or not there exists a polynomial q(t) so that p(t)eλ1 t = q(t)eλ2 t . We note that this is only possible when p(t) = 0 and q(t) = 0 by Lemma 1.3. Lemma 1.5. If p(t) is any polynomial of degree less than or equal (n − 1) then
(D − λ1 )m p(t)eλ2 t = q(t)eλ2 t where q(t) is a polynomial of degree at most the degree of p(t). Proof. Consider the case m = 1. We have
(D − λ1 ) p(t)eλ2 t = q(t)eλ2 t where q(t) = p (t) + (λ2 − λ1 )p(t) which is a polynomial of degree p(t). You can now iterate this result for general < > 0. 36
Lemma 1.6. If L(y) = y (n) +a1 y (n−1) +· · ·+an y has real coefficients and p(r) = rn +a1 r(n−1) + · · · + an . Then p(z) = p(z) for all z ∈ C. Therefore if p(α + βi) = 0 then p(α − βi) = 0. Proof. For every z1 , z2 ∈ C, we have z1 + z2 = z1 + z2 and z1 z2 = z1 z2 which also implies z1n = z1 n . From Lemma 1.6 we know that for a differential operator L with real coefficients, all complex roots must occur in complex conjugate pairs (counting multiplicity) and from Lemma 1.2 we know that for a pair of complex roots λ = α + βi each of multiplicity k, a set of 2k linearly independent solutions is given for j = 0, · · · , (k − 1) by j λt j αt t e =t e cos(βt) + i sin(βt) . j λt
j αt
t e =t e
cos(βt) − i sin(βt) .
From this we see that there is a set of real solutions given as a linear combination of these solutions by 1 j λt j αt λt t e cos(βt) = t e + e , 2 and 1 j λt j αt λt x e sin(βt) = t e − e . 2i We already know from Lemma 1.2 that tj eλt and tj eλt are linearly independent. Suppose we have a linear combination cj tj eαt cos(βt) + dj tj eαt sin(βt) = 0. This would imply that
(cj − dj i) j λt (cj + dj i) j λt te + t e = 0, 2 2 but since these functions are independent this implies (cj − dj i) = 0, (cj + dj i) = 0, which implies cj = dj = 0. Combining these results we have the main theorem:
Theorem 1.3. If L = y (n) + a1 y (n−1) + · · · + an y has real coefficients and we assume that the polynomial p(r) = rn + a1 r(n−1) + · · · + an has zeros given by r1 , r1 , r2 , r2 , · · · , r , r , r2 +1 , · · · , rs where rj = αj + βj i, j = 1, · · · ,
