Ordinary Differential Equations Lecture Notes

Ordinary Differential Equations Objectives These notes introduce the analytical solution of ordinary differential equations. Emphasis is placed on simple equations of first and second order, with emphasis on equations with constant coefficients. Brief treatment is given to nonhomogeneous equations of second and higher orders. There is a brief discussion on using MATLAB to obtain symbolic solutions of differential equations. Background and basic definitions A differential equation is an equation, which contains a derivative. The simplest kind of a differential equation is shown below: dy  f ( x) dx with y  y0 at x  x0 [1] In general, differential equations have an infinite number of solutions. In order to obtain a unique solution one or more initial conditions (or boundary conditions) must be specified. In the above example, the statement that y = y0 at x = x0 is an initial condition. (The difference between initial and boundary conditions, which is really one of naming, is discussed below.) The differential equation in [1] can be “solved” as a definite integral. x  f ( x)dx y - y0  [2] x0 The definite integral can be either found from a table of integrals or solved numerically, depending on f(x). The initial (or boundary) condition (y = y0 at x = x0) enters the solution directly. Changes in the values of y0 or x0 affect the ultimate solution for y. A simple change – making the right hand side a function of x and y, f(x,y), instead of a function of x alone – gives a much more complicated problem. dy  f ( x, y) dx with y  y0 at x  x0 [3] We can formally write the solution to this equation just as we wrote equation [2] for the solution to equation [1]. y - y0  x  f ( x, y)dx x0 [4] Here the definite integral can no longer be evaluated simply. Thus, alternative approaches are needed. Equation [4] is used in the derivation of some numerical algorithms. The (unknown) exact value of f(x,y) is replaced by an approximate interpolation polynomial which is a function of x only. In the theory of differential equations, several approaches are used to provide analytical solutions to the differential equations. Regardless of the approach used, one can always check to see a proposed solution is correct by substituting a proposed solution into the original differential equation and determining if the solution satisfies the differential equation and the initial or boundary conditions. For example, without knowing how to solve the following differential equation and initial conditions, d2y dy 2   101y  10.4e x 2 dx dx with y  1.1 and dy  0.9 at x  0 dx [5] You can verify that the equation below satisfies the differential equation and boundary conditions. y  e x cos10 x  0.1e x [6] To show this we first set x = 0 to find that y(0) = 1 + 0.1 = 1.1 as required by the first initial condition for y. Taking the first derivative of the proposed solution gives. dy  e x cos10 x  10e x sin 10 x  0.1e x dx [7] Evaluating the first derivative at x = 0 gives dy/dx|x=0 = –1 + 0 + 0.1 = –0.9 as required by the second initial condition. We need the second derivative to show that the solution satisfies the differential equation. This is found by taking the first derivative of equation [7]. d2y  e x cos10 x  10e x sin 10 x  10e x sin 10 x  100e x cos10 x  0.1e x 2 dx  e x cos10 x  20e x sin 10 x  100e x cos10 x  0.1e x [8] Substituting equations [6], [7], and [8] into the original differential equation in [5], shows, after some cancellation, that the original differential equation is satisfied by the solution in equation [6]. d2y dy  2  101y  e  x cos10 x  20e  x sin 10 x  100e  x cos10 x  0.1e x 2 dx dx  2  e  x cos10 x  10e  x sin 10 x  0.1e x  101 e  x cos10 x  0.1e x      1  100  2  101e  x cos10 x  20  20e  x sin 10 x  0.1  0.2  10.1e x  10.4e x In the remainder of these notes we will be interested in showing how analytical solutions to differential equations can be obtained in some simple cases. However, you should recognize that the proof of a solution is the demonstration that the solution satisfies the differential equation and the boundary or initial conditions. [9] Ordinary differential equations involve functions, which have only one independent variable. Thus, they contain only ordinary derivatives. Partial differential equations involve functions with more than one independent variable. Thus, they contain partial derivatives. The abbreviations ODE and PDE are used for ordinary and partial differential equations, respectively. In an ordinary differential equation, we use the general notation that we are trying to solve for a function, y(x), where the equation involves derivatives of y with respect to x. We call y the dependent variable and x the independent variable. Of course, engineering problems will use different symbols, appropriate to the physical problem, in place of x and y. The most common example is the use of the symbol t, suggesting time, as the independent variable. The order of the differential equation is the order of the highest derivative in the equation. Equations [1] and [3] are first-order differential equations. A differential equation with first, second and third order derivatives only would be a third order differential equation. In a linear differential equation, the terms involving the dependent variable and its derivatives are all linear terms. The independent variable may have nonlinear terms. Thus x3d2y/dx2 + y = 0 is a linear, second-order differential equation. ydy/dx + sin(y) = 0 is a nonlinear first-order differential equation. (Either term in this equation – ydy/dx or sin(y) – would make the differential equation nonlinear.) Differential equations need to be accompanied by initial or boundary conditions. An nth order differential equation must have n initial (or boundary) conditions in order to have a unique solution. Although initial and boundary conditions both describe equations that give specific values to the dependent variable (or its derivatives) at one or more points, the term “initial conditions” is usually used when all the conditions are specified at one initial point. The term “boundary conditions” is used when the conditions are specified at two different values of the independent variable. For example, in a second order differential equation for y(x), the specification that y(0) = a and y’(0) = b, would be called two initial conditions. The specification that y(0) = c and y(L) = d, would be called two boundary conditions. The initial or boundary conditions can involve a value of the variable itself, lower-order derivatives of the variable, or equations containing both values of the dependent variable and its lower-order derivatives. An equation that only has terms with the dependent variable and its derivatives may be arranged so that all such terms are on the left-hand side and the right-hand side is zero. Such equations are called homogenous differential equations. Differential equations that contain one or more terms in the independent variable only are called non-homogeneous equations. Some simple ordinary differential equations From previous courses, you should be familiar with the following differential equations and their solutions. If you are not sure about the solutions, just substitute them into the original differential equation. dy k dt with y  y0 at t  t0  d2y dy  k with y  y and  v0 at t  0  0 dt dt 2 y  y0  k (t  t0 ) kt 2 y  v0t  y0 2 [10] [11] dy  ky dt with d2y  k 2 y 2 dx d2y  k2y 2 dx  y  y0 at t  t0   y  y0e k (t t0 ) y  A sin( kx)  B cos(kx) y  A sinh( kx)  B cosh(kx)  A' e kx  B' e  kx [12] [13] [14] In equations [13] and [14] the constants A and B (or A’ and B’) are determined by the initial or boundary conditions. We call solutions such as these general solutions. A general solution of the differential equation has constants that can be modified to represent any initial or boundary conditions. A particular solution is sometimes defined as a solution that satisfies the differential equation and boundary conditions. However, the term particular solution has a slightly different meaning in some cases that we will consider later – the solution of non-homogenous equations. Note that we have used t as the independent variable in equations [10] to [12] and x as the independent variable in equations [13] and [14]. There are four possible functions that can be a solution to equation [13]: sin(kx), cos(kx), eikx, and e-ikx, where i2 = -1. Similarly, there are four possible functions that can be a solution to equation [14]: sinh(kx), cosh(kx), ekx, and e-kx. In each of these cases the four possible solutions are not linearly independent.1 The minimum number of functions with which all solutions to the differential equation can be expressed is called a basis set for the solutions. The solutions shown above for equations [13] and [14] each contain a pair of functions that are basis sets for the solutions to those equations. One final solution that is useful is the solution to general linear first-order differential equation. This equation can be written as follows. dy  f ( x) y  g ( x) dx [15] This equation has the following solution, where the constant, C, is determined from the initial condition.  f ( x ) dx   g ( x)e  f ( x ) dx dx ye  C         [16] First-order equations First order differential equations are often used to model rate processes. For example, radioactive decay where the content of radioactive nuclei is denoted by the symbol, n, is modeled by the following first-order differential equation. 1 We have the following equations among these various functions: e x  ex sinh( x )  2 e x  ex cosh( x )  2 e ix  e ix sin( x )  2 e ix  e ix cos( x )  2 dn  kn dt [17] For a positive constant, k, this equation tells us that the rate dn/dt is negative and proportional to the amount of radioactive nuclei, n, present. If the initial content at t = 0 is n0 we can multiply this equation by dt/n to obtain the following form that can be integrated directly. dn  kdt  n t n dn  kt  n  k  kdt  ln  n0   kt  n  n0e n0 0 n [18] The half-life for radioactive decay, t1/2, is defined as the time required for the initial radioactivity, n0, to decrease to half its original value. Equation [18] shows us how to compute this half-life.  n0  ln  2  n0    ln( 2) 1   ln  2    ln( 2)  kt1 / 2  t1 / 2  k     [19] With equation [19], we can rewrite the final version of equation [18] to introduce the half life. n  n0e  t t1/ 2 ln(2 ) [20] Similar equations apply to first-order chemical reactions and Newton’s law of cooling. Separable equations – Equation [17] is an example of a separable differential equation. This is a differential equation that can be separated into two sides, each of which is a function of one variable only. Some examples of separable equations and their general solutions are shown below. In each case, C represents a constant that can be determined by specifying an initial condition. (This comes from the usual constant in the general result for an indefinite integral.) dy  f ( x)  y   f ( x)dx  C dx dy dy  f ( x) g ( y )    f ( x)dx  C dx g ( y)  dy  y  h   dx x  [21] dx du  C x h(u )  u Exact forms – describe a set of relations that provide an analytical solution for a first order differential equation. Their use lies mainly in subsequent analytical derivations. Here we consider a first-order differential equation of the following form, where P(x,y) and Q(x,y) are arbitrary functions of x and y. dy P ( x, y )  dx Q ( x, y ) [22] We see that we can rewrite this equation as P( x, y)dx  Q( x, y)dy  0 [23] The title “exact form” comes from the following relationship for a differential df, which depends on dx and dy. df  P( x, y)dx  Q( x, y)dy [24] In this equation we may or may not be able to say that f is a function of x and y, written as f(x,y). If f is a function of x and y we can write the following equation for the total differential df. df  f f dx  dy x y [25] If f = f(x,y), the integral of df between any two points (x1,y1) and (x2,y2) will be simply f(x2,y2) – f(x1,y1) and will not depend on the path chosen for the integral. Such a differential is called an exact differential. This is the origin of the term exact in the discussion of exact forms for first order differential equations. By comparing equations [24] and [25] we see that, for exact differentials, Exact df  P( x, y)dx  Q( x, y)dy  P( x, y)  f f and Q( x, y)  x y [26] Since the order of differentiation can be reversed for mixed second-order partial derivatives, we can write the following formula for P(x,y) and Q(x,y) in exact differentials only. For exact differentials :  2 f P( x, y )  2 f Q( x, y )    y xy x yx [27] Consider the following two examples dy x2  y2 (1)  dx 2 xy dy x 2  y 2 (2)  dx 2 xy [28] If we compare equation [22] with the first example, we see that P(x,y) = x2 + y2 and Q(x,y) = 2xy for this example. In this case P( x, y ) Q( x, y)  2y   2y y x [29] So, the first example is an exact form. In the second example we can define P(x,y) = x2 + y2 and Q(x,y) = -2xy. In this case P( x, y ) Q( x, y)  2y   2 y y x [30] This is not an exact differential. We would have reached the same conclusion if we had used the minus sign in the definition of P rather than the definition of Q. Having defined an exact form and seen how to determine if we have such a form, we can now see how to use this fact to integrate equation [22], dy P ( x, y ) P Q . If in this   y x dx Q ( x, y ) equation, we know that there is some function f such that df = Pdx + Qdy, where P = ∂f/∂x and Q = ∂f/∂y. Furthermore, our differential equation, written in the form of equation [23] tells us that this unknown function f, has the following equation: df = Pdx + Qdy = 0. That is, df is zero so f is a constant. Although we are not interested in this mysterious function, f, we can use the results in this paragraph to integrate equation [22] and get a relationship between x and y. To integrate the exact form, we start by considering an integration of the differential equation holding y constant (so that dy = 0). f   df   P( x, y)dx  g ( y)  C1 [31] y  const In this integration we have an unknown function of y, g(y), in place of the usual constant of integration, because we are ignoring the y dependence of P(x,y) in the integration. In addition, we have the result that f = C1, a constant, since df = 0. We next take the partial derivative of equation [31] with respect to y to obtain the following expression for Q. Q  g ( y ) f      P( x, y )dx   y y  y  const y  [32] Solving this equation for dg(y)/dy gives.  dg ( y )    Q( x, y )    P( x, y )dx   h( y ) dy y  y  const  [33] Although P(x,y) and Q(x,y) involve x, our assumption of an exact form and the resulting derivation tell us that the result for dg(y)/dy must be a function of y only. We expect that when we evaluate this expression the terms containing x will cancel. We will be left with an expression that we can integrate to find g(y). g( y)     dg ( y )     dy  C2   h( y )dy  C2   Q( x, y )    P( x, y )dx  dy  C2 [34] dy y  y  const      We can now substitute this expression for g(y) into equation [31] for f.       f   df   P( x, y )dx   Q( x, y )    P( x, y )dx  dy  C2  C1 y  y  const    y  const   [35] We can combine the two constants into a single constant, C, and remove references to the unwanted function, f, to obtain the following relationship between x and y.       P ( x , y ) dx  Q ( x , y )  P ( x , y ) dx    dy  C    y  y  const   y  const   [36] To see how we could apply this, consider the exact example in equation [28]. That example had P(x,y) = x2 + y2 and Q(x,y) = 2xy. We first find the integral of P(x,y)dx at constant y.  P( x, y)dx  y  const  x y  const x  y dx   y x 3 3 2 2 2 [37] To evaluate the function g(y) we have to take the partial derivative of equation [37], with respect to y, and subtract the result from Q(x,y). h( y )    dg ( y )     x3  Q( x, y )    P( x, y )dx   2 xy    y 2 x   2 xy  2 yx  0 dy y  y  const y  3   [38] Having h(y) = 0 gives a particularly simple result for our solution.     x3   2  P( x, y)dx   Q( x, y)  y   P( x, y)dx dy  3  y x  0  C y  const  y  const   [39] We can solve this equation for y as follows. C x2 y  x 3 We can verify that this is a solution of our original differential equation, [40] dy x2  y2  as dx 2 xy follows. Taking the first derivative of equation [40] and using equation [40] to introduce y into the result gives. dy 1  C x 2      dx 2  x 3  1 / 2 1  C 2 x2   C 2x  1  C 2x     [41]        2    3  2 y  x 2 3  2 yx  x 3   x We can rearrange this equation, using equation [39] to replace C, to obtain the result below; this result shows that our solution satisfies the original differential equation. dy 1  1  x3 2 x2  1  x2 2 x2  x2  y 2 2  2     y x     y   dx 2 yx  x  3 2 yx  3 3  2 yx  3  [42] In summary, the steps used for solving equations using the exact form are outlined below. To solve dy P ( x, y ) P Q , first check to see if P and Q satisfy the equation: . If this   y x dx Q ( x, y ) equation is not satisfied you cannot use this method. If it is satisfied, proceed with the steps below. 1. Integrate P(x,y)dx holding y constant. Call this result u(x,y) =  P( x, y )dx . y const 2. Obtain ∂u/∂y by taking the y derivative of the result from step 1. 3. Find h(y) defined as Q(x,y) – ∂u/∂y. If the result for h(y) contains the variable x, there is some error in the calculations so far. 4. Integrate the expression just found for h(y) to obtain ∫h(y)dy. 5. The solution to the differential equation is u(x,y) + ∫h(y)dy = C, where C is found from the initial condition. Integrating factors – It is sometimes possible to obtain a factor, F, which can convert an inexact differential equation in the form of equation [22] into an exact form. This F factor is used to multiply the entire equation when the values of P and Q in the equation do not satisfy the requirements for using the exact-form method. If we multiply both P and Q by a factor, F, we obtain the following differential equation in place of equation [22]. dy FP( x, y )  dx FQ( x, y) [43] Equation [43] will have the exact form if the following equation holds. ( FP) ( FQ)  y x [44] We want to find a factor F that will satisfy equation [44]. Once we find such a factor, we can use the method outlined above for solving the exact form where we replace P and Q by FP and FQ. Integrating factors can be found by trial and error. In certain cases, the integrating factor will be a function of only one variable (x or y). In the derivations below, we show how to compute the integrating factor when this occurs. Consider first the case where the result will be a function of x only. If we apply the product rule to equation [44], divide the result by FQ, set ∂F/∂y = 0 and rearrange, we obtain the following result. P F P Q F P F F Q y y y x x 1 F 1 dF d ln F 1  P Q         R( x) F x F dx dx Q  y x  F [45] If (∂P/∂y – ∂Q/∂x)/Q is a function of x only, we can integrate equation [45] to obtain the integrating factor as ln F = ∫R(x)dx or F  e R ( x ) dx [46] In a similar fashion, we can obtain the following result if we assume F is a function of y only. P F P Q F P F F Q y y y x x 1 F 1 dF d ln F 1  Q P         S ( y) F y F dy dy P  x y  F [47] If (∂Q/∂x – ∂P/∂y)/P is a function of y only, we can integrate equation [47] to obtain the integrating factor as ln F = ∫S(y)dy or F  e S ( y ) dy [48] To apply the results of equations [46] or [48] we must first assume that the integrating factor is a function of x only or y only. To do this we compute (∂P/∂y – ∂Q/∂x)/Q; if this is a function of x only then we can apply equation [46] to get F. If this does not work, we can compute (∂Q/∂x – ∂P/∂y)/P. If this is a function of y only then we can apply equation [48]. Once we find the integrating factor, we then have to apply the solution process for the exact form. If neither approach works, then any possible integrating factor is a function of both x and y. To illustrate this, we will derive the solution for equation [15], dy  f ( x) y  g ( x) . dx We can write this in the form of equation [23] as follows. dy   f ( x) y  g ( x)dx  0 [49] Comparing this with equation [23] we see that we have defined P(x,y) = f(x)y-g(x) and Q(x,y) = 1, giving ∂P/∂y = f(x) and ∂Q/∂x = 0. So, the equation as proposed is not an exact form. Let’s see if we have an integrating factor that is a function of x only. To do this, we see if R(x) as defined in equation [45] is truly a function of x only. R( x)   1  P Q  1       f ( x) 0   f ( x) Q  y x  1   Since R(x) is a function of x only, we can apply equation [46] to obtain the integrating factor. [50] F  e R ( x ) dx  e f ( x ) dx [51] Applying this integrating factor to P and Q, we have the following values for FP and FQ. FP   f ( x) y  g ( x)e FQ  e f ( x ) dx f ( x ) dx [52] We can now apply the steps for integrating the first-order differential equation corresponding to the exact form, following the steps outlined on page 9, replace P and Q in those steps by the values of FP and FQ shown in equation [52].  FP( x, y)dx  1. Integrate FP(x,y)dx holding y constant. Here we obtain u(x,y) = y const   f ( x) y  g ( x)e y  const  f ( x ) dx dx  y f ( x)e  f ( x ) dx dx  g ( x)e  f ( x ) dx dx   We can simplify the integral with f(x) by noting that that  f ( x )e  f ( x ) dx u ( x, y)  ye  dx  de  f ( x ) dx f ( x ) dx   g ( x )e   e f ( x ) dx d  e   f ( x ) dx f ( x ) dx  dx 3. Find h(y) defined as FQ(x,y) – ∂u/∂y. In this example, we have f ( x ) dx  e f ( x ) dx so . This gives the following expression: 2. Obtain ∂u/∂y by taking the y derivative of the result from step 1. Here e  f ( x ) dx f ( x)dx e  f ( x ) dx u  e y FQ  u  y  0  h( y) 4. Integrate the expression just found for h(y) to obtain ∫h(y)dy. Since we found h(y) = 0 in the previous step, its indefinite integral will also be zero. 5. The solution to the differential equation is u(x,y) + ∫h(y)dy = C, where C is found from the initial condition. With the result for u(x,y) from part one and ∫h(y)dy = 0, we have u ( x, y)   h( y)dy  C  Dividing this equation by e f ( x ) dx ye  f ( x ) dx   g ( x )e  f ( x ) dx dx  C and rearranging gives the solution for y.  f ( x ) dx  f ( x ) dx ye  C   g ( x )e  dx    This is the result that was given previously in equation [16]. Existence and uniqueness of solutions to first order equations Linear or nonlinear first-order differential equations may be written in the form dy/dx = f(x,y) with the initial condition that y(x0) = y0. Although we have solved such equations above, more complex first-order differential equations may have to be solved numerically. Regardless of the complexity of the equation we are interested in knowing if the equation has a solution and if the solution is unique. The existence and uniqueness of the solutions to this differential equation depend on the properties of f(x,y) and its partial derivative f/y. • In order for solutions to exist, f(x,y) must be continuous and |f(x,y)| must be less than some number, say K. • For the solutions to be unique, the partial derivative f/y must be continuous and |f/y| must be less than some other number, say M. • The function f(x,y) and derivative f/y must be continuous in a rectangular region, R, about the initial point, (x0, y0); the rectangular region is defined by the area |x – x0| < a and |y – y0| < b. • The solution to the differential equation exists for at least all x in the interval |x – x0| < , where  is the minimum of a and b/K. We are generally interested in having a solution in some range |x – x0| around the initial condition. We see that we will have a solution so long as f(x,y) is continuous for |x – x0| < a and the maximum absolute value of the derivative is less than the ratio b/K. In this ratio b represents the maximum expected value of y and K is the maximum expected value of |f(x,y)|. The solution we obtain will be unique so long as |f/y| remains bounded in the region of the solution. Second order equations Second-order differential equations are among the most common in mechanical engineering applications. Many of these equations arise from Newton’s second law of motion, F = ma, where the acceleration is the second derivative of the displacement. We start by considering linear, second-order differential equations. The most general such equation has the following form. d2y dy  p ( x)  q ( x) y  r ( x) 2 dx dx [53] Here we assume that, if the physical model has a factor multiplying the second derivative, we can divide the entire equation by that factor. The resulting equation has no factor multiplying the second derivative. A simpler class of differential equations results if the right hand side term, r(x) is zero. This is called the linear, second-order, homogenous differential equation. d2y dy  p ( x )  q( x) y  0 dx dx 2 [54] For this linear homogenous differential equation, we have the general result that any linear combination of linearly independent solutions to the equation is also a solution. For example, if y1 and y2 are solutions, then the linear combination y = c 1y1 + c2y2 is also a solution. This is true for the linear homogenous equation only. For the solution of second-order, linear, homogenous equations, we will generally have two linearly independent solutions that provide a basis for all solutions of the differential equation. These two independent solutions can be added in the form shown above, y = c 1y1 + c2y2, to give a general solution to the differential equation. The values of c 1 and c2 are then found by fitting initial conditions or boundary conditions for the problem. Two such conditions are required. Initial conditions typically specify the value of y and its first derivative at some (initial) value of x. Boundary conditions specify the value of y at two different x locations. Constant coefficients – The easiest case to consider is the equation with constant coefficients shown below. d2y dy  a  by  0 2 dx dx [55] Two solutions to this equation are shown below. y1  e1x y2  e  2 x [56] Where 1 and2 are the roots to the following quadratic equation.  a  a 2  4b  a a2 1 ,  2    b 2 2 4 [57] The general solution to equation [55] is a linear combination of the two solutions in equation [56]. y  C1 y1  C2 y2  C1e1x  C2 2 x [58] The two solutions in equation [56] will not be linearly independent if 1 =2; this will occur if a2 = 4b so that 1 =2 = -a/2 and y1 = e-ax/2. In this case, we use a method known as reduction of order to find the second solution. We start by writing the second solution, y2, in terms of the first solution, y1, and an unknown function, u(x). The derivation shown below finds an expression for u(x) that gives us our second solution. y2  uy1  ue1x  ue  ax / 2 [59] Substituting this equation into equation [55] gives. d 2 y2 dy2 d 2 y1 d 2u dy1 du dy1 du a  by2  u 2  y1 2  2 y1  ay1  buy1  0 2 dx dx dx dx dx dx dx dx [60] We can combine the three terms multiplied by u to get a factor which is the same as the original differential equation.  d 2 y1  dy1 d 2u du  dy1  u 2  a  by1   y1 2   2  ay1   0 dx dx  dx dx   dx  [61] Since y1 is a solution of the differential equation in [56], the term in brackets is zero. Setting this term to zero and substituting y1 = e-ax/2 and dy1/dx = (–a/2) e-ax/2 gives the result, shown below, that e-ax/2d2u/dx2 = 0. e  ax / 2 2 d 2u du   a  ax / 2   ax / 2   ax / 2 d u  2  e 0   ae e dx 2 dx   2 dx 2   [62] Equation [62] can only be satisfied if d 2u  0  u  Ax  B dx 2 [63] Since y2 = uy1, we have the following solution for y2. y2  uy1   Ax  B e ax / 2 [64] The general solution to equation [55] when a2 = 4b is given by a linear combination of the solution in equation [64] and y1 = e-ax/2. Since the solution for y1 is contained as a linear factor in the solution for y2, we can use the following pair of solutions for equation [55] in the double root case, when a2 = 4b. Both solutions are then used to give the general solution. y1  e  ax / 2 y2  xe  ax / 2 y  C1 y1  C2 y2  C1e  ax / 2  C2 xe  ax / 2  C1  C2 x e  ax / 2 [65] When a2 < 4b, the argument of the square root in equation [57] is negative and we will have complex values for 1 and2. In this case, we can define the argument of the square root as –2, and use this definition to write the values for 1 and2 as shown below, where i2 = –1. a2  b 4 2 1 ,  2  a  i 2 [66] We can get a modified form of the solution in equation [56] in this case, that gives a better indication of the actual behavior. To do this we use the Euler relationship for complex exponentials. ei  cos   i sin  Substituting equation [66] into equation [56] gives the following result.2 2 In the final step, we use the trigonometric relations that cos(–x) = cos(x) and sin(–x) = –sin(x) [67] y1  e 1x  e  ax / 2eix  e  ax / 2 cos   i sin  y2  e 2 x  e  ax / 2ei  e  ax / 2 cos x   i sin  x   e  ax / 2 cos x  i sin x [68] The general solution is a linear combination of the two solutions in equation [68]. y  C1 y1  C2 y2  C1e  ax / 2 cos x  i sin x  C2e  ax / 2 cos x  i sin x y  e  ax / 2 A cos x  B sin x [69] In the final step of this derivation, we have defined A = C1 + C2 and B = i(C1 - C2). However, in practice, we can use the final form of equation [69] as the general solution when a2 < 4b and use initial or boundary conditions to determine the constants A and B. If we are given the initial values of y and dy/dx as y0 and v0, respectively, then we can find the constants in the general solution for each of the three cases considered above: (1) two distinct real roots, (2) the double root, and (3) two complex roots. For two distinct, real roots, equation [58] gives the following equations for the initial conditions.  ( 0) y0  y (0)  C1e 1 ( 0)  C2 2  C1  C2 dy  ( 0) v0   1C1e 1 ( 0)   2C2 2  1C1   2C2 dx x  0 [70] Solving for C1 and C2 gives. C1   2 y 0  v0  2  1 C2  v0  1 y 0  2  1 [71] Substituting these results into equation [58] gives the general solution for two distinct real roots in terms of the initial conditions on y and dy/dx. y  C1e1x  C2 2 x   2 y0  v0 1x v0  1 y0 2 x e  e  2  1  2  1 [72] When there is a double root, the solution is given by equation [65]. Using that equation for the initial conditions on y and dy/dx gives the following result. v0  dy dx x 0 y 0  y (0)  C1  C 2 (0) e  a ( 0) / 2  C1 a  a  C1  C 2 (0)   e a ( 0) / 2  C 2 e a ( 0) / 2   C1  C 2 2  2 Here, C`1 = y0, and C2 = v0 + ay0/2, and the solution for y is. [73] ay     y   y0   v0  0  x e  ax / 2 2     [74] Finally, in the case of complex solutions, equation [69] gives the following equations for the initial conditions on y0 and v0. y0  e  a ( 0) / 2 A cos (0)  B sin (0)  A aA  a v0    e  a ( 0) / 2 A cos (0)  B sin (0)  e  a ( 0) / 2B cos (0)  A sin (0)    B 2  2 [75] This gives A = y0 and B = (v0 + ay0/2)/ so that the general solution for the specified initial conditions is. ay  1   y  e  ax / 2  y0 cos x   v0  0  sin x   2    [76] Linear combinations of sine and cosine of x can be written in terms of cos(x – ) by using the trigonometric formula for the cosine of the difference of two angles. C cosx    C cos  cos x  C sin  sin x C cosx    A cos x  B sin x [77] The two expressions for Ccos(x – ) in the above equations are equivalent if the following two equations hold. A  C cos  B  C sin  [78] The relationships between the constants A and B for the sine and cosine expression and the constants C and  for the cos(x – ) expression are shown below. A2  B 2  C 2 cos2   C 2 sin 2   C 2  C 2  A2  B 2 B C sin  B   tan     tan 1 A C cos  A [79] [80] We can rewrite equation [76] as follows y  Ce  ax / 2 cosx   where the values of C and  can be found by comparing equation [76] with equations [80] and [81]. [81] C y02 ay  1   2  v0  0  2     1  ay    tan 1   v0  0  2   y0  2 [82] The solutions can be cast into dimensionless forms. In the case where the initial displacement, y 0 is nonzero, we can divide by this initial condition to obtain a solution in terms of y/y0. If y0 = 0, then v0 must be nonzero to have a solution other than y = 0. In this case, we can divide the solution for y by av0/b to get a solution in terms of the dimensionless quantity by/av0. Here we consider only the case where y0 is nonzero and divide equation [74] for the double root by y0 to obtain y   v0 a    ax / 2   2v0  ax   ax / 2  1     x  e  1    1  e y0   y0 2   ay   0  2 [83] This equation gives the dimensionless form for y/y0 as a function of ax/2 with 2v0/ay0 as a parameter. Dividing equation [76] for the trigonometric solution by y0 gives.    v y a   sin x   e  ax / 2 cos x   0  y0  y0 2    [84] From the definition of w in equation [66], we can write a2  b 4 2   a 4b 1 2 a2 [85] Substituting this result into equation [84] gives. 1 / 2   ax 4b   v   ax 4b  y 4b    ax / 2 0  e  1    1  1  [86] cos   sin  2 2   y  a 2 y0 2 2 a a      0    In this case the dimensionless form gives y/y0 as a function of ax/2 with two dimensionless parameters: v0/y0 and 4b/a2 – 1. Non-homogenous equations The solution to a linear non-homogenous equation such as the second-order equation shown below, d2y dy  p( x)  q( x) y  r ( x) 2 dx dx can be written in terms of the solution, yH, to the homogenous equation [87] dy H d 2 yH  p ( x )  q ( x) y H  0 dx dx 2 [88] The total solution is the sum of the homogenous solution and a particular solution, yP, that satisfies equation [87]. y  yH  yP [89] We can see that this is the solution to equation [87] by substituting equation [89] into equation [87]. This gives d 2  yH  yP  d  yH  yP    q ( x ) y H  y P   p ( x ) dx dx 2 d 2 yH dy H d 2 yP dy P  p ( x )  q ( x ) y   p ( x )  q ( x) y P  r ( x) H dx dx dx 2 dx 2 [90] Since yH satisfies equation [88], the first three terms in the second row of equation [90] are zero and the remaining terms give the result defined for yP: yP satisfies equation [87]. The solution to the non-homogenous equation, then, proceeds by first finding the solution to the homogenous equation then by finding the particular solution. An important part of this process is that the arbitrary constants in the homogenous solution should not be determined until the final solution, the sum of the homogenous and particular solution is found. One method for finding the particular solution is known as the method of undetermined coefficients. In this method, a trial solution for yP is proposed using the trial solutions shown in the table below. For these r(x) Start with this yP as a trial solution ax r(x) = Ae n r(x) = Ax r(x) = Asin x ax yP = Be n yP = a0 + a1x + … + anx r(x) = Acos x ax sin x ax r(x) = Ae cos x r(x) = Ae yP = B sin x + C cos x yP = e ax (B sin x + C cos x) If r(x) contains more than one of the terms shown on the left, include the corresponding yP terms in the general solution for yP. If r(x) contains an nth order polynomial in x, yP should include a polynomial with all possible powers of x from x 0 to xn. If any term in r(x) is proportional to part of the solutions for yH multiply the proposed yP in the table above by x. E.g., if both r(x) and yH have a term in eax, with the same value of a, yP should contain a term in xeax instead of eax. The solutions proposed for yP in the table above have undetermined coefficients. These coefficients are found by substituting the proposed solution in to the differential equation and equating coefficients of like terms on both sides of the resulting equation. For example, we can apply the method of undetermined coefficients to the solution of the following equation. d2y dy  3  2 y  x2 2 dx dx [91] First, we find the solution to the homogenous equation. d 2 yH dy H  3  2y  0 dx dx 2 [92] We have previously show that the solution to this equation is given by equations [55] and [56], where we have to find the roots of the characteristic equation from equation [57]. For this problem, those roots are found as follows.  a  a 2  4b  3  32  4(2) 1 ,  2    1,2 2 2 [93] Thus, the solution to the homogenous equation is given by the following equation. yH  C1e t  C2e 2t [94] Since r(x) is a second order polynomial in x, we have to use the following equation for yP(x). yP  a0  a1 x  a2 x 2 [95] Substituting this particular solution into the original differential equation in [91] gives   d 2 yP dy  3 P  2 yP  2a2  3a1  2a2 x   2 a0  a1 x  a2 x 2  x 2 2 dx dx [96] Setting terms in like powers of on both sides of the equation equal to each other gives. x 0 terms : 2a2  3a1  2a0  0 x1 terms : x 2 terms : 6a2  2a1  0 2 a2  1 [97] We can easily solve these equations, starting with the last one and working backwards, to find a2 = ½, a1 = -3/2, and a0 = 7/4. This gives the particular solution shown below. yP  7 3 1  x  x2 4 2 2 [98] You can substitute this solution into equation [91] to verify that it satisfies the differential equation. The solution to the differential equation is the sum of the particular solution just found the homogenous solution from equation [94]. 7 3 1 y  yH  yP  C1e  x  C2e  2 x   x  x 2 4 2 2 [99] Only after we have the combined solution can we match the initial or boundary conditions. If we have an initial condition on both y and dy/dx as y0 and g0, we have to satisfy the following equations. 7 3 1 7  (0)  (0)  C1  C2  4 2 2 4 dy 3 3  g 0  C1e  0  2C2e  2( 0)   (0)  C1  2C2  2 2 dx 0 y (0)  y0  C1e  0  C2e  2( 0)  [100] Solving this pair of equations for the two unknowns gives C1 = 2y0 + g0 – 2, and C2 = ¼ - y0 – g0. Substituting these values into equation [99] gives the solution to equation [91] that satisfies the initial conditions that y(0) = y0 and dy/dx|0 = g0 as follows. 7 3 1 1  y  2 y0  g 0  2e  x    y0  g 0 e  2 x   x  x 2 4 2 2 4  [101] Higher order equations with constant coefficients Higher order differential equations, with constant coefficients, can be written in the following form. d n y n 1 d m y   am  r ( x) dx n m  0 dx m [102] Here we use the notation that the zeroth derivative of a function is the function itself. The solution to this equation is given, as before, as the sum of a homogenous and particular solution. The homogenous solution is written as follows. n y H   Cm e  m x [103] m 1 Where m are the roots of the following characteristic equation. n  n 1  amm  0 m0 [104] The particular solution can be found by the method of undetermined coefficients as described previously. Symbolic solutions of differential equations in MATLAB MATLAB uses the command dsolve to obtain symbolic solutions of differential equations. The equation and the initial conditions are entered as string arguments to this function. By default, the independent variable is t and the dependent variable is y. The symbol D is used to indicate a derivative. For example, Dy means dy/dt; D2y means d2y/dt2 and so forth. The differential equation d2y/dx2 + 3 dy/dx + 2y = ex cos x is specified in MATLAB as the following formula: ‘D2y + 3 * Dy + 2 * y = exp(x) * cos(x)’. Note the use of the single quotation marks to delimit a string. The dsolve function can be called with a single argument giving the differential equation. In this case a general solution is returned with constants that have to be evaluated. You can also place initial and boundary conditions in the arguments to dsolve. Such conditions are specified by the notation y(a) = b to indicate that y has a value of b at t = a. For example, the initial condition that y is 3 at t = 0 would be indicated as ‘y(0) = 3’. Initial or boundary conditions on derivatives are indicated using the D notation. An initial condition that dy/dt = 1 at t = 0 would be entered as ‘Dy(0) = 1’. Other variables can be used as the independent and dependent variable by entering them in the strings used to define the differential equation and boundary conditions. In the first example below, the independent variable is set to x by using this variable in the definition of the differential equation. Solving the differential equation d2y/dx2 + 3 dy/dx + 2y = ex cos x with initial conditions that y(0) = 3 and dy/dx = 1 at x = 0 would be done by the following command: dsolve(‘D2y + 3 * Dy + 2 * y = exp(x) * cos(x)’, ‘y(0) = 3’, ‘Dy(0) = 1’) MATLAB produces the following answer in response to this dsolve command. 1/2*exp(x)*cos(x)-exp(-2*t)*(-1/2*cos(x)*cosh(x)1/2*cos(x)*sinh(x)+4)+exp(-t)*(-cos(x)*cosh(x)-cos(x)*sinh(x)+7) The following text was copied from a MATLAB session, using MATLAB, to solve some differential equation problems from an earlier edition of Kreyszig. The % symbol used below is the MATLAB code to indicate that the following characters are a comment. >> dsolve('D2y+y=2*t','y(0)=-1','Dy(0)=8') ans = 6*sin(t)-cos(t)+2*t >> dsolve('D2y-4*y=exp(-2*t)-2*t','y(0)=0','Dy(0)=0') ans = -1/16*exp(2*t)+1/8*exp(-2*t)+1/16*(-1+8*t*exp(2*t)-4*t)*exp(-2*t >> dsolve('D2y+25*y=24*sin(t)','y(0)=1','Dy(0)=1' ans = cos(5*t)+sin(t) >> t = 0:0.1:20; %Set t values from 0 to 20 in increments of 0.1 >> plot(t,y) %Plot y = cos(5*t)+sin(t) for t values defined above %Plot not shown here >> dsolve('D3y+3*D2y+3*Dy+y=exp(-t)*sin(t)','y(0)=2','Dy(0)=0','D2y(0)=-1') ans = exp(-t)*cos(t)+exp(-t)+exp(-t)*t^2+2*exp(-t)*t