Advanced Engineering Mathematics with MATLAB
986 Pages

Advanced Engineering Mathematics with MATLAB

Course Number: 341 235, Spring 2011

College/University: Universidad Anáhuac

Word Count: 141886


Document Preview

a general result that is extremely important. Because the integrand contains nonanalytic points along and inside the region enclosed by our two curves, as shown by the CauchyRiemann equations, the results depend upon the path taken. Since complex integrations often involve integrands that have nonanalytic points, many line integrations depend upon the contour taken. Example 1.4.2 Let us integrate the entire...

Unformatted Document Excerpt
Coursehero >> Mexico >> Universidad Anáhuac >> 341 235

Course Hero has millions of student submitted documents similar to the one
below including study guides, practice problems, reference materials, practice exams, textbook help and tutor support.

Course Hero has millions of student submitted documents similar to the one below including study guides, practice problems, reference materials, practice exams, textbook help and tutor support.

general a result that is extremely important. Because the integrand contains nonanalytic points along and inside the region enclosed by our two curves, as shown by the CauchyRiemann equations, the results depend upon the path taken. Since complex integrations often involve integrands that have nonanalytic points, many line integrations depend upon the contour taken. Example 1.4.2 Let us integrate the entire function f(z)=z2 along the two paths from z=0 to z=2+i shown in Figure 1.4.2. For the first integration, x=2y Figure 1.4.2: Contour used in Example 1.4.2. 26 Advanced Engineering Mathematics with MATLAB while along the second path we have two straight paths: z=0 to z=2 and z=2 to z=2+i. For the first contour integration, (1.4.9) (1.4.10) (1.4.11) (1.4.12) (1.4.13) For our second integration, (1.4.14) Along C2a we find that y=dy=0 so that (1.4.15) and (1.4.16) Complex Variables because x=2 and dx=0. Consequently, 27 (1.4.17) Figure 1.4.3: Contour used in Example 1.4.3. In this problem we obtained the same results from two different contours of integration. Exploring other contours, we would find that the results are always the same; the integration is path-independent. But what makes these results path-independent while the integration in Example 1.4.1 was not? Perhaps it is the fact that the integrand is analytic everywhere on the complex plane and there are no nonanalytic points. We will explore this later. Finally, an important class of line integrals involves closed contours. We denote this special subclass of line integrals by placing a circle on the integral sign: the following examples: Example 1.4.3 Let us integrate f(z)=z around the closed contour shown in Figure 1.4.3. From Figure 1.4.3, (1.4.18) Consider now Now (1.4.19) 28 Advanced Engineering Mathematics with MATLAB (1.4.20) and (1.4.21) where we used z=e integral equals zero. i around the portion of the unit circle. Therefore, the closed line Example 1.4.4 Let us integrate f(z)=1/(z a) around any circle centered on z=a. The Cauchy-Riemann equations show that f(z) is a meromorphic function. It is analytic everywhere except at the isolated singularity z=a. If we introduce polar coordinates by letting z a=re i and dz=ire id , (1.4.22) Note that the integrand becomes undefined at z=a. Furthermore, the answer is independent of the size of the circle. Our example suggests that when we have a closed contour integration it is the behavior of the function within the contour rather than the exact shape of the closed contour that is of importance. We will return to this point in later sections. Problems 1. Evaluate around the circle |z|=1 taken in the counterclockwise direction. around the square with vertices at (0, 0), (1, 0), (1, 1), and (0, 1) 2. Evaluate taken in the counterclockwise direction. 3. Evaluate 4. Evaluate 5. Evaluate C |z|dz along the right half of the circle |z|=1 from z= i to z=i. along the line y=x from ( 1, 1) to (1, 1). along the line y=x2 from (0, 0) to (1, 1). 6. Evaluate where C is (a) the upper semicircle |z|=1 and (b) the lower semicircle |z|=1. If z=re i, restrict < < . Take both contours in the counterclockwise direction. Complex Variables 1.5 THE CAUCHY-GOURSAT THEOREM 29 In the previous section we showed how to evaluate line integrations by brute-force reduction to real-valued integrals. In general, this direct approach is quite difficult and we would like to apply some of the deeper properties of complex analysis to work smarter. In the remaining portions of this chapter we introduce several theorems that will do just that. Figure 1.5.1: Diagram used in proving the Cauchy-Goursat theorem. If we scan over the examples worked in the previous section, we see considerable differences when the function was analytic inside and on the contour and when it was not. We may formalize this anecdotal evidence into the following theorem: Cauchy-Goursat theorem2: Let f(z) be analytic in a domain D and let C be a simple Jordan cur e3 inside D so that f(z) is analytic on and inside of C. Then Proof: Let C denote theut those that we can only solve numerically. Problems Find the order and state whether the following ordinary differential equations are linear or nonlinear: 1. y /y=x2+x 2. y2y =x+3 3. sin(y )=5y 4. y =y 5. y =3x2 6. (y3) =1 3y 7. y =y3 8. y 4y +5y=sin(x) 9. y +xy=cos(y ) 10. (2x+y)dx+(x 3y) dy=0 11. (1+x2)y =(1+y)2 12. yy =x(y2+1) 13. y +y+y2=x+ex 14. y +cos(x)y +y=0 15. x2y +x1/2(y )3+y=ex 16. y +xy +ey=x2 2.2 SEPARATION OF VARIABLES The simplest method of solving a first-order ordinary differential equation, if it works, is separation of variables. It has the advantage of handling both linear and nonlinear problems, especially autonomous equations.1 From integral calculus, we already met this technique when we solved the first-order differential equation First-Order Ordinary Differential Equations 75 (2.2.1) By multiplying both sides of (2.2.1) by dx, we obtain dy=f(x) dx. (2.2.2) At this point we note that the left side of (2.2.2) contains only y while the right side is purely a function of x. Hence, we can integrate directly and find that (2.2.3) For this technique to work, we must be able to rewrite the differential equation so that all of the y dependence appears on one side of the equation while the x dependence is on the other. Finally we must be able to carry out the integration on both sides of the equation. One of the interesting aspects of our analysis is the appearance of the arbitrary constant C in (2.2.3). To evaluate this constant we need more information. The most common method is to require that the dependent variable give a particular value for a particular value of x. Because the independent variable x often denotes time, this condition is usually called an initial condition, even in cases when the independent variable is not time. Example 2.2.1 Let us solve the ordinary differential equation (2.2.4) Because we can separate variables by rewriting (2.2.4) as (2.2.5) its solution is simply 1 An autonomous equation is a differential equation where the independent variable does not explicitly appear in the equation, such as y =f(y). 76 Advanced Engineering Mathematics with MATLAB ye y e y=ln |x|+C (2.2.6) by direct integration. Example 2.2.2 Let us solve (2.2.7) subject to the initial condition y(0)=1. Multiplying (2.2.7) by dx, we find that dy+y dx=xexy dx, (2.2.8) or (2.2.9) A quick check shows that the left side of (2.2.9) contains only the dependent variable y while the right side depends solely on x and we have separated the variables onto one side or the other. Finally, integrating both sides of (2.2.9), we have ln(y)=xex ex x+C. (2.2.10) Since y(0)=1, C=1 and y(x)=exp[(x 1)ex+1 x]. (2.2.11) In addition to the tried-and-true method of solving ordinary differential equation by hand, scientific computational packages such as MATLAB provide symbolic toolboxes that are First-Order Ordinary Differential Equations 77 designed to do the work for you. In the present case, typing dsolve('Dy+y=x*exp(x)*y','y(0)=1','x') yields ans = 1/exp( 1)*exp( x+x*exp(x) exp(x)) which is equivalent to (2.2.11). Our success here should not be overly generalized. Sometimes these toolboxes give the answer in a rather obscure form or they fail completely. For example, in the previous example, MATLAB gives the answer ans = lambertw((1og (x) +C1)*exp( 1)) 1 The MATLAB function lambertw is Lambert's W function, where w=lambertw(x) is the solution to wew=x. Using this definition, we can construct the solution as expressed in (2.2.6). Example 2.2.3 Consider the nonlinear differential equation x2y +y2=0. (2.2.12) Separating variables, we find that (2.2.13) Equation (2.2.13) shows the wide variety of solutions possible for an ordinary differential equation. For example, if we require that y(0)=0, then there are infinitely many different solutions satisfying this initial condition because C can take on any value. On the other hand, if we require that y(0)=1, there is no solution because we cannot choose any constant C such that y(0)=1. Finally, if we have the initial condition that y(1)=2, then there is only one possible solution corresponding to Consider now the trial solution y=0. Does it satisfy (2.2.12)? Yes, it does. On the other hand, there is no choice of C which yields this solution. The solution y=0 is called a singular solution to (2.2.12). Singular solutions are solutions to a differential equation which cannot be obtained from a solution with arbitrary constants. 78 Advanced Engineering Mathematics with MATLAB Figure 2.2.1: The solution to (2.2.13) when C= 2, 0, 2, 4. Finally, we illustrate (2.2.13) using MATLAB. This is one of MATLAB's strengths--the ability to convert an abstract equation into a concrete picture. Here the MATLAB script clear hold on x = 5:0.5:5; for c = 2:2:4 y = x ./ (c*x 1); if (c== 2) subplot(2,2,1), plot(x,y,'*') axis tight; title('c = 2'); ylabel ('y','Fontsize',20); end if (c== 0) subplot(2,2,2), plot(x,y,'^') axis tight; title('c = 0'); end if (c== 2) subplot(2,2,3), plot(x,y,'s') axis tight; title ('c = 2'); xlabel('x' ,'Fontsize',20); ylabel ('y','Fontsize',20); end if (c== 4) subplot(2,2,4), plot(x,y,'h') axis tight; title ('c = 4'); xlabel ('x','Fontsize',20); end end yields Figure 2.2.1 which illustrates (2.2.13) when C= 2, 0, 2, and 4. The previous example showed that first-order ordinary differential equations may have a unique solution, no solution, or many solutions. From a complete study2 of these equations, we have the following theorem: 2 The proof of the existence and uniqueness of first-order ordinary differential equations is beyond the scop of this book. See Ince, E.L., 1956: Ordinary Differential Equations. Dover Publications, Inc., Chapter 3. First-Order Ordinary Differential Equations 79 Theorem: Existence and Uniqueness Suppose some real-valued function f(x, y) is continuous on some rectangle in the xy-plane containing the point (a, b) in its interior. Then the initial alue problem (2.2.14) has at least one solution on the same open interval I containing the point x=a. Furthermore, if the partial derivati e f/ y is continuous on that rectangle, then the solution is unique on some (perhaps smaller) open inter al I0 containing the point x=a. Example 2.2.4 Consider the initial-value problem y =3y1/3/2 with y(0)=1. Here f(x, y)=3y1/3/2 and fy=y 2/3/2. Because fy is continuous over a small rectangle containing the point (0, 1), there is a unique solution around x=0, namely y=(x+1)3/2, which satisfies the differential equation and the initial condition. On the other hand, if the initial condition reads y(0)=0, then fy is not continuous on any rectangle containing the point (0, 0) and there is no unique solution. For example, two solutions to this initial-value problem, valid on any open interval that includes x=0, are y1(x)=x3/2 and (2.2.15) Example 2.2.5: Hydrostatic equation Consider an atmosphere where its density varies only in the vertical direction. The pressure at the surface equals the weight per unit horizontal area of all of the air from sea level to outer space. As you move upward, the amount of air remaining above decreases and so does the pressure. This is why we experience pressure sensations in our ears when ascending or descending in an elevator or airplane. If we rise the small distance dz, there must be a corresponding small decrease in the pressure, dp. This pressure drop must equal the loss of weight in the column per unit area, g dz. Therefore, the pressure is governed by the differential equation dp= g dz, commonly called the hydrostatic equation. (2.2.16) 80 Advanced Engineering Mathematics with MATLAB To solve (2.2.16), we must express in terms of pressure. For example, in an isothermal atmosphere at constant temperature Ts, the ideal gas law gives p= RTs, where R is the gas constant. Substituting this into (2.2.16) and separating variables yields (2.2.17) Integrating (2.2.17) gives (2.2.18) Thus, the pressure (and density) of an isothermal atmosphere decreases exponentially with height. In particular, it decreases by e 1 over the distance RTs/g, the so-called "scale height." Example 2.2.6: Terminal velocity As an object moves through a fluid, its viscosity resists the motion. Let us find the motion of a mass m as it falls toward the earth under the force of gravity when the drag varies as the square of the velocity. From Newton's second law, the equation of motion is (2.2.19) where denotes the velocity, g is the gravitational acceleration, and CD is the drag coefficient. We choose the coordinate system so that a downward velocity is positive. Equation (2.2.19) can be solved using the technique of separation of variables if we change from time t as the independent variable to the distance traveled x from the point of release. This modification yields the differential equation (2.2.20) First-Order Ordinary Differential Equations 81 since =dx/dt. Separating the variables leads to (2.2.21) or (2.2.22) where k=CD/m and =0 for x=0. Taking the inverse of the natural logarithm, we finally obtain (2.2.23) Thus, as the distance that the object falls increases, so does the velocity and it eventually approaches a constant value commonly known as the terminal velocity. Because the drag coefficient CD varies with the superficial area of the object while the mass depends on the volume, k increases as an object becomes smaller, resulting in a smaller terminal velocity. Consequently, although a human being of normal size will acquire a terminal velocity of approximately 120 mph, a mouse, on the other hand, can fall any distance without injury. Example 2.2.7: Interest rate Consider a bank account that has been set up to pay out a constant rate of P dollars per year for the purchase of a car. This account has the special feature that it pays an annual interest rate of r on the current balance. We would like to know the balance in the account at any time t. Although financial transactions occur at regularly spaced intervals, an excellent approximation can be obtained by treating the amount in the account x(t) as a continuous function of time governed by the equation x(t+ t) x(t)+rx(t) t P t, (2.2.24) where we have assumed that both the payment and interest are paid in time increments of t. As the time between payments tends to zero, we obtain the first-order ordinary differential equation 82 Advanced Engineering Mathematics with MATLAB (2.2.25) If we denote the initial deposit into this account by x(0), then at any subsequent time x(t)=x(0)ert P (ert 1)/r. (2.2.26) Although we could compute x(t) as a function of P, r, and x(0), there are only three separate cases that merit our close attention. If P/r>x(0), then the account will eventually equal zero at rt=ln{P/[P rx(0)]}. On the other hand, if P/r<x(0), the amount of money in the account will grow without bound. Finally, the case x(0)=P/r is the equilibrium case where the amount of money paid out balances the growth of money due to interest so that the account always has the balance of P/r. Example 2.2.8: Steady-state flow of heat When the inner and outer walls of a body, for example the inner and outer walls of a house, are maintained at different constant temperatures, heat will flow from the warmer wall to the colder one. When each surface parallel to a wall has attained a constant temperature, the flow of heat has reached a steady state. In a steady state flow of heat, each surface parallel to a wall, because its temperature is now constant, is referred to as an isothermal surface. Isothermal surfaces at different distances from an interior wall will have different temperatures. In many cases the temperature of an isothermal surface is only a function of its distance x from the interior wall, and the rate of flow of heat Q in a unit time across such a surface is proportional both to the area A of the surface and to dT/dx, where T is the temperature of the isothermal surface. Hence, (2.2.27) where is called the thermal conductivity of the material between the walls. In place of a flat wall, let us consider a hollow cylinder whose inner and outer surfaces are located at r=r1 and r=r2, respectively. At steady state, (2.2.27) becomes (2.2.28) assuming no heat generation within the cyli means: M(tx, ty)=tnM(x, y) and N(tx, ty)=tnN(x, y). For example, the ordinary differential equation (x2+y2) dx+(x2 xy) dy=0 (2.3.2) is a homogeneous equation because both coefficients are homogeneous functions of degree 2: M(tx, ty)=t2x2+t2y2=t2(x2+y2)=t2M(x, y), and N(tx, ty)=t2x2 t2xy=t2(x2 xy)=t2N(x, y). (2.3.4) (2.3.3) Why is it useful to recognize homogeneous ordinary differential equations? Let us set y=ux so that (2.3.2) becomes (x2+u2x2) dx+(x2 ux2)(u dx+x du)=0. (2.3.5) First-Order Ordinary Differential Equations 91 Then, x2(1+u) dx+x3(1 u) du=0, (2.3.6) (2.3.7) or (2.3.8) Integrating (2.3.8), u+2ln|1+u|+ln|x|=ln|c|, (2.3.9) (2.3.10) (2.3.11) or (x+y)2=cxey/x. (2.3.12) Problems First show that the following differential equations are homogeneous and then find their solution. Then use MATLAB to plot your solution. Try and find the symbolic solution using MATLAB's dsolve. 92 Advanced Engineering Mathematics with MATLAB 1. 2. 3. 4. 5. 6. 7. 8. y =sec(y/x)+y/x y =ey/x+y/x. 2.4 EXACT EQUATIONS Consider the multivariable function z=f(x, y). Then the total derivative is (2.4.1) If the solution to a first-order ordinary differential equation can be written as f(x, y)=c, then the corresponding differential equation is M(x, y)dx+N(x, y) dy=0. (2.4.2) How do we know if we have an exact equation (2.4.2)? From the definition of M(x, y) and N(x, y), (2.4.3) if M(x, y) and N(x, y) and their first-order partial derivatives are continuous. Consequently, if we can show that our ordinary differential equation is exact, we can integrate First-Order Ordinary Differential Equations 93 (2.4.4) to find the solution f(x, y)=c. Example 2.4.1 Let us check and see if [y2 cos(x) 3x2y 2x] dx+[2y sin(x) x3+ln(y)] dy=0 (2.4.5) is exact. Since M (x, y)=y2 cos(x) 3x2y 2x, and N(x, y)=2y sin(x) x3+ln(y), we find that (2.4.6) and (2.4.7) Because Nx=My, (2.4.5) is an exact equation. Example 2.4.2 Because (2.4.5) is an exact equation, let us find its solution. Starting with (2.4.8) direct integration gives f(x, y)=y2sin(x) x3y x2+g(y). Substituting (2.4.9) into the equation fy=N, we obtain (2.4.9) 94 Advanced Engineering Mathematics with MATLAB (2.4.10) Thus, g (y)=ln(y), or g(y)=y ln(y) y+C. Therefore, the solution to the ordinary differential equation (2.4.5) is y2 sin(x) x3y x2+y ln(y) y=c. (2.4.11) Example 2.4.3 Consider the differential equation (x+y)dx+x ln(x) dy=0 (2.4.12) on the interval (0, ). A quick check shows that (2.4.12) is not exact since (2.4.13) However, if we multiply (2.4.12) by 1/x so that it becomes (2.4.14) then this modified differential equation is exact because (2.4.15) Therefore, the solution to (2.4.12) is x+y ln(x)=C. (2.4.16) This mysterious function that converts an inexact differential equation into an exact one is called an integrating factor. Unfortunately there is no general rule for finding one unless the equation is linear. First-Order Ordinary Differential Equations 95 Problems Show that the following equations are exact. Then solve them, using MATLAB to plot them. Finally try and find the symbolic solution using MATLAB's dsolve. 1. 2xyy =x2 y2 2. (x+y)y +y=x 3. (y2 1) dx+[2xy sin(y)] dy=0 4. [sin(y) 2xy+x2] dx+[x cos(y) x2] dy=0 5. y dx/x2+(1/x+1/y) dy=0 6. (3x2 6xy) dx (3x2+2y) dy=0 7. y sin(xy) dx+x sin(xy) dy=0 8. (2xy2+3x2) dx+2x2y dy=0 9. (2xy3+5x4y) dx+(3x2y2+x5+1) dy=0 10. (x3+y/x) dx+[y2+ln(x)] dy=0 11. [x+e y+x ln(y)] dy+[y ln(y)+ex] dx=0 12. cos(4y2) dx 8xy sin(4y2) dy=0 13. sin2(x+y) dx cos2(x+y) dy=0 14. Show that the integrating factor for (x y)y + y(1 y)=0 is (y)= ya/(1 y)a+2, a+1=1/ . Then show that the solution is 2.5 LINEAR EQUATIONS In the case of first-order ordinary differential equations, any differential equation of the form (2.5.1) is said to be linear. Consider now the linear ordinary differential equation (2.5.2) 96 Advanced Engineering Mathematics with MATLAB or (2.5.3) Let us now multiply (2.5.3) by x 4. (How we knew that it should be x 4 and not something else will be addressed shortly.) This magical factor is called an integrating factor because (2.5.3) can be rewritten (2.5.4) or (2.5.5) Thus, our introduction of the integrating factor x 4 allows us to use the differentiation product rule in reverse and collapse the right side of (2.5.4) into a single x derivative of a function of x times y. If we had selected the incorrect integrating factor, the right side would not have collapsed into this useful form. With (2.5.5), we may integrate both sides and find that (2.5.6) or (2.5.7) or y=x4(x 1)ex+Cx4. (2.5.8) From this example, it is clear that finding the integrating factor is crucial to solving firstorder, linear, ordinary differential equations. To do this, let us first rewrite (2.5.1) by dividing through by a1(x) so that it becomes First-Order Ordinary Differential Equations 97 (2.5.9) or dy+[P(x)y Q(x)] dx=0. (2.5.10) If we denote the integrating factor by (x), then (x)dy+ (x)[P(x)y Q(x)] dx=0. (2.5.11) Clearly, we can solve (2.5.11) by direct integration if it is an exact equation. If this is true, then (2.5.12) or (2.5.13) Integrating (2.5.13), (2.5.14) Note that we do not need a constant of integration in (2.5.14) because (2.5.11) is unaffected by a constant multiple. It is also interesting that the integrating factor only depends on P(x) and not Q(x). We can summarize our findings in the following theorem. Theorem: Linear First-Order Equation If the functions P(x) and Q(x) are continuous on the open interval I containing the point x0, then the initial-value problem 98 Advanced Engineering Mathematics with MATLAB has a unique solution y(x) on I, given by with an appropriate value of C and (x) is defined by (2.5.14). The procedure for implementing this theorem is as follows: Step 1: If necessary, divide the differential equation by the coefficient of dy/dx. This gives an equation of the form (2.5.9) and we can find P(x) by inspection. Step 2: Find the integrating factor by (2.5.14). Step 3: Multiply the equation created in Step 1 by the integrating factor. Step 4:Run the derivative product rule in reverse, collapsing the left side of the differential equation into the form d[ (x)y]/dx. If you are unable to do this, you have made a mistake. Step 5: Integrate both sides of the differential equation to find the solution. The following examples illustrate the technique. Example 2.5.1 Let us solve the linear, first-order ordinary differential equation xy y=4x ln(x). (2.5.15) We begin by dividing through by x to convert (2.5.15) into its canonical form. This yields (2.5.16) From (2.5.16), we see that P(x)=1/x. Consequently, from (2.5.14), we have that (2.5.17) Multiplying (2.5.16) by the integrating factor, we find that (2.5.18) First-Order Ordinary Differential Equations 99 or (2.5.19) Integrating both sides of (2.5.19), (2.5.20) Multiplying (2.5.20) through by x yields the general solution y=2x ln2(x)+Cx. (2.5.21) Although it is nice to have a closed form solution, considerable insight can be gained by graphing the solution for a wide variety of initial conditions. To illustrate this, consider the MATLAB script clear % use symbolic toolbox to solve (2.5.15) y = dsolve('x*Dy y=4*x*log(x)' , 'y(1) = c','x'); % take the symbolic version of the solution % and convert it into executable code solution = inline(vectorize(y),'x','c'); close all; axes; hold on % now plot the solution for a wide variety of initial conditions x = 0.1:0.1:2; for c = 2:4 if (c= = 2) plot(x, solution(x,c),'.'); end if (c= = 1) plot(x, solution(x, c),'o'); end if (c= = 0) plot(x, solution(x, c),'x'); end if (c= = 1) plot(x, solution(x, c),'+'); end if (c= = 2) plot(x, solution(x, c),'*'); end if (c= = 3) plot(x, solution(x, c),'s'); end if (c= = 4) plot(x, solution(x, c),'d'); end end axis tight xlabel ('x','Fontsize ',20); ylabel ('y','Fontsize',20) legend ('c = 2','c = 1','c = 0','c = 1',... 'c = 2','c = 3','c = 4'); legend boxoff This script does two things. First, it uses MATLAB's symbolic toolbox to solve (2.5.15). Alternatively we could have used (2.5.21) and introduced it as a function. The second portion of this script plots this solution for e slope equals f(x, y). This graphical representation is known as the direction field or slope field of (2.6.1). Starting with the initial point (x0, y0), we can then construct the solution curve by extending the initial line segment in such a manner that the tangent of the solution curve parallels the direction field at each point through which the curve passes. Before the days of computers, it was common to first draw lines of constant slope (isoclines) or f(x, y)=c. Because along any isocline all of the line segments had the same 112 Advanced Engineering Mathematics with MATLAB slope, considerable computational savings were realized. Today, computer software exists which perform these graphical computations with great speed. To illustrate this technique, consider the ordinary differential equation (2.6.2) Its exact solution is x(t)=Cet+t2+2t+2, (2.6.3) where C is an arbitrary constant. Using the MATLAB script clear % create grid points in t and x [t , x] = meshgrid( 2:0.2:3, 1:0.2:2); % load in the slope slope = x t.*t; % find the length of the vector (1, slope) length = sqrt(1 + slope .* slope) ; % create and plot the vector arrows quiver(t, x, 1./length,slope./length, 0.5) axis equal tight hold on % plot the exact solution for various initial conditions tt = [ 2:0.2:3]; for cva1 = 10:1:10 x_exact = cval * exp(tt) + tt.*tt + 2*tt + 2; plot(tt, x_exact) xlabel('t' , 'Fontsize',20) ylabel('x' , 'Fontsize',20) end we show in Figure 2.6.1 the directional field associated with (2.6.2) along with some of the particular solutions. Clearly the vectors are parallel to the various particular solutions. Therefore, without knowing the solution, we could choose an arbitrary initial condition and sketch its behavior at subsequent times. The same holds true for nonlinear equations. Rest points and autonomous equations In the case of autonomous differential equations (equations where the independent variable does not explicitly appear in the equation), considerable information can be gleaned from a graphical analysis of the equation. First-Order Ordinary Differential Equations 113 Consider the nonlinear ordinary differential equation (2.6.4) Figure 2.6.1: The direction field for (2.6.2). The solid lines are plots of the solution with various initial conditions. The time derivative x vanishes at x= 1, 0, 1. Consequently, if x(0)=0, x(t) will remain zero forever. Similarly, if x(0)=1 or x(0)= 1, then x(t) will equal 1 or 1 for all time. For this reason, values of x for which the derivative x is zero are called rest points, equilibrium points, or critical points of the differential equation. The behavior of solutions near rest points is often of considerable interest. For example, what happens to the solution when x is near one of the rest points x= 1, 0, 1? Consider the point x=0. For x slightly greater than zero, x <0. For x slightly less than 0, x >0. Therefore, for any initial value of x near x=0, x will tend to zero. In this case, the point x=0 is an asymptotically stable critical point because whenever x is perturbed away from the critical point, it tends to return there again. Turning to the point x=1, for x slightly greater than 1, x >0; for x slightly less than 1, x <0. Because any x near x=1, but not equal to 1, will move away from x=1, the point x=1 is called an unstable critical point. A similar analysis applies at the point x= 1. This procedure of determining the behavior of an ordinary differential equation near its critical points is called a graphical stability analysis. Phase line A graphical representation of the results of our graphical stability analysis is the phase line. On a phase line, the equilibrium points are denoted by circles. See Figure 2.6.2. Also on the phase line we identify the sign of x for all values of x. From the sign of x , we then indicate whether x is increasing or deceasing by an appropriate arrow. If the arrow points toward the right, x is increasing; toward the left x decreases. Then, by knowing the sign of the derivative for all values of x, together with the starting value of x,(x)x3 ln(x)+u(x) [3x2 ln(x)+x2], and y =u (x)x3 ln(x)+2u (x) [3x2 ln(x)+x2]+u(x) [6x ln(x)+5x]. (3.0.9) (3.0.10) Substitution of y(x), y (x), and y (x) into (3.0.8) yields x5 ln(x)u +[x4 ln(x)+2x4] u =0. (3.0.11) Setting u = , separation of variables leads to (3.0.12) Note how our replacement of u (x) with (x) has reduced the second-order ordinary differential equation to a first-order one. Solving (3.0.12), we find that 128 Advanced Engineering Mathematics with MATLAB (3.0.13) and (3.0.14) Because y(x)=u(x)x3 ln(x), the complete solution is y(x)=C1x3+C2x3 ln(x). (3.0.15) Substitution of (3.0.15) into (3.0.8) confirms that we have the correct solution. We can verify our answer by using the symbolic toolbox in MATLAB. Typing the command: dsolve('x*x*D2y-5*x*Dy+9*y=0','x') yields ans = C1*x^3+C2*x^3*log(x) In summary, we can reduce (in principle) any higher-order linear ordinary differential , equations into a system of first-order ordinary differential equations. This system of differential equations can then be solved using techniques from the previous chapter. In Chapter 14 we will pursue this idea further. Right now, however, we will introduce methods that allow us to find the solution in a more direct manner. Example 3.0.3 An autonomous differential equation is one where the independent variable does not appear explicitly. In certain cases we can reduce the order of the differential equation and then solve it. Consider the autonomous ordinary differential equation y =2y3. (3.0.16) The trick here is to note that Higher-Order Ordinary Differential Equations 129 (3.0.17) where v=dy/dx. Integrating both sides of (3.0.17), we find that 2 =y4+C1. (3.0.18) Solving for v, (3.0.19) Integrating once more, we have the final result that (3.0.20) Problems For the following differential equations, use reduction of order to find a second solution. Can you obtain the general solution using dsolve in MATLAB? 1. xy +2y =0, y1(x)=1 2. y +y 2y=0, y1(x)=ex 3. x2y +4xy 4y=0, y1(x)=x 4. xy (x+1)y +y=0, y1 (x)=ex 5. (2x x2)y +2(x 1)y 2y=0, y1(x)=x 1 6. y +tan(x)y 6 cot2(x)y=0, y1(x)=sin3(x) 7. 8. y +ay +b(1+ax bx2)y=0, y1(x)=e bx2 /2 Solve the following autonomous ordinary differential equations: 9. yy =y 2 10. y =2yy , y(0)=y (0)=1 11. yy =y +y 2 12. 2yy =1+y 2 13. y =e2y, y(0)=0, y (0)=1 14. 15. Solve the nonlinear second-order ordinary differential equation 130 Advanced Engineering Mathematics with MATLAB by (1) reducing it to the Bernoulli equation (2) solving for (x), and finally (3) integrating u =v to find u(x). 16. Consider the differential equation a2(x)y + 1(x)y + 0(x)y=0, a2(x) 0. Show that this ordinary differential equation can be rewritten using the substitution 3.1 HOMOGENEOUS LINEAR EQUATIONS WITH CONSTANT COEFFICIENTS In our drive for more efficient methods to solve higher-order, linear, ordinary differential equations, let us examine the simplest possible case of a homogeneous differential equation with constant coefficients: (3.1.1) Although we could explore (3.1.1) in its most general form, we will begin by studying the second-order version, namely y +by +cy=0, (3.1.2) since it is the next step up the ladder in complexity from first-order ordinary differential equations. Higher-Order Ordinary Differential Equations 131 Motivated by the fact that mx the solution ax to the first-order ordinary differential equation y +ay=0 is y(x)=C1e , we make the educated guess that the solution to (3.1.2) is y(x)=Ae . Direct substitution into (3.1.2) yields (am2 +bm+c) Ae mx =0. (3.1.3) Because A 0 or we would have a trivial solution and since emx 0 for arbitrary x, (3.1.3) simplifies to am2 +bm+c=0. (3.1.4) Equation (3.1.4) is called the auxiliary or characteristic equation. At this point we must consider three separate cases. Distinct real roots In this case the roots to (3.1.4) are real and unequal. Let us denote these roots by m=m1, and m=m2. Thus, we have the two solutions: (3.1.5) We will now show that the most general solution to (3.1.2) is (3.1.6) This result follows from the p linearly dependent on the real line? To find out, we compute the Wronskian or (3.1.42) (3.1.43) Therefore, x, xex, and x2ex are linearly independent. Having introduced this concept of linear independence, we are now ready to address the question of how many linearly independent solutions a homogeneous linear equation has. Theorem: On any interval I over which an n-th order homogeneous linear differential equation is normal, the equation has n linearly independent solutions y1(x), y2(x),..., yn(x) and any particular solution of the equation on I can be expressed as a linear combination of these linearly independent solutions. Proof : Again for convenience and clarity we prove this theorem for the special case of n=2. Let y1(x) and y2(x) denote solutions on I of (3.1.2). We know that these solutions exist by the existence theorem and have the following values: (3.1.44) at some point a on I. To establish the linear independence of y1 and y2 we note that, if C1y1(x)+C2y2(x)=0 holds identically on I, then Because x=a lies in I, we have that C1y1(a)+C2y2(a)=0, there too. (3.1.45) and 140 Advanced Engineering Mathematics with MATLAB (3.1.46) which yields C1=C2=0 after substituting (3.1.44). Hence, the solutions y1 and y2 are linearly independent. To complete the proof we must now show that any particular solution of (3.1.2) can be expressed as a linear combination of y1 and y2. Because y, y1, and y2 are all solutions of (3.1.2) on I, so is the function Y(x)=y(x) y( )y1(a) y (a)y2(x), (3.1.47) where y( ) and y (a) are the values of the solution y and its derivative at x=a. Evaluating Y and Y at x=a, we have that Y (a)=y(a) y( )y1( ) y (a)y2(a)=y(a) y(a)=0, (3.1.48) and (3.1.49) Thus, Y is the trivial solution to (3.1.2). Hence, for every x in I, y(x) y(a)y1(x) y (a)y2(x)=0. (3.1.50) Solving (3.1.50) for y(x), we see that y is expressible as the linear combination y(x)=y(a)y1(x)+y (a)y2(x) (3.1.51) of y1 and y2, and the proof is complete for n=2. Problems Find the general solution to the following differential equations. Check your general solution by using dsolve in MATLAB. 1. y +6y +5y=0 Higher-Order Ordinary Differential Equations 141 2. y 6y +10y=0 3. y 2y +y=0 4. y 3y +2y=0 5. y 4y +8y=0 6. y +6y +9y=0 7. y +6y 40y=0 8. y +4y +5y=0 9. y +8y +25y=0 10. 4y 12y +9y=0 11. y +8y +16y=0 12. y +4y =0 13. y +4y =0 14. y +2y +y =0 15. y 8y=0 16. y 3y +3y y =0 17. The simplest differential equation with "memory"--its past behavior affects the present--is Solve this integro-differential equation by differentiating it with respect to t to eliminate the integral. 3.2 SIMPLE HARMONIC MOTION Second-order, linear, ordinary differential equations often arise in mechanical or electrical problems. The purpose of this section is to illustrate how the techniques that we just derived may be applied to these problems. We begin by considering the mass-spring system illustrated in Figure 3.2.1 where a mass m is attached to a flexible spring suspended from a rigid support. If there were no spring, then the mass would simply fall downward due to the gravitational force mg. Because there is no motion, the gravitational force must be balanced by an upward force due to the presence of the spring. This upward force is usually assumed to obey Hooke's law which states that the restoring force is opposite to the direction of elongation and proportional to 142 Advanced Engineering Mathematics with MATLAB Figure 3.2.1: Various configurations of a mass/spring system. The spring alone has a length L which increases to L+s when the mass is attached. During simple harmonic motion, the length of the mass/spring system varies as L+s+x. the amount of elongation. Mathematically the equilibrium condition can be expressed mg=ks. Consider now what happens when we disturb this equilibrium. This may occur in one of two ways: We could move the mass either upward or downward and then release it. Another method would be to impart an initial velocity to the mass. In either case, the motion of the mass/spring system would be governed by Newton's second lalution is (3.4.13) and the general solution is (3.4.14) We can verify our result by using the symbolic toolbox in MATLAB. Typing the command: dsolve('D2y-2*Dy+y=x+sin(x)','x') yields ans = x+2+1/2*cos(x)+C1*exp(x)+C2*exp(x)*x Example 3.4.3 Let us find the particular solution to y +y 2y=xex (3.4.15) by the method of undetermined coefficients. From the form of the right side of (3.4.15), we guess the particular solution 158 Advanced Engineering Mathematics with MATLAB yp(x)=Axex+Bex. (3.4.16) Therefore, (3.4.17) and (3.4.18) Substituting into (3.4.15), we find that 3Aex=xex. (3.4.19) Clearly we cannot choose a constant A such that (3.4.19) is satisfied. What went wrong? To understand why, let us find the homogeneous or complementary solution to (3.4.15); it is yH(x)=C1e 2x +C2ex. (3.4.20) Therefore, one of the assumed particular solutions, Bex, is also a homogeneous solution and cannot possibly give a nonzero left side when substituted into the differential equation. Consequently, it would appear that the method of undetermined coefficients does not work when one of the terms on the right side is also a homogeneous solution. Before we give up, let us recall that we had a similar situation in the case of linear homogeneous second-order ordinary differential equations when the roots from the auxiliary equation were equal. There we found one of the homogeneous solution was We eventually found that the second solution was Could such a solution work here? Let us try. We begin by modifying (3.4.16) by multiplying it by x. Thus, our new guess for the particular solution reads yp(x)=Ax2ex+Bxex. (3.4.21) Higher-Order Ordinary Differential Equations 159 Then, (3.4.22) and (3.4.23) Substituting (3.4.21) into (3.4.15) gives (3.4.24) Grouping together terms that vary as xex, we find that 6A=1. Similarly, terms that vary as ex yield 2A+3B=0. Therefore, (3.4.25) so that the general solution is (3.4.26) In summary, the method of finding particular solutions to higher-order ordinary differential equations by the method of undetermined coefficients is as follows: G Step 1: Find the homogeneous solution to the differential equation. G Step 2: Make an initial guess at the particular solution. The form of yp(x) is a linear combination of all linearly independent functions that are generated by repeated differentiations of f(x). G Step 3: If any of the terms in yp(x) given in Step 2 duplicate any of the homogeneous solutions, then that particular term in yp(x) must be multiplied by xn, where n is the smallest positive integer that eliminates the duplication. 160 Advanced Engineering Mathematics with MATLAB Example 3.4.4 Let us apply the method of undetermined coefficients to solve y +y=sin(x) e3x cos(5x). (3.4.27) We begin by first finding the solution to the homogeneous version of (3.4.27): (3.4.28) Its solution is yH(x)=A cos(x)+B sin(x). (3.4.29) To find the particular solution we examine the right side of (3.4.27) or f(x)=sin(x) e3x cos(5x). (3.4.30) Taking a few derivatives of f(x), we find that f (x)=cos(x) 3e3x cos(5x)+5e3x sin(5x), (3.4.31) f (x)= sin(x) 9e3x cos(5x)+30e3x sin(5x)+25e3x cos(5x), (3.4.32) and so forth. Therefore, our guess at the particular solution is yp(x)=Cx sin(x)+Dx cos(x)+Ee3x cos(5x)+Fe3x sin(5x). (3.4.33) Why have we chosen x sin (x) and x cos (x) rather than sin (x) and cos (x)? Because sin(x) and cos(x) are homogeneous solutions to (3.4.27), we must multiply them by a power of x. Higher-Order Ordinary Differential Equations 161 Since (3.4.34) (3.4.35) =sin(x) e3x cos(5x). (3.4.36) Therefore, 2C=0, 2D=1, 30F 15E= 1, and 30E+15F=0. Solving this system of and Thus, the general solution is (3.4.37) equations yields C=0, Problems Use the method of undetermined coefficients to find the general solution of the following differential equations. Verify your solution by using dsolve in MATLAB. 1. y +4y +3y=x+1 2. y y=ex 2e 2x 3. y +2y +2y=2x2+2x+4 4. y +y =x2+x 5. y +2y =2x+5 e 2x 6. y 4y +4y=(x+1)e2x 7. y + 4y +4y=xex 8. y 4y=4sinh(2x) 9. y + 9y=x cos(3x) 10. y +y=sin(ave yielded closed form solutions. To show that this is not necessarily so, let us solve y 4y=e2x/x (3.6.53) by variation of parameters. Higher-Order Ordinary Differential Equations 181 Again we begin by solving the homogeneous differential equation (3.6.54) which has the solution yH(x)=Ae2x+Be 2x . (3.6.55) Thus, our two independent solutions are y1(x)=e2x and y2(x)=e Therefore, the particular solution equals yp(x)=e2xu1(x)+e 2x 2x . u2(x). (3.6.56) From (3.6.14), we have that (3.6.57) while (3.6.58) Solving for and we find that (3.6.59) or (3.6.60) and 182 Advanced Engineering Mathematics with MATLAB (3.6.61) or (3.6.62) Therefore, the general solution is y(x)=Ae2x+Be 2x +e2xu1(x)+e 2x u2(x) (3.6.63) (3.6.64) Problems Use variation of parameters to find the general solution for the following differential equations. Then see if you can obtain your solution by using dsolve in MATLAB. 1. y 4y +3y=e x 2. y y 2y=x 3. y 4y=xex 4. y +9y=2sec(x) 5. y +4y +4y=xe 2x 6. y +2ay =sin2( x) 7. y 4y +4y=(x+1)e2x 8. y 4y=sin2(x) 9. y 2y +y=ex/x 10. y +y=tan(x) 3.7 EULER-CAUCHY EQUATION The Euler-Cauchy or equidimensional equation is a linear differential equation of the form (3.7.1) Higher-Order Ordinary Differential Equations 183 where an, n 1 ,..., a0 are constants. The important point here is that in each term the power to which x is raised equals the order of differentiation. To illustrate this equation, we will focus on the homogeneous, secondorder, ordinary differential equation (3.7.2) The solution of higher-order ordinary differential equations follows by analog. If we wish to solve the nonhomogeneous equation (3.7.3) we can do so by applying variation of parameters using the complementary solutions that satisfy (3.7.2). Our analysis starts by trying a solution of the form y=xm, where m is presently undetermined. The first and second derivatives are (3.7.4) respectively. Consequently, substitution yields the differential equation (3.7.5) =am(m 1)xm+bmxm+cxm (3.7.6) =[ m(m 1)+bm+c]xm. (3.7.7) Thus, y=xm is a solution of the differential equation whenever m is a solution of the auxiliary equation am(m 1)+bm+c=0, or am2+(b a)m+c=0. (3.7.8) 184 Advanced Engineering Mathematics with MATLAB At this point we must consider three different cases which depend upon the values of a, b, and c. Distinct real roots Let m1 and m2 denote the real roots of (3.7.8) such that m1 m2. Then, (3.7.9) are homogeneous solutions to (3.7.2). Therefore, the general solution is (3.7.10) Repeated real roots If the roots of (3.7.8) are repeated [m1=m2= (b a)/2], then we presently have only one solution, To construct the second solution y2, we use reduction in order. We begin by first rewriting the Euler-Cauchy equation as (3.7.11) Letting P(x)=b/(ax), we have (3.7.12) (3.7.13) (3.7.14) The general solution is then Higher-Order Ordinary Differential Equations 185 (3.7.15) For higher-order equations, if m1 is a root of multiplicity k, then it can be shown that are the k linearly independent solutions. Therefore, the general solution of the differential equation equals a linear combination of these k solutions. Conjugate complex roots If the roots of (3.7.8) are the complex conjugate pair m1= +i , and m2= and are real and >0, then a solution is y(x)=C1x +i i , where +C2xa i . (3.7.16) However, because xi =[eln(x)]i =ei ln(x) , we have by Euler's formula (3.7.17) xi =cos[ ln(x)]+i sin [ ln (x)], and x i =cos [ ln(x)] i sin [ ln(x)], (3.7.18) Substitution into (3.7.16) leads to y(x)=C3x cos [ ln(x)]+C4x sin [ ln(x)], (3.7.19) where C3=C1+C2, and C4=iC1 iC2. Example 3.7.1 Let us find the general solution to 186 Advanced Engineering Mathematics with MATLAB x2y +5xy 12y=ln(x) (3.7.20) by the method of undetermined coefficients and variation of parameters. In the case of undetermined coefficients, we begin by letting t=ln(x) and y(x)=Y(t). Substituting these variables into (3.7.20), we find that Y +4Y 12Y=t. (3.7.21) The homogeneous solution to (3.7.21) is YH(t)=A e 6t +B e2t, (3.7.22)ly placed tire treads has its largest amplitude at small n. This produces one loud tone plus strong harmonic overtones because the fundamental and its overtones are the dominant terms in the Fourier series representation. 212 Advanced Engineering Mathematics with MATLAB Clearly this loud, monotone whine is undesirable. How might we avoid it? Just as soldiers mardhing in step produce a loud uniform sound, we suspect that our uniform tread pattern is the problem. Therefore, let us now vary the interval between the treads so that the distance between any tread and its nearest neighbor is not equal as illustrated in Figure 4.1.6. Again we perform its Fourier analysis and obtain that Figure 4.1.6: Temporal spacing and frequency spectrum of nonuniformly spaced snow tire treads. (4.1.40) (4.1.41) (4.1.42) Fourier Series 213 (4.1.43) (4.1.44) and (4.1.45) (4.1.46) (4.1.47) The MATLAB script epsilon = pi/12; % set up parameter for fs coefficient n = 1:20; % number of harmonics arg1 = (pi/2)*n; arg2 = (pi/4)*n; arg3 = epsi1on*n; % compute the fourier coefficient a_n an = (cos(arg1) + cos(arg2)).*sin (arg3); an = (2/pi) *an./n; % compute the fourier coefficient b_n bn = (sin(arg2) - sin(arg1)).*sin(arg3); bn = (2/pi) *bn. /n; % compute the magnitude cn = 0.5 *sqrt (an.*an + bn.*bn); % add in the a_0 term cn = [2*epsilon/pi,cn]; n = [0,n]; clf % clear any figures axes ('FontSize',20) % set font size stem (n,cn,'filled') % plot spectrum set(gca,'PlotBoxAspectRatio', [8 4 1]) % set aspect ratio xlabel('n') % label x-axis ylabel('(a_n^2 + b_n^2)^{1/2}/2') % label y-axis, 214 Advanced Engineering Mathematics with MATLAB was used to compute the amplitude of each harmonic as a function of n and the results were plotted. See Figure 4.1.6. The important point is that our new choice for the spacing of the treads has reduced or eliminated some of the harmonics compared to the case of equally spaced treads. On the negative side we have excited some of the harmonics that were previously absent. However, the net effect is advantageous because the treads produce less noise at more frequencies rather than a lot of noise at a few select frequencies. If we were to extend this technique so that the treads occurred at completely random positions, then the treads would produce very little noise at many frequencies and the total noise would be comparable to that generated by other sources within the car. To find the distribution of treads with the whitest noise8 is a process of trial and error. Assuming a distribution, we can perform a Fourier analysis to obtain its frequency spectrum. If annoying peaks are present in the spectrum, we can then adjust the elements in the tread distribution that may contribute to the peak and analyze the revised distribution. You are finished when no peaks appear. Problems Find the Fourier series for the following functions. Using MATLAB, plot the Fourier spectrum. Then plot various partial sums and compare them against the exact function. 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. f(t)=eat, L < t < L f(t)=t+t2, L < t < L 8 White noise is sound that is analogous to white light in that it is uniformly distributed throughout the complete audible sound spectrum. Fourier Series 11. 12. 13. 14. 15. 16. 215 4.2 PROPERTIES OF FOURIER SERIES In the previous section we introduced the Fourier series and showed how to compute one given the function f(t). In this section we examine some particular properties of these series. Differentiation of a Fourier series In certain instances we only have the Fourier series representation of a function f(t). Can we find the derivative or the integral of f(t) merely by differentiating or integrating the Fourier series term by term? Is this permitted? Let us consider the case of differentiation first. Consider a function f(t) of period 2L which has the derivative f (t). Let us assume that we can expand f (t) as a Fourier series. This implies that f (t) is continuous except for a finite number of discontinuities and f(t) is continuous over an interval that starts at t= and ends at t= +2nto account, we can also write (5.1.19) as (5.1.20) Finally, from the law of the mean of integrals, we have the sifting property that (5.1.21) if a<t0<b. This property is given its name because (t t0) acts as a sieve, selecting from all possible values of f(t) its value at t=t0. We can also use several other functions with equal validity to represent the delta function. These include the limiting case of the following rectangular or triangular distributions: (5.1.22) or (5.1.23) The Fourier Transform and the Gaussian function: 275 (5.1.24) Note that the delta function is an even function. Example 5.1.3: Multiple Fourier transforms The concept of Fourier transforms can be extended to multivariable functions. Consider a two-dimensional function f(x, y). Then, holding y constant, (5.1.25) Then, holding constant, (5.1.26) Therefore, the double Fourier transform of f(x, y) is (5.1.27) assuming that the integral exists. In a similar manner, we can compute f(x, y) given F( , ) by reversing the process. Starting with (5.1.28) followed by (5.1.29) 276 Advanced Engineering Mathematics with MATLAB we find that (5.1.30) Example 5.1.4: Computation of Fourier transforms using MATLAB The Heaviside (unit) step function is a piecewise continuous function defined by (5.1.31) where a 0. We will have much to say about this very useful function in the chapter on Laplace transforms. Presently we will use it to express functions whose definition changes over different ranges of t. For example, the "top hat" function (5.1.9) can be rewritten f(t)=H(t+a) H(t a). We can see that this is correct by considering various ranges of t. For example, if t< a, both step functions equal zero and f(t)=0. On the other hand, if t>a, both step functions equal one and again f(t)=0. Finally, for a<t<a, the first step function equals one while the second one equals zero. In this case, f(t)=1. Therefore, f(t)=H(t+a) H(t a) is equivalent to (5.1.9). This ability to rewrite functions in terms of the step function is crucial if you want to use MATLAB to compute Fourier transform via the MATLAB routine fourier. For example, how would we compute the Fourier transform of the signum function? The MATLAB commands Problems 1. (a) Show that the Fourier transform of f(t)=e is |t| , a>0, Using MATLAB, plot the amplitude and phase spectra for this transform. (b) Use MATLAB's fourier to find F( ). The Fourier Transform 2. (a) Show that the Fourier transform of f(t)=te is a|t| 277 , a>0, Using MATLAB, plot the amplitude and phase spectra for this transform. (b) Use MATLAB's fourier to find F( ). 3. (a) Show that the Fourier transform of is Using MATLAB, plot the amplitude and phase spectra for this transform. (b) Use MATLAB's fourier to find F( ). 4. (a) Show that the Fourier transform of is Using MATLAB, plot the amplitude and phase spectra for this transform. (b) Rewrite f(t) in terms of step functions. Then use MATLAB's fourier to find F( ). 5. (a) Show that the Fourier transform of is 278 Advanced Engineering Mathematics with MATLAB Using MATLAB, plot the amplitude and phase spectra for this transform. (b) Rewrite f(t) in terms of step functions. Then use MATLAB's fourier to find F( ). 6. (a) Show that the Fourier transform of is Using MATLAB, plot the amplitude and phase spectra for this transform. (b) Rewrite f(t) in terms of step functions. Then use MATLAB's fourier to find F( ). 7. (a) Show that the Fourier transform of is Using MATLAB, plot the amplitude and phase spectra for this transform. (b) Rewrite f(t) in terms of step functions. Then use MATLAB's fourier to find F( ). 8. (a) Show that the Fourier transform of is Using MATLAB, plot the amplitude and phase spectra for this transform. (b) Rewrite f(t) in terms of step functions. Then use MATLAB's fourier to find F( ). The Fourier Transform 9. (a) Show that the Fourier transform of 279 is Using MATLAB, plot the amplitude and phase spectra for this transform. (b) Rewrite f(t) in terms of step functions. Then use MATLAB's fourier find to F( ). 10. (a) Show that the FouriNS BY FOURIER TRANSFORMS 321 As with Laplace transforms, we may use Fourier transforms to solve ordinary differential equations. However, this method gives only the particular solution and we must find the complementary solution separately. Consider the differential equation (5.6.1) Taking the Fourier transform of both sides of (5.6.1), (5.6.2) where we used the derivative rule (5.3.19) to obtain the transform of y and Y( )=F[y(t)]. Therefore, (5.6.3) Applying the inversion integral to (5.6.3), (5.6.4) We evaluate (5.6.4) by contour integration. For t>0 we close the line integral with an infinite semicircle in the upper half of the -plane. The integration along this arc equals zero by Jordan's lemma. Within this closed contour we have a second-order pole at z=i. Therefore, (5.6.5) (5.6.6) 322 Advanced Engineering Mathematics with MATLAB and (5.6.7) For t<0, we again close the line integral with an infinite semicircle but this time it is in the lower half of the -plane. The contribution from the line integral along the arc vanishes by Jordan's lemma. Within the contour, we have a simple pole at z= i. Therefore, (5.6.8) and (5.6.9) The minus sign in front of the 2 i results from the contour being taken in the clockwise direction or negative sense. Using the step function, we can combine (5.6.7) and (5.6.9) into the single expression (5.6.10) Note that we only found the particular or forced solution to (5.6.1). The most general solution therefore requires that we add the complementary solution Ae t, yielding (5.6.11) The arbitrary constant A would be determined by the initial condition which we have not specified. We could also have solved this problem using MATLAB. The MATLAB script The Fourier Transform 323 which is equivalent to (5.6.10). Consider now a more general problem of y +y=f(t), <t< (5.6.12) where we assume that f(t) has the Fourier transform F( ). Then the Fouriertransformed solution to (5.6.12) is (5.6.13) or y(t)=g(t)*f(t), (5.6.14) where g(t)=F l[1/(1+ i)]=e tH(t). Thus, we can obtain our solution in one of two ways. First, we can take the Fourier transform of f(t), multiply this transform by G( ), and finally compute the inverse. The second method requires a convolution of f(t) with g(t). Which method is easiest depends upon f(t) and g(t). The function g(t) can also be viewed as the particular solution of (5.6.12) resulting from the forcing function (t), the Dirac delta function, because F[ (t)]=1. Traditionally this forced solution g(t) is called the Green's function and G( ) is called the frequency response or steady-state transfer function of our system. Engineers often extensively study the frequency response in their analysis rather than the Green's function because the frequency response is easier to obtain experimentally and the output from a linear system is just the product of two transforms [see (5.6.13)] rather than an integration. 324 Advanced Engineering Mathematics with MATLAB In summary, we can use Fourier transforms to find particular solutions to differential equations. The complete solution consists of this particular solution plus any homogeneous solution that we need to satisfy the initial conditions. Convolution of the Green's function with the forcing function also gives the particular solution. Example 5.6.1: Spectrum of a damped harmonic oscillator Second-order differential equations are ubiquitous in engineering. In electrical engineering many electrical circuits are governed by second-order, linear ordinary differential equations. In mechanical engineering they arise during the application of Newton's second law. For example, in mechanics the damped oscillations of a mass m attached to a spring with a spring constant k and damped with a velocity dependent resistance are governed by the equation my +cy +ky=f(t), (5.6.15) where y(t) denotes the displacement of the oscillator from its equilibrium position, c denotes the damping coefficient, and f(t) denotes the forcing. Assuming that both f(t) and y(t) have Fourier transforms, let us analyze this systetransform of the system of ordinary differential equations results in an algebraic set of equations containing Y1(s), Y2(s),..., Yn(s). By some method we solve this set of equations and invert each transform Y1(s), Y2(s),..., Yn(s) in turn to give y1(t), y2(t),..., yn(t). The following examples will illustrate the details of the process. Example 6.8.1 Let us solve the ordinary differential equation y +2y =8t, (6.8.2) subject to the initial conditions that y (0)=y(0)=0. Taking the Laplace transform of both sides of (6.8.2), (6.8.3) or (6.8.4) where Y(s), Substituting the initial conditions into (6.8.4) and solving for (6.8.5) (6.8.6) The Laplace Transform 385 Matching powers of s in the numerators of (6.8.6), C+D=0, B+2C=0, A+2B=0, and 2A=8 or A=4, B= 2, C=1 and D= 1. Therefore, (6.8.7) Finally, performing term-by-term inversion of (6.8.7), the final solution is y(t)=2t2 2t+1 e 2t . (6.8.8) We could have performed the same operations using the symbolic toolbox with MATLAB. The MATLAB script clear % define symbolic variables syms s t Y % take Laplace transform of left side of differential equation LHS = laplace (diff (diff (sym ( 'y (t) ' ) ) ) +2*diff (sym('y(t)'))); % take Laplace transform of right side of differential equation RHS = laplace(8*t); % set Y for Laplace transform of y % and introduce initial conditions newLHS = subs(LHS,'laplace(y(t),t,s)','y(0)','D(y)(0)',Y,0,0); % solve for Y Y = solve (newLHS-RHS, Y) ; % invert Laplace transform and find y(t) y = ilaplace(Y,s,t) yields the result y = 1 exp( 2*t) 2*t+2*t^2 which agrees with (6.8.8). Example 6.8.2 Let us solve the ordinary differential equation y +y=H(t) H(t 1) (6.8.9) 386 Advanced Engineering Mathematics with MATLAB with the initial conditions that y (0)=y(0)=0. Taking the Laplace transform of both sides of (6.8.9), (6.8.10) where Y(s), Substituting the initial conditions into (6.8.10) and solving for (6.8.11) Using the second shifting theorem, the final solution is y(t)=1 cos(t) [1 cos(t 1)]H(t 1). We can check our results using the MATLAB script clear % define symbolic variables syms s t Y % take Laplace transform of left side of differential equation LHS = laplace(diff(diff(sym('y(t)'))) +sym('y(t)')); % take Laplace transform of right side of differential equation RHS = laplace('Heaviside(t) Heaviside(t 1)',t,s); % set Y for Laplace transform of y % and introduce initial conditions newLHS = subs(LHS,'laplace(y(t),t,s)','y(0)','D(y)(0)',Y,0,0); % solve for Y Y = solve(newLHS-RHS,Y); % invert Laplace transform and find y(t) y = ilaplace (Y,s,t) (6.8.12) which yields y = 1 cos(t) Heaviside(t 1)+Heaviside(t 1)*cos(t 1) The Laplace Transform Example 6.8.3 Let us solve the ordinary differential equation y +2y +y=f(t) 387 (6.8.13) with the initial conditions that y (0)=y(0)=0, where f(t) is an unknown function whose Laplace transform exists. Taking the Laplace transform of both sides of (6.8.13), s2Y(s) sy(0) y (0)+2sY(s) 2y(0)+Y(s)=F(s), (6.8.14) where Y(s), Substituting the initial conditions into (6.8.14) and solving for (6.8.15) We wrote (6.8.15) in this form because the transform Y(s) equals the product of two transforms 1/(s+1)2 and F(s). Therefore, by the convolution theorem we can immediately write (6.8.16) Without knowing f(t), this is as far as we can go. Example 6.8.4: Forced harmonic oscillator Let us solve the simple harmonic oscillator forced by a harmonic forcing y + 2 y=cos( t), (6.8.17) subject to the initial conditions that y (0)=y(0)=0. Although the complete solution could be found by summing the complementary solution and a particular solution obtained, say, from the method of undetermined coefficients, we now illustrate how we can use Laplace transforms to solve this problem. Taking the Laplace transform of both sides of (6.8.17), substituting in the initial conditions, and solving for Y(s), 388 Advanced Engineering Mathematics with MATLAB (6.8.18) and (6.8.19) Equation (6.8.19) gives an oscillation that grows linearly with time although the forcing function is simply periodic. Why does this occur? e response for a discrete system with a transform function given by (7.5.21). num = [1 0 0]; % enter the coefficients of the denominator % of the transfer function (7.5.21) den = [1 -1.5 0.5]; % create the Kronecker delta sequence x = [1 zeros(1, 20)]; % find the impulse response 472 Advanced Engineering Mathematics with MATLAB y = filter(num,den,x); % plot impulse response plot(y,'o'),axis([0 20 0.5 2.5]) xlabel('n+1','Fontsize',20) ylabel('impulse response', 'Fontsize',20) Figure 7.5.2 shows the computed impulse response. The asymptotic limit is two, so the system is marginally stable as we found before. We note in closing that the same procedure can be used to find the inverse of any ztransform which consists of a ratio of two polynomials. Here we simply set G(z) equal to the given z-transform and perform the same analysis. Problems For the following time-discrete systems, find the transfer function and determine whether the systems are unstable, marginally stable, or stable. Check your answer by graphing the impulse response using MATLAB. 1. yn=yn 1Hn 1+xn 2. yn=2yn 1Hn 1 yn 2Hn 2+xn 3. yn=3yn 1Hn 1+xn 4. Chapter 8 The Hilbert Transform In addition to the Fourier, Laplace, and z-transforms, there are many other linear transforms which have their own special niche in engineering. Examples include Hankel, Walsh, Radon, and Hartley transforms. In this chapter we consider the Hilbert transform which is a commonly used technique for relating the real and imaginary parts of a spectral response, particularly in communication theory. We begin our study of Hilbert transforms by first defining them and then exploring their properties. Next, we develop the concept of the analytic signal. Finally, we explore a property of Hilbert transforms that is frequently applied to data analysis: the Kramers-Kronig relationship. 8.1 DEFINITION In Chapter 7 we motivated the development of z-transforms by exploring the concept of the ideal sampler. In the case of Hilbert transforms, we introduce another fundamental operation, namely quadrature phase shifting or the ideal Hilbert transformer. This procedure does nothing more than shift the phase of all input frequency components by /2. Hilbert transformers are frequently used in communication systems and signal processing; examples include the generation of single-sideband modulated signals and radar and speech signal processing. Because a /2 phase shift is equivalent to multiplying the Fourier transform of a signal by e i /2= i, and because phase shifting must be an odd function of frequency,1 the transfer function of the phase shifter is G( )= i sgn( ), where sgn( ) is defined by (5.2.11). In other words, if X( ) denotes the input spectrum to the phase shifter, the output spectrum must be i sgn( )X( ). If the process is repeated, the total phase shift is , a complete phase reversal of all frequency components. The output spectrum then equals [ i sgn( )]2X( )= X( ). This agrees with the notion of phase reversal because the output function is x(t). Consider now the impulse response of the quadrature phase shifter, g(t)= F 1[G( )]. From the definition of Fourier transforms, (8.1.1) and (8.1.2) 1 For a real function the phase of its Fourier transform must be an odd function of . 474 Advanced Engineering Mathematics with MATLAB Since G ( )= 2i ( ), the corresponding impulse response is (8.1.3) Consequently, if x(t) is the input to a quadrature phase shifter, the superposition integral gives the output time function as (8.1.4) as the Hilbert transform of x(t), although some authors use the We shall define negative of (8.1.4) corresponding to a+ /2 phase shift. The transform is also called the harmonic conjugate of x(t). In similar fashion, is the Hilbert transform of the Hilbert transform of x(t) and corresponds to the output of two cascaded phase shifters. However, this output is known to be x(t), so = x(t), and we arrive at the inverse Hilbert transform relationship that (8.1.5) The Hilbert Transform 475 Figure 8.1.1: Descended from a Prussian middle-class family, David Hito the right at the speed c if we wish to keep in the same position relative to the nearest crest and trough. The quantities k, , and c are the wavenumber, frequency, and phase speed or wave-velocity, respectively. The relationship =kc holds between the frequency and phase speed. It may seem paradoxical that we are talking about traveling waves in a problem dealing with waves confined on a string of length L. Actually we are dealing with standing waves because at the same time that a wave is propagating to the right its mirror image is running to the left so that there is no resultant progressive wave motion. Figures 10.3.1 and 10.3.2 illustrate our solution. Figure 10.3.1 gives various cross sections. The single large peak at t=0 breaks into two smaller peaks which race towards the two ends. At each end, they reflect and turn upside down as they propagate back towards x=L/2 at ct/L=1. This large, negative peak at x=L/2 again breaks apart, with the two smaller peaks propagating towards the endpoints. They reflect and again become positive peaks as they propagate back to x=L/2 at ct/L=2. After that time, the whole process repeats itself. MATLAB can used to examine the solution in its totality. The script % set parameters for the calculation clear ; M = 50 ; dx = 0.02 ; dt = 0.02 ; % compute Fourier coefficients sign = 32 ; for m = 1 : M temp1 = (2*m 1)*pi ; temp2 = sin (temp1/8) ; a (m) = sign * temp2 * temp2 / (temp1 * temp1) ; sign = sign ; end The Wave Equation 581 Figure 10.3.1: The vibration of a string u(x, t)/h at various positions x/L at the times ct/L=0, 0.2, 0.4, 0.6, 0.8, and 1. For times 1<ct/L<2 the pictures appear in reverse time order. % compute grid and initialize solution X = [0 : dx : 1] ; T = [0 : dt : 2] ; u = zeros (length (T) , length (X) ) ; XX = repmat (X,[length (T) 1] ) ; TT = repmat (T',[1 length (X)] ) ; % compute solution from (10.3.33) for m = 1:M temp1 = (2*m 1)*pi; u = u + a(m) .* sin(templ*XX) .* cos(temp1*TT); end % plot space/time picture of the solution surf (XX, TT, u) xlabel ( ' DISTANCE' , ' Fontsize' , 20) ; ylabel ( 'TIME','Fontsize',20) zlabel('SOLUTION','Fontsize',20) gives a three-dimensional view of (10.3.33). The solution can be viewed in many different prospects using the interactive capacity of MATLAB. An important dimension to the vibrating string problem is the fact that the wavenumber kn is not a free parameter but has been restricted to the values of n /L. This restriction on wavenumber is common in wave problems dealing with limited domains (for example, a building, ship, lake, or planet) and these oscillations are given the special name of normal modes or natural vibrations. In our problem of the vibrating string, all of the components propagate with the same phase speed. That is, all of the waves, regardless of wavenumber 582 Advanced Engineering Mathematics with MATLAB Figure 10.3.2: Two-dimensional plot of the vibration of a string u(x, t)/h at various times ct/L and positions x/L. kn, move the characteristic distance c t or c t after the time interval t elapsed. In the next example we will see that this is not always true. Example 10.3.2: Dispersion In the preceding example, the solution to the vibrating string problem consisted of two simple waves, each propagating with a phase speed c to the right and left. In many problems where the equations of motion are a little more complicated than (10.3.1), all of the harmonics no longer propagate with the same phase speed but at a speed that depends upon the wavenumber. In such systems the phase relation varies between the harmonics and these systems are referred to as dispersive. A modification of the vibrating string problem provides a simple illustration. We now subject each element of the string to an additional applied force which is proportional to its displacement: (10.3.35)) where h>0 is constant. For example, if we embed the string in a thin sheet of rubber, then in addition to the restoring force due to tension, there is a restoring force due to the rubber on each portion of the string. From its use inditions defined for |x|< . Then illustrate your solution using MATLAB. 1. u(x, 0)=2 sin(x) cos(x) 2. u(x, 0)=x sin(x) 3. u(x, 0)=1/(x2+1) 4. u(x, 0)=e x 5. u(x, 0)=cos( x/2) 6. u(x, 0)=sin(3x) 7. Assuming that the functions F and G that ut(x, 0)=cos(x) ut(x, 0)=cos(2x) ut(x, 0)=ex ut(x, 0)=1/(x2+1) ut(x, 0)=sinh(ax) ut(x, 0)=sin(2x) sin(x) are differentiable, show by direct substitution 606 Advanced Engineering Mathematics with MATLAB and are the D'Alembert solutions to the hyperbolic system where c2=E/ and E, k, and are constants. 8. D'Alembert's solution can also be used in problems over the limited domain 0<x<L. To illustrate this, let us solve the wave equation (10.4.1) with the initial conditions u(x, 0)=0, ut(x, 0)=Vmax(1 x/L), 0<x<L, and the boundary conditions u(0, t)=u(L, t) 0, 0<t. Step 1: Show that the solution to this problem is where along with the periodicity conditions V0( )=V0( into the boundary conditions. ), and V0(L+ )= V0(L ) to take care of those cases when the argument of V0( ) is outside of (0, L). Hint: Substitute the solution Step 2: Show that at any point x within the interval (0, L), the solution repeats with a period of 2L/c if ct>2L. Therefore, if we know the behavior of the solution for the time interval 0<ct<2L, we know the behavior for any other time. Step 3: Show that the solution at any point x within the interval (0, L) and time t+L/c, where 0<ct<L, is the mirror image (about u=0) of the solution at the point L x and time t, where 0<ct<L. Step 4: Show that the maximum value of u(x, t) occurs at x=ct, where 0<x<L and when 0<ct<L. At that point, The Wave Equation 607 where umax equals the largest magnitude of u(x, t) for any time t. Plot umax as a function x and show that it a parabola. Hint: Find the maximum value of u(x, t) when 0<x ct and ct x<L with 0<x+ct<L or L<x+ct<2L. 10.5 THE LAPLACE TRANSFORM METHOD The solution of linear partial differential equations by Laplace transforms is the most commonly employed analytic technique after separation of variables. Because the transform consists solely of an integration with respect to time, the transform U(x, s) of the solution of the wave equation u(x, t) is (10.5.1) assuming that the wave equation only varies in a single spatial variable x and time t. Partial derivatives involving time have transforms similar to those that we encountered in the case of functions of a single variable. They include (10.5.2) and (10.5.3) These transforms introduce the initial conditions via u(x, 0) and ut(x, 0). On the other hand, derivatives involving x become (10.5.4) and (10.5.5) Because the transformation eliminates the time variable, only U(x, s) and its derivatives remain in the equation. Consequently, we transform the partial differential equation into a 608 Advanced Engineering Mathematics with MATLAB boundary-value problem involving an ordinary differential equation. Because this equation is often easier to solve than a partial differential equation, the use of Laplace transforms considerably simplifies the original problem. Of course, the Laplace transforms must exist for this technique to work. The following schematic summarizes the Laplace transform method: In the following examples, we illustrate transform methods by solving the classic equation of telegraphy as it applies to a uniform transmission line. The line has a resistance R, an inductance L, a capacitance C, and a leakage conductance G per unit length. We denote the current in the direction of positive x by I; V is the voltage drop across the transmission line at the point x. The dependent variables I and V are functions of both distance x along the line and time t. To derive the differential equations that govern the current and voltage in the line, consider the points A at x and B at x+ x in Figure 10.5.1. The current and voltage at A are I(x, t) and V(x, t); at B, to B is and the and Therefore, the voltage drop from A Figure 10.5.1: Schematic of an uniform transmission line. The Wave Equation 609 current in the line is Neglecting terms that are proportionalDlirpfyb',31)A 820 Advanced Engineering Mathematics with MATLAB end if (iter == 256) subplot (2,2,4), [cs,h] = contourf (X,Y,u); clabel(cs,h,'Fontsize',16) axis tight; title('after 256 iterations','Fontsize',20); xlabel ('X/L','Fontsize',20); ylabel ('Y/H','Fontsize',20); end end The initial guess everywhere except along the top boundary was zero. In Figure 12.7.1 we illustrate the numerical solution after 4, 16, 64, and 256 iterations where we have taken 11 grid points in the x and y directions. Project: Successive Over-Relaxation The fundamental difficulty with relaxation methods used in solving Laplace's equation is the rate of convergence. Assuming x= y, the most popular method for accelerating convergence of these techniques is successive over-relaxation (SOR): where Laplace's Equation 821 Figure 12.7.1: The solution to Laplace's equation by the Gauss-Seidel method. The boundary conditions are ux(0, y)=ux(L, y)=uy(x, 0)=0, and u(x, H)=1+x/L. Figure 12.7.2: The number of iterations required so that |Rm, n| 10 3 as a function of during the iterative solution of the problem posed in the project. We used x= y=0.01, and L=z0=1. The iteration count for the boundary conditions stated in Step 1 is given by the solid line while the iteration count for the boundary conditions given in Step 2 is shown by the dotted line. The initial guess equaled zero. 822 Advanced Engineering Mathematics with MATLAB Most numerical methods books dealing with partial differential equations discuss the theoretical reasons behind this technique;27 the optimum value always lies between one and two. Step 1: Write a MATLAB script that uses the Gauss-Seidel method to numerically solve Laplace's equation for 0 x L, 0 y z0 with the following boundary conditions: u(x, 0)=0, u(x, z0)=1+x/L, u(0, y)=y/z0, and u(L, y)=2y/z0. Because this solution will act as "truth" in this project, you should iterate until the solution does not change. Step 2: Now redo the calculation using successive over-relaxation. Count the number of iterations until |Rm,n| 10 3 for all m and n. Plot the number of iterations as a function of . How does the curve change with resolution x? Step 3: Redo Steps 1 and 2 with the exception of u(0, y)=u(L, y)=0. How has the convergence rate changed? Can you explain why? How sensitive are your results to the first guess? 27 For example, Young, D.M., 1971: Iterative Solution of Large Linear Systems. Academic Press, 570 pp. Chapter 13 Vector Calculus Physicists invented vectors and vector operations to facilitate their mathematical expression of such diverse topics as mechanics and electromagnetism. In this chapter we focus on multivariable differentiations and integrations of vector fields, such as the velocity of a fluid, where the vector field is solely a function of its position. 13.1 REVIEW The physical sciences and engineering abound with vectors and scalars. Scalars are physical quantities which only possess magnitude. Examples include mass, temperature, density, and pressure. Vectors are physical quantities that possess both magnitude and direction. Examples include velocity, acceleration, and force. We shall denote vectors by boldfaced letters. Two vectors are equal if they have the same magnitude and direction. From the limitless number of possible vectors, two special cases are the zero vector 0 which has no magnitude and unspecified direction and the unit vector which has unit magnitude. The most convenient method for expressing a vector analytically is in terms of its components. A vector a in three-dimensional real space is any order triplet of real numbers (components) a1, a2, and a3 such that a=a1i+ a2j+a3k, where a1i, a2j, and a3k are vectors which lie along the coordinate axes and have their origin at a common initial point. The magnitude, length, or norm of a vector a, |a|, equals A particularly important vector is the position vector, defined by r=xi+yj+zk. As in the case of scalars, certain arithmetic rules hold. Addition and subtraction are very similar to their scalar counterparts: a+b=(a1+b1)i+(a2+b2)j+(a3+2i z2j+yzk and the surface is the exterior of the hemispheric surface of x2+y2+z2=16 above the plane z=2. 9. F=yi+xj+yk and the surface is the top of the surface z=x+1, where 1 x 1 and 1 y 1. 10. F=zi+xj 3zk and the surface is the top side of the plane x+y+z=2a that lies above the square 0 x<a, 0<y<a in the xy-plane. 11. F=(y2+z2)i+(x2+z2)j+(x2+y2)k and the surface is the top side of the surface z=1 x2 with 1 x 1 and 2 y 2. 12. F=y2i+xzj k and the surface is the cone normal pointing away from the z-axis. with the Vector Calculus 855 13. F=y2i+x2j+5zk and the surface is the top side of the plane z=y+1, where 1 x 1 and 1 y 1. 14. F= yi+xj+zk and the surface is the exterior or bottom side of the paraboloid z=x2+y2, where 0 z 1. 15. F= yi+xj+6z2k and the surface is the exterior of the paraboloids z=4 x2 y2 and z=x2+y2. 13.6 GREEN'S LEMMA Consider a rectangle in the xy-plane which is bounded by the lines x=a, x=b, y=c, and y=d. We assume that the boundary of the rectangle is a piece-wise smooth curve which we denote by C. If we have a continuously differentiable vector function F=P(x, y)i+Q(x, y)j at each point of enclosed region R, then (13.6.1) (13.6.2) (13.6.3) where the last integral is a closed line integral counterclockwise around the rectangle because the horizontal sides vanish since dy=0. By similar arguments, (13.6.4) so that (13.6.5) 856 Advanced Engineering Mathematics with MATLAB This result, often known as Green's lemma, may be expressed in vector form as (13.6.6) Although this proof was for a rectangular area, it can be generalized to any simply closed region on the xy-plaue as follows. Consider an area which is surrounded by simply closed curves. Within the closed contour we can divide the area into an infinite number of infinitesimally small rectangles and apply (13.6.6) to each rectangle. When we sum up all of these rectangles, we find where the integration is over the entire surface area. On the other hand, away from the boundary, the line integral along any one edge of a rectangle cancels the line integral along the same edge in a contiguous Figure 13.6.1: Diagram for the verification of Green's lemma in Example 13.6.1. rectangle. Thus, the only nonvanishing contribution from the line integrals arises from the outside boundary of the domain Example 13.6.1 Let us verify Green's lemma using the vector field F=(3x2 8y2)i+ (4y 6xy)j and the enclosed area lies between the curves The two curves intersect at x=0 and x=1. See Figure 13.6.1. We begin with the line integral: (13.6.7) Vector Calculus 857 (13.6.8) In (13.6.7) we used y=x2 in the first integral and the areal integration, in our return integration. For (13.6.9) (13.6.10 and Green's lemma is verified in this particular case. Example 13.6.2 Let us redo Example 13.6.1 except that the closed contour is the triangular region defined by the lines x=0, y=0, and x+y=1. The line integral is (13.6.11) (13.6.12) (13.6.13) On the other hand, the areal integration is (13.6.14) 858 Advanced Engineering Mathematics with MATLAB (13.6.15 and Green's lemma is verified in this particular case. Example 13.6.3 Let us verify Green's lemma using the vector field F=(3x+4y)i+(2x 3y)j and the closed contour is a circle of radius two centered at the origin of the xy-plane. See Figure 13.6.2. Beginning with the line integration, (13.6.16) (13.6.17) (13.6.18) = 8 . (13.6.19) Figure 13.6.2: Diagram for the verification of Green's lemma in Example 13.6.3. Vector Calculus 859 For the areal integration, (13.6.20 and Green's lemma is verified in the special case. Problems Verify Green's lemma for the following two-dimensional vector fields and contours: 1. F=(x2+4y)i+(y x)j and the contour is the square bounded by the lines x=0, y=0, x=1, and y=1. 2. F=(x y)i+xyj and the contour is the square bounded by the lines x=0, y=0, x=1, and y 1. 3. F= y2i+x2j and the contour is the triangle bounded by the lines x=1, y=0, and y=x. 4. F=(xy x2)i+x2yj and the contour is the triangle bounded by the lines y=0, x=1, and y=x. 5. F=sin(y)i+xcos(y)j and the contourries with x=0 and z=3. Hence, (13.7.12) Hence, (13.7.13) Turning to the other side of the equation, (13.7.14) Our line integral has been such that the normal vector must be n=k. Therefore, (13.7.15) 864 Advanced Engineering Mathematics with MATLAB Figure 13.7.3: Diagram for the verification of Stokes' theorem in Example 13.7.2. and Stokes' theorem is verified for this special case. Example 13.7.2 Let us verify Stokes' theorem using the vector field F=(x2 y)i+4zj+ x2k, where the closed contour consists of the x and y coordinate axes and that portion of the circle x2+y2=a2 that lies in the first quadrant with z=1. See Figure 13.7.3. The line integral consists of three parts: (13.7.16) Along C1, x varies while y=0 and z=1. Therefore, (13.7.17) Along the circle C2, we use polar coordinates with x=a cos(t), y=a sin(t), and z=1. Therefore, (13.7.18) (13.7.19) Vector Calculus 865 (13.7.20) (13.7.21) Figure 13.7.4: Diagram for the verification of Stokes' theorem in Example 13.7.3. because dx= a sin(t)dt, and dy=a cos(t)dt. Finally, along C3, y varies with x=0 and z=1 Therefore, (13.7.22) so that (13.7.23) Turning to the other side of the equation, (13.7.24) 866 Advanced Engineering Mathematics with MATLAB From the path of our line integral, our unit normal vector must be n=k. Then, (13.7.25) and Stokes' theorem is verified for this case. Example 13.7.3 Let us verify Stokes' theorem using the vector field F=2yzi (x+3y 2)j+(x2+z)k, where the closed triangular region is that portion of the plane x+y+z=1 that lies in the first octant. As shown in Figure 13.7.4, the closed line integration consists of three line integrals: (13.7.26) Along C1, z=0 and y=1 x. Therefore, using x as the independent variable, (13.7.27) Along C2, x=0 and y=1 z. Thus, (13.7.28) Finally, along C3, y=0 and z=1 x. Hence, (13.7.29) Vector Calculus 867 Thus, (13.7.30) On the other hand, (13.7.31) To find n d , we use the general coordinate system x=u, y= , and z= 1 u r=ui+ j+(1 u )k and . Therefore, (13.7.32) Thus, (13.7.33) (13.7.34) (13.7.35) (13.7.36) and Stokes' theorem is verified for this case. Problems Verify Stokes' theorem using the following vector fields and surfaces: 1. F=5yi 5xj+3zk and the surface S is that portion of the plane z=1 with the square at the vertices (0, 0, 1), (1, 0, 1), (1, 1, 1), and (0, 1, 1). 2. F=x2i+y2j+z2k and the surface S is the rectangular portion of the plane z=2 defined by the corners (0, 0, 2), (2, 0, 2), (2, 1, 2), and (0, 1, 2). 3. F=zi+xj+yk and the surface S is the triangular portion of the plane z=1 defined by the vertices (0, 0, 1), (2, 0, 1), and (0, 2, 1). 868 Advanced Engineering Mathematics with MATLAB 4. F=2zi 3xj+4yk and the surface S is that portion of the plane z=5 within the cylinder x2+y2=4. 5. F=zi+xj+yk and the surface S is that portion of the plane z=3 bounded by the lines y=0, x=0, and x2+y2=4. 6. F=(2z+x)i+(y z)j+(x+y)k and the surface S is the interior of the triangularly shaped plane with vertices at (1, 0, 0), (0, 1, 0), and (0, 0, 1). 7. F=zi+xj+yk and the surface S is that portion of the plane 2x+y+2z=6 in the first octant. 8. F=xi+xzj+yk and the surface S is that portion of the paraboloid z=9 x2 y2 within the cylinder x2+y2=4. 13.8 DIVERGENCE THEOREM Although Stokes' theorem is useful in computing closed line integrals, it is usually very difficult to go the other way and convert a surface integral into a closed line integral because the integrand must have a very special form, namely In this section we introduce a theorem that allows with equal facility the conversion of a closed surface integral into a volume integral and ice ersa. Furthermore, if we can convert a given surface integral into a closed one by the introduction of a simple surface (for example, closing a hemispheric surface by adding an equatorial plate), it may be easier to use the divergence theorem and subtract off the contribution from the new surface integral rather than do the original problem. This relationship between a closed surface integral and a volume integral involving the d remaining portion of this section, we show some operations that may be performed on a determinant to introduce the desired zeros. Most of the properties follow from the expansion of determinants by cofactors. G Rule 1: For every square matrix A, det(AT)=det(A). The proof is left as an exercise. G Rule 2: If any two rows or columns of A are identical, det(A)=0. To see that this is true, consider the following 3 3 matrix: (14.2.10) G Rule 3: The determinant of a triangular matrix is equal to the product of its diagonal elements. Linear Algebra If A is lower triangular, successive expansions by elements in the first column give 891 (14.2.11) =...=a11a22...ann. (14.2.12) If A is upper triangular, successive expansions by elements of the first row prove the property. Rule 4: If a square matrix A has either a row or a column of all zeros, then det(A)=0. The proof is left as an exercise. Rule 5: If each element in one row (column) of a determinant is multiplied by a number c, the value of the determinant is multiplied by c. Suppose that |B| has been obtained from |A| by multiplying row i (column j) of |A| by c. Upon expanding |B| in terms of row i (column j) each term in the expansion contains c as a factor. Factor out the common c, the result is just c times the expansion |A| by the same row (column). Rule 6: If each element of a row (or a column) of a determinant can be expressed as a binomial, the determinant can be written as the sum of two determinants. To understand this property, consider the following 3 3 determinant: (14.2.13) The proof follows by expanding the determinant by the row (or column) that contains the binomials. Rute 7: If B is a matrix obtained by interchanging any two rows (columns) of a square matrix A, then det(B)= det(A). The proof is by induction. It is easily shown for any 2 2 matrix. Assume that this rule holds of any (n 1) (n 1) matrix. If A is n n, then let B be a matrix formed by interchanging rows i and j. Expanding |B| and |A| by a different row, say k, we have that (14.2.14) where Mks and Nks are the minors formed by deleting row k, column s from |B| and |A|, respectively. For s=1, 2,..., n, we obtain Nks and Mks by interchanging rows i and j. By the induction hypothesis and recalling that Nks and Mks are (n 1) (n 1) determinants, 892 Advanced Engineering Mathematics with MATLAB Nks= Mks for s=1, 2,..., n. Hence, |B|= |A|. Similar arguments hold if two columns are interchanged. G Rule 8: If one row (column) of a square matrix A equals to a number c times some other row (column), then det(A)=0. Suppose one row of a square matrix A is equal to c times some other row. If c=0, then |A|=0. If c 0, then |A|=c|B|, where |B|=0 because |B| has two identical rows. A similar argument holds for two columns. G Rule 9: The value of det(A) is unchanged if any arbitrary multiple of any line (row or column) is added to any other line. To see that this is true, consider the simple example: (14.2.15) where c 0. The first determinant on the left side is our original determinant. In the second determinant, we again expand the first column and find that (14.2.16) Example 14.2.2 Let us evaluate using a combination of the properties stated above and expansion by cofactors. By adding or subtracting the first row to the other rows, we have that (14.2.17) (14.2.18) Linear Algebra 893 (14.2.19) Problems Evaluate the following determinants. Check your answer using MATLAB. 1. 2. 3. 4. 5. 6. 7. 8. 9. Using the properties of determinants, show that This determinant is called Vandermonde's determinant. 10. Show that 894 Advanced Engineering Mathematics with MATLAB 11. Show that if all of the elements of a row or column are zero, then det(A)= 0. 12. Prove that det(AT)=det(A). 14.3 CRAMER'S RULE One of the most popular methods for solving simple systems of linear equations is Cramer's rule.4 It is very useful for 2 2 systems, acceptable for 3 3 systems, and of doubtful use for 4 4 or larger systems. Let us have n equations with n unknowns, Ax=b. Cramer's rule states that (14.3.1) where Ai. 5. 7. 9. yn=2n n 1 11. xn=2+( 1)n; yn=1+( 1)n 13. xn=1 2( 6)n; yn= 7( 6)n n+1 0 T)+1] )/(1 ) Answers To the Odd-Numbered Problems Section 7.5 1. marginally stable 3. unstable Section 8.1 7. Section 8.2 5. w(t)=u(t)* (t)= e Section 8.3 1. z(t)=ei Section 8.4 3. Section 9.1 1. n=(2n 1)2 2/(4L2), yn(x)=cos[(2n 1) x/(2L)] 3. 5. 7. 9. 11. (a) (b) (c) 13. n n 0 n t 1 941 sin(t) = 1, y0(x)=e x and n =n2, yn(x)=sin(nx) n cos(nx) = n4 4/L4, yn(x)=sin(n x/L) =n2 2, yn(x) sin[n ln(x)] =(2n 1)2 2/4, yn(x)=sin[(2n 1) ln(x)/2] n n 0 =0, y0(x)=1; =n2 2, yn(x)=cos[n ln(x)] =n2+1, yn(x)=sin[n ln(x)]/x n 15. =0, y0(x)=1; yn(x)=cosh( nx)+cos( nx) tanh( n)[sinh( nx)+ sin( nx)], where n=1, 2, 3,..., and Section 9.3 1. is the nth root of tanh( )= tan( ). 3. 942 Advanced Engineering Mathematics with MATLAB Section 9.4 1. 3. 5. Section 10.3 1. 3. 5. 7. 9. where kn is the nth solution of J0(2k)=0. 11. Section 10.4 1. u(x, t)=sin(2x) cos(2ct)+cos(x) sin(ct)/c 3. 5. Section 10.5 1. 3. u(x,t)=sin( x) cos( t) sin( x) sin( t)/ Answers To the Odd-Numbered Problems 5. 943 7. 9. u(x,t)=xt te x+sinh(t)e x+[1 e (t x)+t x sinh(t x)] H(t x) 11. 13. 15. 17. Section 11.3 1. 3. 5. 7. 9. 11. 944 Advanced Engineering Mathematics with MATLAB 13. 15. 17. 19. 21. 23. 25. 27. where 29. n is the nth root of tan( )=hL/k. where kn denotes the nth root of k tan(kL)=k2/a2 31. 33. where kn is the nth root of k cot(k)=1 A, n <kn <(n+1) . 35. where kn is the nth root of J0(k)=0. Answers To the Odd-Numbered Problems 37. where kn is the nth root of J0(k)=0. 39. where kn is the nth root of kJ1(kL)=hJ0(kL). Section 11.4 1. 3. 5. 7. 9. 945 11. 13. 15. 17. 946 Advanced Engineering Mathematics with MATLAB where 19. where 21. where 23. where kn is the nth root of J0(k)=0. 25. where kn is the nth root of J0(k)=0. Section 11.5 1. 3. Section 12.3 1. n n is the nth root of is the nth root of cot( )=(3a+ 2)/3a. 3. 5. 7. u(x, y)=1 9. 11. u(x, y)=1 13. u(x,y)=T0+ T cos(2 x/ )e 2 y/ Answers To the Odd-Numbered Problems 15. where kn is the nth root of J0(k)=0. 17. where kn is the nth root of J1(kb)=0. 19. where kn is the nth root of J1(k)=0. 21. 23. where kn is the nth root of kJ0(k)=J1(k). 25. where kn is the nth root of kJ1(k)=BJ0(k). 29. 31. u(r, )=T0 Section 12.4 1. 3. 947 5. 948 Advanced Engineering Mathematics with MATLAB Section 12.5 1. Section 12.6 1. Section 13.1 1. a b= 3i+19j+10k 3. a b=i 8j+7k 5. a b= 3i 2j 5k 9. 11. 13. Plane parallel to the xy plane at height of z=3, n=k 15. Paraboloid, 17. A plane, 19. A parabola of infinite extent along the y-axis, 21. y=2/(x+1); z=exp[(y 1)/y] 23. y=x; z2=y/(3y 2) Section 13.2 1. 3. 5. 7. 9. Answers To the Odd-Numbered Problems 11. 949 13. Section 13.3 1. 16/7 + 2/(3 ) 3. e2+2e8/3+e64/2 13/6 5. 4 7. 0 9. 2 Section 13.4 1. (x, y, z)=x2y+y2z+4z+constant 3. (x, y, z)=xyz+constant 5. (x, y, z)=x2 sin(y)+xe3z+4z+constant 7. (x, y, z)=xe2z+y3+constant 9. (x, y, z)=xy+xz+constant Section 13.5 1. 1/2 3. 0 5. 27/2 7. 5 9. 0 11. 40/3 13. 86/3 15. 96 Section 13.6 1. 5 3. 1 5. 0 7. 0 9. 16 11. 2 Section 13.7 1. 10 3. 2 5. 7. 45/2 950 Advanced Engineering Mathematics with MATLAB Section 13.8 1. 3. 5. 7. 3 16 4 5/12 Section 14.1 1. 3. 5. 7. 9. 11. 13. yes 15. 17. 19. yes no 21. 23. 25. 27. Answers To the Odd-Numbered Problems 29. 951 Section 14.2 1. 7 3. 1 5. 24 7. 3 Section 14.3 1. 3. x1=0, x2=0, x3= 2 Section 14.4 1. x2=2, x1=1 3. x3= , x2 , x1= 5. x3= , x2=2 , x1= 1 7. x3=2.2, x2=2.6, x1=1 9. 11. Section 14.5 1. 3. 5. 7. 952 Advanced Engineering Mathematics with MATLAB Section 14.6 1. 3. 5. 7. 9. 11. 13. 15. 17. 19. Index abscissa of convergence, 284 absolute value of a complex number, 2ff addition of complex numbers, 2 of matrices, 745 of vectors, 694 age of the earth, 608 aliasing, 219ff amplitude modulation, 246 of a complex number, 4 spectrum, 229 analytic complex function, 11 derivative of, 11 analytic signal, 417ff Archimedes' principle, 740ff argument of a complex number, 4 autonomous ordinary differential eq, 65, 110 auxiliary eq, 112 back substitution, 749, 761 band-pass functions, 415 Bernoulli equation, 92 Bessel eq of order n, 460ff function of the first kind, 462 expansion in, 466ff function of the second kind, 462 function, modified, 464 recurrence formulas, 466 Bessel, Friedrich Wilhelm, 460 Biot number, 563 boundary conditions Cauchy, 482 Dirichlet, 548 Neumann, 548 Robin, 548 boundary-value problem, 108 branches of a complex function, 11 principal, 4 Bromwich contour, 351 Bromwich integral, 351 Bromwich, Thomas J.I'A., 352 954 Index carrier frequency, 246 Cauchy boundary condition, 482 data, 482 integral formula, 28ff principle value, 54ff problem, 482 Cauchy, Augustin-Louis, 13 Cauchy-Goursat theorem, 25 Cauchy-Riemann eqs, 13 centered finite differences, 535 characteristic eq, 112 functions, 425 polynomial, 770 value, 424, 770 vector, 770 characteristics, 504 chemical reaction, 73ff circular frequency, 121 circulation, 710 closed contour integral, 23, 707 surface integral, 714 cofactor, 752 column in a matrix, 744 column vector, 744 complementary error function, 289 complementary solution of an ordinary differential eq, 132 complex conjugate, 1 envelope, 419 Fourier coefficients, 199, 439 Fourier integral, 228ff Fourier series, 198ff number, 1 plane, 4 -valued function, 9ff variable, 1 components of a vector, 693 compound interest, 71, 389ff conformable for addition of matrices, 745 for multiplication of matrices, 745 conservative field, 712 consistency in finite differencing for the heat eq, 629 for the wave eq, 537 Index consistent system of linear eqs, 760 contour integral, 20ff convergence of finite difference solution for heat eq, 630 for wave eq, 538 of a Fourier integral, 230 of Fourier series, 169 convolution theorem for Fourier transform, 269ff for Hilbert transform, 412 for Laplace transform, 317ff for z-transform, 371ff Coriolis force, 695 Cramer's rule, 756 Crank-Nicolson method, 629 critically damped, 127 critical points, 97, 158 stable, 97, 158 stable node, 160 unstable, 97, 159 cross product, 694 curl, 703 curve simply closed, 25 space, 694 cutoff frequency, 524 d'Alembert, Jean Le Rond, 504 d'Alembert's solution, 503ff damped constant, 126 damped harmonic motion, 126ff deformation principle, 26 degenerate eigenvalue problem, 433 del operator, 696 delay differential eq, 338 (Dirac) delta function, 231ff, 294ff de Moivre's theorem, 3 design of film projectors, 313ff design of wind vane, 129 determinant, 751ff diagonal, principal, 744 difference eq, 359 differential eqs, 61 first-order, 61ff linear, 63 linear, first-order, 83 linear nth-order, 107 nonlinear, 63 order, 62 955 956 Index ordinary, 61ff partial, 61 type, 61 differentiation of a Fourier series, 181 diffusivity, 547 direction fields, 95 Dirichlet conditions, 169 Dirichlet, Peter G.Lejeune-, 171 Dirichlet problem, 548 dispersion, 490 divergence theorem, 733ff of a vector, 702 division of complex numbers, 2 dot product, 694 double Fourier series, 682 dual Fourier-Bessel series, 659 dual integral eq, 654 Duhamel's theorem for the heat eq, 614ff for ordinary differential eqs, 349 eigenfunctions, 425ff expansion in, 437 orthogonality of, 434 eigenvalue of a matrix, 770 of a Sturm-Liouville problem, 424ff eigenvalue problem for matrices, 770ff for ordinary differential eqs, 423ff singular, 424 eigenvectors, 770ff orthogonality of, 778 electrical circuit, 87, 143, 333ff electrostatic potential, 644 element of a matrix, 744 elementary row operations, 760 elliptic partial differential eq, 635 entire complex function, 11 equilibrium points, 97, 158 equivalent systems of linear eqs, 760 error function, 289 essential singularity, 36 Euler-Cauchy equation, 152ff Euler's formula, 3 Euler method, 98 exact ordinary differential eq, 79 explicit numerical method for heat eq, 627 Index for wave eq, 535 exponential order, 284 existence of ordinary differential eq first-order, 69 n-order, 108 evaluation of partial sums using z-transforms, 383 fast Fourier transform, 219 filter, 222 final-value theorem for Laplace transforms, 302 for z-transforms, 369 finite differences approximation to derivatives, 535ff finite Fourier series, 212ff first shifting theorem, 297 first-order ordinary differential eq, 61 linear, 83ff flux lines, 698 folding frequency, 220 forced harmonic motion, 137ff Fourier coefficients, 168 cosine series, 175 cosine transform, 612 number, 557 series in amplitude/phase form, 195ff series of an even function, 189 series of an odd function, 189 series for a multivariable function, 202 series on [-L, L], 167ff sine series, 175 sine transform, 612 Fourier, Joseph, 170 Fourier-Bessel coefficients, 467 expansion, 466 Fourier-Legendre coefficients, 451 expansion, 451 Fourier transform, 227ff basic properties of, 241ff convolution for, 269ff inverse of, 228, 254ff method of solving the heat eq, 607 of a constant, 238 of derivatives, 245 957 958 Index of multivariable functions, 233 of sign function, 239 of step function, 240 free undamped motion, 121 frequency convolution, 271 frequency modulation, 247 frequency response, 275 frequency spectrum, 230 for a damped harmonic oscillator, 275ff for low frequency filter, 279 function even extension of, 189 generalized, 296 multivalued complex, 10 odd extension of, 189 single-valued complex, 10 vector-valued, 696ff fundamental of a Fourier series, 168 Gauss, Carl Friedrich, 734 Gauss's divergence theorem, 733ff Gaussian elimination, 763 Gauss-Seidel method, 687 general solution to ordinary differential eq, 64 generalized functions, 296 generating function for Legendre polynomials, 448 Gibbs phenomenon, 185ff gradient, 697 graphical stability analysis, 97 Green's function, 275ff, 345ff for a damped harmonic oscillator, 277 for the Klein-Gordon eq, 524 for low frequency filter, 278 Green's lemma, 722ff grid point, 534 groundwater flow, 639ff half-range expansions, 189ff Hankel transform, 653 harmonic function, 18, 636 conjugate, 18 harmonic of a Fourier series, 168 heat conduction in a rotating satellite, 207ff within a metallic sphere, 663ff heat dissipation in disc brakes, 595ff Index heat eq, 545ff for an infinite cylinder, 569ff for a semi-infinite bar, 607ff nonhomogeneous, 547 one-dimensional, 549ff within a solid sphere, 567ff Heaviside expansion theorem, 309ff step function, 291ff Heaviside, Oliver, 292 Hilbert, David, 401 Hilbert pair, 401 Hilbert transform, 399ff and convolution, 412 and derivatives, 412 and shifting, 411 and time scaling, 412 discrete, 408 linearity of, 410 product theorem, 414 holomorphic complex function, 11 homogeneous ordinary differential eq, 78, 107 solution to ordinary differential eq, 132 system of linear eqs, 748 hydraulic potential, 639 hydrostatic eq, 69 hyperbolic partial differential eq, 479 ideal Hilbert transformer, 399 ideal sampler, 361 imaginary part of a complex number, 1 impulse function, see (Dirac) delta function impulse response, 345 inconsistent system of linear eqs, 759 indicial admittance for heat eq, 617 for ordinary differential eqs, 345 inertia supercharging, 191 initial conditions, 482ff -value problem, 107, 326ff initial-boundary-value problem, 548 initial-value theorem for Laplace transforms, 301 for z-transforms, 369 inner product, 745 959 960 Index integral curve, 157 integral eq of convolution type, 321 integrals complex contour, 19 Fourier type, evaluation of, 260 line, 707ff real, evaluation of, 45 integration of a Fourier series, 182ff interest rate, 71, 389 integrating factor, 82 inverse discrete Fourier transform, 213 Fourier transform, 228, 254ff Hilbert transform, 400 Laplace transform, 309ff, 317, 350ff z-transform, 375ff inversion formula for the Fourier transform, 228 for the Laplace transform, 350ff for the z-transform, 380ff inversion of Fourier transform, 254ff by contour integration, 256 by direct integration, 254 by partial fractions, 255 inversion of Laplace transform by contour integration, 350ff by convolution, 317 by partial fractions, 309ff in amplitude/phase form, 312ff inversion of z-transform by contour integration, 380ff by partial fractions, 378ff by power series, 375ff by recursion, 376ff irrotational, 704 isoclines, 95 isolated singularities, 15 iterative meth 694 positively oriented curve, 28 potential function, 712ff power content, 184 power spectrum, 251 principal branch, 4 principal diagonal, 744 principle of linear superposition, 113, 486 quadrature phase shifting, 399 quieting snow tires, 176ff radiation condition, 483 radius of convergence, 33 rank of a matrix, 763 real part of a complex number, 1 real definite integrals, Index evaluation of, 45 recurrence relation in finite differencing, 162 for Bessel functions, 466 for Legendre polynomial, 449ff reduction in order, 109 regular complex function, 11 regular Sturm-Liouville problem, 424 relaxation methods, 687ff removable singularity, 36 residue, 35 residue theorem, 39ff resonance, 142, 205, 329 rest points, 97 Riemann, G.F.B., 14 Robin problem, 548 Rodrigues' formula, 448 root locus method, 277 roots of a complex number, 6ff row echelon form, 763 row of a matrix, 744 row vector, 744 Runge, Carl, 102 Runge-Kutta method, 101, 163 scalar, 693 Schouten-van der Pol theorem for Laplace transforms, 357 Schwarz's integral formula, 679 second shifting theorem, 298 separation of variables for heat eq, 548ff for Laplace's eq, 639ff for ordinary differential eq, 65 for Poisson's eq, 680ff for wave eq, 483ff shifting in the s variable, 297 in the t variable, 242, 298 in the variable, 246 sifting property, 233 simple closed curve, 25 eigenvalue, 426 pole, 36 simple harmonic oscillator, 120ff, 329 simply close curve, 25 sinc function, 229 single sideband signal, 419 single-valued complex function, 10 965 966 Index singular solutions to ordinary differential eq, 67 Sturm-Liouville problem, 424 singularity essential, 36 isolated, 36 pole of order n, 36 removable, 36 slope field, 95 solenoidal, 703 solution curve, 95 solution of ordinary differential eq by Fourier series, 203ff by Fourier transform, 273ff space curve, 694 spectral radius, 770 spectrum of a matrix, 770 stability of numerical methods by Fourier method for heat eq, 629 for wave eq, 538 by matrix method for wave eq, 775 steady-state heat eq, 71, 556 steady-state output, 98 steady-state solution to ordinary differential eq, 139 steady-state transfer function, 275 step function, 291ff step response, 345 Stokes, Sir George Gabriel, 727 Stokes' theorem, 726ff streamlines, 698 Sturm, Charles, 424 Sturm-Liouville eq, 423 problem, 423ff subtraction of complex numbers, 2 of matrices, 745 of vectors, 694 successive over-relaxation, 691 superposition integral for heat eq, 616ff for ordinary differential eqs, 349 superposition principle, 486 surface conductance, 562 surface integral, 714ff system of linear differential eqs, 779ff homogeneous eqs, 748 Index nonhomogeneous eqs, 748 tangent vector, 694 Taylor expansion, 33 terminal velocity, 71, 90 telegraph eq, 493, 513ff thermal conductivity, 546 threadline eq, 509ff time shifting, 242, 297 trajectories, 157 transfer function, 344 transform Fourier, 227ff Hilbert, 399ff Laplace, 283ff z-transform, 359ff transient solution to ordinary differential eq, 139 transmission line, 513ff transpose of a matrix, 747 tridiagonal matrix, solution of, 748ff underdamped, 127 underdetermined system of linear eqs, 767 uniqueness of ordinary differential eq first-order, 69 nth-order, 108 unit normal, 698 step function, 291ff vector, 693 Vandermonde's determinant, 756 variation of parameters, 145ff vector, 693, 744 vector element of area, 717 vibrating string, 489 vibrating threadline, 509 vibration of floating body, 124 Volterra eq of the second kind, 321 volume integral, 733ff wave eq, 479ff damped, 492ff for a circular membrane, 496ff for an infinite domain, 503ff one-dimensional, 481 weight function, 434 Wronskian, 118 967 968 Index zero vector, 693 z-transform, 359ff basic properties of, 367 convolution for, 371 final-value theorem for, 369 initial-value theorem for, 369 inverse of, 375ff of periodic sequences, 370 of a sequence multiplied by an exponential sequence, 367 of a shifted sequence, 367 solving of difference eqs, 386ff