### Create a StudySoup account

#### Be part of our community, it's free to join!

Already have a StudySoup account? Login here

# Ordinary Differential Equations I MATH 640

CSU

GPA 3.78

### View Full Document

## 46

## 0

## Popular in Course

## Popular in Mathematics (M)

This 230 page Class Notes was uploaded by Melvina Keeling on Monday September 21, 2015. The Class Notes belongs to MATH 640 at Colorado State University taught by Oleg Emanouilov in Fall. Since its upload, it has received 46 views. For similar materials see /class/210083/math-640-colorado-state-university in Mathematics (M) at Colorado State University.

## Similar to MATH 640 at CSU

## Reviews for Ordinary Differential Equations I

### What is Karma?

#### Karma is the currency of StudySoup.

#### You can buy or earn more Karma at anytime and redeem it for class notes, study guides, flashcards, and more!

Date Created: 09/21/15

Undergraduate Texts in Mathematics Editars S Axler KA Ribet Undergraduate Texts in Mathematics Abbott Understanding Analysis Anglin Mathematics A Concise History and Philosophy Readings in Mathematics AnglinLambek The Heritage of Thales Readings in Mathematics Apostol Introduction to Analytic Number Theory Second edition Armstrong Basic Topology Armstrong Groups and Symmetry Axler Linear Algebra Done Right Second edition Beardon Limits A New Approach to Real Analysis BakNewman Complex Analysis Second edition BanchoffWermer Linear Algebra Through Geometry Second edition Berberian A First Course in Real Analysis Bix Conics and Cubics A Concrete Introduction to Algebraic Curves Br maud An Introduction to Probabilistic Modeling Bressoud Factorization and Primality Testing Bressoud Second Year Calculus Readings in Mathematics Brickman Mathematical Introduction to Linear Programming and Game Theory Browder Mathematical Analysis An Introduction Buchmann Introduction to Cryptography Buskesvan Rooij Topological Spaces From Distance to Neighborhood Callahan The Geometry of Spacetime An Introduction to Special and General Relavitity Cartervan Brunt The Lebesgue Stieltjes Integral A Practical Introduction Cederberg A Course in Modern Geometries Second edition ChambertLoir A Field Guide to Algebra Childs A Concrete Introduction to Higher Algebra Second edition ChungAitSahlia Elementary Probability Theory With Stochastic Processes and an Introduction to Mathematical Finance Fourth edition CoxLittleO Shea Ideals Varieties and Algorithms Second edition Croom Basic Concepts of Algebraic Topology Curtis Linear Algebra An Introductory Approach Fourth edition DaeppGorkin Reading Writing and Proving A Closer Look at Mathematics Devlin The Joy of Sets Fundamentals of Contemporary Set Theory Second edition Dixmier General Topology Driver Why Math EbbinghausFlumThomas Mathematical Logic Second edition Edgar Measure Topology and Fractal Geometry Elaydi An Introduction to Difference Equations Third edition Erdt isSuranyi Topics in the Theory of Numbers Estep Practical Analysis in One Variable Exner An Accompaniment to Higher Mathematics Exner Inside Calculus FineRosenberger The Fundamental Theory of Algebra Fischer Intermediate Real Analysis FlaniganKazdan Calculus Two Linear and Nonlinear Functions Second edition Fleming Functions of Several Variables Second edition F oulds Combinatorial Optimization for Undergraduates Foulds Optimization Techniques An Introduction Franklin Methods of Mathematical Economics continued after index J David Logan A First Course in Differential Equations With 55 Figures Springer J David Logan Willa Cather Professor of Mathematics Department of Mathematics University of Nebraska at Lincoln Lincoln NE 68 5880130 USA dloganmathunledu Editorial Board S AXler KA Ribet Mathematics Department Department of Mathematics San Francisco State University University of California at Berkeley San Francisco CA 94132 Berkeley CA 947203840 USA US axlerSFSUedu ribetmathberkeleyedu Mathematics Subject Classification 2000 34xx 15xx Library of Congress Control Number 2005926697 hardcover Library of Congress Control Number 2005926698 softcover ISBN10 0387259635 hardcover ISBN13 9780387259635 ISBN10 0387259643 softcover ISBN13 978038725964 2 2006 Springer ScienceBusiness Media Inc Al rights reserve T is work may not be translated or copied in whole or in part without the written permission ofthe publisher Springer ScienceBusiness Media Inc 233 Spring Street New York NY 10013 USA except for brief excerpts in connection wit reviews or scholarly analysis Use in connection wit any orm of information storage and retrieval electronic adaptation computer software or by similar or dissimilar methodology now wn or ereafter developed is forbi en The use in this publication oftrade names trademarks service marks and similar terms even if they are not identi ie as such is not to be taken as an expression of opinion as to whether or not they are subject to proprietary rig ts Printed in the United Smtes of America SEA 9 8 7 6 5 4 3 2 1 sp ringeronline com Dedicated 107 Reece Charles Logan Jaren Logan Golightly Con ten ts Preface i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i xi To the Student i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i xiii 1 Differential Equations and Models i i i i i i i i i i i i i i i i i i i i i i i i H 1 1 1 Differential Equations i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i 2 111 Equations and Solutions i i i i i i i i i i i i i i i i i i i i i i i i i i i i i 2 112 Geometrical Interpretation i i i i i i i i i i i i i i i i i i i i i i i i i i 9 12 Pure Time Equations i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i 13 13 Mathematical Models i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i 19 1 31 Particle Dynamics i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i 21 1 32 Autonomous Differential Equations i i i i i i i i i i i i i i i i i i i 28 1 33 Stability and Bifurcation i i i i i i i i i i i i i i i i i i i i i i i i i i i i 41 1 34 Heat Transfer i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i 45 135 Chemical Reactors i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i 48 136 Electric Circuits i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i 51 2 Analytic Solutions and Approximations i i i i i i i i i i i i i i i i i i i i i 55 2 1 Separation of Variables i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i 55 2 2 FirstOrder Linear Equations i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i 61 2 3 Approximation i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i 70 2 31 Picard lteration i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i 71 2 32 Numerical Methods i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i 74 233 Error Analysis i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i H 78 viii Contents 3 SecondOrder Differential Equations l l l l l l l l l l l l l l l l l l l l l l l l 83 3 1 Particle Mechanics i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i 84 32 Linear Equations With Constant Coef cients i i i i i i i i i i i i i i i i i 87 313 The Nonhomogeneous Equation i i i i i i i i i i i i i i i i i i i i i i i i i i i i 95 31311 Undetermined Coef cients i i i i i i i i i i i i i i i i i i i i i i i i i i i 96 31312 Resonance i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i 102 34 Variable Coef cients i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i 105 31411 CauchyiEuler Equation i i i i i i i i i i i i i i i i i i i i i i i i i i i i i 106 31412 Power Series Solutions i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i 109 3143 Reduction of Order i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i 111 3144 Variation of Parameters i i i i i i i i i i i i i i i i i i i i i i i i i i i i i 112 3 5 Boundary Value Problems and Heat Flow i i i i i i i i i i i i i i i i i i i 117 3 6 HigherOrder Equations i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i 124 317 Summary and Review i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i 127 4 Laplace Transforms l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l 133 411 De nition and Basic Properties i i i i i i i i i i i i i i i i i i i i i i i i i i i i 133 412 Initial Value Problems i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i 140 413 The Convolution Property i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i 145 4 4 Discontinuous Sources i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i 149 415 Point Sources i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i 152 416 Table of Laplace Transforms i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i 157 5 Linear Systems l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l 159 511 Introduction i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i 159 5 2 Matrices i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i 165 5 3 TwoDimensional Systems i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i 179 5311 Solutions and Linear Orbits i i i i i i i i i i i i i i i i i i i i i i i i i i 179 51312 The Eigenvalue Problem i i i i i i i i i i i i i i i i i i i i i i i i i i i i 185 51313 Real Unequal Eigenvalues i i i i i i i i i i i i i i i i i i i i i i i i i i i 187 5314 Complex Eigenvalues i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i 189 51315 Real Repeated Eigenvalues i i i i i i i i i i i i i i i i i i i i i i i i i i 191 5136 Stability l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l 194 5 4 Nonhomogeneous Systems i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i 198 5 5 ThreeDimensional Systems i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i 204 6 Nonlinear Systems l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l 209 6 1 Nonlinear Models i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i 209 611 Phase Plane Phenomena i i i i i i i i i i i i i i i i i i i i i i i i i i i i 209 61112 The Lotka7Volterra Model i i i i i i i i i i i i i i i i i i i i i i i i i i i 217 61113 Holling Functional Responses i i i i i i i i i i i i i i i i i i i i i i i i 221 61114 An Epidemic Model i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i 223 Contents ix 62 Numerical Methods 229 63 Linearization and Stability 233 64 Periodic Solutions 246 641 The Poincar iBendixson Theorem 249 Appendix A References 255 Appendix B Computer Algebra Systems 257 El Maple 258 B2 MATLAB 260 Appendix C Sample Examinations 265 Appendix D Solutions and Hints to Selected Exercises 271 Index 287 Preface There are many excellent texts on elementary differential equations designed for the standard sophomore course However in spite of the fact that most courses are one semester in length the texts have evolved into calculuslike presen tations that include a large collection of methods and applications packaged with student manuals and Webbased notes projects and supplements All of this comes in several hundred pages of text with busy formats Most students do not have the time or desire to read voluminous texts and explore internet supplements The format of this differential equations book is different397 it is a onesemester brief treatment of the basic ideas models and solution methods Its limited coverage places it somewhere between an outline and a detailed text book I have tried to write concisely to the point and in plain language Many worked examples and exercises are included A student who works through this primer will have the tools to go to the next level in applying differential equa tions to problems in engineering science and applied mathematics It can give some instructors who want more concise coverage an alternative to existing texts The numerical solution of differential equations is a central activity in sci ence and engineering and it is absolutely necessary to teach students some aspects of scienti c computation as early as possible I tried to build in ex ibility regarding a computer environment The text allows students to use a calculator or a computer algebra system to solve some problems numerically and symbolically and templates of MATLAB and Maple programs and com mands are given in an appendix The instructor can include as much of this or as little of this as he or she desires For many years I have taught this material to students who have had a standard threesemester calculus sequencer It was well received by those who xii Preface appreciated having a small de nitive parcel of material to learn Moreover this text gives students the opportunity to start reading mathematics at a slightly higher level than experienced in precalculus and calculus Therefore the book can be a bridge in their progress to study more advanced material at the junior enior level where books leave a lot to the reader and are not packaged in elementary formats Chapters 1 2 3 5 and 6 should be covered in orderi They provide a route to geometric understanding the phase plane and the qualitative ideas that are important in differential equations Included are the usual treatments of separable and linear rstorder equations along with secondorder linear ho mogeneous and nonhomogeneous equations There are many applications to ecology physics engineering and other areas These topics will give students key skills in the subject Chapter 4 on Laplace transforms can be covered at any time after Chapter 3 or even omittedi Always an issue in teaching differ ential equations is how much linear algebra to cover In two extended sections in Chapter 5 we introduce a moderate amount of matrix theory including solv ing linear systems determinants and the eigenvalue problemi In spite of the book s brevity it still contains slightly more material than can be comfortably covered in a single threehour semester coursei Generally I assign most of the exercises397 hints and solutions for selected problems are given in Appendix Dr I welcome suggestions comments and corrections Contact information is on my Web site http wwwmathun1edu dlogan where additional items may be found I would like to thank John Polking at Rice University for permitting me to use his MATLAB program pplane7 to draw some of the phase plane diagrams and Mark Spencer at Springer for his enthusiastic support of this project Fi nally I would like to thank Tess for her continual encouragement and support for my work David Logan Lincoln Nebraska To the Student What is a course in differential equations about Here are some informal preparatory remarks to give you some sense of the subject before we take it up seriouslyi You are familiar with algebra problems and solving algebraic equations For example the solutions to the quadratic equation 12 7 z 0 are easily found to be I 0 and z l which are numbers A differential equation sometimes abbreviated DE is another type of equation where the unknown is not a number but a function We will call it ut and think of it as a function of time A DE also contains derivatives of the unknown function which are also not known So a DE is an equation that relates an unknown function to some of its derivatives A simple example of a DE is ut ut7 where u t denotes the derivative of We ask what function ut solves this equation That is what function ut has a derivative that is equal to itself From calculus you know that one such function is ut e the exponential function We say this function is a solution of the DE or it solves the DE Is it the only one If you try ut Ce where C is any constant whatsoever you will also nd it is a solution So differential equations have lots of solutions fortunately we will see they are quite similar and the fact that there are many allows some exibility in imposing other desired conditions This DE was very simple and we could guess the answer from our calcu lus knowledgei But unfortunately or fortunatelyl differential equations are usually more complicated Consider for example the DE u t 21 t 2ut 0i xiv To the Student This equation involves the unknown function and both its rst and second derivatives We seek a function for which its second derivative plus twice its rst derivative plus twice the function itself is zero Now can you quickly guess a function ut that solves this equation It is not likely An answer is ut 6quot cos ti And ut 6quot sint works as well Let s check this last one by using the product rule and calculating its derivatives ut 6quot sin t ut 6quot cos t 7 6quot sin t u t 76quot sint 7 26quot cos t Git sinti Then u t 21 t 2ut 76quot sint 7 26quot cos t 6quot sint 26quot cost 7 6quot sin t 26quot sint 0 So it works The function ut 6quot sint solves the equation u t 2u t 2ut 0 In fact ut A64 sint B6quot cost is a solution regardless of the values of the constants A and B So again differential equations have lots of solutions Partly the subject of differential equations is about developing methods for nding solutions Why differential equations Why are they so important to deserve a course of study Well differential equations arise naturally as models in areas of sci ence engineering economics and lots of other subjects Physical systems bi ological systems economic systems7all these are marked by change Differen tial equations model realworld systems by describing how they change The unknown function ut could be the current in an electrical circuit the concen tration of a chemical undergoing reaction the population of an animal species in an ecosystem or the demand for a commodity in a microeconomy Differ ential equations are laws that dictate change and the unknown ut for which we solve describes exactly how the changes occur In fact much of the reason that the calculus was developed by Isaac Newton was to describe motion and to solve differential equations To the Student xv For example suppose a particle of mass m moves along a line with constant velocity VI Suddenly say at time t 07 there is imposed an external resistive force F on the particle that is proportional to its velocity v vt for times t gt 0 Notice that the particle will slow down and its velocity will change From this information can we predict the velocity vt of the particle at any time t gt 0 We learned in calculus that Newton s second law of motion states that the mass of the particle times its acceleration equals the force or ma Fl We also learned that the derivative of velocity is acceleration so a vt Therefore if we write the force as F ikvt where k is the proportionality constant and the minus sign indicates the force opposes the motion then mvt ikvti This is a differential equation for the unknown velocity vti If we can nd a function vt that works in the equation and also satis es 110 V0 then we will have determined the velocity of the particle Can you guess a solution After a little practice in Chapter 1 we will be able to solve the equation and nd that the velocity decays exponentially397 it is given by w Vedam t2 0 Let s check that it works k mvt mVI aimm ikl e ktm ikvti Moreover 110 Vol So it does check The differential equation itself is a model that governs the dynamics of the particle We set it up using Newton s second law and it contains the unknown function vt along with its derivative 1 The solution vt dictates how the system evolvesi In this text we study differential equations and their applications We ad dress two principal questions 1 How do we nd an appropriate DE to model a physical problem 2 How do we understand or solve the DE after we obtain it We learn modeling by examining models that others have studied such as Newton s second law and we try to create some of our own through exercisesi We gain understanding and learn solution techniques by practice Now we are ready Read the text carefully with pencil and paper in hand and work through all the examples Make a commitment to solve most of the exercises You will be rewarded with a knowledge of one of the monuments of mathematics and science 1 Differential Equations and Models In science engineering economics and in most areas where there is a quanti tative component we are greatly interested in describing how systems evolve in time that is in describing a system s dynamics In the simplest one dimensional case the state of a system at any time t is denoted by a function which we generically write as u We think of the dependent variable u as the state variable of a system that is varying with time t which is the independent variable Thus knowing u is tantamount to knowing what state the system is in at time t For example ut could be the population of an animal species in an ecosystem the concentration of a chemical substance in the blood the number of infected individuals in a flu epidemic the current in an electrical circuit the speed of a spacecraft the mass of a decaying isotope or the monthly sales of an advertised itemi Knowledge of ut for a given sys tem tells us exactly how the state of the system is changing in time Figure 11 shows a time series plot of a generic state function We always use the variable u for a generic state397 but if the state is population then we may use p or N if the state is voltage we may use Vi For mechanical systems we often use I zt for the position One way to obtain the state ut for a given system is to take measurements at different times and t the data to obtain a nice formula for Or we might read ut off an oscilloscope or some other gauge or monitor Such curves or formulas may tell us how a system behaves in time but they do not give us insight into why a system behaves in the way we observe Therefore we try to formulate explanatory models that underpin the understanding we seek Often these models are dynamic equations that relate the state ut to its rates of 2 1 Differential Equations and Models state time Figure 11 Time series plot of a generic state function u ut for a system change as expressed by its derivatives ut7 u t7 m and so on Such equations are called differential equations and many laws of nature take the form of such equations For example Newton s second law for the motion for a mass acted upon by external forces can be expressed as a differential equation for the unknown position I zt of the mass In summary a differential equation is an equation that describes how a state ut changes A common strategy in science engineering economics etc is to formulate a basic principle in terms of a differential equation for an unknown state that characterizes a system and then solve the equation to determine the state thereby determining how the system evolves in time 11 Differential Equations 111 Equations and Solutions A differential equation abbreviated DE is simply an equation for an un known state function u ut that connects the state function and some of its derivatives Several notations are used for the derivative including du u7E7u7M The overdot notation is common in physics and engineering397 mostly we use the simple mime notation The reader should be familiar with the de nition of the derivative ut h 7 ut h 7 mm 11 Differential Equations 3 For small h the difference quotient on the right side is often taken as an approximation for the derivative Similarly the second derivative is denoted by 2 du u 7W7u7m and so forth397 the nth derivative is denoted by um The rst derivative of a quantity is the rate of change of the quantity measuring how fast the quantity is changing and the second derivative measures how fast the rate is changing For example if the state of a mechanical system is position then its rst derivative is velocity and its second derivative is acceleration or the rate of change of velocity Differential equations are equations that relate states to their rates of change and many natural laws are expressed in this manner The order of the highest derivative that occurs in the DE is called the order of the equation Example 11 Three examples of differential equations are t9 sint 07 LqHRqI q sinth7 P 0 7 T100 fl The rst equation models the angular de ections t9 9t of a pendulum of length l39 the second models the charge 4 qt on a capacitor in an electri cal circuit containing an inductor resistor and a capacitor where the current is driven by a sinusoidal electromotive force operating at frequency w in the last equation called the logistics equation the state function p pt rep resents the population of an animal species in a closed ecosystem397 39r is the population growth rate and K represents the capacity of the ecosystem to sup port the population The unspeci ed constants in the various equations 1 L7 R C7 w 7 and K are called parameters and they can take any value we choose Most differential equations that model physical processes contain such parameters The constant 9 in the pendulum equation is a xed parameter representing the acceleration of gravity on earth In mks units 9 98 meters per secondsquaredi The unknown in each equation 9t qt and pt is the state function The rst two equations are secondorder and the third equation is rstorderi Note that all the state variables in all these equations depend on time t Because time dependence is understood we often save space and drop 4 1 Differential Equations and Models that dependence when writing differential equations So for example in the rst equation 9 means 9t and 9 means t9 ti In this chapter we focus on rstorder differential equations and their ori gins We write a generic rstorder equation for an unknown state u ut in the form 1 ft 11 When we have solved for the derivative we say the equation is in normal form There are several words we use to classify DEs and the reader should learn theml If f does not depend explicitly on t ie the DE has the form 1 then we call the DE autonomous Otherwise it is nonautonomous For example the equation 1 73112 2 is autonomous but u 73112 cost is nonautonomous If f is a linear function in the variable u then we say 11 is linear else it is nonlinear For example the equation 1 73112 2 is nonlinear because ft u 73112 l 2 is a quadratic function of u not a linear one The general form of a rstorder linear equation is u plttgtu m where p and q are known functions Note that in a linear equation both u and 1 occur alone and to the rst power but the time variable t can occur in any manner Linear equations occur often in theory and applications and their study forms a signi cant part of the subject of differential equations A function u ut is a solution1 of the DE 11 on an interval I z a lt t lt b if it is differentiable on I and when substituted into the equation it satis es the equation identically for all t 6 I39 that is W mum te 1 Therefore a function is a solution if when substituted into the equation every term cancels out In a differential equation the solution is an unknown state function to be found For example in u 7ue the unknown is a function u ut39 we ask what function ut has the property that its derivative is the same as the negative of the function plus 6quot Example 12 This example illustrates what we might expect from a rstorder linear DEl Consider the DE 1 7u equot 1 We are overburdening the notation by using the same symbol u to denote both a Variable and a function It would be more precise to write u 40t is a solution77 but we choose to stick to the common use and abuse of a single letter 11 Differential Equations 5 The state function ut te t is a solution to the DE on the interval I 700 lt t lt 00 Later we learn how to find this solution In fact for any constant C the function ut t Cequot is a solution We can verify this by direct substitution of u and u into the DE397 using the product rule for differentiation u t C7equot equot 7u 6 Therefore ut satis es the DE regardless of the value of C We say that this expression ut t Cequot represents a oneparameter family of solutions one solution for each value of C This example illustrates the usual state of affairs for any rstorder linear DEithere is a oneparameter family of solutions depending upon an arbitrary constant C This family of solutions is called a general solution The fact that there are many solutions to rst order differential equations turns out to be fortunate because we can adjust the constant C to obtain a speci c solution that satis es other conditions that might apply in a physical problem eg a requirement that the system be in some known state at time t 0 For example if we require u0 1 then C 1 and we obtain a particular solution ut t1equoti Figure 12 shows a plot of the oneparameter family of solutions for several values of Cl Here we are using the word parameter in a different way from that in Example lil397 there the word parameter refers to a physical number in the equation itself that is xed yet arbitrary like resistance in a circuit 2 02 Ag 072 i as 05 1 15 2 25 tiu on the interval Figure 12 Time series plots of several solutions to u e 71 g t 3 The solution curves or the oneparameter family of solutions are ut t Cequot7 where C is an arbitrary constant here taking several values between 72 and 2 6 1 Differential Equations and Models An initial value problem abbreviated IVP for a rstorder DE is the problem of nding a solution u ut to 1 1 that satis es an initial condi tion ut0 uo7 Where to is some xed value of time and uo is a xed state We Write the IV concisely as u m u lVP 7 7 1 2 ut0 uoi The initial condition usually picks out a speci c value of the arbitrary constant C that appears in the general solution of the equation So it selects one of the many possible states that satisfy the differential equation The accompanying graph gure 13 depicts a solution to an IV Figure 13 Solution to an initial value problemi The fundamental questions are a is there a solution curve passing through the given point b is the curve the only one and c What is the interval 15 on Which the solution exists Geometrically solving an initial value problem means to nd a solution to the DE that passes through a speci ed point to7 uo in the plane Referring to Example 12 the IV 1 711 equot7 u0 1 has solution ut t Dequot Which is valid for all times t The solution curve passes through the point 07 17 corresponding to the initial condition 110 1 Again the initial condition selects one of the many solutions of the DE397 it xes the value of the arbitrary constant C There are many interesting mathematical questions about initial value prob lems 11 Differential Equations 7 1 Existence Given an initial value problem is there a solution This is the question of existence Note that there may be a solution even if we cannot nd a formula for it 2 Uniqueness If there is a solution is the solution unique That is is it the only solution This is the question of uniqueness 3i Interval of existence For which times t does the solution to the initial value problem exist Obtaining resolution of these theoretical issues is an interesting and worth while endeavor and it is the subject of advanced courses and books in differ ential equations In this text we only brie y discuss these matters The next three examples illustrate why these are reasonable questions Example 13 Consider the initial value problem 11 ut 7 37 u1 2 This problem has no solution because the derivative of u is not de ned in an interval containing the initial time t 1 There cannot be a solution curve passing through the point 12 Example 14 Consider the initial value problem 11 21112 u0 0 The reader should verify that both ut 0 and ut t2 are solutions to this initial value problem on t gt 0 Thus it does not have a unique solution More than one state evolves from the initial state Example 15 Consider the two similar initial value problems u 17112 11007 u 1 112 u0 0 The rst has solution 6271 7175 m7 8 1 Differential Equations and Models which exists for every value of ti Yet the second has solution ut tan t7 which exists only on the interval 7 lt t lt So the solution to the rst initial value problem is de ned for all times but the solution to the second blows up in nite time These two problems are quite similar yet the times for which their solutions exist are quite different The following theorem which is proved in advanced books provides partial answers to the questions raised above The theorem basically states that if the right side ft7 u of the differential equation is nice enough then there is a unique solution in a neighborhood the initial value Theorem 16 Let the function f be continuous on the open rectangle R z a lt t lt b c lt u lt d in the tu plane and consider the initial value problem ft 7 1 3 W750 um where to7 uo lies in the rectangle Ri Then the IV 13 has a solution u ut on some interval 15 containing to where 15 C a7bl If in addition the partial derivative2 fu t7 u is continuous on R then 13 has a unique solution The interval of existence is the set of time values for which the solution to the initial value problem exists Theorem 16 is called a local existence theorem because it guarantees a solution only in a small neighborhood of the initial time to the theorem does not state how large the interval of existence is Observe that the rectangle R mentioned in the theorem is open and hence the initial point cannot lie on its boundary In Example 15 both right sides of the equations ftu 1 7 u2 and ft7 u 1 u2 are continuous in the plane and their partial derivatives fu 7211 and fu 2o are continuous in in the plane So the initial value problem for each would have a unique solution regardless of the initial condition In addition to theoretical questions there are central issues from the view point of modeling and applications these are the questions we mentioned in the To the Student section 1 How do we determine a differential equation that models or governs a given physical observation or phenomenon 2 We use subscripts to denote partial derivatives7 and s0 fu 11 Differential Equations 9 2 How do we nd a solution either analytically approximately graphically or numerically u ut of a differential equation The rst question is addressed throughout this book by formulating model equations for systems in particles dynamics circuit theory biology and in other areas We learn some basic principles that sharpen our ability to invent explanatory models given by differential equations The second question is one of developing methods and our approach is to illustrate some standard ana lytic techniques that have become part of the subject By an analytic method we mean manipulations that lead to a formula for the solution39 such formulas are called analytic solutions or Closedform solutions For most realworld problems it is dif cult or impossible to obtain an analytic solution By a nu merical solution we mean an approximate solution that is obtained by some computer algorithm a numerical solution can be represented by a data set ta ble of numbers or by a graph In real physical problems numerical methods are the ones most often used Approximate solutions can be formulas that approximate the actual solution eg a polynomial formula or they can be numerical solutions Almost always we are interested in obtaining a graphical representation of the solution Often we apply qualitative methods These are methods designed to obtain important information from the DE without actually solving it either numerically or analytically For a simple example consider the DE u u2 it Because u gt 0 we know that all solution curves are increasing Or for the DE u u2 7 t2 we know solution curves have a hor izontal tangent as they cross the straight lines u it Quantitative methods emphasize understanding the underlying model recognizing properties of the DE interpreting the various terms and using graphical properties to our ben e t in interpreting the equation and plotting the solutions39 often these aspects are more important than actually learning specialized methods for obtaining a solution formulai Many of the methods both analytic and numerical can be performed easily on computer algebra systems such as Maple Mathematica or MATLAB and some can be performed on advanced calculators that have a builtin computer algebra systemi Although we often use a computer algebra system to our ad vantage especially to perform tedious calculations our goal is to understand concepts and develop techniquei Appendix B contains information on using MATLAB and Maple 112 Geometrical Interpretation What does a differential equation u ftu tell us geometrically At each point tu of the tuplane the value of ft u is the slope u of the solution 10 1 Differential Equations and Models curve u ut that goes through that point This is because ut ft7 1W This fact suggests a simple graphical method for constructing approximate solution curves for a differential equation Through each point of a selected set of points t7 u in some rectangular region or window of the tu plane we draw a short line segment with slope ft7 The collection of all these line segments or minitangents form the direction eld or slope eld for the equation We may then sketch solution curves that t this direction eld397 the curves must have the property that at each point the tangent line has the same slope as the slope of the direction eld For example the slope eld for the differential equation u 711 2t is de ned by the right side of the differential equation ft7 u 711 2t The slope eld at the point 24 is 1623 732 4 5 This means the solution curve that passes through the point 2 4 has slope 5 Because it is tedious to calculate several minitangents simple programs have been developed for calculators and computer algebra systems that perform this task automatically for us Figure 14 shows a slope eld and several solution curves that have been t into the eld Figure 14 The slope eld in the window 72 g t g 4 74 g u g 8 with several approximate solution curves for the DE 1 7u 2t Notice that a problem in differential equations is just opposite of that in differential calculusi ln calculus we know the function curve and are asked to 11 Differential Equations 11 nd the derivative slope397 in differential equations we know the slopes and try to nd the state function that ts themi Also observe that the simplicity of autonomous equations no time t depen dence on the right side u u shows itself in the slope eld In this case the slope eld is independent of time so on each horizontal line in the tu plane where u has the same value the slope eld is the same For example the DE u 3u5 7 u is autonomous and along the horizontal line u 2 the slope eld has value 18 This means solution curves cross the line u 2 with a relatively steep slope u 18 EXERCISES H Verify that the two differential equations in Example 15 have solutions as stated 2 From the set of commonly encountered functions guess a nonzero solution u ut to the DE 11 u2i 9 Show that ut lnt C is a oneparameter family of solutions of the E u 6 where C is an arbitrary constant Plot several members of this family Find and plot a particular solution that satis es the initial condition u0 0 F Find a solution u ut of u 2u t2 4t 7 in the form of a quadratic function of ti 5 Find values of m such that u tm is a solution to Ztu u 533 Plot the oneparameter family of curves ut t2 7 C63 and nd a differential equation whose solution is this family 7 Show that the oneparameter family of straight lines u Ct is a solution to the differential equation tu 7 u 0 for any value of the constant Ci 5 0 Classify the rstorder equations as linear or nonlinear autonomous or nonautonomousi a 1M 2t3u 7 6 b cos tu 7 2u sinu 0 c 11 d 71 7 Su 0i 1 Differential Equations and Models to H O H H H N H 9 HH we H on H H H 00 i Explain Example 14 in the context of Theorem 16 In particular explain why the theorem does not apply to this initial value problemi Which hy pothesis fails i Verify that the initial value problem 1 u0 0 has in nitely many solutions of the form u l film where a gt 0 Sketch these solutions for different values of a What hypoth esis fails in Theorem 16 tga tgta7 i Consider the linear differential equation 1 ptu 10 ls it true that the sum of two solutions is again a solution ls a constant times a solution again a solution Answer these same questions if qt 0 Show that if ul is a solution to u ptu and ug is a solution to u ptu 40 then ul ug is a solution to u ptu 475 By hand sketch the slope eld for the DE 1 ul 7 in the window 0 g t g 87 0 g u g 8 at integer points What is the value of the slope eld along the lines u 0 and u 4 Show that ut 0 and ut 4 are constant solutions to the DE On your slope eld plot draw in several solution curvesi Using a software package sketch the slope eld in the window 74 g t g 4 72 g u g 2 for the equation 1 l 7 u2 and draw several approximate solution curvesi Lines and curves in the tu plane where the slope eld is zero are called nullclines For the given DE nd the nullclinesi Graph the locus of points where the slope eld is equal to 73 Repeat Exercise 13 for the equation 1 t7 u2i i In the tu plane plot the nullclines of the differential equation 1 2112 u 7 4V i Using concavity show that the secondorder DE u 7 u 0 cannot have a solution other than the u 0 solution that takes the value zero more than once Hint construct a contradiction argument7if it takes the value zero twice it must have a negative minimum or positive maximuml i For any solution u ut of the DE u 7 u 0 show that ii2 7 u2 C where C is a constant Plot this one parameterfamily of curves on a 1111 set of axes i Show that if ul and ug are both solutions to the DE 1 ptu 07 then ulug is constant 12 Pure Time Equations 13 19 Show that the linear initial value problem 2 7 1 11 ut 110l7 has a continuously differentiable solution ie a solution whose rst deriva tive is continuous given by at2l tlt0 ut7bt21 tgt0 for any constants a and 12 Yet there is no solution if u0 y 1 Do these facts contradict Theorem 16 12 Pure Time Equations In this section we solve the simplest type of differential equation First we need to recall the fundamental theorem of calculus which is basic and used regularly in differential equations For reference we state the two standard forms of the theoremi They show that differentiation and integration are inverse processesi Fundamental Theorem of Calculus I If u is a differentiable function the integral of its derivative is b a 1404 ub 7 my Fundamental Theorem of Calculus II If g is a continuous function the derivative of an integral with variable upper limit is g gltsgtdsglttgt where the lower limit a is any number This last expression states that the function fgsds is an antiderivative of 9 ie a function whose derivative is 9 Notice that fgsds C is also an antiderivative for any value of Or The simplest differential equation is one of the form u 9t7 14 where the right side of the differential equation is a given known function gti This equation is called a pure time equation Thus we seek a state function u whose derivative is gti The fundamental theorem of calculus H u must be 14 1 Differential Equations and Models an antiderivative of g We can write this fact as ut fgsds C or using the inde nite integral notation as ut gtdtC 15 where C is an arbitrary constant called the constant of integration Recall that antiderivatives of a function differ by an additive constanti Thus all so lutions of 14 are given by 15 and 15 is called the general solution The fact that 15 solves 1 4 follows from the fundamental theorem of calculus Hi Example 17 Find the general solution to the differential equation u t2 7 1 Because the right side depends only on t the solution u is an antiderivative of the right side or 1 ut gt3 7 tC where C is an arbitrary constant This is the general solution and it graphs as a family of cubic curves in the tu plane one curve for each value of Or A particular antiderivative or solution can be determined by imposing an initial condition that picks out a speci c value of the constant C and hence a speci c curve For example if ul 2 then l3 7 l C 2 giving C The solution to the initial value problem is then ut t3 7 t Example 18 For equations of the form u gt we can take two successive antiderivatives to nd the general solution The following sequence of calculations shows howl Consider the DE u t2i Then 1 u t22tC1 u t3t201t02i Here C1 and C2 are two arbitrary constantsi For secondorder equations we always expect two arbitrary constants or a twoparameter family of solutions It takes two auxiliary conditions to determine the arbitrary constantsi In this example if u0 l and if u 0 0 then cl l and Cg l and we obtain the particular solution u g t2 1 12 Pure Time Equations 15 Example 19 The autonomous equation u u cannot be solved by direct integration because the right side is not a known function of t it depends on u which is the unknown in the problemi Equations with the unknown u on the right side are not pure time equationsi Often it is not possible to nd a simple expression for the antiderivative7 sint z analytic expressions for their antiderivativesi In these cases we must represent 2 and equot have no s1mple or inde nite integral For example the functions the antiderivative of g as ut tgsd8 C with a variable upper limit Here a is any xed value of time and C is an arbitrary constant We have used the dummy variable of integration 8 to avoid confusion with the upper limit of integration the independent time variable t It is really not advisable to write ut fgtdti Example 110 Solve the initial value problem u 72 u0 2 The right side of the differential equation has no simple expression for its an tiderivativei Therefore we write the antiderivative in the form L ut eiszds Ci 0 The common strategy is to take the lower limit of integration to be the initial value of t here zeroi Then u0 2 implies C 2 and we obtain the solution to the initial value problem in the form of an integral 2 ut eisgds2i 1 6 0 If we had written the solution of the differential equation as ut equot2dt C7 1 Differential Equations and Models in terms of an inde nite integral then there would be no way to use the initial condition to evaluate the constant of integration or evaluate the solution at a particular value of t We emphasize that integrals With avariable upper limit of integration de ne a function Referring to Example 110 we can de ne the special function erf called the error function by m 2 f 424 er e 5 w 0 Thefactorl W in front ofthe integral normalizes the function to force erfltxgt 1 Up to this constant multiple the erf function gives the area under a bell shaped curve exp7s2 from 0 to t In terms of this special function the solution 16 can be Written ut 2 g erfti The erf function Which is plotted in gure 15 is an important function in probability and statistics and in diffusion processes lts values are tabulated in computer algebra systems and mathematical handbooks I Figure 15 Graph of the erf functioni Functions de ned by integrals are common in the applied sciences and are equally important as functions de ned by simple algebraic formulas To the 12 Pure Time Equations 17 point the reader should recall that the natural logarithm can be de ned by the integral t1 lnt ds tgt0i 13 One important viewpoint is that differential equations often de ne special func tions For example the initial value problem u ult1gt 07 can be used to de ne the natural logarithm function ln ti Other special func tions of mathematical physics and engineering for example Bessel functions Legendre polynomials and so on are usually de ned as solutions to differen tial equations By solving the differential equation numerically we can obtain values of the special functions more ef ciently than looking those values up in tabulated formi We end this section with the observation that one can nd solution formu las using computer algebra systems like Maple MATLAB Mathematica etc and calculators equipped with computer algebra systemsi Computer algebra systems do symbolic computationi Below we show the basic syntax in Maple Mathematica and on a Tl SQ that returns the general solution to a differential equation or the solution to an initial value problemi MATLAB has a special addon symbolic package that has similar commandsi Our interest in this text is to use MATLAB for scienti c computation rather than symbolic calculationi Additional information on computing environments is in Appendix B The general solution of the rstorder differential equation u ftu can be obtained as follows deSolve u f t u t u Tl89 dsolvediffut tft ut ut Maple DSolve u t f t ut ut t Mathematica To solve the initial value problem u ftu7 ua b the syntax is deSolveu ftu and uab t u Tl89 dsolvediffutt ftut uab ut Maple DSolve u t f tut uEa b ut t Mathematica EXERCISES 1 Using antiderivatives nd the general solution to the pure time equation u tcost2 and then nd the particular solution satisfying the initial condition uO 1 Graph the particular solution on the interval 75 5 18 1 Differential Equations and Models 2 Solve the initial value problem 1 a271 ul 4i 3 Find a function ut that satis es the initial value problem u 73W ul lu 1 2i 4 Find all state functions that solve the differential equation 1 te mi 5 Find the solution to the initial value problem 1 Lit ul 07 in terms of an integral Graph the solution on the interval 14 by using numerical integration to calculate values of the integral 6 The differential equation 1 3u 6quot can be converted into a pure time equation for a new dependent variable y using the transformation u yegti Find the pure time equation for y solve it and then determine the general solution u of the original equation 7 Generalize the method of Exercise 6 by devising an algorithm to solve u au 4t where a is any constant and q is a given function In fact show that t ut 06 quot eat e asqsdsi 0 Using the fundamental theorem of calculus verify that this function does solve u au qti 8 Use the chain rule and the fundamental theorem of calculus to compute the derivative of erfsin t 9 Exact equationsi Consider a differential equation written in the non normal form ftu gtuu 0 If there is a function h ht7 u for which ht ht huu 07 or by the chain rule just h u 0 Such equations are f and hu y then the differential equation becomes called exact equations because the left side is exactly a total derivative of the function h ht7 The general solution to the equation is therefore given implicitly by htu C where C is an arbitrary constant a Show that ftu gtuu 0 is exact if and only if fu g b Use part a to check if the following equations are exact If the equa tion is exact7 nd the general solution by solving ht f and hu g for h you may want to review the method of nding potential func tions associated with a conservative force eld from your multivariable calculus course it u3 StuQu 0 iii t3 112 lntu 0i isinuiusint 111 u tcosucost H t Mathematical Models 19 H O H H H N H CA3 i An integral equation is an equation where the unknown ut appears under an integral signi Use the fundamental theorem of calculus to show that the integral equation 2 ut e pquot5usds A 1014 constants 0 can be transformed into an initial value problem for i Show that the integral equation t ut 6721 susds 0 can be transformed into an initial value problem for Show by integration that the initial value problem 13 can be trans formed into the integral equation ut uo 0 fsusdsl i From the de nition of the derivative a difference quotient approximation to the rst derivative is u t E Use Taylor s theorem to show that an approximation for the second derivative is N ut h 7 2ut ut 7 h u t Recall that Taylor s expansion for a function u about the point t with increment h is ut h ut uth git ah a a t Use this and and a similar formula for ut 7 13 Mathematical Models By a mathematical model we mean an equation or set of equations that describes some physical problem or phenomenon that has its origin in science engineering or some other area Here we are interested in differential equation models By mathematical modeling we mean the process by which we obtain and analyze the model This process includes introducing the important and relevant quantities or variables involved in the model making modelspeci c as sumptions about those quantities solving the model equations by some method 20 1 Differential Equations and Models and then comparing the solutions to real data and interpreting the results Of ten the solution method involves computer simulation This comparison may lead to revision and re nement until we are satis ed that the model accurately describes the physical situation and is predictive of other similar observations Therefore the subject of mathematical modeling involves physical intuition formulation of equations solution methods and analysis Overall in mathe matical modeling the overarching objective is to make sense of the world as we observe it often by inventing caricatures of reality Scienti c exactness is t t quotM Model A quot 39 depend strongly on the assumptions and changing the assumptions changes the modeli u i sometimes sacri ced for If some assumptions are less critical than others we say the model is robust to those assumptions The best strategy to learn modeling is to begin with simple examples and then graduate to more dif cult ones The reader is already familiar with some models In an elementary science or calculus course we learn that Newton s second law force equals mass times acceleration governs mechanical systems such as falling bodies397 Newton s inversesquare law of gravitation describes the motion of the planets 7 Ohm s law in circuit theory dictates the voltage drop across a resistor in terms of the current397 or the law of mass action in chemistry describes how fast chemical reactions occur In this course we learn new models based on differential equations The importance of differential equations as a subject matter lies in the fact that differential equations describe many physical phenomena and laws in many areas of application In this section we introduce some simple problems and develop differential equations that model the physical processes involved The rst step in modeling is to select the relevant variables independent and dependent and parameters that describe the probleml Physical quantities have dimensions such as time distance degrees and so on or corresponding units such as seconds meters and degrees Celsius The equations we write down as models must be dimensionally correctl Apples cannot equal orangesi Verifying that each term in our model has the same dimensions is the rst task in obtaining a correct equationi Also checking dimensions can often give us insight into what a term in the model might be We always should be aware of the dimensions of the quantities both variables and parameters in a model and we should always try to identify the physical meaning of the terms in the equations we obtain All of these comments about modeling are perhaps best summarized in a quote attributed to the famous psychologist Carl Jung Science is the art of creating suitable illusions which the fool believes or argues against but the wise man enjoys their beauty and ingenuity without being blind to the fact they are human veils and curtains concealing the abysmal darkness of the unknowablel 13 Mathematical Models 21 When one begins to feel too con dent in the correctness of the model he or she should recall this quote 131 Particle Dynamics In the late 16th and early 17th centuries scientists were beginning to quanti tatively understand the basic laws of motion Galileo for example rolled balls down inclined planes and dropped them from different heights in an effort to understand dynamical laws But it was Isaac Newton in the mid1600s who developed calculus and the theory of gravitation who nally wrote down a basic law of motion known now as Newton s second law that is in reality a differential equation for the state of the dynamical systemi For a particle of mass m moving along a straight line under the in uence of a speci ed external force F the law dictates that mass times acceleration equals the force on the particle or m1 Ft7 171 Newton s second law This is a secondorder differential equation for the unknown location or position I zt of the particle The force F may depend on time t position I 10 or velocity 1 0 This DE is called the equation of motion or the dynamical equation for the systemi For secondorder differential equations we impose two initial conditions 10 IO and zO vo which x the initial position and initial velocity of the particle respectively We expect that if the initial position and velocity are known then the equation of motion should determine the state for all times t gt 0 Example 111 Suppose a particle of mass m is falling downward through a viscous uid and the uid exerts a resistive force on the particle proportional to the square of its velocity We measure positive distance downward from the top of the uid surface There are two forces on the particle gravity and uid resistance The gravitational force is my and is positive because it tends to move the mass in a positive downward direction397 the resistive force is fax2 and it is negative because it opposes positive downward motion The net force is then F mgi aza and the equation of motion is m1 mgi az22i This second order equation can immediately be reformulated as a rstorder differential equation for the velocity v 1quot Clearly vgi vi m 22 1 Differential Equations and Models v equiliblium solution vt VT ll t Figure 16 Generic solution curves or time series plots for the model 1 g 7 am12l For 1 lt UT the solution curves are increasing because 1 gt 0397 for 1 gt UT the solution curves are decreasing because 1 lt 0 All the solution curves approach the constant terminal velocity solution vt up If we impose an initial velocity 110 v0 then this equation and the initial condition gives an initial value problem for 1 11t Without solving the DE we can obtain important qualitative information from the DE itselfl Over a long time if the uid were deep we would observe that the falling mass would approach a constant terminal velocity le Physically the terminal velocity occurs when the two forces the gravitational force and resistive force balancei Thus 0 g 7 avgnm or my UT i a By direct substitution we note that 1t UT is a constant solution of the differential equation with initial condition 110 up We call such a constant solution an equilibrium or steadystate solution It is clear that regard less of the initial velocity the system approaches this equilibrium state This supposition is supported by the observation that 1 gt 0 when 1 lt UT and 1 lt 0 when 1 gt up Figure 16 shows what we expect illustrating several generic solution curves time series plots for different initial velocitiesl To nd the position zt of the object we would integrate the velocity vt once it is determined397 that is zt fgvsdsl Example 112 A ball of mass m is tossed upward from a building of height h with initial velocity vol If we ignore air resistance then the only force is that due to grav 13 Mathematical Models 23 ity having magnitude mg directed downwardi Taking the positive direction upward with z 0 at the ground the model that governs the motion ie the height I zt of the ball is the initial value problem m1 7mg 10 h7 zO vol Note that the force is negative because the positive direction is upwardl Because the right side is a known function a constant in this case the differential equation is a pure time equation and can be solved directly by integration antiderivativesi 1f z t 9 ie the second derivative is the constant g then the rst derivative must be z t 7975 517 where cl is some constant the constant of integration We can evaluate cl using the initial condition zO vol We have zO 9 X 0 cl 1107 giving cl vol Therefore at any time the velocity is given by zt 7975 vol Repeating we take another antiderivativei Then 1 2 zt 75975 vot 52 where Cg is some constant Using 10 h we nd that Cg hi Therefore the height of the ball at any time t is given by the familiar physics formula 1 1t 7597 vot hi Example 113 Imagine a mass m lying on a table and connected to a spring which is in turn attached to a rigid wall gure 17 At time t 0 we displace the mass a positive distance 10 to the right of equilibrium and then release it If we ignore friction on the table then the mass executes simple harmonic motion397 that is it oscillates back and forth at a xed frequency To set up a model for the motion we follow the doctrine of mechanics and write down Newton s second law of motion m1 F7 where the state function I zt is the position of the mass at time t we take I 0 to be the equilibrium position and z gt 0 to the right and F is the external force All that is required is to impose the form of the force Experiments con rm that if the displacement is not too large which we assume then the force exerted by the spring is proportional to its displacement from equilibriumi That is F ikzi 1 7 24 1 Differential Equations and Models The minus sign appears because the force opposes positive motion The pro portionality constant k having dimensions of force per unit distance is called the spring constant or stiffness of the spring and equation 17 is called Hooke s laWi Not every spring behaves in this manner but Hooke s law is used as a model for some springs397 it is an example of what in engineering is called a constitutive relation It is an empirical result rather than a law of nature To give a little more justi cation for Hooke s law suppose the force F depends on the displacement I through F Fz with F0 0 Then by Taylor s theorem Fz F0F 0z F O12 ikzF 012 where we have de ned FO iki So Hooke s law has a general validity if the displacement is small allowing the higherorder terms in the series to be neglected We can measure the stiffness k of a spring by letting it hang from a ceiling without the mass attached397 then attach the mass m and measure the elongation L after it comes to rest The force of gravity mg must balance the restoring force 161 of the spring so k mgLi Therefore assuming a Hookean spring we have the equation of motion m1 7km 18 which is the springmass equation The initial conditions released at time zero at position 10 are 10 Io zO 0 We expect oscillatory motioni If we attempt a solution of 18 of the form zt Acos wt for some frequency w and amplitude A we nd upon substitu tion that w xkm and A 10 Therefore the displacement of the mass is given by zt 10 cos kmti This solution represents an oscillation of amplitude 10 frequency km and period ZWxkmi Example 114 Continuing with Example 113 if there is damping caused for example by friction or submerging the system in a liquid then the springmass equation 13 Mathematical Models 25 equilibrium 2620 I xz i l H I I Figure 17 Springmass oscillatori must be modi ed to account for the damping force The simplest assumption again a constitutive relation is to take the resistive force Fr to be proportional to the velocity of the mass Thus also assuming Hooke s law for the spring force F5 we have the damped springmass equation m1 Fr FS 7017 km The positive constant c is the damping constanti Both forces have negative signs because both oppose positive to the right motioni For this case we expect some sort of oscillatory behavior with the amplitude decreasing during each oscillationi 1n Exercise 1 you will show that solutions representing decaying oscillations do in fact occuri Example 115 For conservative mechanical systems another technique for obtaining the equa tion of motion is to apply the conservation of energy law the kinetic energy plus the potential energy remain constant We illustrate this method by nding the equation governing a frictionless pendulum of length 1 whose bob has mass mi See gure 18 As a state variable we choose the angle 9 that it makes with the verticali As time passes the bob traces out an arc on a circle of radius 139 we let 8 denote the arclength measured from rest 9 0 along the arc By geometry 8 9 As the bob moves its kinetic energy is onehalf its mass times the velocitysquared397 its potential energy is mgh where h is the height above zeropotential energy level taken where the pendulum is at rest There fore ms2 mgl17 cos 9 E7 where E is the constant energy In terms of the angle 9 10932 g17 cos 9 C7 19 26 1 Differential Equations and Models where C Emli The initial conditions are 90 90 and 90 we where 90 and we are the initial angular displacement and angular velocity respectively As it stands the differential equation 19 is firstorder397 the constant C can be determined by evaluating the differential equation at t 0 We get C lwg 91 7 cos 90 By differentiation with respect to t we can write 19 as 9 3an90i 110 z This is a secondorder nonlinear DE in 6t called the pendulum equation It can also be derived directly from Newton s second law by determining the forces which we leave as an exercise Exercise 6 We summarize by stating that for a conservative mechanical system the equation of motion can be found either by determining the energies and applying the conservation of energy law or by nding the forces and using Newton s second law of motion I s16 9 I m fa kl r mg Figure 18 A pendulum consisting of a mass m attached to a rigid weightless rod of length ll The force of gravity is my directed downward The potential energy is mgh where h is the height of the mass above the equilibrium positioni 13 Mathematical Models EXERCISES 1 0 9 F 5 533 I 00 When a mass of 03 kg is placed on a spring hanging from the ceiling it elongates the spring 15 cm What is the stiffness k of the spring l Consider a damped springmass system Whose position zt is governed by the equation m1 701 7 km Show that this equation can have a decayingoscillation solution of the form zt 6 cos wt Hint By substituting into the differential equations show that the decay constant A and frequency M can be determined in terms of the given parameters m c and kl A car of mass m is moving at speed V When it has to brake The brakes apply a constant force F until the car comes to rest How long does it take the car to stop How far does the car go before stopping Now With speci c data compare the time and distance it takes to stop if you are going 30 mph vs 35 mph Take m 1000 kg and F 6500 N Write a short paragraph on recommended speed limits in a residential areas Derive the pendulum equation 110 from the conservation of energy law 19 by differentiation A pendulum of length 05 meters has a bob of mass 01 kg If the pendulum is released from rest at an angle of 15 degrees nd the total energy in the system Derive the pendulum equation 110 by resolving the gravitational force on the bob in the tangential and normal directions along the arc of motion and then applying Newton s second law Note that only the tangential component affects the motion i If the amplitude of the oscillations of a pendulum is small then sint9 is nearly equal to 9 Why and the nonlinear equation 110 is approximated by the linear equation 9 909 0 a Show that the approximate linear equation has a solution of the form 9 t Acoswt for some value of m which also satis es the initial conditions 90 A 90 0 What is the period of the oscillation b A 650 lb Wrecking ball is suspended on a 20 m cord from the top of a crane The ball hanging vertically at rest against the building is pulled back a small distance and then released How soon does it strike the building l An enemy cannon at distance L from a fort can re a cannon ball from the top of a hill at height H above the ground level With a muzzle velocity 1 How high should the wall of the fort be to guarantee that a cannon ball 28 1 Differential Equations and Models will not go over the wall Observe that the enemy can adjust the angle of its shoti Hint lgnoring air resistance the governing equations follow from resolving Newton s second law for the horizontal and vertical components of the force m1 0 and my 7mg 132 Autonomous Di erential Equations In this section we introduce some simple qualitative methods to understand the dynamics of an autonomous differential equation 1 We introduce the methods in the context of population ecology as well as in some other areas in the life sciences Models in biology often have a different character from fundamental laws in the physical sciences such as Newton s second law of motion in mechanics or Maxwell s equations in electrodynamicsl Ecological systems are highly complex and it is often impossible to include every possible factor in a model39 the chore of modeling often comes in knowing what effects are important and what effects are minor Many models in ecology are often not based on physical law but rather on observation experiment and reasoning Ecology is the study of how organisms interact with their environment A fundamental problem in population ecology is to determine what mechanisms operate to regulate animal populations Let p pt denote the population of an animal species at time t For the human population Tl Malthus in the late 1700s proposed the model which states that the per capita growth rate is constant where the constant 7 gt 0 is the growth rate given in dimensions of time li We can regard 39r as the birth rate minus the death rate or 7 b 7 d This per capita law is same as 7 P T177 which says that the growth rate is proportional to the population It is easily veri ed check this that a oneparameter family of solutions is given by 100 GET where C is any constant If there is an initial condition imposed that is 170 p0 then C p0 and we have picked out a particular solution pt p06 of the DE that is the one that satis es the initial condition Therefore the Malthus model predicts exponential population growth gure 19 13 Mathematical Models 29 p Izexprt m t Figure 19 The Malthus model for population growth pt poeml The reader should note a difference between the phrases per capita growth rate and growth rate To say that the per capita growth rate is 2 per time is to say that p p 002 which gives exponential growth397 to say that the growth rate is 2 animals per time is to say 10 002 which forces pt to be of the form pt 0027M I7 K constant which is linear growth In animal populations for fairly obvious reasons we do not expect expo nential growth over long times Environmental factors and competition for re sources limit the population when it gets larger Therefore we might expect the per capita growth rate 7 which is constant in the Malthus model to decrease as the population increases The simplest assumption is a linearly decreasing per capita growth rate where the rate becomes zero at some maximum carry ing capacity Kl See gure lllOl This gives the logistics model of population growth developed by P Verhulst in the 1800s by P p p 7 li or 7 li l llll p lt K p plt K lt gt Clearly we may write this autonomous equation in the form T 2 7 7 l p p K The rst term is a positive growth term which is just the Malthus terml The second term which is quadratic in p decreases the population growth rate and is the competition terml Note that if there were p animals then there would be about p2 encounters among theml So the competition term is proportional to the number of possible encounters which is a reasonable modell Exercise 11 presents an alternate derivation of the logistics model based on food supply 3O 1 Differential Equations and Models rapK VPHPK 0 K 391 390 K 39P Figure 110 Plots of the logistics model of population growth The left plot shows the per capita growth rate vs population and the right plot shows the growth rate vs populationl Both plots give important interpretations of the model For any initial condition 100 p0 we can nd the formula for the solution to the logistics equation 111 You will solve the logistics equation in Exercise 8 But there are qualitative properties of solutions that can be exposed without actually nding the solution Often all we may want are qualitative features of a model First we note that there are two constant solutions to 111 pt 0 and pt K corresponding to no animals extinction and to the number of animals represented by the carrying capacity respectively These constant solutions are found by setting the right side of the equation equal to zero because that forces 10 07 or p constant The constant solutions are called steadystate or equilibrium solutions If the population is between p 0 and p K the right side of 111 is positive giving 10 gt 0397 for these population numbers the population is increasing If the population is larger than the carrying capacity K then the right side of 111 is negative and the population is decreasing These facts can also be observed from the growth rate plot in gure 110 These observations can be represented conveniently on a phase line plot as shown in gure lllll We rst plot the growth rate 17 vs p which in this case is a parabola opening downward The points of intersection on the p axis are the equilibrium solutions 0 and Kl We then indicate by a directional arrow on the p axis those values of p where the solution pt is increasing where p gt 0 or decreasing p lt 0 Thus the arrow points to the right when the graph of the growth rate is above the axis and it points to the left when the graph is below the axis In this context we call the p axis a phase line We can regard the phase line as a onedimensional parametric solution space with the population p pt tracing out points on that line as t increases In the range 0 lt p lt K the arrow points right because p gt 0 So 13 Mathematical Models 31 pt increases in this range For p gt K the arrow points left because 10 lt 0 The population pt decreases in this range These qualitative features can be easily transferred to time series plots gure 1 12 showing pt vsl t for different initial conditions Both the phase line and the time series plots imply that regardless of the initial population if nonzero the population approaches the carrying capac ity Kl This equilibrium population p K is called an attractorl The zero population is also an equilibrium population But near zero we have 10 gt 0 and so the population diverges away from zero We say the equilibrium popu lation p 0 is a repellerl We are considering only positive populations so we ignore the fact that p 0 could be approached on the left side In summary our analysis has determined the complete qualitative behavior of the logistics population modell we growth rate Figure 111 The p axis is the phase line on which arrows indicate an in creasing or decreasing population for certain ranges of p This qualitative method used to analyze the logistics model is applicable to any autonomous equation 11 112 The equilibrium solutions are the constant solutions which are roots of the algebraic equation 0 Thus if uquot is an equilibrium then 0i 32 1 Differential Equations and Models K 0 equilibria p gt 0 increasing 0 t Figure 1 12 Time series plots of solutions to the logistics equation for various initial conditions For 0 lt p lt K the population increases and approaches K whereas for p gt K the population decreases to Kl If 100 K then pt K for all times t gt 0397 this is the equilibrium solution These are the values where the graph of vs u intersects the u axis We always assume the equilibria are isolated397 that is if u is an equilibrium then there is an open interval containing u that contains no other equilibrial Figure 113 shows a generic plot where the equilibria are u a 12 ct In between the equilibria we can observe the values of u for which the population is increasing gt 0 or decreasing lt 0 We can then place arrows on the phase line or the u axis in between the equilibria showing direction of the movement increasing or decreasing as time increases If desired the information from the phase line can be translated into time series plots of ut vs t gure 114 In between the constant equilibrium solutions the other solution curves increase or decrease397 oscillations are not possible Moreover assuming f is a wellbehaved function is continuous solution curves actually approach the equilibria getting closer and closer as time increases By uniqueness the curves never intersect the constant equilibrium solutions On the phase line if arrows on both sides of an equilibrium point toward that equilibrium point then we say the equilibrium point is an attractor If both of the arrows point away the equilibrium is called a repellerl Attractors are called asymptotically stable because if the system is in that constant equilibrium state and then it is given a small perturbation iiel a change or bump to a nearby state then it just returns to that state as t4 ltxgtl It is clear that real systems will seek out the stable states Repellers are unstable because a small perturbation can cause the system to go to a different equilib 13 Mathematical Models 33 graph of u 0 a b 0 Stable unstable semistable Figure 113 A generic plot showing which is u vsl u The points of intersection a b c on the u axis are the equilibrial The arrows on the uaxis or phase line show how the state u changes with time between the equilibrial The direction of the arrows is read from the plot of They are to the right when gt 0 and to the left when lt 0 The phase line can either be drawn as a separate line with arrows as in gure 111 or the arrows can be drawn directly on the uaxis of the plot as is done here gt Figure 114 Time series plots of 112 for different initial conditions The constant solutions are the equilibrial 34 1 Differential Equations and Models rium or even go off to in nity In the logistics model for population growth we observe gure 111 that the equilibrium u K is an asymptotically stable attractor and the zero population u 0 is unstable397 all solutions approach the carrying capacity u K at t e 00 Finally if one of the arrows points toward the equilibrium and one points away we say the equilibrium is semistable Semistable equilibria are not stable We emphasize that when we say an equilibrium if is asymptotically stable our understanding is that this is with respect to small perturbations To x the idea consider a population of sh in a lake that is in an asymptotically stable state u A small death event say caused by some toxic chemical that is dumped into the lake will cause the population to drop Asymptotic stability means that the system will return the original state u over time We call this local asymptotic stability If many sh are killed by the pollution then the perturbation is not small and there is no guarantee that the sh population will return to the original state u For example a catastrophe or bonanza could cause the population to jump beyond some other equilibriuml If the population returns to the state uquot for all perturbations no matter how large then the state u is called globally asymptotically stable A more precise de nition of local asymptotic stability can be given as follows An isolated equilibrium state uquot of 112 is locally asymptotically stable if there is an open interval I containing uquot with limtnaroo ut uquot for any solution u ut of 112 with u0 in 1 That is each solution starting in I converges to u Note that a semistable point is not asymptotically stable397 such points are in fact not stable Example 116 Dimensionless Models When we formulate a mathematical model we some times trade in the dimensioned quantities in our equations for dimensionless ones In doing so we obtain a dimensionless model often containing fewer pa rameters than the original model The idea is simple If for example time t is the independent variable in a model of population growth and a constant 7 with dimension time l7 representing the per capita growth rate appears in the model then the variable 739 tT l it has no dimensions that is it is dimensionless time divided by time It can serve as a new independent vari able in the model representing dimensionless time or time measured relative 1 is a time scale in the probleml Every to the inverse growth rate We say T variable in a model has a natural scale with which we can measure its relative value these scales are found from the parameters in the probleml The pop ulation p of an animal species in a geographical region can be scaled by the carrying capacity K of the region which is the number of animals the region 13 Mathematical Models 35 can support Then the variable P pK is dimensionless animals divided by animals and represents the fraction of the region s capacity that is lled If the carrying capacity is large the actual population p could be large requiring us to work with and plot big numbers However the dimensionless population P is represented by smaller numbers which are easier to deal with and plot For some models selecting dimensionless dependent and independent variables can pay off in great bene tsiit can help us understand the magnitude of var ious terms in the equations and it can reduce the number of parameters in a problem thus giving simpli cationi We illustrate this procedure for the initial value problem for the logistics model mm 7 a pm poA 113 There are two variables in the problem the independent variable t measured in time and the dependent variable p measured in animals There are three parameters in the problem the carrying capacity K and initial population p0 both measured in animals and the growth rate 39r measured in ltimei Let us de ne new dimensionless variables 739 rt tT 1 and P pKi These represent a dimensionless time and a dimensionless population P is measured relative to the carrying capacity and t is measured relative to the 1 growth rate39 the values K and T are called scalesi Now we transform the DE into the new dimensionless variablesi First we transform the derivative d d KP dP P 7 7 K dt dwT T E Then the logistics DE in 113 becomes dP KP TKE TKP1 7 Y dP P 1 7 P i d7 ln dimensionless variables 739 and P the parameters in the DE disappeared Next the initial condition becomes KP0 gig or 130 0 7 where a pOK is a dimensionless parameter animals divided by animals In summary the dimensioned model 113 with three parameters can be replaced by the dimensionless model with only a single dimensionless parameter a dP 7 E P1P7 130 a 114 36 1 Differential Equations and Models What this tells us is that although three parameters appear in the original problem only a single combination of those parameters is relevant We may as well work with the simpler equivalent dimensionless model 114 where populations are measured relative to the carrying capacity and time is measured relative to how fast the population is growing For example if the carrying capacity is K 3007 000 and the dimensioned p varies between 0 lt p lt 3007 000 it is much simpler to have dimensionless populations P with 0 lt P lt 1 Furthermore in the simpli ed form 114 it is easy to see that the equilibria are P 0 and P l the latter corresponding to the carrying capacity p Ki We have pointed out that an autonomous model can be easily analyzed qualitatively without ever nding the solution In this paragraph we introduce a simple method for solving a general autonomous equation 1 u 115 The method is called separation of variables If we divide both sides of the equation by we get 1 u 1 f 00 Now remembering that u is a function of t we integrate both sides with respect to t to obtain 1 u dt1dt0t0 K 7 where C is an arbitrary constant A substitution u ut du utdt reduces the integral on the left and we obtain 1 d t C 1 16 u i i f W This equation once the integral is calculated de nes the general solution u ut of 115 implicitly We may or may not be able to actually calculate the integral and solve for u in terms of t to determine an explicit solution u This method of separating the variables putting all the terms with u on the left side is a basic technique in differential equations397 it is adapted to more general equations in Chapter 2 Example 117 Consider the growthidecay model 1 Tu 1 17 13 Mathematical Models 37 where T is a given constant If T lt 0 then the equation models exponential decay if T gt 0 then the equation models exponential growth eigi popu lation growth as in the Malthus model We apply the separation of variables methodi Dividing by u we could divide by Tu but we choose to leave the constant the right side and taking antiderivatives gives ludtTdtCi u Because u dt du we can write lduTtCi u lnlul Tt C or eTHC ecemi Integrating gives This means either u ace or u fecemi Therefore the general solution of the growthidecay equation can be written compactly as ut Clem7 where Cl has been written for iec and is an arbitrary constant If an initial condition u0 uo 118 is prescribed on 117 it is straightforward to show that C1 uo and the solution to the initial value problem 1177118 is ut uoerti The growth iecay equation and its solution given in Example 116 occur often enough in applications that they are worthy of memorization The equa tion models processes like growth of a population mortality death growth of principal in a money account where the interest is compounded continuously at rate T and radioactive decay like the decay of Carbonl4 used in carbon datingl EXERCISES 1 The Allee effect At low population densities it may be dif cult for an animal to reproduce because of a limited number of suitable mates A population model that predicts this behavior is the Allee model Ci Allee 188571955 pTplt 71gtlt17gt7 0ltaltKi 1 Differential Equations and Models 0 CAD F 01 Find the per capita growth rate and plot the per capita rate vs pi Graph 17 vs p determine the equilibrium populations and draw the phase line Which equilibria are attractors and which are repellers Which are asymp totically stable From the phase line plot describe the long time behavior of the system for different initial populations and sketch generic time series plots for different initial conditions Modify the logistics model to include harvesting That is assume that the animal population grows logistically while at the same time animals are being removed by hunting shing or whatever at a constant rate of h animals per unit time What is the governing DE Determine the equilibriai Which are asymptotically stable Explain how the system will behave for different initial conditions Does the population ever become extinct i The Ricker population law is 70 0 Tpe 1 7 where 7 and a are constants Determine the dimensions of 7 and a At what population is the growth rate maximum Make a generic sketch of the per capita growth rate and write a brief explanation of how a population behaves under this law ls it possible to use the separation of variables method to nd a simple formula for pt In this exercise we introduce a simple model of growth of an individual organism over time For simplicity we assume it is shaped like a cube having sides equal to L Lti Organisms grow because they assimilate nutrients and then use those nutrients in their energy budget for maintenance and to build structure It is conjectured that the organism s growth rate in volume equals the assimilation rate minus the rate food is used Food is assimilated at a rate proportional to its surface area because food must ultimately pass across the cell walls397 food is used at a rate proportional to its volume because ultimately cells are threedimensional Show that the differential equation governing its size Lt can be written L t 07 bL where a and b are positive parameters What is the maximum length the organism can reach Using separation of variables show that if the length of the organism at time t 0 is L0 0 it is very small then the length is given by LO l 7 6 1 Does this function seem like a reasonable model for growth i In a classical ecological study of budworm outbreaks in Canadian r forests researchers proposed that the budworm population N was governed by the 13 Mathematical Models 39 533 90K law K where the rst term on the right represents logistics growth and where PN is a birdpredation rate given by N TN lt17 5 7 PN 7 aN2 7 N 2 122 Sketch a graph of the birdpredation rate vsi N and discuss its meaning PW What are the dimensions of all the constants and variables in the model Select new dimensionless independent and dependent variables by t 7N 7 bT b and reformulate the model in dimensionless variables and dimensionless constantsi Working with the dimensionless model show that there is at least one and at most three positive equilibrium populations What can be said about their stability Use the method of separation of variables to nd the general solution to the following autonomous differential equations a 1M b 11 6 29 c u l u2i d u Su 7 a where a is a constant e u 7 H uugi f 11 eugi i In Exercises 6 a7f nd the solution to the resulting lVP when u0 1 Find the general solution to the logistics equation u Tu1 7 using the separation of variables methodi Hint use the partial fractions decom position 1 7 1K 1K uK7u 77 Kiu Show that a solution curve that crosses the line u K2 has an in ection point at that position 1 Differential Equations and Models to H 53 H H H N H 9 H g H 5 i In the usual Malthus growth law N i Carbon dating The halflife of Carbonl4 is 5730 years That is7 it takes this many years for half of a sample of Carbonl4 to decay If the decay of Carbonl4 is modeled by the DE 1 iku where u is the amount of Carbonl4 nd the decay constant kl Answer 0000121 yr 1i In an artifact the percentage of the original Carbon14 remaining at the present day was measured to be 20 percent How old is the artifact In 1950 charcoal from the Lascaux Cave in France gave an average count of 009 disintegrations of 014 per minute per gram Living wood gives 6 68 disintegrations Estimate the date that individuals lived in the cave TN for a population of size N assume the growth rate is a linear function of food availability F that is 7 12F where b is the conversion factor of food into newborns Assume that FT is the total constant food in the system with FT F 5N where cN is amount of food already consumed Write down a differential equation for the population Ni What is the carrying capacity What is the population as t gets large One model of tumor growth is the Gompertz equation 13 iaRln 7 where R Rt is the tumor radius and a and k are positive constants Find the equilibria and analyze their stability Can you solve this differen tial equation for Rt A population model is given by p TPP7m where 39r and m are positive constants State reasons for calling this the explosioniextinction model i In a xed population of N individuals let I be the number of individuals infected by a certain disease and let S be the number susceptible to the disease with 15 Ni Assume that the rate that individuals are becoming infected is proportional to the number of infectives times the number of susceptibles7 or I aSL where the positive constant a is the transmission coef cient Assume no individual gets over the disease once it is contracted lf 0 Io is a small number of individuals infected at t 07 nd an initial value problem for the number infected at time t Explain how the disease evolves Over a long time7 how many contract the disease In Example lill we modeled the velocity of an object falling in a uid by the equation mv mgi av lf 110 07 nd an analytic formula for Mt 13 Mathematical Models 41 133 Stability and Bifurcation Differential equations coming from modeling physical phenomena almost al ways contain one or more parameters It is of great interest to determine how equilibrium solutions depend upon those parameters For example the logistics growth equation p 0 T100 f has two parameters the growth rate 7 and the carrying capacity Ki Let us add harvesting7 that is we remove animals at a constant rate H gt 0 We can think of a sh population where sh are caught at a given rate Hi Then we have the model 10 T100 7 7 Ht 119 We now ask how possible equilibrium solutions and their stability depend upon the rate of harvesting Hi Because there are three parameters in the problem we can nondimensionalize to simplify it We introduce new dimensionless variables by p F7 That is we measure populations relative to the carrying capacity and time u TTti relative to the inverse growth rate In terms of these dimensionless variables 1 19 simpli es to check this 11uliu7h7 where h HTK is a single dimensionless parameter representing the ratio of the harvesting rate to the product of the growth rate and carrying capacity We can now study the effects of changing h to see how harvesting in uences the steadystate sh populations in the model In dimensionless form we think of h as the harvesting parameter information about changing h will give us information about changing Hi The equilibrium solutions of the dimensionless model are roots of the quadratic equation fltugt ult17ugt7ho which are 1 1 u i 1i4hl The growth rate is plotted in gure 115 for different values of h For h lt 14 there are two positive equilibrium populations The graph of in this case is concave down and the phase line shows that the smaller one is unstable and the larger one is asymptotically stable As h increases these pop ulations begin to come together and at h 14 there is only a single unstable 42 1 Differential Equations and Models equilibriumi For h gt 14 the equilibrium populations cease to exist So when harvesting is small there are two equilibria one being stable39 as harvesting increases the equilibrium disappears We say that a bifurcation bifurcation means dividing occurs at the value h 14 This is the value where there is a signi cant change in the character of the equilibriai For h 2 14 the popula tion will become extinct regardless of the initial condition because lt 0 for all All these facts can be conveniently represented on a bifurcation diagrami See gure 116 In a bifurcation diagram we plot the equilibrium solutions uquot vs the parameter h In this context h is called the bifurcation parameter The plot is a parabola opening to the left We observe that the upper branch of the parabola corresponds to the larger equilibrium and all so lutions represented by that branch are asymptotically stable39 the lower branch corresponding to the smaller solution is unstable u h390 A h18 h14 h12 phase line A unstable stable Figure 115 Plots of ul 7 u 7 h for different values of hr The phase line is plotted in the case h lSl Finally we give an analytic criterion that allows us to determine stability of an equilibrium solution by simple calculusi Let 1 120 be a given autonomous systems and uquot an isolated equilibrium solution so that 0 We observe from figure 113 that when the slope of the graph of at the equilibrium point is negative the graph falls from left to right and both arrows on the phase line point toward the equilibrium point Therefore a condition that guarantees the equilibrium point u is asymptotically stable is f u lt 0 Similarly if the graph of strictly increases as it passes through the equilibrium then f u gt 0 and the equilibrium is unstable If the slope of is zero at the equilibrium then any pattern of arrows is possible and there is no information about stabilityi lf 0 then uquot is a critical point of 13 Mathematical Models 43 stable Figure 116 Bifurcation diagram plot of the equilibrium solution as a func tion of the bifurcation parameter h uquot i vl 7 4h For h gt i there are no equilibria and for h lt there are two with the larger one being stablel A M MH bifurcation occurs at h f and could be a local maximum local minimum or have an inflection point If there is a local maximum or local minimum then uquot is semistable which is not stable If there is an inflection point then f changes sign at uquot and we obtain either a repeller or an attractor depending on how the concavity changes negative to positive or positive to negative Theorem 118 Let uquot be an isolated equilibrium for the autonomous system 120 If f lt 07 then uquot is asymptotically stable397 if f 11 gt 07 then u is unstable If f 07 then there is no information about stabilityl Example 119 Consider the logistics equation u Tul 7 The equilibria are uquot 0 and u Kl The derivative of is 7 7 2TuKl Evaluating the derivative at the equilibria gives fO 39r gt 07 7T lt 0 Therefore u 0 is unstable and u K is asymptotically stable EXERCISES 1 A sh population in a lake is harvested at a constant rate and it grows logisticallyl The growth rate is 02 per month the carrying capacity is 40 0 9 F 1 Differential Equations and Models thousand and the harvesting rate is 15 thousand per month Write down the model equation nd the equilibria and classify them as stable or unstable Will the sh population ever become extinct What is the most likely longterm sh population i For the following equations nd the equilibria and sketch the phase linei Determine the type of stability of all the equilibria Use Theorem 118 to con rm stability or instabilityi a 1 u237 b 1 2ul 7 u 7 c 1 4 7 u2 7 u For the following models which contain a parameter h nd the equilibria in terms of h and determine their stability Construct a bifurcation diagram showing how the equilibria depend upon h and label the branches of the curves in the diagram as unstable or stable a 1 hu7 u2l b 117 uu2 7hi Consider the model 1 A 7 bu 7 aug where a and b are xed positive constants and A is a parameter that may vary a If A lt b show that there is a single equilibrium and that it is asymp totically stablel b If A gt I nd all the equilibria and determine their stability c Sketch the bifurcation diagram showing how the equilibria vary with Al Label each branch of the curves shown in the bifurcation diagram as stable or unstable l The biomass P of a plant grows logistically with intrinsic growth rate 7 and carrying capacity Ki At the same time it is consumed by herbivores at a rate aP b P per herbivore where a and b are positive constants The model is P aPH 7 7 7 P 7 TP 1 K b P where H is the density of herbivores Assume aH gt M and assume 7 K a and b are xed Plot as a function of P the growth rate and the consumption rate for several values of H on the same set of axes and 13 Mathematical Models 45 identify the values of P that give equilibriai What happens to the equilibria as the herbivory H is steadily increased from a small value to a large value Draw a bifurcation diagram showing this effect That is plot equilibrium solutions vs the parameter Hi If herbivory is slowly increased so that the plants become extinct and then it is decreased slowly back to a low level do the plants return 533 A deer population grows logistically and is harvested at a rate proportional to its population size The dynamics of population growth is modeled by P P 7 P1 7 7 AP where A is the per capita harvesting rater Use a bifurcation diagram to explain the effects on the equilibrium deer population when A is slowly increased from a small value to a large value gt7 Draw a bifurcation diagram for the model u u3 7 u h7 where h is the bifurcation parameteri Label branches of the curves as stable or unstable 90 A where A is a parameterl Draw the Consider the model u uu 7 e bifurcation diagram plotting the equilibrium solutions uquot vsi Al Label each curve on the diagram as stable or unstable 134 Heat Transfer An object of uniform temperature To eg a potato is placed in an oven of temperature Tel It is observed that over time the potato heats up and eventually its temperature becomes that of the oven environment Tel We want a model that governs the temperature Tt of the potato at any time t Newton s law of cooling heating a constitutive model inferred from experiment dictates that the rate of change of the temperature of the object is proportional to the difference between the temperature of the object and the environmental temperature That is T 7hT 7 T8 121 The positive proportionality constant h is the heat loss coef cient There is a fundamental assumption here that the heat is instantly and uniformly distributed throughout the body and there are no temperature gradients or spatial variations in the body itself From the DE we observe that T T8 is an equilibrium solutioni lfT gt T8 then T lt 07 and the temperature decreases 7 if T lt Te then T gt 07 and the temperature increases Plotting the phase line easily shows that this equilibrium is stable Exercisel 46 1 Differential Equations and Models We can nd a formula for the temperature Tt satisfying 121 using the separation of variables method introduced in the last section Here for variety we illustrate another simple method that uses a Change of variables Let u T 7 Tel Then 1 T and 121 may be written 1 7hui This is the decay equation and we have memorized its general solution u Ce hti Therefore T 7 Te Ce h 7 or Tt Te Ce hti This is the general solution of 121 If we impose an initial condition T0 To then one nds C To 7 Te7 giving Tt Te T0 7 Tad We can now see clearly that Tt 7gt T8 as t 7gt 00 A plot of the solution showing how an object heats up is given in gure 117 o39 t Figure 117 Temperature history in Newton s law of cooling If the environmental or ambient temperature fluctuates then T8 is not constant but rather a function of time T8 The governing equation becomes T 7hT 7 Tea In this case there are no constant or equilibrium solutionsi Writing this model in a different way T 7hT hTe The rst term on the right is internal to the system the body being heated and considered alone with zero ambient temperature leads to an exponentially 7hT has solution T Ce htl There fore there is a transient governed by the natural system that decays away The decaying temperature recall that T 13 Mathematical Models 47 external environmental temperature Te t gives rise to timedependent dynam ics and eventually takes over to drive the system39 we say the system is driven or forced by the environmental temperature In Chapter 2 we develop methods to solve this equation with time J r J in the 39 function EXERCISES 1 A small solid initially of temperature 22 C is placed in an ice bath of 0 C It is found experimentally by measuring temperatures at different times that the natural logarithm of the temperature Tt of the solid plots as a linear function of time t39 that is lnT 7atb Show that this equation is consistent with Newton s law of cooling 1f the temperature of the object is 8 C degrees after two hours what is the heat loss coef cient When will the solid be 2 C N A small turkey at room temperature 70 F is places into an oven at 350013 If h 042 per hour is the heat loss coef cient for turkey meat how long should you cook the turkey so that it is uniformly 200 F Comment on the validity of the assumptions being made in this model CA3 The body of a murder victim was discovered at 1100 AM The medical examiner arrived at 1130 AM and found the temperature of the body was 946013 The temperature of the room was 70 F One hour later in the same room he took the body temperature again and found that it was 934013 Estimate the time of death g Suppose the temperature inside your winter home is 68 F at 100 PM and your furnace then fails 1f the outside temperature is 10 F and you notice that by 1000 PM the inside temperature is 57 F what will be the temperature in your home the next morning at 600 AM 5 The temperature TO of an exothermic chemically reacting sample placed in a furnace is governed by the initial value problem T 7kT 7 Te qe eT T0 To where the term qe eT is the rate heat is generated by the reaction What are the dimensions of all the constants 16 Te 4 To and 9 in the problem Scale time by k 1 and temperature by T8 to obtain the dimensionless model 4w 7 1 ae b p7 we v 48 1 Differential Equations and Models for appropriately chosen dimensionless parameters a b and cl Fix a 1 HOW many positive equilibria are possible depending upon the value of l Hint Graph the heat loss term and the heat generation term vs 1 on the same set of axes for different values of 1 135 Chemical Reactors A continuously stirred tank reactor also called a chemostat or compart ment is a basic unit of many physical chemical and biological processes A continuously stirred tank reactor is a wellde ned geometric volume or entity Where substances enter react and are then discharged A chemostat could be an organ in our body a polluted lake an industrial chemical reactor or even an ecosystemi See gure 118 V CW 139 Cin 100 Figure 118 A chemostat or continuously stirred tank reactori We illustrate a reactor model With a speci c example Consider an industrial pond With constant volume V cubic meters Suppose that polluted water con taining a toxic chemical of concentration Cir grams per cubic meter is dumped into the pond at a constant volumetric ow rate of 4 cubic meters per day At the same time the continuously mixed solution in the pond is drained off at the same ow rate qr If the pond is initially at concentration C0 What is the concentration Ct of the chemical in the pond at any time t The key idea in all chemical mixture problems is to obtain a model by conserving mass the rate of change of mass in the pond must equal the rate mass ows in minus the rate mass ows out The total mass in the pond at any time is VC and the mass ow rate is the volumetric ow rate times the mass concentration397 thus mass balance dictates VC 4C 7 40 Hence the initial value problem for the chemical concentration is VC inn 7 4107 00 Co 122 Where C0 is the initial concentration in the tank This initial value problem can be solved by the separation of variables method or the change of variables method Section 134 See Exercise 1 13 Mathematical Models 49 Now suppose we add degradation of the chemical while it is in the pond assuming that it degrades to inert products at a rate proportional to the amount present We represent this decay rate as T0 gm per cubic meter per day where 7 is constant Then the model equation becomes VC 4C 7 qC 7 TVCi Notice that we include a factor V in the last term to make the model dimen sionally correct A similar model holds when the volumetric ow rates in and out are different which gives a changing volume Vti Letting gin and gout denote those ow rates respectively we have Vt0 7 qum 7 game 7 New where Vt V0 gm 7 gou t and where V0 is the initial volume Methods developed in Chapter 2 show how this equation is solved EXERCISES H Solve the initial value problem 122 and obtain a formula for the concen tration in the reactor at time t E0 An industrial pond having volume 100 m3 is full of pure water Contami nated water containing a toxic chemical of concentration 00002 kg per m3 is then is pumped into the pond with a volumetric ow rate of 05 m3 per minute The contents are wellmixed and pumped out at the same ow rate Write down an initial value problem for the contaminant concentration Ct in the pond at any time t Determine the equilibrium concentration and its stability Find a formula for the concentration Ct 3 In the preceding problem change the ow rate out of the pond to 06 m3 per minute How long will it take the pond to empty Write down a revised initial value problemi g i A vat of volume 1000 gallons initially contains 5 lbs of salt For t gt 0 a salt brine of concentration 01 lbs per gallon is pumped into the tank at the rate of 2 gallons per minute397 the perfectly stirred mixture is pumped out at the same ow rate Derive a formula for the concentration of salt in the tank at any time t Check your answer on a computer algebra system and sketch a graph of the concentration vs time 5 Consider a chemostat of constant volume where a chemical C is pumped into the reactor at constant concentration and constant ow rate While in the reactor it reacts according to C C 7 products From the law of mass action the rate of the reaction is 39r kCQ where k is the rate constant If the concentration of C in the reactor is given by Ct then 1 Differential Equations and Models 533 to mass balance leads the governing equation VC qciniqcikVCQi Find the equilibrium states and analyze their stability Redo this problem after nondimensionalizing the equation pick time scale Vq and concentration scale Ciquot Work Exercise 5 if the rate of reaction is given by MichaelisiMenten kinetics aC Tb C Where a and b are positive constants A batch reactor is a reactor of volume V Where there are no in and out ow rates Reactants are loaded instantaneously and then allowed to react over a time T called the residence time Then the contents are expelled instantaneously Fermentation reactors and even sacular stomachs of some animals can be modeled as batch reactors If a chemical is loaded in a batch reactor and it degrades With rate 7 O kC given in mass per unit time per unit volume What is the residence time required for 90 percent of the chemical to degrade l The Monod equation for conversion of a chemical substrate of concen tration C into its products is def aC E 7 ibJr C Where a and b are positive constants This equation With MichaelisiMenten kinetics describes how the substrate is being used up through chemical reaction If in addition to reaction the substrate is added to the solution at a constant rate R Write down a differential equation for C Find the equilibrium solution and explain how the substrate concentration evolves for various initial conditions i Consider the chemical reaction AB A C Where one molecule of A reacts With one molecule of B to produce one molecule of C and the rate of the reaction is k the rate constant By the laW of mass action in chemistry the reaction rate is 7 kab Where a and I represent the timedependent concentrations of the reactants A and Bi Thus the rates of change of the reactants and product are governed by the three equations a 71ml 1 71ml 5 kabl If initially a0 a0 120 b0 and 00 0 With co gt 20 nd a single rstorder differential equation that involves only the concentration a atl What is the limiting concentration limtmgt00 at What are the other two limiting concentrations 13 Mathematical Models 51 10 Digestion in the stomach gut in some organisms can be modeled as a chemical reactor of volume V Where food enters and is broken down into nutrient products Which are then absorbed across the gut lining397 the food7 product mixture in the stomach is perfectly stirred and exits at the same rate as it entered Let So be the concentration of a substrate food con sumed at rate 4 volume per time In the gut the rate of substrate break down into the nutrient product 5 7 P7 is given by kVS Where k is the rate constant and S St is the substrate concentration The nutrient product of concentration P Pt7 is then absorbed across the gut bound ary at a rate aVP Where a is the absorption constanti At all times the contents are thoroughly stirred and leave the gut at the ow rate q a Show that the model equations are VS L150 7 L15 7 kVS VP kVS 7 aVP 7 413 b Suppose the organism eats continuously in a steadystate mode Where the concentrations become constanti Find the steadysteady or equi librium concentrations 58 and Pei 0 V Some ecologists believe that animals regulate their consumption rate in order to maximize the absorption rate of nutrients Show that the maximum nutrient concentration P8 occurs When the consumption rate is L a kVi 39 39 39 akSQV d Show that the maximum absorptlon rate 1s therefore J E i 136 Electric Circuits Our modern technological society is lled With electronic devices of all types At the basis of these are electrical circuits The simplest circuit unit is the loop in gure 119 that contains an electromotive force emf Et a battery or generator that supplies energy a resistor an inductor and a capacitor all connected in series A capacitor stores electrical energy on its two plates a resistor dissipates energy usually in the form of heat and an inductor acts as a choke that resists changes in current A basic law in electricity Kirchhoff s law tells us that the sum of the voltage drops across the circuit elements as measured eg by a voltmeter in a loop must equal the applied emfi ln symbols VL VR V0 52 1 Differential Equations and Models Figure 119 An RCL circuit With an electromotive force Et supplying the electrical energy This law comes from conservation of energy in a current loop and it is derived in elementary physics textsl A voltage drop across an element is an energy potential that equals the amount of work required to move a charge across that element Let I t denote the current in amperes or charge per second in the circuit and let 4 qt denote the charge in coulombs on the capacitorl These quantities are related by q 1 There are several choices of state variables to describe the response of the circuit charge on the capacitor 4 current I or voltage V0 across the capacitorl Let us Write Kirchhoff s law in terms of charge By Ohm s law the voltage drop across the resistor is proportional to the current or VR RI Where the proportionality constant R is called the resistance measured in ohmsl The voltage drop across a capacitor is proportional to the charge on the capacitor or 71 047 Where C is the capacitance measured in farads Finally the voltage drop V0 across an inductor is proportional to how fast the current is changing or VL LI7 Where L is the inductance measured in henrysl Substituting these voltage drops into Kirchhoff s law gives 1 LI RI 54 Et 13 Mathematical Models 53 or using 4 I 1 54 This is the RCL Circuit equation which is a secondorder DE for the charge q The initial conditions are Lq Rq q0 qo qO 10 0 These express the initial charge on the capacitor and the initial current in the circuit Here Et may be a given constant eigi Et 12 for a lZvolt battery or may be a oscillating function of time t eigi Et A coswt for an alternating voltage potential of amplitude A and frequency m If there is no inductor then the resulting RC circuit is modeled by the rstorder equation 1 7 54 lf Et is constant this equation can be solved using separation of variables or the change of variables method Exercise 2 We show how to solve secondorder differential equations in Chapter 3 Rq E EXERCISES 1 Write down the equation that governs an RC circuit With a lZvolt battery taking R l and C 5 Determine the equilibrium solution and its stability lf q0 5 nd a formula for qti Find the current 1ti Plot the charge and the current on the same set of axes 0 i In an arbitrary RC circuit With constant emf E use the method of sepa ration of variables to derive the formula 4t KequotRC EC for the charge on the capacitor Where K is an arbitrary constant If 970 40 What is K 9 An RCL circuit With an applied emf given by E t has initial charge 970 go and initial current 10 10 What is I m Write down the circuit equation and the initial conditions in terms of current 1ti F Formulate the governing equation of an RCL circuit in terms of the current 75 When the circuit has an emf given by Et Acoswti What are the appropriate initial conditions 539 Find the DE model for the charge in an LC circuit With no emf Show that the response of the circuit may have the form qt Acoswt for some amplitude A and frequency wt 54 1 Differential Equations and Models 6 Consider a standard RCL circuit With no emf but With a voltage drop across the resistor given by a nonlinear function of current 1 1 3 V37 lt I 71 This replaces Ohm s law If C L 1 nd a differential equation for the current Ht in the circuit 7 Write the RCL circuit equation With the voltage as the unknown state functioni 2 Analytic Solutions and Approximations In the last chapter we studied several rstorder DE models and a few elemen tary techniques to help understand the qualitative behavior of the models In this chapter we introduce analytic solution techniques for rstorder equations and some general methods of approximation including numerical methods 21 Separation of Variables In Section 132 we presented a simple algorithm to obtain an analytic solution to an autonomous equation 1 u called separation of variables Now we show that this method is applicable to a more general class of equations A separable equation is a rstorder differential where the right side can be factored into a product of a function oft and a function of u That is a separable equation has the form 1 gthui 21 To solve separable equations we take the expression involving u to the left side and then integrate with respect to t remembering that u Therefore dividing by hu and taking the antiderivatives of both sides with respect to t gives 1 7 Wu dt 7 gtdt0 56 2 Analytic Solutions and Approximations where C is an arbitrary constant of integration Both antiderivatives generate an arbitrary constant but we have combined them into a single constant C Next we change variables in the integral on the left by letting u ut7 so that du utdti Hence 1 Wdu gtdt C This equation once the integrations are performed yields an equation of the form Ct C7 22 which de nes the general solution u implicitly as a function of t We call 22 the implicit solution To obtain an explicit solution u ut we must solve 22 for u in terms of t397 this may or may not be possible As an aside we recall that if the antiderivatives have no simple expressions then we write the antiderivatives with limits on the integrals Example 21 Solve the initial value problem 7 t 1 7 2u 7 We recognize the differential equation as separable because the right side is 0 1 Bringing the 2u term to the left side and integrating gives u0 ll 2uu dt t1dt C 2udu t2tCi Therefore 1 u2 it t 0 This equation is the general implicit solution We can solve for u to obtain two forms for explicit solutions 1 ui t2t0i Which sign do we take The initial condition requires that u be positive Thus we take the plus sign and apply u0 l to get C 1 The solution to the initial value problem is therefore 1 t2 t 1 u HZ 21 Separation of Variables 57 This solution is valid as long as the expression under the radical is not negative In the present case the solution is de ned for all times t 6 R and so the interval of existence is the entire real line Example 22 Solve the initial value problem 2 it 1 7 ul 4 t Note that we might expect trouble at t 0 because the derivative is unde ned there The equation is separable so we separate variables and integrate with l u equot EWdtTdta We can integrate the left side exactly but the integral on the right cannot be respect to t resolved in closed formi Hence we write it with variable limits and we have t 67 Tdt Cl 1 Judiciously we chose the lower limit as t 1 so that the initial condition would be easy to apply Clearly we get C 2 Therefore teit dt2 1 t utlt1 edt2gt3 This solution is valid on 0 lt t lt 00 In spite of the apparent complicated form of the solution which contains an integral it is not dif cult to plot using a computer algebra systemi The plot is shown in gure 21 The method of separation of variables is a key technique in differential equations Many important models turn out to be separable not the least of which is the autonomous equationi EXERCISES 1 Find the general solution in explicit form of the following equations a u b u 111129 cosu 58 2 Analytic Solutions and Approximations Figure 21 This plot is obtained on the interval 006 using the Maple command plotevalf2intexpss slntA2 t0ii0i6 y0n6i c 1 tlu2li d u u0i E0 Find the solution to the initial value problem 1t2e u7 u0 ln 27 and determine the interval of existence 3 Draw the phase line associated With the DE 1 u4 7 112 and then solve the DE subject to the initial condition u0 ll Hint for the integration you Will need a partial fractions expansion 1 a b c u47u2 2u27u7 Where a b and c are to be determined 4 Find the general solution in implicit form to the equation 4 7 2t u i Bu 7 5 Find the solution When u1 3 and plot the solution What is its interval of existence 21 Separation of Variables 7 5X7 50 H O 2213 Solve the initial value problem 1 my ut0 uo and nd the interval of existence when uo lt 0 when uo gt 0 and when uo 0 Find the general solution of the DE 1 6tu 7 l23l Show that there is no value of the arbitrary constant that gives the solution u 1 A solution to a DE that cannot be obtained from the general solution by xing a value of the arbitrary constant is called a singular solution Find the general solution of the DE T2 7 t2u tu 07 where T is a xed positive parameterl Find the solution to the initial value problem when uT2 1 What is the interval of existence Allometric growth describes temporal relationships between sizes of dif ferent parts of organisms as they grow eg the leaf area and the stem diameter of a plant We say two sizes ul and ug are allometn39cally related if their relative growth rates are proportional or a1 u2 7 ul Show that if ul and ug are allometrically related then ul Cu7 for some constant C A differential equation of the form He where the right depends only on the ratio of u and t is called a homo geneous Show that the substitution u ty transforms a homogeneous equation into a rstorder separable equation for y Use this method to solve the equation 4t2 3112 u l Ztu l Solve the initial value problem d dt 11662 7 u0 3i 6 2 Analytic Solutions and Approximations H H H 0 4 CAD H t H 5 m 533 m 7 i Find the general solution u uT of the DE Tum in where p is a positive constant i A population of uo individuals all has HIV but none has the symptoms of AIDS Let ut denote the number that does not have AIDS at time t gt 0 If Tt is the per capita rate of individuals showing AIDS symptoms the conversion rate from HIV to AIDS then u u 77ft In the simplest case we can take 7 to be a linear function of time or Tt at Find ut and sketch the solution when a 0 2 and uo 100 At what time is the rate of conversion maximum i An arrow of mass m is shot vertically upward with initial velocity 160 ftseci It experiences both the deceleration of gravity and a deceleration of magnitude mv2800 due to air resistance How high does the arrow go In very cold weather the thickness of ice on a pond increases at a rate inversely proportional to its thickness If the ice initially is 005 inches thick and 4 hours later it is 007 inches thick how thick will it be in 10 hours Write the solution to the initial value problem 1 1 7112642 u0 5 in terms of the erf function erft f e sgdsi Use separation of variables to solve the following problems Write the so lution explicitly when possible a u ptu7 where pt is a given continuous function b 1 72m u1 2 Plot the solution on 0 g t 2 7211 0 lt t lt 1 7 7 cu7 7u27 1 t 2 7 u075i Find a continuous solution on the interval 0 g t g 3 and plot the solution A certain patch is populated with a cohort of newly hatched grasshoppers numbering ho As time proceeds they die of natural causes at per capita rate m and they are eaten by spiders at the rate aH1bH per spider where 22 FirstOrder Linear Equations 61 H H is the population of grasshoppers and a and b are positive constants Thus the dynamics is given by aH H imH 7 S 1 bH 7 where S is the spider population and time is given in days a Determine the units on the constants m a and 12 b Choose new dimensionless variables 739 mt and h 12H7 and reformu late the differential equation and initial condition in a dimensionless problem for h h7 i In your differential equation you should have a single dimensionless constant given by A aSmi c Solve the dimensionless initial value problem to obtain a formula for h7 i What is lirnTA00 h7 i Let N0 be the number of individuals in a cohort at time t 0 and N NO be the number of those individuals alive at time t If m is the constant per capita mortality rate then N N 7m which gives Nt Noe mh The survivorship function is de ned by St NtN0 and ST therefore gives the probability of an individual living to age T In the case of a constant per capita mortality the survivorship curve is a decaying exponential a What fraction die before age T Calculate the fraction of individuals that die between age a and age 12 b If the per capita death rate depends on time or m mt nd a formula for the survivorship function your answer will contain an in tegrali c What do you think the human survivorship curve looks like 22 FirstOrder Linear Equations A differential equation of the form a mm m 23gt is called a rstorder linear equation The given functions p and q are as sumed to be continuous lf qt 0 then the equation is called homogeneous7 62 2 Analytic Solutions and Approximations otherwise it is called nonhomogeneous Linear equations have a nice struc ture to their solution set and we are able to derive the general solution The homogeneous equation 1 ptu7 24 without the nonhomogeneous term qt can readily be solved by separation of variables to obtain uhtCeP where Ptptdt 25 where C is an arbitrary constant We have placed a subscript h on this solution to distinguish it from the solution of the nonhomogeneous equation 23 The solution uht to the homogeneous equation is sometimes called the comple mentary solution397 we just refer to it as the homogeneous solutioni Also note that we have used the inde nite integral notation for the function PO in some cases we have to represent Pt in the form with limits of integration We always choose the lower limit a to be the value of time where the initial condition is prescribed We now describe a standard technique to solve the nonhomogeneous equa tion 23 The idea is to try a solution of the form 25 where we let the constant C in the homogeneous solution vary as a function of t397 we then substi tute this form into 23 to determine the C Cti The method is for obvious reasons called variation of parameters1 Thus assume a solution to 23 of the form ut CteP i Then plugging in may CteP P t pt0tePlt gt qa But P p and therefore two of the terms cancel giving OWEN 407 Ct e P qti Integration yields Ct e P qtdtK 1 Another method using integrating factors is presented in the Exercises 22 FirstOrder Linear Equations 63 where K is a constant of integration So we have ut e P qtdtKgt GPO Kepm GPO e Ptqtdt 26 which is the general solution to the general linear nonhomogeneous equation 23 If the antiderivative in the last equation cannot be calculated explicitly then we write the solution as L ut Kelp eP e P7qTdTi We urge the reader not to memorize these formulas397 rather remember the method and apply it to each problem as you solve it Example 23 Find the general solution to l t u 7tu The homogeneous equation is u u and has solution uht Cefumn Calm Ctr Therefore we vary the parameter C and assume a solution of the original non homogeneous equation of the form ut Ctti Substituting into the equation we get u Ct C tt Cttt3 C t t Therefore Ct t2dt t3 K and the general solution to the original equation is ut tg Kgt t Kti The arbitrary constant K can be determined by an initial condition 64 2 Analytic Solutions and Approximations Example 24 Consider the DE 1 Zu t 27 The associated homogeneous equation is u Zu which has solution uh C623 Therefore we assume the solution of 27 is of the form ut Ct62 i Substituting into the original equation gives Ct2e2 Owe 20mgt t C t te Q i Integrating Ct te 2 dt K fie m t 1 K The integral was calculated analytically using integration by parts Therefore the general solution of 27 is ut lt7ie 2t2t 1 Kgt 6 KeQ 7i2t1i Notice that the general solution is composed of two terms uht and upt de ned by 1 uh K6217 up 71275 1 We know uh is the general solution to the homogeneous equation and it is easy to show that up is a particular solution to the nonhomogeneous equation 27 So the general solution to the nonhomogeneous equation 27 is the sum of the general solution to the associated homogeneous equation and any particular solution to the nonhomogeneous equation 27 Example 24 illustrates a general principle that reveals the structure of the solution to a rstorder linear DE The general solution can be written as the sum of the solution to the homogeneous equation and any particular solution of the nonhomogeneous equationi Precisely the basic structure theorem for rstorder linear equations states 22 FirstOrder Linear Equations 65 Theorem 25 Structure Theorem The general solution of the nonhomogeneous equation u mm W is the sum of the general solution uh of the homogeneous equation 1 ptu and a particular solution up to the nonhomogeneous equation In symbols 1W 11W upt7 Where uht Kepm and up GPO fe P qtdt and Where Pt fptdti Example 26 Consider an RC electrical circuit Where the resistance is R 1 and the capac itance is C 05 Initially the charge on the capacitor is 40 5 The current is driven by an emf that generates a variable voltage of sinti How does the circuit respond The governing DE for the charge qt on the capacitor is 1 Rq at sint7 or substituting the given parameters 4 724 sinti 28 The homogeneous equation 4 24 has solution qh Ce mi We assume the solution to the nonhomogeneous equation has the form L Cte 2 i Substi tuting into 28 gives Ct 724672 C tqe 2t 72Cte 2t sint7 Ct 62 sinti Integrating 2t 2t 2 1 Ct e s1ntdtKe gslntigcostK7 Where K is a constant of integration The integral was calculated using software or one can use integration by parts Therefore the general solution of 28 is 72t 2 1 72E qtCte gslntigcostKe l 66 2 Analytic Solutions and Approximations Next we apply the initial condition q0 5 to obtain K 265 Therefore the solution to the initial value problem is 2 1 26 72 qt 7 gslnti Ecost Fe i The solution is consistent with Theorem 25 Also there is an important phys ical interpretation of the solution The homogeneous solution is the transient response qht e 2t that depends upon the initial charge and decays over a time397 what remains over a long time is the particular solution which is regarded as the steadystate response qpt sint 7 cos t The ho mogeneous solution ignores the forcing term the emf whereas the particular solution arises from the forcing termi After a long time the applied emf drives the response of the systemi This behavior is characteristic of forced linear equa tions coming from circuit theory and mechanics The solution is a sum of two terms a contribution due to the internal system and initial data the decay ing transient and a contribution due to the external forcing term the steady responsei Figure 22 shows a plot of the solution 3 iransieni i swamestate response Figure 22 Response of the circuit in Example 26 showing the initial tran sient and the longtime steadystatel 22 FirstOrder Linear Equations 67 Example 27 Sales Response to Advertising The eld of economics has always been a source of interesting phenomena modeled by differential equations In this example we set up a simple model that allows management to assess the effectiveness of an advertising campaign Let S 5t be the monthly sales of an itemi In the absence of advertising it is observed from sales history data that the logarithm of the monthly sales decreases linearly in time or lnS fat 12 Thus 5quot iaS and sales are modeled by exponential decay To keep sales up advertising is required If there is a lot of advertising then sales tend to saturate at some maximum value 5 M this is because there are only nitely many consumers The rate of increase in sales due to advertising is jointly proportional to the advertising rate At and to the degree the market is not saturated397 that is TAO gt i The constant 7 measures the effectiveness of the advertising campaign The term MAZS is a measure of the market share that has still not purchased the producti Then combining both natural sales decay and advertising we obtain the model M S S iaS39 TA t 7 i lt gt lt M gt The rst term on the right is the natural decay rate and the second term is the rate of sales increase due to advertising which drives the sales As it stands because the advertising rate A is not constant there are no equilibria constant solutions We can rearrange the terms and write the equation in the form 5 a gt 5 may 29 Now we recognize that the sales are governed by a rstorder linear DE The Exercises request some solutions for different advertising strategiesi EXERCISES 1 Find the general solution of u 7u ti 2 Find the general solution of u 7u e 3 Show that the general solution to the DE 1 au 1 t is given by z ut Ce a e aquot51 sdsi 0 2 Analytic Solutions and Approximations H H H H N H g 15335 5 0 50 i A decaying battery generating 2006 5 volts is connected in series with a 20 ohm resistor and a 001 farad capacitori Assuming q 0 at t 0 nd the charge and current for all t gt 0 Show that the charge reaches a maximum and nd the time it is reached Solve u u 3t by introducing y uquot Solve u t u2 by letting y t u i Express the general solution of the equation 1 Ztu l in terms of the erf functioni Find the solution to the initial value problem 1 pu 17 u0 uo7 where p and q are constants Find a formula for the general solution to the DE 1 1011 405 where p is constant Find the solution satisfying ut0 uoi A differential equation of the form 1 atu gtun is called a Bernoulli equation and it arises in many applications Show that the Bernoulli equation can be reduced to the linear equation y 1 naty 1 i 7090 by changing the dependent variable from u to y via y ul Solve the Bernoulli equations see Exercise 10 l a au Stuui b 1 ul ue H 1 c u tumi Initially a tank contains 60 gal of pure water Then brine containing 1 lb of salt per gallon enters the tank at 2 galmini The perfectly mixed solution is drained off at 3 galmini Determine the amount in lbs of salt in the tank up until the time it emptiesi A large industrial retention pond of volume V initially free of pollutants was subject to the in ow of a contaminant produced in the factory s pro cessing planti Over a period of 12 days the EPA found that the in ow con centration of the contaminant decreased linearly in time to zero from its initial value of a grams per volume its ow rate 4 volume per day be ing constanti During the 12 days the spillage to the local stream was also 22 FirstOrder Linear Equations H F H 01 H 533 H I H 5 0 H to L What is the concentration in the pond after 12 days Do a numerical ex periment using a computer algebra system where V 6000 cubic meters 12 20 days a 003 grams per cubic meter and q 50 cubic meters per day With this data how long would it take for the concentration in the pond to get below the required EPA level of 000001 grams per cubic meter if fresh water is pumped into the pond at the same ow rate with the same spillover Determine the dimensions of the various quantities in the salesiadvertising model 29 lfA is constant what is the equilibrium i Technology Transfer Suppose a new innovation is introduced at time t 0 in a community of N possible users eg a new pesticide introduced to a community of farmers Let zt be the number of users who have adopted the innovation at time t If the rate of adoption of the innovation is jointly proportional to the number of adoptions and the number of those who have not adopted write down a DE model for Describe qualitatively how zt changes in time Find a formula for A house is initially at 12 degrees Celsius when its heatingicooling system fails The outside temperature varies according to T8 9 10 cos 27d where time is given in days The heat loss coef cient is h 3 degrees per day Find a formula for the temperature variation in the house and plot it along with T8 on the same set of axes What is the time lag between the maximum inside and outside temperature i In the sales response to advertising model 29 assume 50 SO and that advertising is constant over a xed time period T and is then removed That is 0 S t g T t gt T a A t 7 lt gt 07 Find a formula for the sales 50 Hint solve the problem on two intervals and piece together the solutions in a continuous way In a community having a xed population N the rate that people hear a rumor is proportional to the number of people who have not yet heard the rumor Write down a DE for the number of people P who have heard the rumor Over a long time how many will hear the rumor Is this a believable model i An object of mass m l is dropped from rest at a large height and as it falls it experiences the force of gravity mg and a timedependent resistive force of magnitude FT u where v is its velocity Write down an initial value problem that governs its velocity and nd a formula for the solution 2 Analytic Solutions and Approximations 20 The MacArthur7Wilson model of the dynamics of species eg bird 0 H species that inhabit an island located near a mainland was developed in the 1960s Let P be the number of species in the source pool on the mainland and let S 50 be the number of species on the island Assume that the rate of change of the number of species is 5X7 where X is the colonization rate and M is the extinction rate In the MacArthuriWilson model E F37 where I and E are the maximum colonization and extinction rates respec XI17 and M tively a Over a long time what is the expected equilibrium for the number of species inhabiting the island ls this equilibrium stable Given 50 50 nd an analytic formula for 50 O U VV Suppose there are two islands one large and one small with the larger island having the smaller maximum extinction rater Both have the same colonization rater Show that the smaller island will eventually have fewer speciesi Integrating Factor Method There is another popular method called the integrating factor method for solving rstorder linear equations writ ten in the form 1 7 ptu qti If this equation is multiplied by e fm m7 called an integrating factor show that uei Irma 4 HOW Note that the left side is a total derivative Next integrate both sides and show that you obtain 26 Use this method to solve Exercises 1 and 2 above 2 3 Approximation The fact is that most differential equations cannot be solved with simple ana lytic formulasi Therefore we are interested in developing methods to approxi mate solutions to differential equations Approximations can come in the form 23 Approximation 71 of a formula or a data set obtained by computer methods The latter forms the basis of modern scienti c computation 231 Picard Iteration We rst introduce an iterative procedure called Picard iteration Picard 18561941 that is adapted from the classical xed point method to approx imate solutions of nonlinear algebraic equations In Picard iteration we begin with an assumed rst approximation of the solution to an initial value prob lem and then calculate successively better approximations by an iterative or recursive procedure The result is a set of recursive analytic formulas that ap proximate the solution We rst review the standard xed point method for algebraic equations Example 28 Consider the problem of solving the nonlinear algebraic equation I cos 1 Graphically it is clear that there is a unique solution because the curves y z and y cosz cross at a single point Analytically we can approximate the root by making an initial guess IO and then successively calculate better approxi mations via zk1 coszk for k 012 For example if we choose 10 09 then 11 coszo cos0i9 0 622 12 cos 11 cos0i622 0 813 13 cos 12 cos0i8l3 0687 I4 coszg cos0i687 0773 I5 cos I4 cos0i773 0 716nn Thus we have generated a sequence of approximations 09 0622 0813 0687 0773 0716 If we continue the process the sequence converges to 1 0739 which is the solution to z cosz to three decimal places This method called xed point iteration can be applied to general algebraic equations of the form I The iterative procedure zk1 gzk k 012m will converge to a root 1 provided lgzl lt l and the initial guess 10 is suf ciently close to I The conditions stipulate that the graph of g is not too steep its absolute slope at the root must be bounded by one and the initial guess is close to the root 72 2 Analytic Solutions and Approximations We pick up on this iteration idea for algebraic equations to obtain an ap proximation method for solving the initial value problem IVP u 1 uh 2 10 ut0 uol First we turn this initial value problem into an equivalent integral equation by integrating the DE from to to t ut uo fsusdsl Now we de ne a type of xed point iteration called Picard iteration that is based on this integral equation formulationl We de ne the iteration scheme t uk1t u0 fsuksds k 012 211 to where uo t is an initial approximation we often take the initial approximation to be the constant function uo t uol Proceeding in this manner we generate a sequence u1t ug t ug tm of iterates called Picard iterates that under certain conditions converge to the solution of the original initial value problem 210 Example 29 Consider the linear initial value problem 1 2tl u u0 0 Then the iteration scheme is z uk1t 231 uksds k 012M 0 Take uo 0 then r u1t 231 0ds t2l 0 Then 1 t u2t 231 u1sds 231 32ds t2 gt 0 0 Next r r u3t 281u2sds uk1t 2818284d8 t2 t4 tF l 0 0 23 Approximation 73 In this manner we generate a sequence of approximations to the solution to the 1V 1n the present case one can verify that the analytic solution to the 1V is t2 ut e 7 1 The Taylor series expansion of this function is 2 1 1 1 ut e 71 t4 t5 t2 2 6 n 7 and it converges for all t Therefore the successive approximations generated by Picard iteration are the partial sums of this series and they converge to the exact solution The Picard procedure 211 is especially important from a theoretical view point The method forms the basis of an existence proof for the solution to a general nonlinear initial value problem397 the idea is to show that there is a limit to the sequence of approximations and that limit is the solution to the initial value problemi This topic is discussed in advanced texts on differential equa tions Practically however Picard iteration is not especially useful for problems in science and engineering There are other methods based upon numerical al gorithms that give highly accurate approximationsi We discuss these methods in the next section Finally we point out that Picard iteration is guaranteed to converge if the right side of the equation ftu is regular enough397 speci cally the rst partial derivatives of f must be continuous in an open rectangle of the tu plane containing the initial point However convergence is only guaranteed locally in a small interval about to EXERCISES 1 Consider the initial value problem 1 11 u0 0 Apply Picard iteration with no 0 and compute four terms If the process continues to what function will the resulting series converge 0 1 Apply Picard iteration to the initial value problem 1t7u7 11017 to obtain three Picard iterates taking no 1 Plot each iterate and the exact solution on the same set of axes 74 2 Analytic Solutions and Approximations 232 Numerical Methods As already emphasized most differential equations cannot be solved analyti cally by a simple formula In this section we develop a class of methods that solve an initial value problem numerically using a computer algorithmi In in dustry and science differential equations are almost always solved numerically because most realworld problems are too complicated to solve analyticallyi And even if the problem can be solved analytically often the solution is in the form of a complicated integral that has to be resolved by a computer calcula tion anyway So why not just begin with a computational approach in the rst place We study numerical approximations by a method belonging to a class called nite difference methods Here is the basic idea Suppose we want to solve the following initial value problem on the interval 0 g t g T u ft7u7 u0 uo 2 12 Rather than seek a continuous solution de ned at each time t we develop a strategy of discretizing the problem to determine an approximation at discrete times in the interval of interest Therefore the plan is to replace the continuous time model 212 with an approximate discrete time model that is amenable to computer solution To this end we divide the interval 0 g t g T into N segments of constant length h called the stepsizei Thus the stepsize is h TNi This de nes a set of equally spaced discrete times 0 t07t17t27m7tN T where tn nh n 0127M7Ni Now suppose we know the solution utn of the initial value problem at time tni How could we estimate the solution at time n1 Let us integrate the DE 212 from tn to tn1 and use the fundamental theorem of calculus We get the equation zn1 utn1 7 utn ft7 udti 2 13 tn The integral can be approximated using the lefthand rule giving utn1 7 utn z hftn7 If we denote by un the approximation of the solution utn at t tn then this last formula suggests the recursion formula un1 un hftn7 214 If u0 uo then 214 provides an algorithm for calculating approximations ul ug U3 etci recursively at times t1 t2 tgpn This method is called the Eu ler method named after the Swiss mathematician Ll Euler 170771783 The 23 Approximation 75 discrete approximation consisting of the values 110111 ug ug etc is called a numerical solution to the initial value problemi The discrete values approx imate the graph of the exact solution and often they are connected with line segments to obtain a continuous curve It seems evident that the smaller the stepsize h the better the approximation One can show that the cumulative error over an interval 0 g t g T is bounded by the stepsize h thus the Euler method is said to be of order hi Example 210 Consider the initial value problem 1ltu7 u0 025 Here ftu l tu and the Euler difference equation 214 with stepsize h is unl un h1 inun un hl nhun7 n 07 17 2 37 we take h 0i Beginning with uo 025 we have ul uo 0i11 00i1u0 025 0i11 0350 Then in ul 0i1110i1u1 0 35 0i11 10i10i35 0454 Next ug in 0 1120i1u2 0454 0i11 20i10i454 0563 Continuing in this manner we generate a sequence of numbers at all the discrete time points We often connect the approximations by straight line segments to generate a continuous curve In gure 23 we compare the discrete solution to the exact solution ut 61220i25 f 6 522d3 Because it is tedious to do numerical calculations by hand one can program a calculator or write a simple set of instructions for a computer algebra system to do the work for us Most calculators and computer algebra systems have builtin programs that implement the Euler algorithm automatically Below is a MATLAB m le to perform the calculations in Example 210 and plot the approximate solution on the interval 01 We take 10 steps so the stepsize is h 110 Oili 76 2 Analytic Solutions and Approximations function eulerlD Tmaxl397 N10397 hTmaxN397 u0i25397 uhistory0i25397 for nlN397 uuhlnhu397 uhistoryuhistory u397 end T0hTmaX397 plotTuhistory xlabel time t ylabel u Figure 23 The numerical solution t with a continuous curve and exact so lution in Example 210 Here h Oili A closer approximation can be obtained with a smaller stepsizei ln science and engineering we often write simple programs that implement recursive algorithms397 that way we know the skeleton of our calculations7 which is often preferred to plugging into an unknown black box containing a canned programi There is another insightful way to understand the Euler algorithm using the direction eldi Beginning at the initial value7 we take uo To nd ul the approximation at t1 we march from t07u0 along the direction eld segment with slope ft07 uo until we reach the point t17 ul on the vertical line t tli 23 Approximation 77 Then from thul we march along the direction eld segment with slope ft1u1 until we reach t27u2i From t27u2 we march along the direction eld segment with slope ft27 ug until we reach t37u3i We continue in this manner until we reach tN Ti So how do we calculate the un Inductively let us assume we are at tmun and want to calculate un1i We march along the straight line segment with slope ftnun to tn17unli Thus writing the slope of this segment in two different ways M famun tnl tn But tn1 7 tn h7 and therefore we obtain unl un hf m un7 which is again the Euler formula In summary the Euler method computes approximate values by moving in the direction of the slope eld at each point This explains why the the numerical solution in Example 210 gure 23 lags behind the increasing exact solution The Euler algorithm is the simplest method for numerically approximating the solution to a differential equation To obtain a more accurate method we can approximate the integral on the right side of 213 by the trapezoidal rule giving um 7 an gWW mummy 215 This difference equation is not as simple as it may rst appear It does not give the un1 explicitly in terms of the un because the un1 is tied up in a possibly nonlinear term on the right side Such a difference equation is called an implicit equation At each step we would have to solve a nonlinear algebraic equation for the un139 we can do this numerically which would be time consuming Does it pay off in more accuracy The answer is yes The Euler algorithm makes a cumulative error over an interval proportional to the stepsize h whereas the implicit method makes an error of order it Observe that h2 lt h when h is smalli A better approach which avoids having to solve a nonlinear algebraic equa tion at each step is to replace the un1 on the right side of 215 by the un1 calculated by the simple Euler method in 214 That is we compute a predictor 5M1 un hftn7 un7 2 16 and then use that to calculate a corrector um un hwn u flttnhan1gti 2 17 78 2 Analytic Solutions and Approximations This algorithm is an example of a predictoricorrector method and again the cumulative error is proportional to h2 an improvement to the Euler method This method is called the modi ed Euler method also Heun s method and the secondorder RungeiKutta method The Euler and modi ed Euler methods are two of many numerical con structs to solve differential equations Because solving differential equations is so important in science and engineering and because realworld models are usually quite complicated great efforts have gone into developing accurate ef cient methods The most popular algorithm is the highly accurate fourthorder HungeLKutta method Where the cumulative error over a bounded interval is proportional to h4i The RungePKutta update formula is h un1 un 061 k2 k3 k47 Where k1 ftn7un7 k2 ftn un k1gt7 ks fan un 9 k4 ftn h7un him We do not derive the formulas here but they follow by approximating the integral in 213 by Simpson s rule and then averagingi The RungeiKutta method is built in on computer algebra systems and on scienti c calculatorsi Note that the order of the error makes a big difference in the accuracy If h 01 then the cumulative errors over an interval for the Euler modi ed Euler and RungePKutta methods are proportional to 01 0 01 and 00001 respectivelyi 233 Error Analysis Readers Who want a detailed account of the errors involved in numerical al gorithms should consult a text on numerical analysis or on numerical solution of differential equations In this section we give only a brief elaboration of the comments made in the last section on the order of the error involved in Euler s methodi Consider again the initial value problem 1 ft7u7 u0 uo 2 18 23 Approximation 79 on the interval 0 g t g T with solution For our argument we assume u has a continuous second derivative on the interval which implies that the second derivative is bounded The Euler method which gives approximations un at the discrete points tn nh n 12 in N is the recursive algorithm un1 un hftn7 219 We want to calculate the error made in performing one step of the Euler al gorithmi Suppose at the point tn the approximation un is exact397 that is un Then we calculate the error at the next step Let En1 utn17un1 denote the error at the n lst stepi Evaluating the DE at t tn we get utn flttn7utn Recall from calculus Taylor s theorem with remainder that if u has two con tinuous derivatives then um um h um u tnh gitmm utn hftnutn u 7 nh2 un hftn7 un git Twyla 220 where the second derivative is evaluated at some point Tn in the interval tn7tn1i Subtracting 219 from 220 gives 1 En1 Eu 7nh2i So if un is exact the Euler algorithm makes an error proportional to h2 in computing un1i So at each step the Euler algorithm gives an error of order it This is called the local errori Notice that the absolute error is Enl u fnlh2 CVhQ7 where C is an absolute bound for the second derivative of u on the entire interval397 that is litOM S C for 0 g t g T If we apply the Euler method over an entire interval of length T where T Nh and N the number of steps then we expect to mallte a cumulative error of N times the local error or an error bounded by a constant times h This is why we say the cumulative error in Euler s method is order hi An example will con rm this calculation Consider the initial value problem for the growthidecay equation 1 ku u0 uo with exact solution ut uoekti The Euler method is un1 un hkun 1 hkuni 80 2 Analytic Solutions and Approximations We can iterate to nd the exact formula for the sequence of approximations ul 1 hku0 ug 1 hku1 1 hk2uo7 ug 1 hku2 1 hk3uo7 un 1 hknuoi One can calculate the cumulative error in applying the method over an interval 0 g t g T With T Nh Where N is the total number of steps We have EN uT iuN uolekT 7 1 hkNi The exponential term in the parentheses can be expressed in its Taylor series ekT 1 ST kT2 and the second term can be expanded using the binomial theorem 1 hkN 1 Nhk Walk my Using EN u01kTkT2717Nhk7Whk27ihkN iuoTkQ 2 h terms containing at least it So the cumulative error is the order of the stepsize hi EXERCISES 1 Use the Euler method and the modi ed Euler method to numerically solve the initial value problem 1 012511779 u0 2 on the interval 0 g t g 2 using a stepsize h 0125 Find the exact solution and compare it both graphically and in tabular form to the numerical solutions Perform the same calculation With h 011 h 01017 and h 010011 Con rm that the cumulative error at t 2 is roughly order h for the Euler method and order h2 for the modi ed Euler methodi Use the Euler method to solve the initial value problem 1 ucos t7 u0 1 on the interval 0 g t g 20 With 50 100 200 and 400 steps Compare to the exact solution and comment on the accuracy of the nu N merical algorithmi 9 A population of bacteria given in millions of organisms is governed by the law 1 0 6u 17 L u0 012 Km 7 7 Where in a periodically varying environment the carrying capacity is Kt 10 019 sint7 and time is given in days Plot the bacteria population for 40 days 23 Approximation g 01 q 00 Consider the initial value problem for the decay equation 1 77m u0 uo Here 7 is a given positive decay constant Find the exact solution to the initial value problem and the exact solution to the sequence of difference approximations un1 un 7 h39run de ned by the Euler method Does the discrete solution give a good approximation to the exact solution for all stepsizes h What are the constraints on h Suppose the temperature inside your Winter home is 68 degrees at 100 PM and your furnace then fails If the outside temperature has an hourly variation over each day given by 15 10 cos 1 degrees Where t 0 repre sents 200 PM and you notice that by 1000 PM the inside temperature is 57 degrees What Will be the temperature in your home the next morning at 600 AM Sketch a plot showing the temperature inside your home and the outside air temperature Write a program in your computer algebra system that uses the Runge Kutta method for solving the initial value problem 212 and use the program to numerically solve the problem 1 7u22t7 u0 1 Consider the initial value problem 1 511766quot u0 1 Find the exact solution and plot it on the interval 0 g t g 3 Next use the Euler method With h 01 to obtain a numerical solution Explain the results of this numerical experiment Consider the initial value problem 1 7u 15 7 11equot1u17 110l7 Where a is a parameter This model arises in the study of a chemically reacting uid passing through a continuously stirred tank reactor Where the reaction gives off heat The variable u is related to the temperature in the reactor Logan 1997 pp 4307434 Plot the solution for a 52 and for a 53 to show that the model is sensitive to small changes in the parameter a this sensitivity is called structural instability Can you explain Why this occurs Plot the bifurcation diagram With bifurcation parameter a 3 SecondOrder Differential Equations Secondorder differential equations are one of the most widely studied classes of differential equations in mathematics physical science and engineering One sure reason is that Newton s second law of motion is expressed as a law that involves acceleration of a particle which is the second derivative of position Thus onedimensional mechanical systems are governed naturally by a second order equationsi There are really two strategies in dealing with a secondorder differential equation We can always turn a single secondorder differential equation into a system of two simultaneous rstorder equations and study the systemi Or we can deal with the equation itself as it stands For example consider the damped springmass equation m1 7km 7 CIi From Section 13 we recall that this equation models the decaying oscillations of a mass m under the action of two forces a restoring force 7161 caused by the spring and a frictional force 751 caused by the damping mechanismi This equation is nothing more than a statement of Newton s second law of motion We can easily transform this equation into a system of two rstorder equations with two unknowns by selecting a second unknown state function y yt de ned by yt z t 7 thus y is the velocity Then my 7km 7 0y So the secondorder equation is equivalent to i 1797 k y 7 r e y m 84 3 SecondOrder Differential Equations This is a simultaneous system of two equations in two unknowns the position zt and the velocity yt39 both equations are rstorder Why do this Have we gained advantage ls the system easier to solve than the single equation The answers to these questions emerge as we study both types of equations in the sequeli Here we just make some general remarks that you may nd cryptici It is probably easier to nd the solution formula to the secondorder equation di rectly But the rstorder system embodies a geometrical structure that reveals the underlying dynamics in a far superior way And rstorder systems arise just as naturally as secondorder equations in many other areas of application Ultimately it comes down to one s perspective and what information one wants to get from the physical systemi Both viewpoints are important In this chapter we develop some methods for understanding and solving a single secondorder equation In Chapters 5 and 6 we examine systems of rstorder equations 31 Particle Mechanics Some secondorder differential equations can be reduced essentially to a single rstorder equation that can be handled by methods from Chapter 2 We place the discussion in the context of particle mechanics to illustrate some of the standard techniques The general form of Newton s law is m1 Ft7 171 31 where z zt is the displacement from equilibriumi a If the force does not depend on the position I then 31 is m1 Ft7 1 We can make the velocity substitution y z to obtain my F t y7 which is a rstorder differential equation that can be solved with the methods of the preceding chapteri Once the velocity y yt is found then the position zt can be recovered by antidifferentiation or zt fytdt C b If the force does not depend explicitly on time t then 31 becomes m1 Fz7 1 Again we introduce y 1quot Using the chain rule to compute the second deriva tive acceleration a M gay dt dz dt dz 31 Particle Mechanics 85 Then d my F 067 y7 which is a rstorder differential equation for the velocity y in terms of the position It If we solve this equation to obtain y yz then we can recover zt by solving the equation 1 by separation of variables c In the important special case where the force F depends only on the position I we say F is a conservative forcei See also Example 15 Then using the same calculation as in b Newton s law becomes nwj i F 17 which is a separable equation We may integrate both sides with respect to z to get my dz FzdzE dz 7 or my2 Fzdz El Note that the left side is the kinetic energy onehalf the mass times the velocity squaredl We use the symbol E for the constant of integration because it must have dimensions of energy We recall from calculus that the potential energy function Vz is de ned by idVdz Fz or the force is the negative gradient of the potential Then f Fzdz 7Vz and we have 1 Emy2 Vz E7 32 which is the energy conservation theorem the kinetic plus potential energy for a conservative system is constant The constant E which represents the total energy in the system can be computed from knowledge of the initial position 10 IO and initial velocity y0 yo or E yg Vzol We regard the conservation of energy law as a reduction of Newton s second law 7 the latter is a secondorder equation whereas 32 is a rstorder equation if we replace the position y by dzdti It may be recast into 1 igJEilz 33 This equation is separable and its solution would give I The appropri ate sign is taken depending upon whether the velocity is positive or negative during a certain phase of the motion Usually we analyze conservative systems in phase space myspace or the phase plane by plotting y vs I from equation 32 for different values of the parameter El The result is a oneparameter family of curves or orbits in 86 3 SecondOrder Differential Equations the my plane along which the motion occurs The set of these curves forms the phase diagram for the systeml On these orbits we do not know how I and y depend upon time t unless we solve 33 But we do know how velocity relates to position Example 31 Consider a springmass system without dampingl The governing equation is m1 lm where k is the spring constant The force is 7km and the potential energy Vz is given by Vz 77kzdz 312 We have picked the constant of integration to be zero which automatically sets the zero level of potential energy at z 0 ie V0 0 Conservation of energy is expressed by 32 or my2 12 E7 which plots as a family of concentric ellipses in the my phase plane one ellipse for each value of El These curves represent oscillations and the mass tracks on one of these orbits in the phase plane continually cycling as time passes397 the position and velocity cycle back and forth At this point we could attempt to solve 33 to determine how I varies in time but in the next section we nd an easier method to solve secondorder linear equations for zt directlyl EXERCISES 1 Consider a dynamical system governed by the equation I 71 13 Hence m 1 Find the potential energy Vz with V0 0 How much total energy E is in the system if 10 2 and zO 1 Plot the orbit in the my phase plane of a particle having this amount of total energy 2 In a conservative system show that the conservation of energy law 32 can be obtained by multiplying the governing equation m1 Fz by z and noting that 21 CAD i In a conservative system derive the relation d t i I C 2E VI which gives time as an antiderivative of an expression that is a function of position 32 Linear Equations with Constant Coefficients 87 4 A bullet is shot from a gun with muzzle velocity 700 meters per second horizontally at a point target 100 meters away Neglecting air resistance by how much does the bullet miss its target 01 i Solve the following differential equations by reducing them to rstorder equations a I i zi b I IIquot c I 741 d I I e tz z 4t 533 In a nonlinear springmass system the equation governing displacement is I 7213 Show that conservation of energy for the system can be expressed as if C 7 14 where C is a constant Plot this set of orbits in the phase plane for different values of C If 10 10 gt 0 and zO 07 show that the period of oscillations is T 7 4 1 dr 10 0 1 Sketch a graph of the period T vs 10 Hint in 33 separate variables and integrate over onefourth of a period 32 Linear Equations with Constant Coef cients We recall two models from Chapter 1 For a springmass system with damping the displacement zt satis es mm 01 k1 0 The current t in an RCL circuit with no emf satis es LI R1 I 0 The similarity between these two models is called the mechanicalelectrical analogy The spring constant k is analogous to the inverse capacitance lC397 both a spring and a capacitor store energy The damping constant c is analogous to the resistance R both friction in a mechanical system and a resistor in an electrical system dissipate energy The mass m is analogous to the inductance 88 3 SecondOrder Differential Equations L397 both represent inertia in the systemi All of the equations we examine in the next few sections can be regarded as either circuit equations or mechanical problems After dividing by the leading coef cient both equations above have the form u pu qu 07 34 where p and q are constants An equation of the form 34 is called a second order linear equation with constant coef cients Because zero is on the right side physically there is no external force or emf the equation is homo geneous Often the equation is accompanied by initial data of the form u0 A7 uO Bi 35 The problem of solving 34 subject to 35 is called the initial value prob lem IVP Here the initial conditions are given at t 0 but they could be given at any time t tor Fundamental to our discussion is the following existenceuniqueness theorem which we assume to be true It is proved in ad vanced textsi Theorem 32 The initial value problem 343 5 has a unique solution that exists on 700 lt t lt 00 The plan is this We rst note that the DE 3 4 always has two independent solutions u1t and ug t by independent we mean one is not a multiple of the other We prove this fact by actually exhibiting the solutions explicitlyi If we multiply each by an arbitrary constant and form the combination ut clul t Cgu1t where cl and Cg are the arbitrary constants then we can easily check that ut is also a solution to 34 This combination is called the general solution to 34 We prove at the end of this section that all solutions to 34 are contained in this combination Finally to solve the initial value problem we use the initial conditions 35 to uniquely determine the constants cl and 52 We try a solution to 34 of the form u e 7 where A is to be determined We suspect something like this might work because every term in 34 has to be the same type of function in order for cancellation to occur397 thus u u and u must be the same form which suggests an exponential for u Substitution of u 6 into 34 instantly leads to A2pAq 0 3 6 32 Linear Equations with Constant Coefficients 89 which is a quadratic equation for the unknown A Equation 36 is called the Characteristic equation Solving we obtain roots 1 A 5PM P2 7 44 These roots of the characteristic equation are called the Characteristic values or roots corresponding to the differential equation 34 There are three cases depending upon whether the discriminant p2 7 4g is positive zero or negative The reader should memorize these three cases and the forms of the solution Case 1 If p2 7 4g gt 0 then there are two real unequal characteristic values A1 and A2 Hence there are two independent exponentialtype solutions MO 6M7 1120 6M7 and the general solution to 34 is ut 016A ageMt 37 Case 2 If p2 7 4g 0 then there is a double root A 7102 Then one solution is ul e A second independent solution in this case is ug te Therefore the general solution to 34 in this case is ut 016 CgteM 38 Case 3 If p2 7 44 lt 0 then the roots of the characteristic equation are complex conjugates having the form A a i Therefore two complex solutions of 34 are 40411917 40443 To manufacture real solutions we use a fundamental result that holds for all linear homogeneous equations Theorem 33 If u gt iht is a complex solution to the differential equation 34 then its real and imaginary parts gt and ht7 are real solutions 9O 3 SecondOrder Differential Equations The simple proof is requested in the Exercises Let us take the rst of the complex solutions given above and expand it into its real and imaginary parts using Euler s formula elm cos t isin ti We have ew ememt emcos t i sin t em cos t 23960 sin tl Therefore by Theorem 33 ul em cos t and ug em sin t are two real independent solutions to equation 34 If we tallte the second of the complex solutions 6a mt instead of e mt then we get the same two real solutions Consequently in the case that the characteristic values are complex A ozii 7 the general solution to DE 34 is ut Clem cos t 0260 sin t 39 In the case of complex eigenvalues we recall from trigonometry that 39 can be written differently as ut emcl cos t Cg sin t emAcos t 7 go where A is called the amplitude and go is the phase This latter form is called the phase7amplitude form of the general solution Written in this form A and go play the role of the two arbitrary constants instead of cl and 02 One can show that that all these constants are related by 52 A Mc 537 g0 arctan l 51 This is because the cosine of difference expands to Acos t 7 go Acos t cos so Asin t sin gal Comparing this expression to cl cos t Cg sin t gives Acosgo cl Asingo 02 Squaring and adding this last set of equations determines A and dividing the set of equations determines gal Observe that the solution in the complex case is oscillatory in nature with em multiplying the amplitude A If a lt 0 then the solution will be a decaying oscillation and if a gt 0 the solution will be a growing oscillationl If a 0 then the solution is ut cl cos t Cg sin t Acos t 7 go and it oscillates with constant amplitude A and period 27r l The frequency is called the natural frequency of the systemi 32 Linear Equations with Constant Coefficients 91 There is some useful terminology used in engineering to describe the motion of a springmass system with damping governed by the equation m1 cz k1 0 The characteristic equation is mA2 CA k 07 with roots A 7 7ci vs 7 472176 7 2m If the roots are complex 52 lt 4mk then the system is underdamped rep resenting a decaying oscillation if the roots are real and equal 52 4mk then the system is critically damped decay no oscillations and at most one pass through equilibrium 1 0 7 if the roots are real and distinct 52 gt 472116 then the system is overdamped a strong decay toward 1 0 The same terminology can be applied to an RCL circuit Example 34 The differential equation u 7 u 7 12u 0 has characteristic equation A2 7 A 712 0 with roots A 737 4 These are real and distinct and so the general solution to the DE is u 016 l 0264 Over a long time the contribution 6 3 decays and the solution is dominated by the 6 term Thus eventually the solution grows exponentially Example 35 The differential equation u 41 4u 0 has characteristic equation A2 4A l 4 07 with roots A 72 72 Thus the eigenvalues are real and equal and the general solution is u 016 2 Cgte h This solution decays as time gets large recall that a decaying exponential dominates the linear growth term t so that te Q goes to zero Example 36 The differential equation u 2u 2u 0 models a damped springmass system with m 1 c 2 and k 2 It has characteristic equation A2 2A 1 0 The quadratic formula gives complex roots A 71 i 2239 Therefore the general solution is u olequot cos 2t l 026quot sin 2t7 92 3 SecondOrder Differential Equations representing a decaying oscillation Here the natural frequency of the un damped oscillation is 2 In phase7amplitude form we can write u Aequotcos2t 7 go Let us assume that the mass is given an initial velocity of 3 from an initial position of 1 Then the initial conditions are u0 1 u 0 3 We can use these conditions directly to determine either cl and 02 in the rst form of the solution or A and go in the phaseamplitude formi Going the latter route we apply the rst condition to get u0 A670 cos20 7 go Acosgo 1 To apply the other initial condition we need the derivative We get u 72Aequot sin2t 7 go Then uO 7214670 sin20 7 go 2Asingo 3 Therefore we have 3 Acosgo 1 Asingo Squaring both equations and summing gives A2 1347 so the amplitude is A 134i Note that the cosine is positive and the sine is positive so the phase angle lies in the rst quadrant The phase is 3 go arctan i 0983 radiansi Therefore the solution to the initial value problem is 13 u 16420th 7 01983 This solution represents a decaying oscillationi The oscillatory part has natural frequency 2 and the period is 7L See gure 31 The phase has the effect of translating the cos 2t term by 09832 which is called the phase shift To summarize we have observed that the differential equation 34 always has two independent solutions ul t and 11207 and that the combination ut clu1t Cgul t is also a solution called the general solution Now as promised we show that the general solution contains all possible solutions to 314 To see this let ul t and ug t be special solutions that satisfy the initial conditions mm 717 u3lt0gt 70 3 2 Linear Equations with Constant Coefficients Figure 31 Plot of the solution and 11203 07 20 17 respectively Theorem 32 implies these two solutions exist Now let vt be any solution of 34 It will satisfy some conditions at t 0 say 110 a and 10 12 But the function ut au1t bu1t satis es those same initial conditions u0 a and uO 12 Must ut therefore equal vt Yes by the uniqueness theorem Theorem 32 Therefore vt au1tbu1t and the solution vt is contained in the general solution Two equations occur so frequently that it is worthwhile to memorize them along with their solutions The pure oscillatory equation u k2u 0 has characteristic roots A iki and the general solution is u cl cos kt 62 sin ktl On the other hand the equation 17 k2u 0 94 3 SecondOrder Differential Equations has characteristic roots A ik and thus the general solution is u Clem 0267163 This latter equation can also be Written in terms of the the hyperbolic functions cosh and sinh as u Cl cosh kt 02 sinh kt7 Where em 67m kt 7 67kt T7 cosh kt sinh kt 6 Sometimes the hyperbolic form of the general solution is easier to work With EXERCISES H Find the general solution of the following equations a u 7 41 4u 0 b u u 4u 0 c u 7 51 6u 0 d u 9u 0 e u 7 21 0 f u 7 12u 0i 2 Find the solution to the initial value problem u u u 07 u0 u 0 l7 and Write it in phaseamplitude formi 9 A damped springmass system is modeled by the initial value problem u 0i125u u 07 u0 27 uO 0 Find the solution and sketch its graph over the time interval 0 g t g 50 If the solution is Written in the form ut Aequot15 coswt 7 go nd A7 w and gal 4 For Which values of the parameters a and b if any Will the solutions to u 7 2au bu 0 oscillate With no decay iiei7 be periodic Oscillate With decay Decay Without oscillations 5 An RCL circuit has equation LI I I 0i Characterize the types of current responses that are possible depending upon the value of the inductance Li 533 An oscillator With damping is governed by the equation z 3az bz 07 Where a and b are positive parameters Plot the set of points in the ab plane iiei ab parameter space Where the system is critically dampedi 539quot t The Nonhomogeneous Equation 95 7 Find a DE that has general solution ut 0164 Cge B i 8 Find a DE that has solution ut 6 2te 3 i What are the initial conditions 9 Find a DE that has solution ut sin 4t 3 cos 4t 10 Find a DE that has general solution ut A cosh 5t B sinh 5t7 where A and B are arbitrary constants Find the arbitrary constants when u0 2 and u 0 0 11 Find a DE that has solution ut e 2tsin 475 3cos 4t What are the initial conditions 12 Describe the current response Hit of a LC circuit with L 5 henrys C 2 farads with 0 0 IO 1 13 Prove Theorem 33 by substituting u into the equation and separating real and imaginary parts using linearity Then use the fact that a complex quantity is zero if and only if its real and imaginary parts are zero 33 The Nonhomogeneous Equation In the last section we solved the homogeneous equation u pu qu 0 310 Now we consider the nonhomogeneous equation u 0u qu 9t7 311 where a known term gt called a source term or forcing term is included on the right side In mechanics it represents an applied timedependent force397 in a circuit it represents an applied voltage an emf such as a battery or generatori There is a general structure theorem analogous to Theorem 12 for rstorder linear equations that dictates the form of the solution to the nonhomogeneous equationi Theorem 37 All solutions ofthe nonhomogeneous equation 311 are given by the sum of the general solution to the homogeneous equation 3 10 and any particular solution to the nonhomogeneous equationi That is the general solution to 311 is 1W Clu1t C2u1t upt7 96 3 SecondOrder Differential Equations where ul and ug are independent solutions to 310 and up is any solution to 311 This result is very easy to showi lf ut is any solution whatsoever of 311 and up t is a particular solution then ut7up t must satisfy the homogeneous equation 310 Therefore by the results in the last section we must have ut 7 upt clu1t Cgu1ti 331 Undetermined Coe icients We already know how to nd the solution of the homogeneous equation so we need techniques to nd a particular solution up to 311 One method that works for many equations is to simply make a judicious guess depending upon the form of the source term gtl Of cially this method is called the method of undetermined coef cients because we eventually have to nd numerical coef cients in our guessi This works because all the terms on the left side of 311 must add up to give gtl So the particular solution cannot be too wild if gt is not too wild397 in fact it nearly must have the same form as gtl The method is successful for forcing terms that are exponential functions sines and cosines polynomials and sums and products of these common functionsi Here are some basic rules without some caveats which come later The capital letters in the list below denote known constants in the source term gt and the lowercase letters denote coef cients to be determined in the trial form of the particular solution when it is substituted into the differential equation 1 If gt A67 is an exponential then we try an exponential up ae Y i 2i lfgt Asinwt or gt Acoswt then we try a combination u Z7 a sinwt 12 coswti 9 lf gt Ant An71tn 1 A0 is a polynomial of degree n then we try up ant an71tn 1 a0 a polynomial of degree ml 4 If gt Antn 147L717 l A1t Awe 7 then we try up antn an71tn 1 a1t a0e397 l 5 If gt Ae Y sinwt or gt Ae Y cos wt then we try up ae lt sinwt be coswti If the source term gt is a sum of two different types we take the net guess to be a sum of the two individual guessesi For example ifgt Ste 1 76 2 7 a polynomial plus an exponential then a good guess would be up atbce 2ti The following examples show how the method works 33 The Nonhomogeneous Equation 97 Example 38 Find a particular solution to the differential equation uH7u7u5t73 The right side gt 5t 7 3 is a polynomial of degree 1 so we try up at 12 Substituting 7a 7at b 5t 7 3 Equating like terms constant term and terms involving t gives 7a 7b 73 and 70 5 Therefore a 57 and b 71649 A particular solution to the equation is therefore 5 16 upt 7 7 Et Example 39 Consider the equation u 314 Su 6672 The homogeneous equation has characteristic polynomial A23A3 07 Which has roots A 7 i Thus the solution to the homogeneous equation is 7312 JP 7312 uht cle cos 7t 626 s1n 7t To nd a particular solution to the nonhomogeneous equation note that gt 66 2 Therefore we guess up ae m Substituting this trial function into the nonhomogeneous equation gives after canceling 6 2 the equation 4a 7 6a 30 6 Thus a 1 and a particular solution to the nonhomogeneous equation is up 6 2 The general solution to the original nonhomogeneous equation is ut claim2 cos gt waist2 sin gt 672 Example 310 Find a particular solution to the DE u 2u sin 3t Our basic rule above dictates we try a solution ofthe form up a sin 3tb cos 3t Then upon substituting 79a sin 3t 7 91 cos 3t 2a sin 3t 21 cos 3t sin 3t Equating like terms gives 79a 2a 1 and b 0 there are no cosine terms on the right side Hence a 717 and a particular solution is up 7 sin 3t 98 3 SecondOrder Differential Equations For this equation because there is no rst derivative we did not need a cosine term in the guess If there were a rst derivative a cosine would have been required Example 311 Next we modify Example 310 and consider u 9u sin St The rules dictate the trial function up asin 3t bcos 3t Substituting into 315 yields 79a sin 3t 7 91 cos 3t 9a sin 3t 91 cos 3t sin 3t But the terms on the left cancel completely and we get 0 sin 3t an absurdityl The method failedl This is because the homogeneous equation u 9u 0 has eigenvalues A iSi which lead to independent solutions ul sin 3t and ug cos St The forcing term gt sin St is not independent from those two basic solutions397 it duplicates one of them and in this case the method as presented above fails The fact that we get 0 when we substitute our trial function into the equation is no surprise7it is a solution to the homogeneous equationi To remedy this problem we can modify our original guess by multiplying it by t That is we attempt a particular solution of the form u Z7 ta sin3t bcos 3t Calculating the second derivative ug and substituting along with up into the original equation leads to show thisl 6a cos 3t 7 6 sin 3t sin 3t Hence a 0 and b 716 We have found a particular solution u 7tcos3ti 17 Therefore the general solution of the original nonhomogeneous equation is the homogeneous solution plus the particular solution Wt Ci COS 375 62 sin 3t 7 t cos 3t Notice that the solution to the homogeneous equation is oscillatory and re mains bounded the particular solution oscillates without bound because of the increasing time factor t multiplying that terml 33 The Nonhomogeneous Equation 99 The technique for nding the form of the particular solution that we used in the preceding example works in general397 this is the main caveat in the set of rules listed above Caveati If a term in the initial t7ial guess for a particular solution up dupli cates one of the basic solutions for the homogeneous equation then modify the guess by multiplying by the smallest power of t that eliminates the duplication Example 312 Consider the DE u 7 4u u 5te2 i The initial guess for a particular solution is up at be2 i But7 as you can check e2t and te2t are basic solutions to the homogeneous equation u 7 4u u 0i Multiplying the rst guess by t gives up at2 bte2 7 which still does not eliminate the duplication because of the te2t termi So multiply by another t to get up at3 bt2e2 i Now no term in the guess duplicates one of the basic homogeneous solutions and so this is the correct form of the particular solution If desired7 we can substitute this form into the differential equation to determine the exact values of the coef cients a and 12 But without actually nding the coef cients7 the form of the general solution is ut cle2t Cgte2t at3 bt2e2 i The constants cl and Cg could be determined at this point by initial conditions if given Sometimes knowing the form of the solution is enough Example 313 Consider an RCL circuit where R 2 L C l and the current is driven by an electromotive force of 2sin St The circuit equation for the voltage Vt across the capacitor is V 2V V 2sin3ti For initial data we take V0 47 V 0 0 We recognize this as a nonhomogeneous linear equation with constant coef cients So the general solution will be the sum of the general solution to the homogeneous equation V 2V V0 100 3 SecondOrder Differential Equations plus any particular solution to the nonhomogeneous equation The homoge neous equation has characteristic equation A2 2 1 0 with a double root A 71 Thus the homogeneous solution is Vh 6401 027i Notice that this solution regardless of the values of the constants will decay away in time397 this part of the solution is called the transient response of the circuit To nd a particular solution we use undetermined coef cients and assume it has the form Vp asin3t bcos Stl Substituting this into the nonhomogeneous equation gives a pair of linear equa tions for a and b 74a73b17 7a79b0i We nd a 70158 and b 70123 Therefore the general solution is Vt 7 e cl 7 at 7 0 158sin 3t 7 0123 cos 3t Now we apply the initial conditions Easily V0 4 implies cl 4123 Next we nd V t so that we can apply the condition V 0 0 Leaving this as an exercise we nd Cg 4597 Therefore the voltage on the capacitor is Vt 7 e 4i123 45972 7 0158 sin 3t 7 0123 cos 3th As we observed the rst term always decays as time increases Therefore we are left with only the particular solution 70 158 sin 3t 7 0123 cos 3t which takes over in time It is called the steadystate response of the circuit gure 32 The method of undetermined coef cients works for nonhomogeneous rst order linear equations as well provided the equation has constant coef cientsi Example 314 Consider the equation 1 qu gti The homogeneous solution is uht Ce q i Provided gt has the right form a particular solution up t can be found by the method of undetermined coef cients exactly as for secondorder equations make a trial guess and substitute into the equation to determine the coef cients in the guess The general solu tion to the nonhomogeneous equation is then ut uh tup For example consider the equation 1 7 Su t 7 2i 33 The Nonhomogeneous Equation 101 2 4 8 8 10 12 2 Figure 32 A plot of the voltage Vt in Example 313 Initially there is a transient caused by the initial conditions It decays away and is replaced by a steadystate response7 an oscillation7 that is caused by the forcing terml The homogeneous solution is uh C633 To nd a particular solution make the trial guess up at 12 Substituting this into the equation gives a 7 and b Consequently the general solution is ut Ce 7 t EXERCISES 1 Each of the following functions represents gt the right side of a nonho mogeneous equation State the form of an initial trial guess for a particular solution u7 a 373 7 1 b 12 c t263 l d 5sin7tl e 62 cos t t2l 102 3 SecondOrder Differential Equations 9 we 5 0 f te t sin Trti Find the general solution of the following nonhomogeneous equations u 7u tegtl b u 7u 662 i c u u 7 d u 7 317 4u 2t2i u u 964 u u 4673 u 7 4u cos 2ti u u 2u tsin2t Solve the initial value problem u 731 740u 2equot7 u0 07 u 0 1 Find the solution of u 7 21 47 u0 l7 u 0 0i i Find the particular solution to the equation u u 2u sin2 t Hint use a double angle formula to rewrite the right side An RL circuit contains a 2 ohm resistor and a 5 henrys inductor connected in series With a 10 volt battery If the open circuit is suddenly closed at time zero nd the current for all times t gt 0 Plot the current vs time and identify the steadystate response A circuit contains a 10 3 farad capacitor in series With a 20 volt battery and an inductor of 04 henrysi At t 0 both L1 0 and I 0 Find the charge qt on the capacitor and describe the response of the circuit in terms of transients and steadystatesi An RCL circuit contains a battery generating 110 volts The resistance is 16 ohms the inductance is 2 henrys and the capacitance is 002 faradsi lf q0 5 and 10 0 nd the charge qt current response of the circuit Identify the transient solution and the steadystate responsei 332 Resonance The phenomenon of resonance is a key characteristic of oscillating systemsi Resonance occurs When the frequency of a forcing term has the same frequency as the natural oscillations in the system397 resonance gives rise to large amplitude 33 The Nonhomogeneous Equation 103 oscillations To give an example consider a pendulum that is oscillating at its natural frequency What happens when we deliberately force the pendulum say by giving it a tap with our nger in the positive angular direction at a frequency near this natural frequency So every time the bob passes through 9 0 with a positive direction we give it a positive tapi We will clearly increase its amplitude This is the phenomenon of resonance It can occur in circuits where we force by a generator the system at its natural frequency and it can occur in mechanical systems and structures where an external periodic force is applied at the same frequency as the system would naturally oscillatei The results could be disastrous such as a blown circuit or a fallen building397 a building or bridge could have a natural frequency of oscillation and the wind could provide the forcing functioni Another imagined example is a company of soldiers marching in cadence across a suspension bridge at the same frequency as the natural frequency of the structure We consider a model problem illustrating this phenomenon an LC circuit that is forced with a sinusoidal voltage source of frequency 1f L 1 the governing equation for the charge on the capacitor will have the form u w2u sin t7 312 where 0J2 1Ci Assume rst that B a w and take initial conditions u0 07 uO 1 The homogeneous equation has general solution uh cl cos wt 52 sinwt7 which gives natural oscillations of frequency Am A particular solution has the form up asin ti Substituting into the DE gives a 1w2 7 So the general solution of 312 is 1 u clcoswt02 s1nwt msln ti 313 At t 0 we have u 0 and so cl 0 Also uO 1 gives 02 7 l Therefore the solution to the initial value problem is B ww2 7 B2 u 7Ws1nwt wQ sln t 3 14 1 7 52 This solution shows that the charge response is a sum of two oscillations of different frequencies 1f the forcing frequency is close to the natural frequency 0 then the amplitude is bounded but it is obviously large because of the factor 0J2 7 B2 in the denominator Thus the system has large oscillations when is close to am 104 3 SecondOrder Differential Equations What happens if B w Then the general solution in 313 is not valid because there is division by zero and we have to resolve the problemi The circuit equation is u w2u sinwt7 315 where the circuit is forced at the same frequency as its natural frequency The homogeneous solution is the same as before but the particular solution will now have the form up tasinwt bcoswt7 with a factor oft multiplying the terms Therefore the general solution of 315 has the form ut cl cos wt Cg sinwt ta sinwt b coswti Without actually determining the constants we can see the nature of the re sponse Because of the t factor in the particular solution the amplitude of the oscillatory response ut will grow in time This is the phenomenon of pure resonance It occurs when the frequency of the external force is the same as the natural frequency of the systemi What happens if we include damping in the circuit ie a resistor and still force it at the natural frequency Consider u 201 2u sin 2t7 where 20 is a small 0 lt a damping coef cient for example resistance The ho mogeneous equation u 2au 2u 0 has solution u 6quot cl cos 2 7 02t 02 sin 2 7 02ti Now the particular solution has the form up acos 2t 1 sin x2t7 where a and b are constants found by substituting into the DE So the response of the circuit is u e ncl cos 2 7 a2t Cg sin 2 7 02t acos 2t 1 sin 2ti The transient is a decaying oscillation of frequency V2 7 02 and the steady state response is periodic of frequency The solution will remain bounded but its amplitude will be large if a is very small EXERCISES H Graph the solution 314 for several different values of B and mi lnclude values where these two frequencies are closer E0 Find the general solution of the equation u 16u cos 4t 9 Consider a general LC circuit with input voltage Vosin t If B and the capacitance C are known what value of the inductance L would cause resonance 34 Variable Coefficients 105 4 Consider the equation u wQu cos t a Find the solution when the initial conditions are u0 u 0 0 when w a B b Use the trigonometric identity 2 sin A sin B cosA 7 B 7 cosA B to write the solution as a product of sinesi c Take w 55 and 45 and plot the solution in part d Show that the solution in c can be interpreted as a highfrequency response contained in a lowfrequency amplitude envelopei We say the high frequency is modulated by the low frequency This is the phenomenon of beatsi 34 Variable Coef cients Next we consider secondorder linear equations with given variable coef cients W and N u ptu qtu gti 3 16 Except for a few cases these equations cannot be solved in analytic form using familiar functionsi Even the simplest equation of this form u 7tu0 where pt gt 0 and qt it which is called Airy s equation re quires the de nition of a new class of functions Airy functions to characterize the solutions Nevertheless there is a welldeveloped theory for these equations and we list some of the main results We require that the coef cients pt and qt as well as the forcing term gt be continuous functions on the interval I of interest We list some basic properties of these equations397 the reader will observe that these are the same properties shared by secondorder constant coef cient equations studied in Section 32 1 ExistenceUniqueness If I is an open interval and to belongs to I then the initial value problem u 0tu4tu 9W7 317 we a u t0b 318 has a unique solution on L 106 3 SecondOrder Differential Equations 2 Superposition of Solutions lf ul and ug are independent solutions of the associated homogeneous equation u ptu qtu 0 319 on an interval 1 then ut clul Cgug is a solution on the interval I for any constants cl and 02 Moreover all solutions of the homogeneous equation are contained in the general solution 3 Nonhomogeneous Equation All solutions to the nonhomogeneous equation 31 can be represented as the sum of the general solution to the homogeneous equation 3 19 and any particular solution to the non homogeneous equation 317 In symbols 71 Clul WWW 112407 which is called the general solution to 317 The dif culty of course is to nd two independent solutions ul and u to the homogeneous equation and to nd a particular solution As we remarked this task is dif cult for equations with variable coef cients The method of writing down the characteristic polynomial as we did for constant coef cient equations does not work 341 CauchyiEuler Equation One equation that can be solved analytically is an equation of the form I c u gu Jr 771107 t2u btu cu 07 which is called a CauchyiEuler equation In each term the exponent on t coincides with the order of the derivative Observe that we must avoid t 0 in our interval of solution because pt bt and qt ct2 are not continuous at t 0 We try to nd a solution of the form of a power function u tml Think about why this might work Substituting gives the Characteristic equation 72172171bmc07 which is a quadratic equation for ml There are three cases If there are two distinct real roots m1 and m2 then we obtain two independent solutions tm1 and tmil Therefore the general solution is u cltm1 Cgtmzl 34 Variable Coefficients 107 If the characteristic equation has two equal roots m1 m2 m then tm and tm lnt are two independent solutions397 in this case the general solution is u cltm Cgtm ln it When the characteristic equation has complex conjugate roots m a i i6 we note using the properties of logarithms exponentials and Euler s formula that a complex solution is W tat tat fem ew t cos ln t isin 111m The real and imaginary parts of this complex function are therefore real solu tions Theorem 33 So the general solution in the complex case is u cltu cos lnt Cgtu sin ln 7 Figure 33 shows a graph of the function sin5 ln t7 which is a function of the type that appears in this solution Note that this function oscillates less and less as t gets large because lnt grows very slowly As t nears zero it oscillates in nitely many times Because of the scale these oscillations are not apparent on the plot Figure 33 Plot of sin5 hit Example 315 Consider the equation tQuN tu 9u 0i 108 3 SecondOrder Differential Equations The characteristic equation is mm71m9 07 Which has roots m iSii The general solution is therefore u cl cos3 ln t Cg sinSlnti Example 316 Consider the equation um Eu t i We can Write this in Cauchy7Euler form as t2 7 2m 0 Which has characteristic equation mm 7 l 7 2m 0 The roots are m 0 and m 3 Therefore the general solution is ut cl Cgtgl Example 317 Solve the initial value problem t2u 3m u 0 u1 07 u 1 2 The DE is Cauchy7Euler type With characteristic equation mm7 l3ml 0 This has a double root m 71 and so the general solution is ut 1 Cflnti Now ul cl 0 and so ut CTZlnti Taking the derivative u t l 7 mt Then ul Cg 2 Hence the solution to the initial value problem is 2 ut ln t A Cauchy 178971857 and Li Euler 170771783 were great mathematicians Who left an indelible mark on the history of mathematics and science Their names are encountered often in advanced course in mathematics and engineer ing 34 Variable Coefficients 109 342 Power Series Solutions In general how are we to solve variable coef cient equations Some equations can be transformed into the CauchyiEuler equation but that is only a small class If we enter the equation in a computer algebra system such as Maple or Mathematica the system will often return a general solution that is expressed in terms of socalled special functions such as Bessel functions Airy functions and so on We could de ne these special functions by the differential equations that we cannot solve This is much like de ning the natural logarithm function lnt as the solution to the initial value problem 1 ul 07 as in Chapter 1 For example we could de ne functions Ait and Bit the Airy functions as two independent solutions of the DE u 7 tu 0 Many of the properties of these special functions could be derived directly from the differential equation itself But how could we get a formula for those functions One way to get a representation of solutions to equations with variable coef cients is to use power series Let p and L be continuous on an open interval I containing to and also have continuous derivatives of all orders on If Solutions to the secondorder equation with variable coef cients u ptu qtu 07 320 can be approximated near t to by assuming a power series solution of the form ut Zan ito a0 a1t7 to a2t7 t02 a3t 7 to3 i n0 The idea is to simply substitute the series and its derivatives into the differential equation and collect like terms thereby determining the coef cients ani We recall from calculus that a power series converges only at t to for all t or in an interval to 7 R to R where R is the radius of convergence Within its radius of convergence the power series represents a function and the power series may be differentiated term by term to obtain derivatives of the function Example 318 Consider the DE ufi ltu 0 on an interval containing to 0 We have ut a0 a1t a2t2 agts a4t4 x x ut a1 2a2t 3a3t2 4a4t3 r u t 2 6a3t12a4t2 n 110 3 SecondOrder Differential Equations Substituting into the differential equation gives 2a 6a3t 12a4t2 r 7 1 ta0 a1t a2t2 agt3 3 0 Collecting like terms 7 2a2 7 7 a1 6a3t 7a 7 a1 12a4t2 3 0 Therefore 7 2 07 iao a1 6amp3 07 7a27a112a4 07 Notice that all the coef cients can be determined in terms of a0 and all We have 1 1 1 1 1 a2 5am a3 ado a17 a4 lta1 a2 Wi iaob Therefore the power series for the solution ut can be written 1 1 1 1 711 a0 ait t 5amp0 Elta0 a113 5a1 ant4 H 7 12 1 3 1 4 1 3 1 4 7 ao12t 675 2475 39a1t6t 1275 1 7 which gives the general solution as a linear combination of two independent power series solutions 7 1 2 1 3 1 4 u1t 7 12t 6t 24t 1 3 1 4 t t t t r u2 6 12 The two coef cients a0 and a1 can be determined from initial conditions For example if u0 17 710 37 then a0 1 and a1 37 which gives the power series solution 7 12 13 14 13 14 ut 7 12t6t24t 3t6t12t 7 7 12 23 4 713t2t3t24tuu In this example the power series converges for all t We have only calculated ve terms and our truncated power series is an approximation to the actual solution to the initial value problem in a neighborhood of t 0 Figure 314 shows the polynomial approximations by taking the rst term the rst two terms the rst three and so on 34 Variable Coefficients 111 four polynomial approximations Figure 34 Successive polynomial approximations There are many important equations of the form 3120 where the coef cients p and 4 do not satisfy the regularity properties having continuous derivatives of all order mentioned at the beginning of this subsection How ever if p and q are not too illbehaved at to we can still seek a series solution In particular if t 7 t0pt and t 7 t02qt have convergent power series ex pansions in an interval about to then we say to is a regular singular point for 3120 and we attempt a series solution of the form ut tr Zan i to where 7 is some number Substitution of this form into 3120 leads to equations for both 7 and the coef cients anl This technique which is called the Frobenius method is explored in the Exercises 343 Reduction of Order If one solution u1t of the DE u ptu qtu 0 happens to be known then a second linearly independent solution ug t can be found of the form ug t vtu1t for some vt to be determined To nd 112 3 SecondOrder Differential Equations vt we substitute this form for ug t into the differential equation to obtain a rstorder equation for vtl This method is called reduction of order and we illustrate it with an example Example 319 Consider the DE 1 l u 7 u t Qu 0 An obvious solution is u1t ti So let ug vttl Substituting we get 1 1 21 tv 7 v tv t 2vt 07 which simpli es to tv v 0i Letting w 17 we get the rstorder equation tww0l By separating variables and integrating we get w ltl Hence 1 f dt ln t and the second independent solution is ug t tln ti Consequently the general solution of the equation is ut clt Cgt ln ti Note that this example is a CauchyiEuler equation397 but the method works on general linear secondorder equationsl 344 Variation of Parameters There is a general formula for the particular solution to the nonhomogeneous equation u pW u9 620 called the variation of parameters formulal Recall how we attacked the rstorder linear equation wp uy in Chapter 1 We rst found the solution of the associated homogeneous equa tion uptu 07 as C6430 where Pt is an antiderivative of ptl Then we found the solution to the nonhomogeneous equation by varying the constant 34 Variable Coefficients 113 C ie by assuming ut Cte P l Substitution of this into the equation yielded Ct and therefore the solution The same method works for secondorder equations but the calculations are more involved Let ul and ug be independent solutions to the homogeneous equation u ptu qtu 0 Then uh 0111105 0211205 is the general solution of the homogeneous equationi To nd a particular solu tion we vary both parameters cl and Cg and take up t 01tu1t Cg tu2 322 Now we substitute this expression into the nonhomogeneous equation to get expressions for 01t and 52ti This is a tedious task in calculus and algebra and we leave most of the details to the interested readeri But here is how the argument goes We calculate u equationi For notational simplicity we drop the t variable in all of the functions We have and ug so that we can substitute into the u 5114 0212 c lul cgugi There is exibility in our answer so let us set clul 6ng 0 323 Then u 0111 02 0114 021 511qu Substituting these into the nonhomogeneous DE gives 0114 0214 Gail1 6ng ptclu1 0212 qtclu1 02112 gtl Now we observe that ul and ug satisfy the homogeneous equation and this simpli es the last equation to 61114quot cQu2 gtl 3 24 Equations 3 23 and 324 form a system of two linear algebraic equations in the two unknowns 01 and 0 If we solve these equations and integrate we nally obtain readers should ll in the details 516 7 We 52 We 3 25 114 3 SecondOrder Differential Equations Where Wt u1tu 2t 7 u1tu2ti 326 This expression Wt is called the Wronskiani Combining the previous expres sions gives the variation of parameters formula for the particular solution of 321 up ul Wm u2t dt The general solution of 321 is the homogeneous solution uht plus this par ticular solution If the antiderivatives in 325 cannot be computed explicitly then the integrals should be Written With a variable limit of integration Example 320 Find a particular solution to the DE u 9u 3sec3ti Here the homogeneous equation u 9u 0 has two independent solutions ul cos 3t and u sin St The Wronskian is Wt 3 cos2 t 3 sin2 3t 3 Therefore sin 3t 3 sec 3t cos 3t 3 sec 3t 01t 7 Tdt 02t Tdti Simplifying 01t itan Stdt 1ncos 3t7 02t ldt t We do not need constants of integration because we seek only the particular solution Therefore the particular solution is 1 up t g lncos 3t tsin St The general solution is ut cl cos 3t 02 sin 3t lncos 3t tsin St The constants may be determined by initial data if given 34 Variable Coefficients 115 When the secondorder equation has constant coef cients and the forcing term is a polynomial exponential sine or cosine then the method of undeter mined coef cients works more easily than the variation of parameters formula For other cases we use the formula or Laplace transform methods which are the subject of Chapter 4 Of course the easiest method of all is to use a com puter algebra systemi When you have paid your dues by using analytic methods on several problems then you have your license and you may use a computer algebra systemi The variation of parameters formula is important because it is often used in the theoretical analysis of problems in advanced differential equations EXERCISES ll Solve the following initial value problems a t2u Stu781i 07 111 07 u 1 2 b t2u tu 07 ul 011 2 c t2u 7 til211 07 ul 011 1i 0 i For what values of B is u t a solution to the equation 1 7 t2u 7 Ztu 2u 0 9 This exercise presents a transformation method for solving a Cauchy7Euler equationi Show that the transformation 1 lnt to a new independent variable I transforms the Cauchy7Euler equation at2u btu cu 0 into an linear equation with constant coef cients Use this method to solve Exercise la 4 Use the power series method to obtain two independent power series so lutions to u u 0 about to 0 and verify that the series are the expansions of cost and sint about t 0 5 Use the power series method to nd the rst three terms of two independent power series solutions to Airy s equation u 7 tu 07 centered at to 0i 533 Find the rst three terms of two independent power series solutions to the equation 1 t2u u 0 near to 0i 7 Solve the rstorder nonlinear initial value problem 11 1 112 u0 17 using a power series methodi Compare the accuracy of the partial sums to the exact solution Hint you will have to square out a power series 5 0 Consider the equation u 7 Ztu 27m 07 which is Hermite s differential equation an important equation in quantum theory Show that if n is a nonnegative integer then there is a polynomial solution Hn t of degree n 116 3 SecondOrder Differential Equations to H O H H H N H 9 H F H 5 H 533 l Consider the equation u 7 Zau 0211 07 Which has solution u e l Which is called a Hermite polynomial of degree n Find H0 t7 M7 H5 t up to a constant multiplei t Use reduction of order to nd a second independent solutioni l One solution of t2 t2 u HTut2u 0 is u1t ti Find a second independent solutioni l One solution of l t2u tu t2 7 Zu 0 is u1t t cos ti Find a second independent solutionl Let yt be one solution of the equation u ptu qtu 0 Show that the reduction of order method With ut vtyt leads to the rstorder linear equation yz 2y pyz0 2 11h Show that t Ce f 1700It 4 gt W7 and then nd a second linear independent solution of the equation in the form of an integral Use ideas from the last exercise to nd a secondorder linear equation that has independent solutions 6 and cos it Let ul and ug be independent solutions of the linear equation u ptu qtu 0 on an interval I and let Wt be the Wronskian of ul and ugl Show that Wt POWVOL and then prove that Wt 0 for all t 6 I or Wt is never zero on 1 Find the general solution of u tu u 0 given that u equot22 is one solution Use the transformation u exp ytdt to convert the secondorder equation u ptu qtu 0 to a Riccati equation y y2 pty qt 0 Conversely show that the Riccati equation can be reduced to the secondorder equation in u using the transformation y u ul Solve the rstorder nonautonomous equation 3 y W2 gy 35 Boundary Value Problems and Heat Flow 117 17 Use the variation of parameters formula to nd a particular solution to the following equations a u u a where a is a constant Note that l and lnt are two independent solutions to the homogeneous equationl b u u tantl e t2u 7 2u tgl H 5 0 Frobenius method Consider the differential equation Bessel s equation of order k 1 k2 u u lit 2 u0 where k is a real number a Show that to 0 is a regular singular point for the equation 710 b Assuming a solution of the form ut tr 200 ant show that 7 ik 0 V In the case that k nd the rst three terms of two independent series solutions to the DE L V Show that if k 0 then the Frobenius method leads to only one series solution and nd the rst three terms The entire series which con verges for all t is denoted by Jo t and is called a Bessel function of the rst kind of order zerol Finding a second independent solution is beyond the scope of our treatment 35 Boundary Value Problems and Heat Flow Let us consider the following problem in steadystate heat conduction A cylin drical uniform metallic bar of length L and crosssectional area A is insulated on its lateral side We assume the left face at z 0 is maintained at T0 degrees and that the right face at z L is held at TL degreesl What is the tempera ture distribution u in the bar after it comes to equilibrium Here represents the temperature of the entire cross section of the bar at position I where 0 lt z lt L We are assuming that heat ows only in the axial direction along the bar and we are assuming that any transients caused by initial tem peratures in the bar have decayed away In other words we have waited long 118 3 SecondOrder Differential Equations enough for the temperature to reach a steady state One might conjecture that the temperature distribution is a linear function of z along the bar397 that is To LETOzi This is indeed the case which we show belowi But also we want to consider a more complicated problems where the bar has both a variable conductivity and an internal heat source along its length An inter nal heat source for example could be resistive heating produced by a current running through the mediuml The physical law that provides the basic model is conservation of energy If z7z dz is any small section of the bar then the rate that heat ows in at z7 minus the rate that heat ows out at z dz7 plus the rate that heat is generated by sources must equal zero because the system is in a steady state See gure 35 laterally insulated Figure 35 Cylindrical bar laterally insulated through which heat is owing in the zdirectionl The temperature is uniform in a xed crosssectioni If we denote by the rate that heat ows to the right at any section z measured in caloriesarea time and we let denote the rate that heat is internally produced at z measured in caloriesvolume time then 7 A45z dz fzAdz 0i Canceling A dividing by dz7 and rearranging gives 451 d1 451 7 dz 7 m Taking the limit as dz a 0 yields WI f I 327 This is an expression of energy conservation in terms of ux But what about temperature Empirically the ux at a section z is found to be pro portional to the negative temperature gradient iuz which measures the 35 Boundary Value Problems and Heat Flow 119 steepness of the temperature distribution or pro le at that point or 328 This is Fourier s heat conduction lawl The given proportionality factor Kx is called the thermal conductivity in units of energy length degrees time which is a measure of how well the bar conducts heat at location It For a uniform bar K is constant The minus sign in 328 means that heat ows from higher temperatures to lower temperatures Fourier s law seems intuitively correct and it conforms with the second law of thermodynamics397 the larger the temperature gradient the faster heat ows from high to low temperatures Combining 3 27 and 328 leads to the equation 7Kzuzf 0 lt z lt L7 329 which is the steadystate heat conduction equationl When the boundary conditions u0 To7 uL Tl7 330 are appended to 329 we obtain a boundary value problem for the tem perature Boundary conditions are conditions imposed on the unknown state u given at different values of the independent variable I unlike initial conditions that are imposed at a single value For boundary value problems we usually use I as the independent variable because boundary conditions usually refer to the boundary of a spatial domain Note that we could expand the heat conduction equation to 7Kltzgtwltzgt e K mu m m 331 but there is little advantage in doing so Example 321 If there are no sources 0 and if the thermal conductivity K is constant then the boundary value problem reduces to u 07 0ltzltL7 u0 To uL Tll Thus the bar is homogeneous and can be characterized by a constant conduc tivity The general solution of u 0 is 01102 applying the boundary conditions determines the constants cl and Cg and gives the linear temperature distribution To TLETO I as we previously conjecturedl 120 3 SecondOrder Differential Equations ln nonuniform systems the thermal conductivity K depends upon location I in the systeml And K may depend upon the temperature u as well Moreover the heat source term f could depend on location and temperature In these cases the steadystate heat conduction equation 329 takes the more general form KW 101W f17u7 which is a nonlinear secondorder equation for the steady temperature distri bution u Boundary conditions at the ends of the bar may also specify the ux rather than the temperature For example in a homogeneous system if heat is injected at z 0 at a rate of N calories per area per time then the left boundary condition takes the form N7 or iKuO Ni Thus a ux condition at an endpoint imposes a condition on the derivative of the temperature at that endpointl In the case that the end at z L say is insulated so that no heat passes through that end then the boundary condition is uL 07 which is called an insulated boundary condition As the reader can see there are a myriad of interesting boundary value problems associated with heat owl Similar equations arise in diffusion processes in biology and chemistry for example in the diffusion oftoxic substances where the unknown is the chemical concentrationl Boundary value problems are much different from initial value problems in that they may have no solution or they may have in nitely many solutions Consider the following Example 322 When K 1 and the heat source term is 9n and both ends of a bar of length L 2 are held at u 0 degrees the boundary value problem becomes iu 9u7 0ltzlt2l u0 07 u20l The general solution to the DE is cl sin 31 Cg cos 31 where cl and Cg are arbitrary constants Applying the boundary condition at z 0 gives u0 cl sin3 0 Cg cos3 0 02 0 So the solution must have the form cl sinSIl Next apply the boundary condition at z 2 Then u2 cl sin6 07 to obtain cl 0 We have shown that the only solution 35 Boundary Value Problems and Heat Flow 121 is 0 There is no nontrivial steady state But if we make the bar length 7r then we obtain the boundary value problem fu 9u7 0ltzlt7rl u0 07 u7r 0 The reader should check that this boundary value problem has in nitely many solutions cl sin 31 where cl is any number If we change the right boundary condition one can check that the boundary value problem iu 9u7 0ltzlt7rl u0 07 117r17 has no solution at alll Example 323 Find all real values of A for which the boundary value problem iu Au7 0 lt z lt Tn 332 u0 07 uTr 07 333 has a nontrivial solution These values are called the eigenvalues and the corresponding nontrivial solutions are called the eigenfunctions lnterpreted in the heat ow context the left boundary is held at zero degrees and the right end is insulated The heat source is Aul We are trying to nd which linear heat sources lead to nontrivial steady states To solve this problem we consider different cases because the form of the solution will be different for A 07 A lt 07 A gt 0 HA 0 then the general solution of u 0 is azbl Then u z al The boundary condition u0 0 implies b 0 and the boundary condition u 7r 0 implies a 0 Therefore when A 07 we get only a trivial solutionl Next consider the case A lt 0 so that the general solution will have the form ut a sinh 7A1 b cosh iAzl The condition u0 0 forces 12 0 Then u t aTAcosh The right boundary condition becomes u 7r a iAcosh7 A 0 07 giving a 0 Recall that coshO 1 Again there is only the trivial solutionl Finally assume A gt 0 Then the general solution takes the form ut asin Az bcos 122 3 SecondOrder Differential Equations The boundary condition u0 0 forces 12 0 Then ut asin AI and u a A cos Applying the right boundary condition gives uTr a AcosA7r 0 Now we do not have to choose a 0 which would again give the trivial solution because we can satisfy this last condition with cos A7r 0 The cosine function is zero at the values 7r2 i mr n 0 123 mTherefore A7r7r2mr n0l23w Solving for A yields Alt2n 21gt2 n0123im Consequently the values of A for which the original boundary value problem has a nontrivial solution are 4 1 E These are the eigenvalues The corre sp onding solutions are um asinlt 2nl gt1 n0123w These are the eigenfunctions Notice that the eigenfunctions are unique only up to a constant multiplei In terms of heat flow the eigenfunctions represent possible steadystate temperature pro les in the bar The eigenvalues are those values A for which the boundary value problem will have steadystate pro lesi Boundary value problems are of great interest in applied mathematics sci ence and engineering They arise in many contexts other than heat flow includ ing wave motion quantum mechanics and the solution of partial differential equations EXERCISES 1 A homogeneous bar of length 40 cm has its left and right ends held at 30 C and 10 C respectively If the temperature in the bar is in steady state what is the temperature in the cross section 12 cm from the left end If the thermal conductivity is K what is the rate that heat is leaving the bar at its right face 35 Boundary Value Problems and Heat Flow 123 2 The thermal conductivity of a bar of length L 20 and crosssectional area A 2 is 1 and an internal heat source is given by 051L 7 If both ends of the bar are maintained at zero degrees What is the steady state temperature distribution in the bar Sketch a graph of What is the rate that heat is leaving the bar at z 20 CAD i For a metal bar of length L With no heat source and thermal conductivity Kz show that the steady temperature in the bar has the form x dy 7 510 m027 Where cl and Cg are constants What is the temperature distribution ifboth ends of the bar are held at zero degrees Find an analytic formula and plot the temperature distribution in the case that l zl If the left end is held at zero degrees and the right end is insulated nd the temperature distribution and plot it 4 Determine the values of A for Which the boundary value problem iu A117 0ltzlt17 u0 07 11107 has a nontrivial solution 5 Consider the nonlinear heat ow problem uu 07 0 lt z lt 7r7 u0 07 W 17 Where the thermal conductivity depends on temperature and is given by ul Find the steadystate temperature distribution 6 Show that if there is a solution u to the boundary value problem 3297330 then the following condition must hold 7KLuL K0u0 0 lnterpret this condition physicallyl 7 Consider the boundary value problem uH w2u 07 u0 a7 uL 12 When does a unique solution exist 124 3 SecondOrder Differential Equations 8 Find all values of A for which the boundary value problem iu72u Au 0ltzltl MO 07 1 0 07 has a nontrivial solution 9 Show that the eigenvalues of the boundary value problem fu Au 0ltzltl we 0 W an 07 are given by the numbers An p3 n 123 m where the 17 are roots of the equation tanp lpi Plot graphs of tanp and 1p and indicate graphically the locations of the values pni Numerically calculate the rst four eigenvaluesl 10 Find the values of A eigenvalues for which the boundary value problem 71211 7 11 Au 1 lt z lt e uG 07 ue7r 07 has a nontrivial solution 3 6 HigherOrder Equations So far we have dealt with rst and secondorder equationsi Higherorder equa tions occur in some applications For example in solid mechanics the vertical de ection y of a beam from its equilibrium satis es a fourthorder equa tion However the applications of higherorder equations are not as extensive as those for their rst and secondorder counterparts Here we outline the basic results for a homogeneous nthorder linear DE with constant coef cients u Pn71unil 1011 pou 0 334 The 17 239 01mn 7 1 are speci ed constants The general solution of 334 has the form ut clul t Cgu2t enunt where u1tu2twunt are independent solutions and where 5152Mcn are arbitrary constantsi In different words the general solution is a linear com bination of n different basic solutions To nd these basic solutions we try 36 HigherOrder Equations 125 the same strategy that worked for a secondorder equation namely assume a solution of the form of an exponential function ut em where A is to be determined Substituting into the equation gives A pHA H r m 00 0 335 which is an nth degree polynomial equation for Al Equation 335 is the Characteristic equation From algebra we know that there are n roots A1A27 in AW Here we are counting multiple roots and complex roots the lat ter will always occur in complex conjugate pairs ai bi A root A a has multiplicity K if A 7 aK appears in the factorization of the characteristic polynomial If the characteristic roots are all real and distinct we will obtain it different basic solutions u1t a u2t 6A2 7iii7unt ek i In this case the general solution of 334 will be a linear combination of these ut cleAll 62 A21Cnekquotti 336 If the roots of 335 are not real and distinct then we proceed as might be expected from our study of secondorder equations A complex conjugate pair A aiib gives rise to two real solutions 6 quot cos bt and 6 quot sin bti A double root A multiplicity 2 leads to two solutions 6 and teMi A triple root A multiplicity 3 leads to three independent solutions 6 t6 296 and so on In this way we can build up from the factorization of the characteristic polynomial a set of n independent basic solutions of 334 The hardest part of the problem is to nd the characteristic roots397 computer algebra systems are often useful for this task As may be expected from our study of secondorder equations an nthorder nonhomogeneous equation of the form um Pn71unil 1011 pou 9t7 337 has a general solution that is the sum of the general solution 336 of the homogeneous equation and a particular solution to the equation 3339 This result is true even if the coef cients pl are functions of ti Example 324 If the characteristic equation for a 6th order equation has roots A 72 i 3i7474747 71 the general solution will be ut 616721 cos 375 626721 sin 375 0364 C4t 4t C5i2 4t 0564 126 3 SecondOrder Differential Equations Initial conditions for an nth order equation 334 at t 0 take the form u0 a1 u 0 a2 mult 1gto an1 Where the ai are given constantsl Thus for an nthorder initial value problem we specify the value of the function and all of its derivatives up to the n 7 1st order at the initial time These initial conditions determine the n arbitrary constants in the general solution and select out a unique solution to the initial value problemi Example 325 Consider um 7 2n 7 Su 5641 The characteristic equation for the homogeneous equation is A372A273A0 A7 31 0 The characteristic roots are A 07 717 37 and therefore the homogeneous equa tion has solution uht cl 0267 03633 The particular solution Will have the form upt ae4 i Substituting into the original nonhomogeneous equation gives a l4l Therefore the general solution to the equation is l ut cl 026quot 0363 1641 The three constants can be determined from initial conditions For example for a thirdorder equation the initial conditions at time t 0 have the form u0 a7 u 0 B7 u 0 77 for some given constants 015 7 Of course initial conditions can be prescribed at any other time to EXERCISES 1 Find the general solution of the following differential equations a um u 0 b u u 1i 539quot N Summary and Review 127 C u u 0 d u 7 1 7Su 0 e u u 2 t St N Solve the initial value problem um7u 74u74u 07 u0 2 u 0 71 u 0 5 9 Write down a linear fthorder differential equation whose general solution is u cl Cgt age 4 e5t54 cos 2t C5 sin 5t 4 Show that the thirdorder equation u Zu 7 51 7 u 0 can be written as an equivalent system of three rstorder equations in the variables u v and w where v u and w Mk 5 What is the general solution of a fourthorder differential equation if the four characteristic roots are A 3 i i 3 i i What is the differential equation 37 Summary and Review One way to think about learning and solving differential equations is in terms of pattern recognition Although this is a very compartmentalized way of thinking it does help our learning process When faced with a differential equation what do we do The rst step is to recognize what type it is It is like a pianist recognizing a certain set of notes in a complicated musical piece and then playing those notes easily because of long hours of practice In differential equations we must practice to recognize an equation and learn the solution technique that works for that equation At this point in your study what kinds of equations should you surely be able to recognize and solve The simplest is the pure time equation 1 gt Here u is the antiderivative ofgt and we sometimes have to write the solution as an integral when we cannot nd a simple form for the antiderivativei The next simplest equation is the separable equation u gtfu7 128 3 SecondOrder Differential Equations where the right side is a product of functions of the dependent and independent variables These are easy just separate variables and integrate Autonomous equations have the form u fu7 where the right side depends only on u These equations are separable should we want to attempt a solution But often for autonomous equations we ap ply qualitative methods to understand the behavior of solutions This includes graphing vs u nding the equilibrium solutions and then drawing arrows on the phase line to determine stability ofthe equilibrium solutions and whether u is increasing or decreasing Nearly always these qualitative methods are su perior to having an actual solution formulai Firstorder autonomous equations cannot have oscillatory solutions Finally the rstorder linear equation is u ptu gti factors a 39 an equa Here we use variation of A t o 39 t quot tion can be solved by multiple methods397 for example u Zu 7 7 is separable linear and autonomous There are other rstorder nonlinear equations that can be solved and some of these were introduced in the Exercises The Bernoulli equation a plttgtu glttgtun can be transformed into a linear equation for the variable y ul 7 and the homogeneous equation 7 u u 7 f can be transformed into a separable equation for the variable y utl Solutions to special and unusual equations can sometimes be found in mathematical handbooks or in computer algebra systems There are really only two secondorder linear equations that can be solved simply These are the equation with constant coefficients cu 121 cu 07 where we have solutions of the form u em with A satisfying the characteristic equation W 12A c 07 and the CauchyiEuler equation at2u btu cu 07 where we have solutions of the form u tm7 where m satis es the charac teristic equation amm 7 l 12m c 0 For these two problems we must distinguish when the roots of the characteristic equation are real and unequal real and equal or complex When the right side of either of these equations is 37 Summary and Review 129 nonzero then the equation is nonhomogeneous Then we can nd particular so lutions using the variation of parameters method which works for all linear equations or use undetermined coef cients which works only for constant coef cient equations with special right sides Nonhomogeneous linear equations with constant coef cients can also be handled by Laplace transforms which are discussed in the next chapter All these methods extend to higherorder equa tions Generally we cannot easily solve homogeneous secondorder linear equa tions with variable coef cients or equations having the form u Ptu WW 0 Many of these equations have solutions that can be written as power series These power series solutions de ne special functions in mathematics such as Bessel functions Hermite polynomials and so forth In any case you can not solve these variable coef cient equations using the characteristic polynomial and nonhomogeneous equations are not amenable to the methods of undeter mined coef cients If you are fortunate enough to nd one solution you can determine a second by reduction of order If you are lucky enough to nd two independent solutions to the homogeneous equation the method of variation of parameters gives a particular solution The basic structure theorem holds for all linear nonhomogeneous equations the general solution is the sum of the general solution to the homogeneous equation and a particular solution This result is fundamental Secondorder equations coming from Newton s second law have the form I Ft7 z 1 These can be reduced to rstorder equations when t or z is missing from the force F or when F Fz7 which is the conservative case The Exercises give you review and practice in identifying and solving dif ferential equationsi EXERCISES 1 Identify each of the differential equations and nd the general solution Some of the solutions may contain an integrali a Zu 5u 7 Su 0 b u 7 Ru 07 where R is a parameter c 1 cost 7 ucosti d u 7 6u e e u 727ui f u 61 9u 5sinti 130 N 9 5 533 H H H q 5 0 i A population is governed by the law 10 Tplt 3 SecondOrder Differential Equations u 78t6i u ut272t2 u u7tu30i 2u u 3u0i k I I l tu 7i u t2u2i m u 73u2i n tu u7 cos2 um 5un 7 61 9631 p Gtu 7 ug 4u 372 7 3tu2u 0i Solve the initial value problem 1 u2 cos t7 u0 27 and nd the interval of existence Solve the initial value problem u u t7 u1 27 and nd the interval of existence i Use the power series method to nd the rst three terms of two independent solutions to u 7i tu 7i tu 0 valid near t 0 For all cases nd the equilibrium solutions for u u 7 a 112 7 a where a is a real parameter and determine their stability Summarize the infor mation on a bifurcation diagrami A spherical water droplet loses volume by evaporation at a rate propor tional to its surface area Find its radius 7 t in terms of the propor tionality constant and its initial radius Tor K717 gt Kap are positive constants Find the equilibria and their stability Describe in words the dynamics of the population where 7 K and a Use the variation of parameters method to nd a particular solution to u 7 u 7 2u coshti If e42 is one solution to the differential equation u 4tu 22t21u 0 nd the solution satisfying the conditions u0 3 u 0 1 Sketch the slope eld for the differential equation 1 7t2 sinu in the window 73 g t g 3 73 g t g 3 and then superimpose on the eld the two solution curves that satisfy u72 1 and u7l 17 respectively 2 i Solve u 4m 7 Tu lnu by making the subst1tutlon y ln ui 37 Summary and Review 12 Adapt your knowledge about solution methods for CauchyiEuler equations to solve the thirdorder initial value problem tsum 7 t2u 2tu 7211 0 With u1 311 2 u 1li Lapace Transforms The Laplace method for solving linear differential equations with constant co ef cients is based upon transforming the differential equation into an algebraic equationi It is especially applicable to models containing a nonhomogeneous forcing term such as the electrical generator in a circuit that is either discon tinuous or is applied only at a single instant of time an impulse This method can be regarded as another tool in addition to variation of parameters and undetermined coef cients for solving nonhomogeneous equa tions It is often a key topic in engineering where the stability properties of linear systems are addressed The material in this chapter is not needed for the remaining chapters so it may be read at any time 41 De nition and Basic Properties A successful strategy for many problems is to transform them into simpler ones that can be solved more easily For example some problems in rectangular co ordinates are better understood and handled in polar coordinates so we make the usual coordinate transformation 1 Tcost and y Tsin 9 After solving the problem in polar coordinates we can return to rectangular coordinates by the inverse transformation 7 12 y2 t9 arctan A similar technique holds true for many differential equations using integral transform meth ods In this chapter we introduce the Laplace transformation which has the 134 4 Laplace Transforms effect of turning a differential equation with state function ut into an algebra problem for an associated transformed function Us39 we can easily solve the algebra problem for US and then return to ut via an inverse transformation The technique is applicable to both homogeneous and nonhomogeneous linear differential equations with constant coef cients and it is a standard method for engineers and applied mathematicians It is particularly useful for differential equations that contain piecewise continuous forcing functions or functions that act as an impulse The transform goes back to the late 1700s and is named for the great French mathematician and scientist Pierre de Laplace although the basic integral goes back earlier to L Eulerl The English engineer 0 Heaviside developed much of the operational calculus for transform methods in the early 1900s Let u ut be a given function de ned on 0 g t lt 00 The Laplace transform of ut is the function Us de ned by Us 000 ute 5tdt 41 provided the improper integral exists The integrand is a function oft and s and we integrate on t leaving a function of 8 Often we represent the Laplace transform in function notation uts Us or just u Usl represents a functionlike operation called an operator or transform whose domain and range are sets of functions39 takes a function ut and transforms it into a new function US see gure 41 In the context of Laplace transfor mations t and u are called the time domain variables and s and U are called the transform domain variables In summary the Laplace transform maps functions ut to functions 13 and is somewhat like mappings we consider in calculus such as y 12 which maps numbers I to numbers y We can compute the Laplace transform of many common functions directly from the de nition 41 Example 41 Let ut eatl Then oo oo 1 1 Us eateistdt eclis dt 6a 5 l8 s gt a 0 0 a 7 s s 7 a 1 In different notation e 1 Sial Observe that this transform exists only for s gt a otherwise the integral does not exist Sometimes we indicate the values of s for which the transformed function Us is de ned 41 Definition and Basic Properties 135 Laplace Transform Us E Figure 41 The Laplace transform as a machine that transforms functions ut to functions Usl Example 42 Let ut t Then using integration by parts 00 675 1 oo 1 UltSgt te stdt t 7 1 eistdt 2 s gt 0 0 5 20 0 5 8 Example 43 The unit switching function ha t is de ned by ha t 0 ift lt a and ha t 1 if t 2 al The switch is off ift lt a and it is on when t 2 al Therefore the function ha t is a step function where the step from 0 to 1 occurs at t al The switching function is also called the Heaviside function The Laplace transform of hat is lhatl Awhate dt Oahate dt 0 hate 5 dt 0e dt1e dt 0 a 1 7 1 i eistligzo equot 7 s gt 0 s s 136 4 Laplace Transforms Example 44 The Heaviside function is useful for expressing multi lined functions in a single formulal For example let 1 0 tlt2 2 t71 2 t 3 fa 57 3ltt 6 07 tgt6 The reader should plot this function This can be written in one line as 1 1 ft h0t t 7 1 7 Ehgt 5 7 t2 7 t71h3t 7 5 7 t2h5tl The rst term switches on the function 12 at t 0397 the second term switches off 12 and switches on til at time t 239 the third term switches off t 71 and switches on 5 7 t2 at t 339 nally the last term switches off 5 7 t2 at t 6 Later we show how to nd Laplace transforms of such functions As you may have already concluded calculating Laplace transforms may be tedious businessl Fortunately generations of mathematicians scientists and engineers have computed the Laplace transforms of many many functions and the results have been catalogued in tables and in software systems Some of the tables are extensive but here we require only a short table which is given at the end of the chapter The table lists a function ut in the rst column and its transform Us or u in the second The various functions in the rst column are discussed in the sequel Computer algebra systems also have commands that calculate the Laplace transform see Appendix B Therefore given ut the Laplace transform US can be computed by the de nition given in formula 41 We can also think of the opposite problem given Us nd a function ut whose Laplace transform is Usl This is the inverse probleml Unfortunately there is no elementary formula that we can write down that computes ut in terms of Us there is a formula but it involves a contour integration in the complex plane In elementary treatments we are satis ed with using tablesl For example if Us 2 then the table gives ut 62 as the function that has Us as its transforml When we think of it this way we say ut 2t is the inverse transform of Us SEQ and l 2t7 71 e 7L Sigll we write In general we use the notation U Zu7 u 1Ul 41 Definition and Basic Properties 137 We think of as an operator transform and 1 as the inverse operation inverse transform The functions ut and Us form a transform pair and they are listed together in two columns of a table Computer algebra systems also supply inverse transformsi One question that should be addressed concerns the existence of the trans formi That is which functions have Laplace transforms Clearly if a function grows too quickly as it gets large then the improper integral will not exist and there will be no transformi There are two conditions that guarantee ex istence and these are reasonable conditions for most problems in science and engineering First we require that ut not grow too fast397 a way of stating this mathematically is to require that there exist constants M gt 0 and a for which iulttgti Mam is valid for all t gt to7 where to is some value of time That is beyond the value to the function is bounded above and below by an exponential function Such functions are said to be of exponential orderi Second we require that ut be piecewise continuous on 0 g t lt 00 In other words the interval 0 g t lt 00 can be divided into intervals on which u is continuous and at any point of discontinuity u has nite left and right limits except possibly at t ltxgti One can prove that if u is piecewise continuous on 0 g t lt gt0 and of exponential order then the Laplace transform Us exists for all s gt Oh What makes the Laplace transform so useful for differential equations is that it turns derivative operations in the time domain into multiplication operations in the transform domain The following theorem gives the crucial operational formulas stating how the derivatives transformi Theorem 45 Let ut be a function and 13 its transformi Then ul sUs 7 110 42 uH 82US 7 su0 7 uO 43 Proof These facts are easily proved using integration by parts We have cm warm ute o 7 7sute dt 0 0 7u0 3Us7 s gt 0 The second operational formula 43 is derived using two successive integration by parts and we leave the calculation to the reader 138 4 Laplace Transforms These formulas allow us to transform a differential equation with unknown ut into an algebraic problem with unknown Usl We solve for Us and then nd ut using the inverse transform u L71 U i We elaborate on this method in the next section Before tackling the solution of differential equations we present additional important and useful properties a Linearity The Laplace transform is a linear operation397 that is the Laplace transform of a sum of two functions is the sum of the Laplace transforms of each and the Laplace transform of a constant times a func tion is the constant times the transform of the function We can express these rules in symbols by a single formula clu 0211 51 u52 vl 44 Here u and v are functions and cl and 02 are any constantsl Similarly the inverse Laplace transform is a linear operation 1clu 521 cl illu Cg ilml 45 b Shift Property The Laplace transform of a function times an exponen tial ute quot shifts the transform of U that is ute quot Us 7 a 46 A 8 Switching Property The Laplace transform of a function that switches on at t a is given by ha tut 7 a Usequot i 439 Proofs of some of these relations follow directly from the de nition of the Laplace transform and they are requested in the Exercises EXERCISES H Use the de nition of the Laplace transform to compute the transform of the square pulse function ut l 1 g t g 239 ut 0 otherwise 2i Derive the operational formula 4 3 3 Sketch the graphs of sin t7 sint 7 7r27 and hW2t sint 7 7r2l Find the Laplace transform of each 4 Find the Laplace transform of t2e 3 l 5 Find sinh kt and cosh kt using the fact that em P H Definition and Basic Properties 139 on 7 90 3 HH HO H N H 9 H F H 5 H 533 i Find 6 3 4 sin kt using the table Find e gt sin 2t using the shift property 46 Using the switching property 47 nd the Laplace transform of the func W From the de nition 41 nd using the integral substitution st 7 2 and then noting few exp739r2d39r tion tlt2 tgt2i i Does the function ut 6 2 have a Laplace transform What about ut lt Commenti Derive the operational formulas 4 6 and 47 i Plot the squarewave function W zltilgtnhnlttgt on the interval t gt 0 and nd its transform Fsi Hint use the geometric serieslzx2ii Show that n Alana 1U s 7mm Derive the formulas 5 WW USL Use these formulas to nd the inverse transform of arctan W mm and use the result to nd Show that sinh t 75 Show that 5 lfthatl ei ls l t al7 and use this formula to compute t2h1ti Find the Laplace transform of the function in Example 44 140 4 Laplace Transforms L ODEin ut gt Algebraic equation in Us solve 1 L lt Us S quotquotquotquot39 Figure 42 A DE for an unknown function ut is transformed to an algebraic equation for its transform Usi The algebraic problem is solved for US in the transform domain and the solution is returned to the original time domain via the inverse transformi 17 The gamma function is de ned by Fz 000 eittxildt z gt 71 a Show that Fnl nFn and Fnl ml for nonnegative integers TL Show that b Show that u a sgt 0 42 Initial Value Problems The following examples illustrate how Laplace transforms are used to solve ini tial value problems for linear differential equations with constant coef cients The method works on equations of all orders and on systems of several equa tions in several unknowns We assume ut is the unknown state function The idea is to take the transform of each term in the equation using the linearity property Then using Theorem 45 reduce all of the derivative terms to alge braic expressions and solve for the transformed state function Usl Finally invert Us to recover the solution Figure 42 illustrates this threestep methodl Example 46 Consider the secondorder initial value problem u w2u 07 u0 07 uO ll 42 Initial Value Problems 141 Taking transforms of both sides and using the linearity property gives uN l w2 lu lo Then Theorem 45 gives 82US 7 su0 7 uO w2Us 07 which is an algebraic equation for the transformed state Usl Using the initial conditions we get 82US 7 l w2Us 0 Solving for the transform function Us gives 1 l w w2s2 ww2327 Us which is the solution in the transform domain Therefore from the table the inverse transform is l ut sinwt7 which is the solution to the original initial value probleml Example 47 Solve the rstorder nonhomogeneous equation u Zu equot u0 0 Taking Laplace transforms of each term ul 2ul leitl 1 31 sUs 7 u0 2U3 Solving for the transformed function Us gives 1 U 51gtlts2gt Now we can look up the inverse transform in the table We nd ut K1 72t 1 l sls2 142 4 Laplace Transforms Example 48 Partial Fractions 1 Sometimes the table may not include an entry for the inverse transform that we seek and so we may have to algebraically manipulate or simplify our expression so that it can be reduced to table entries A standard technique is to expand complex fractions into their partial fraction expansion In the last example we had 1 U s s 1 s 2 We can expand Us as 1 7 a b 31s2 81 327 for some constants a and b to be determined Combining terms on the right side gives 1 7 as2bs1 31s2 7 31s2 7 abs2ab 7 s182 Comparing numerators on the left and right force a b 0 and 2a b 1 Hence a 7b 1 and we have 1 1 71 U s ltswlts2gt 81 2 We have reduced the complex fraction to the sum of two simple easily identi able fractions that are easily found in the table Using the linearity property of the inverse transform 1 1 l lwl l s 1 s 2 it 7 6721 lms 6 Example 49 Partial Fractions 11 A common expression in the transform domain that requires inversion is a fraction of the form 1 U s 32 118 c If the denominator has two distinct real roots then it factors and we can proceed as in the previous example If the denominator has complex roots 42 Initial Value Problems 143 then the following complete the square technique may be used For example consider U5 32 33 6 Then completing the square in the denominator 1 W W 8235 75 6 1 2A 8 3 g This entry is in the table up to a factor of Therefore we multiply and divide by this factor and locate the inverse transform in the table as 2 VB 14 VB 2 e 3W2 sin ti Example 410 In this example we calculate the response of an RC circuit when the emf is a discontinuous function These types of problems occur frequently in engineer ing especially electrical engineering where discontinuous inputs to circuits are commonplace Therefore consider an RC circuit containing a 1 volt battery and with zero initial charge on the capacitor Take R 1 and C 1Sl Assume the switch is turned on from 1 g t g 2 and is otherwise switched off giving a square pulse The governing equation for the charge on the capacitor is 4 3t h1t h2t7 40 0 We apply the basic technique Taking the Laplace transform gives 1 8Q8 we 3Q8 gen 7 Solving for yields 628 1 ss3 1 is 725 ss3e ss3e Now we have to invert which is usually the hardest part Each term on the e4 7 e45 right has the form Usequot 7 and therefore we can apply the switching property 417 From the table we have 144 4 Laplace Transforms 1 2 3 t Figure 43 The charge response is zero up to time t 1 when the switch is closed The charge increases until t 2 when the switch is again opened The charge then decays to zero Therefore7 by the shift property 1 1 71 is 7 73271 383e 31 6 h1tl Similarly7 Putting these two results together gives qlttgt lt17 e3ltHgtgth1lttgt 7 lt17 e3lt 2gtgth2lttgt We can use software to plot the charge responsel See gure 43 Because there are extensive tables and computer algebra systems containing large numbers of inverse transforms the partial fractions technique for inversion is not used as often as in the past EXERCISES 1 Find A B and C for which 1 7A3B C 828717 32 371 Then nd the inverse Laplace transform of 2 Find the inverse transform of the following functions 43 The Convolution Property 145 a US mi b U8 c d ge 45i 9 Solve the following initial value problems using Laplace transforms a u 5u h2t7 W 1 b u u sin2t7 u0 0 c u 7 u 7 6u 7 0 u0 2710 1 d u 7 21 2u 07 710 07 uO 1 e u 7 21 2u equot u0 07 MO 1 f u 7 u 07 110 17 uO 0 g u 041 2u 17 h5t u0 07 u 0 0 h u 9n 7 sin 3t u0 07 u 0 0 i u 7 Zu 17 u0 17 uO 0 4 Use Laplace transforms to solve the two simultaneous differential equations 1 z 7 2y 7 t y 31 y7 with 10 y0 0i Hint use what you know about solving single equations letting z Xs and y Ysi 5 Show that thnutl 1 Uns for n 1237im 43 The Convolution Property The additiVity property of Laplace transforms is stated in 44 the Laplace transform of a sum is the sum of the transforms But what can we say about the Laplace transform of a product of two functions It is not the product of the two Laplace transforms That is if u ut and v vt with u US and v Vs then uv y U3Vsl If this is not true then what is true We 146 4 Laplace Transforms ask it this way What function has transform U3Vs Or7 differently7 what is the inverse transform of U3Vsl The answer may surprise you because it is nothing one would easily guess The function whose transform is U3Vs is the convolution of the two functions ut and vtl It is de ned as follows If u and v are two functions de ned on 07 0 the convolution of u and v denoted by u 96 v is the function de ned by u vt Otu739vt 7 Td7 Sometimes it is convenient to write the convolution as ut vtl The convo lution property of Laplace transforms states that u v U3Vsl It can be stated in terms of the inverse transform as well 1U8Vsl u vi This property is useful because when solving a DE we often end up with a product of transforms 7 we may use this last expression to invert the product The convolution property is straightforward to verify using a multivariable calculus technique interchanging the order of integration The reader should check the following steps ltAtu7vt77d7gt 7 0 lt lt u739vt 7 Td7 gt e Stdt t t u739vt 7 re5 drgt dt 0 u739vt 7 re 5 dtgt d7 0 lt o lt vTe5ltTTgtdTgt u739d739 0 vre dTgt WWW 000 Emma 0 mysw This last expression is U3Vsl vt 7 Te 5 dtgt u7 d7 A 0 A l 0 0 43 The Convolution Property 147 Example 411 Find the convolution of 1 and t2 We have t z 1t2 lt772d7t272t772d7 0 0 t2 t3 t3 t2t72t 2H3 Notice also that the convolution of t2 and l is t is t21 72gt1dT 0 3 In the exercises you are asked to show that u 96 v v u7 so the order of the two functions under convolution does not matter Example 412 Find the inverse of 13 We can do this by partial fractions but here we use convolution We have 3 1 3 71 71 lslt829gtl lslt329gtl z 1sin3t sinSTdT 0 l l7cos3t Example 413 Solve the nonhomogeneous DE u k2u ft7 where f is any given input function and where u0 and u 0 are speci ed initial conditions Taking the Laplace transform 82US 7 su0 7 uO k2Us Fs Then 32 i k2 m Now we can invert each term using the convolution property on the last term Us u0 to get the solution formula WU k 1 z us u0 cos kt sin kt E sin Ht 7 Td39r 0 148 4 Laplace Transforms EXERCISES 1 Compute the convolution of sint and cos ti 2 Compute the convolution oft and It 3 Use the convolution property to nd the general solution of the differential equation u au qt using Laplace transforms 4 Use a change of variables to show that the order of the functions used in the de nition of the convolution does not matter That is u vt v 9f Solve the initial value problem u 7 w2u ft7 u0 uO 0i 5 6 Use Exercise 5 to nd the solution to u 7 4u 17 h1t7 u0 uO 0i 7 Write an integral expression for the inverse transform of Us e 35Fs7 where Fl 8 Find a formula for the solution to the initial value problem u 71 t7 u0 uO 0 9 An integral equation is an equation where the unknown function ut ap pears under an integral sign see also the exercises in Section 12 Consider the integral equation ulttgt ft w 7 Tgtultrgtdn where f and k are given functionsl Using convolution nd a formula for 13 in terms of the transforms of F and K of f and k respectively H 53 Using the idea in the preceding exercise solve the following integral equa tions a ut t7 gai ww b ut f etilu7d7l H H l Solve the integral equation for ut 7 1 u739 ft 7 70 Edi Hint use the gamma function from the Exercise 17 in Section All 44 lquot39 39 Sources 149 44 Discontinuous Sources The problems we are solving have the general form u bu cu ft tgt0 u0 u1u0u2i If f is a continuous function then we can use variation of parameters to nd the particular solution39 if f has the special form of a polynomial exponential sine or cosine or sums and products of these forms we can use the method of undetermined coef cients judicious guessing to nd the particular solution If however f is a piecewise continuous source with different forms on different intervals then we would have to nd the general solution on each interval and determine the arbitrary constants to match up the solutions at the endpoints of the intervals This is an algebraically dif cult task However using Laplace transforms the task is not so tediousi In this section we present additional examples on how to deal with discontinuous forcing functionsi Example 414 As we noted earlier the Heaviside function is very useful for writing piecewise or multilined functions in a single line For example t 0lttlt1 ft 2 1 t 3 0 tgt3 t 2 7 th1t 7 2h3ti The rst term switches on the function t at t 039 the second term switches on the function 2 and switches off the function t at t 139 and the last term switches off the function 2 at t 3 By linearity the Laplace transform of is given by 478 ltl 2 lh1tl lthlw 2 lhstl The second and fourth terms are straightforward from Example 43 and t 13 The third term can be calculated using ftha e 5 ft ai With t we have 75 1 7s 1 75 th1t e t1 3 26 E i Putting all these results together gives 1 2 l 1 2 Fs g geis 7 3 2675 geisgt 7 geigsi 150 4 Laplace Transforms Example 415 Solve the initial value problem u 9u 6 0395 h4t7 u0 uO 07 where the forcing term is an exponential decaying term that switches on at time t 4 The Laplace transform of the forcing term is 1 L 705 t 74s 705z4 72 745 e M e e l e use Then taking the transform of the the equation l 2 7 72 745 s Us 9Us 7 6 8 056 i Whence 1 U 72 74s 3 e s 0 532 9 6 Now we need the shift theoremi But rst we nd the inverse transform of Z W W 5 9 i Here we leave it as an exercise partial fractions to show 71 l 7 36 0395 73cos3t0i5sin3t s 0i5s2 9 7 27 75 Therefore by the shift property 745 t 72 71 e u e f Hoisxsug 3e 050 4 7 3cos 3t 7 4 05 sin3t 7 4 277562 7 which is the solution Notice that the solution does not switch on until t 4 At that time the forcing term turns on producing a transient397 eventually its effects decay away and an oscillating steadystate takes over hm EXERCISES 1 Sketch the function 2h3t 7 2h4 t and nd its Laplace transformi 2 Find the Laplace transform of t2h3ti 3i lnvert Fs s 7 2 4i 4 Sketch the following function write it as a single expression and then nd its transform 37 0 g t lt 2 2 2 lt t lt 7r t 7 7 f 67 7r g t g 7 07 t gt 7 i lquot39 39 Sources 5i 533 7 H 5350 H H H 0 H 9 Find the inverse transform of 1 7 6745 ms 82 Solve the initial value problem 7 cos 2t 0 g t 2 u 1 4 1 0 t gt 2 where u0 u 0 0 Sketch the solution Consider the initial value problem 1 u ft u0 1 where is 0 lt t g 1 given by 1 07 7 72 t gt 1 Solve this problem in two ways a by solving the problem on two intervals and pasting together the solutions in a continuous way and b by Laplace transforms An LC circuit with L C 1 is rampedup with an applied voltage lt lt 6t t 07t79 97 t gt 9 Initially there is no charge on the capacitor and no current Find and sketch a graph of the voltage response on the capacitori Solve u 7u h1t 7 h2t u0 1i Solve the initial value problem 2 2 7 7r 0lttlt1 0 tgt1 where u0 1 and u 0 0 i Let be a periodic function with period pi That is ft 10 for all t gt 0 Show that the Laplace transform of f is given by 1 17 77 5 Fs mO f39re dTi Hint break up the interval foooo into subintervals mp n Up calculate the transform on each subinterval and nally use the geometric series1z12 i 1 Show that the Laplace transform of the periodic squarewave function that takes the value 1 on intervals 0a 2a3a 4a5aw and the value 71 on the intervals a2a 3a4a 5a6aw is tanh 1 Write a single line formula for the function that is 2 between 2n and 2n1 and 1 between 2n 7 1 and Zn where n 0 1 234 152 4 Laplace Transforms 45 Point Sources Many physical and biological processes have source terms that act at a single instant of time For example we can idealize an injection of medicine a shot into the blood stream as occurring at a single instant39 a mechanical system for example a damped springmass system in a shock absorber on a car can be given an impulsive force by hitting a bump in the road39 an electrical circuit can be closed only for an instant which leads to an impulsive applied voltage To x the idea let us consider an RC circuit with a given emf et and with no initial charge on the capacitor In terms of the charge qt on the capacitor the governing circuit equation is qu 4 ea 40 0 48 This is a linear rstorder equation and if the emf is a continuous function or piecewise continuous function the problem can be solved by the methods presented in Chapter 2 or by transform methods We use the latter Taking Laplace transforms and solving for Qs the Laplace transform of qt gives 1 1 Q6 RWEQL where Es is the transform of the emf etl Using the convolution property we have the solution t qt l gainWaf e 49 R 0 But presently we want to consider a special type of electromotive force 6t one given by an voltage impulse that acts only for a single instant ie a quick surge of voltage To x the idea suppose the source is a 1 volt battery Imagine that the circuit is open and we just touch the leads together at a single instant at time t a How does the circuit respond We denote this unit voltage input by et L t which is called a unit impulse at t a The question is how to de ne 6at an energy source that acts at a single instant of time At rst it appears that we should take 6at l ift a and 6at 0 otherwise But this is not correct To illustrate we can substitute into 49 and write I 4t 0 elt TgtRcda7dr 410 If L t 0 at all values oft except t a the integral must be zero because the integrand is zero except at a single point Hence the response of the circuit is qt 0 which is incorrectl Something is clearly wrong with this argument and our tentative de nition of 6a 45 Point Sources 153 The dif culty is with the function Z We must come to terms with the idea of an impulse Actually having the source act at a single instant of time is an idealizationi Rather such a short impulse must occur over a very small interval a 7 52 a 52 where 5 is a small positive number We do not know the actual form of the applied voltage over this interval but we want its average value over the interval to be 1 volt Therefore let us de ne an idealized applied voltage by Ema g a752lttlta52 0 otherwise gazes2m 7 hug2m These idealized impulses are rectangular voltage inputs that get taller and narrower of height 18 and width 5 as 5 gets smalli But their average value over the small interval a 7 52 lt t lt a 52 is always 139 that is HE2 eaE tdt 1 11752 This property should hold for all 5 regardless of how smalli It seems reasonable therefore to de ne the unit impulse L t at t a in a limiting sense having the property aE2 6atdt 1 for all 5 gt 0 11752 Engineers and scientists used this condition along with L t 0 t f a for decades to de ne a unit point source at time t a called the delta function and they developed a calculus that was successful in obtaining solutions to equations having point sources But actually the unit impulse is not a function at all and it was shown in the mid20th century that the unit impulse belongs to a class of so called generalized functions whose actions are not de ned pointwise but rather by how they act when integrated against other functions Mathematically the unit impulse is de ned by the sifting property aalttgt lttgtdt 45w 0 That is when integrated against any nice function 45t the delta function picks out the value of at t at We check that this works in our problemi If we use this sifting property back in 410 then for t gt a we have I LN 67t77R06a7d7 eitiaRC7 t gt a 154 4 Laplace Transforms which is the correct solution Note that qt 0 up until t a because there is no source Furthermore qa lRl Therefore the charge is zero up to time a at which it jumps to the value lR and then decays away To deal with differential equations involving impulses we can use Laplace transforms in a formal way Using the sifting property with 6 obtain we 00 6at gamma e 0 which is a formula for the Laplace transform of the unit impulse function This gives of course the inverse formula 1equot 6a The previous discussion is highly intuitive and lacks a careful mathematical basel However the ideas can be made precise and rigorous We refer to advanced texts for a thorough treatment of generalized functionsl Another common no tation for the unit impulse 6at is Kit 7 a If an impulse has magnitude f0 instead of 1 then we denote it by f06a For example an impulse given a circuit by a 12 volt battery at time t a is 126a Example 416 Solve the initial value problem u u 52t7 u0 u 0 07 with a unit impulse applied at time t 2 Taking the transform 82US SUS 6725 Thus 725 US Using the table it is simple to nd Therefore by the shift property the solution is 71 ut L l 17 e quot2h2tgt The initial conditions are zero and so the solution is zero up until time t 2 when the impulse occursl At that time the solution increases with limit 1 as t 4 00 See gure 44 45 Point Sources 155 Figure 44 Solution in Example 416 EXERCISES 1 0 9 F 5335 Compute few 6 20 3gt264 tdtl l Solve the initial value problem u Su 61t h4t7 u0 1 Sketch the solution Solve the initial value problem u 7 u 5507 u0 uO 0 Sketch the solution Solve the initial value problem u u 5207 u0 uO 0 Sketch the solution lnvert the transform Fs 825 6 35 Solve the initial value problem u 4u u0 52W 5507 uO 0i 156 4 Laplace Transforms I 5 0 l Consider an LC circuit with L C 1 and 110 v 0 07 containing a 1 volt battery where v is the voltage across the capacitor At each of the times t 07 7r7 2w 3W7 M7 mr7 the circuit is closed for a single instantl Determine the resulting voltage response vt on the capacitor Compute the Laplace transform of the unit impulse in a different way from that in this section by calculating the transform of eaE t and then taking the limit as 5 4gt 0 Speci cally show new aloha20gt a has2a 71 5 is 5 Then use l Hospital s rule to compute the limit 2 sinh lim 87 5A0 5 thereby showing 46 Table of Laplace Transforms 46 Table of Laplace Transforms W 118 t 7 n0717273w ta 1quot a1 5a sin kt g cos kt 7 sinh kt W cosh kt 6 quot sin kt W 6 quot cos kt W eat 7 ebt 1 171 57a 57b tneat 7L ut u u t uat ha t ute quot Mt ht tut 7 a 0 u7 vt 7 Td7 f WM t W sUs 7 u0 8 7 sum 7 we 3 Us 7 sn 1u0 7 t 7 un710 gm I355 7 a Usequot U3Vs 7a5 ft a 157 5 Linear Systems Up until now we have focused upon a single differential equation with one unknown state function Yet most physical systems require several states for their characterizationi Therefore we are naturally led to study several differ ential equations for several unknownsi Typically we expect that if there are n unknown states then there will be n differential equations and each DE will contain many of the unknown state functionsi Thus the equations are coupled together in the same way as simultaneous systems of algebraic equations If there are n simultaneous differential equations in n unknowns we call the set of equations an ndimensional systemi 51 Introduction A twodimensional linear homogeneous system of differential equations has the form 1 azby 51 y or dy 52 where a b c and d are constants and where z and y are the unknown states A solution consists of a pair of functions I 1t y yt that when substituted into the equations reduce the equations to identitiesi We can interpret the solution geometrically in two ways First we can plot I zt and y yt 160 5 Linear Systems WWO W W Phase plane Time series Figure 51 Plots showing the two representations of a solution to a system for t Z 0 The plot to the left shows the time series plots m mt y yt and the plot to the right shows the corresponding orbit in the my phase planer vsl t on the same set of axes as shown in gure Sill These are the time series plots and they tell us how the states m and y vary in time Or second we can think of m mt y yt as parametric equations of a curve in an my plane with time t as the parameter along the curve See gure Sill In this latter context the parametric solution representation is called an orbit and the my plane is called the phase planer Other words used to describe a solution curve in the phase plane in addition to orbit are solution curve path and trajectory These words are often used interchangeably ln multi variable calculus the reader probably used the position vector xt Iti to represent this orbit where i and j are the unit vectors but here we use the W lt 83 gt To save vertical space in typesetting we often write this column vector as column vector notation mtytT where T denotes transpose39 transpose means turn the row into a column Mostly we use the phase plane representation of a solution rather than the time series representation The linear system 515 2 has in nitely many orbits each de ned for all times 700 lt t lt 00 When we impose initial conditions which take the form 1050 107 M750 yo then a single orbit is selected out That is the initial value problem consist ing of the system 515 2 and the initial conditions has a unique solution Equations 515 2 also give geometrical information about the direction of the solution curves in the phase plane in much the same way as the slope eld of a single differential equation gives information about the slopes of a solution curve see Section 112 At any point m y in the my plane the right 51 Introduction sides of 515 2 de ne a vector lt gtltwbygt v 7 x 7 7 7 y or dy which is the tangent vector to the solution curve that goes through that point Recall from multivariable calculus that a curve 1t7 ytT has tangent vector V z ty tTl We can plot or have software plot for us this vector at a large set of points in the plane to obtain a vector eld a eld of vectors that indicates the flow or direction of the solution curves as shown in gure 52 The orbits t in so that their tangent vectors coincide with the vector Figure 52 In the phase plane the vector eld v z 7 y7 I yT associated with the system 1 z 7 y y z y and several solution curves 1 1t y yt which spiral out from the origin The vector eld is tangent to the solution curves The orbits approach in nity as time goes forward ie t 7gt 00 and they approach the origin but never reach it as time goes backward ie t 7 fool eld A diagram showing several key orbits is called a phase diagram or phase portrait of the system 50752 The phase portrait may or may not include the vector eld Example 51 We observed in Chapter 3 that a secondorder differential equation can be reformulated as a system of two rstorder equations For example the damped 162 5 Linear Systems spring mass equation m1 7km 7 01 can be rewritten as I 97 k c y 7 1 7 y7 m m where z is position or displacement of the mass from equilibrium and y is its velocity This system has the form of a twodimensional linear systemi In this manner7 mechanical problems can be studied as linear systems With speci c physical parameters k m 1 and c 05 we obtain 1 97 y 71 7 0i5yi The response of this damped springmass system is a decaying oscillationi Fig ure 53 shows a phase diagrami Figure 53 Phase plane diagram and vector eld for the system 1 y y 7 71 70i5y showing several orbits which are spirals approaching the origin These spirals correspond to time series plots of z and y vs t that oscillate and decay 51 Introduction 163 1 Figure 54 Two compartments with arrows indicating the ow rates between them and the decay rate Example 52 In this example we generalize ideas presented in Section 135 on chemical reactors and the reader should review that material before continuing The idea here is to consider linked chemical reactors or several different compart ments Compartmental models play an important role in many areas of science and technology and they often lead to linear systems The compartments may be reservoirs organs the blood system industrial chemical reactors or even classes of individuals Consider the two compartments in gure 54 where a chemical in compartment 1 having volume V liters ows into compartment 2 of volume W liters at the rate of 4 liters per minute In compartment 2 it decays at a rate proportional to its concentration and it ows back into compartment 1 at the same rate 4 At all times both compartments are stirred thoroughly to guarantee perfect mixing Here we could think of two lakes or the blood system and an organ Let I and y denote the concentrations of the chemical in compartments 1 and 2 respectively We measure concentrations in grams per liter The technique for nding model equations that govern the concentrations is the same as that noted in Section 1315 Namely use mass balance Regardless of the number of compartments mass must be balanced in each one the rate of change of mass must equal the rate that mass ows in or is created minus the rate that mass ows out or is consumed The mass in compartment 1 is VI and the mass in compartment 2 is Wyl A mass ow rate mass per time equals the volumetric ow rate volume per time times the concentration mass per volume So the mass ow rate from compartment 1 to 2 is 41 and the mass ow rate from compartment 2 to 1 is qyl Therefore 164 5 Linear Systems balancing rates gives VI iqzqy7 Wy 41 7 qyi WW where k is the decay rate grams per volume per time in compartment 2 The volume W must appear as a factor in the decay rate term to make the dimensions correct We can write this system as 2 i I 7 VzVy7 7 i 7 i y WI ltWkgty7 which is the form of a two dimensional linear systemi In any compartment model the key idea is to account for all sources and sinks in each compartment lfthe volumetric ow rates are not the same then the problem is complicated by variable volumes7 making the problem nonhomogeneous and time dependent But the technique for obtaining the model equations is the same EXERCISES 1 Verify that xt cos 2t7 72 sin 2tT is a solution to the system Sketch a time series plot of the solution and the corresponding orbit in the phase planer Indicate the direction of the orbit as time increases By hand plot several vectors in the vector eld to show the direction of solution curves 2 Verify that 26 t xlt gt lt 36 gt is a solution to the linear system 1 AlerZy7 y7317yi Plot this solution in the my plane for t 6 700700 and nd the tangent vectors along the solution curvei Plot the time series on 700 lt t lt 00 3 Consider the linear system z7zy y4174y7 with initial conditions 10 10 y0 0 Find formulas for the solution zt7yt and plot their time series Hint multiply the rst equation by 4 and add the two equations 01 M Matrices 165 g i In Example 52 let V 1 liter W 05 liters q 005 liters per minute with no decay in compartment 2 If 10 10 grams per liter and y0 0 nd the concentrations zt and yt in the two compartments Hint add appropriate multiples of the equations Plot the time series and the corresponding orbit in the phase planer What are the concentrations in the compartments after a long time 5 ln Exercise 4 assume there is decay in compartment 2 with k 02 grams per liter per minute Find the concentrations and plot the time series graphs and the phase ploti Hint transform the system into a single secondorder equation for 533 In the damped springmass system in Example 51 take m l k 2 and with initial conditions 10 4 and y0 0 Find formulas for the position zt and velocity Show the time series and the orbit in the phase plane I i Let L and I be the charge and the current in an RCL circuit with no electromotive forcei Write down a linear system of rstorder equations that govern the two variables 4 and I Take L l R 0 and C i If q0 8 and 10 0 nd qt and Iti Show a time series plot of the solution and the corresponding orbit in the q phase planer 52 Matrices The study of simultaneous differential equations is greatly facilitated by matri cesi Matrix theory provides a convenient language and notation to express many of the ideas conciselyi Complicated formulas are simpli ed considerably in this framework and matrix notation is more or less independent of dimension In this extended section we present a brief introduction to square matricesi Some of the de nitions and properties are given for general n by n matrices but our focus is on the two and threedimensional cases This section does not rep resent a thorough treatment of matrix theory but rather a limited discussion centered on ideas necessary to discuss solutions of differential equations A square array A of numbers having n rows and n columns is called a square matrix of size n or an n X n matrix we say n by n matrix The number in the ith row and jth column is denoted by aiji General 2 X 2 and 3 X 3 matrices have the form all an all a12 dis Alt gt7 A a21 a22 a23 0421 0422 166 5 Linear Systems The numbers aij are called the entries in the matrix 7 the rst subscript i denotes the row and the second subscript j denotes the column The main diagonal of a square matrix A is the set of elements 0411704227 Wannl We often Write matrices using the brief notation A aij An n Vector x is a list of n numbers 11127m71n Written as a column so vector means column list The numbers 1112 m In in the list are called its components For example Xlt gt I2 is a 2 vectorl Vectors are denoted by lowercase boldface letters like X y etc and matrices are denoted by capital letters like A B etc To minimize space in typesetting we often Write for example a 2vector x as 11 12T Where the T denotes transpose meaning turn the row into a column Two square matrices having the same size can be added entryWiser That is if A clj and B bij are both n X n matrices then the sum A B is an n X n matrix de ned by A B clj bij A square matrix A clj of any size can be multiplied by a constant c by multiplying all the elements of A by the constant in symbols this scalar multiplication is de ned by CA calj Thus 7A iaij and it is clear that A 7A 07 Where 0 is the zero matrix having all entries zero If A and B have the same size then subtraction is de ned by A 7 B A 7B Also A 0 A if 0 has the same size as Al Addition When de ned is both commutative and associativel Therefore the arithmetic rules of addition for n X n matrices are the same as the usual rules for addition of numbers Similar rules hold for addition of column vectors of the same length and multiplication of column vectors by scalars397 these are the de nitions you en countered in multivariable calculus Where n vectors are regarded as elements of Rnl Vectors add componentWise and multiplication of a vector by a scalar multiplies each component of that vector by that scalarl Example 53 l 2 0 72 74 5 Alt344gt Blt7 4 x a gt ylt1gt 1 0 0 6 AB lt10 78gt7 73Blt72 12gt7 6 8 lt30 gty 52 Matrices 167 The product of two square matrices of the same size is not found by mul tiplying entrywiser Rather matrix multiplication is de ned as follows Let A and B be two n X n matrices Then the matrix AB is de ned to be the n X n matrix C cij where the ij entry in the ith row and jth column of the product C is found by taking the product dot product as with vectors of the ith row of A and the jth column of Bi In symbols AB C7 where Cij ai 39 bj aiibij ai2b2j ainbnj7 where ai denotes the 2th row of A7 and bj denotes the jth column of Bi Generally matrix multiplication is not commutative iiei AB y BA so the order in which matrices are multiplied is important However the associative law ABC ABC does hold so you can regroup products as you wish The distributive law connecting addition and multiplication AB C ABAC7 also holds The powers of a square matrix are de ned by A2 AA7 A3 AA27 and so on Example 54 Let A4212 H gt ABi 2135 2432 7 17 14 71105 71402 7174 A2 lt31 333312 i i ogtlt 63 Next we de ne multiplication of an n X n matrix A times an n vector X The Then Also product AX7 with the matrix on the left is de ned to be the n vector whose ith component is ai Xi In other words the ith element in the list Ax is found by taking the product of the ith row of A and the vector x The product XA is not de ned Example 55 When n 2 we have 168 5 Linear Systems For a numerical example take 2 3 5 4 ogtv H71 7 2537 7 31 AX7lt711507gt7lt75gt The special square matrix having ones on the main diagonal and zeros Then elsewhere else is called the identity matrix and is denoted by 11 For example the 2 X 2 and 3 X 3 identities are 1 0 1 0 0 Ilt01gt and I 010 0 0 1 It is easy to see that if A is any square matrix and I is the identity matrix of the same size then AI IA Al Therefore multiplication by the identity matrix does not change the result a situation similar to multiplying real numbers by the unit number 11 If A is an n gtlt n matrix and there exists a matrix B for which AB BA I then B is called the inverse of A and we denote it by B A li If A 1 exists we say A is a nonsingular matrix otherwise it is called singulari One can show that the inverse of a matrix if it exists is unique We never write 1A for the inverse of A1 A useful number associated with a square matrix A is its determinanti The determinant of a square matrix A denoted by detA also by is a number found by combining the elements of the matrix is a special way The determinant of a 1 gtlt 1 matrix is just the single number in the matrix For a 2 X 2 matrix we de ne detAdetlt a b gt ad7cb c d and for a 3 X 3 matrix we de ne a b c det d e f aez39 bfg th 7 069 7 bdz 7 ahfi 513 g h 239 Example 56 We have 2 6 detlt 72 0 gt 20772 6121 52 Matrices 169 There is a general inductive formula that de nes the determinant of an n X n matrix as a sum of n71 X n7 1 matrices Let A aij be an n X n matrix and let Mij denote the n 7 l X n 7 1 matrix found by deleting the 2th row and jth column of A397 the matrix Mij is called the ij minor of Al Then det A is de ned by choosing any xed column J of A and summing the elements a in that column times the determinants of their minors Mi with an associated sign i7 depending upon location in the column That is for any xed J det A EHVHW dam i1 This is called the expansion by minors formula One can show that you get the same value regardless of which column J you user In fact one can expand on any xed row I instead of a column and still obtain the same value detA 271Ija1j detMIji 11 So the determinant is well de ned by these equations The reader should check that these formulas give the values for the 2gtlt 2 and 3gtlt 3 determinants presented above A few comments are in order First the expansion by minors formulas are useful only for small matricesi For an n X n matrix it takes roughly nl arithmetic calculations to compute the determinant using expansion by mi nors which is enormous when n is larger Ef cient computational algorithms to calculate determinants use row reduction methodsi Both computer algebra systems and calculators have routines for calculating determinants Using the determinant we can give a simple formula for the inverse of a 2 gtlt 2 matrix Al Let a b A and suppose det A a 0 Then 1 d 7b A 1 i 54 det A lt 70 a gt So the inverse of a 2 X 2 matrix is found by interchanging the main diagonal elements putting minus signs on the offdiagonal elements and dividing by the determinant There is a similar formula for the inverse of larger matrices397 for completeness we will write the formula down but for the record we com ment that there are more ef cient ways to calculate the inverse With that said the inverse of an n X n matrix A is the n X n matrix whose ij entry is ily v det Mji divided by the determinant of A which is assumed nonzeroi ln symbols A71 1 i j m il T detMl 5A5 170 5 Linear Systems Note that the ij entry of A 1 is computed from the ji minor with indices transposedi In the 3 X 3 case the formula is det M11 idet M21 detMgl A l d A idet M12 detMgg idet M32 6 det M13 idetMgg detMgg Example 57 If 1 2 lt4 3gt7 then 71713727L37277 idetA 741 75 74 1 7 The reader can easily check that AA l L on I wlw l mum U IH Equations 54 and 55 are revealing because they seem to indicate the inverse matrix exists only when the determinant is nonzero you can t divide by zero In fact these two statements are equivalent for any square matrix regardless of its size A 1 exists if and only if detA y 0 This is a major theoretical result in matrix theory and it is a convenient test for invertibility of small matricesi Again for larger matrices it is more ef cient to use row reduction methods to calculate determinants and inversesi The reader should remember the equivalences Ailexists ltgt A is nonsingular ltgt det A f 0 Matrices were developed to represent and study linear algebraic systems n linear algebraic equations in n unknowns in a concise way For example consider two equations in two unknowns 11 12 given in standard form by a1111a1212 b1 a2111a2212 112 Using matrix notation we can write this as a11a12gtlt11gtltbigt a21 a22 12 b2 7 Ax b7 56 or just simply as 52 Matrices 171 Altan a12gt7 Xlt11gt7 bltb1gt a21 a22 12 122 A is the coef cient matrix x is a column vector containing the unknowns where and b is a column vector representing the right side If b 0 the zero vector then the system 56 is called homogeneous Otherwise it is called nonho mogeneous In a twodimensional system each equation represents a line in the plane When b 0 the two lines pass through the origin A solution vector x is represented by a point that lies on both lines There is a unique solution when both lines intersect at a single point39 there are in nitely many solutions when both lines coincide39 there is no solution if the lines are parallel and dif ferent In the case of three equations in three unknowns each equation in the system has the form azl 612 713 d and represents a plane in space If d 0 then the plane passes through the origin The three planes represented by the three equations can intersect in many ways giving no solution no com mon intersection points a unique solution when they intersect at a single point a line of solutions when they intersect in a common line and a plane of solutions when all the equations represent the same plane The following theorem tells us when a linear system Ax b of n equations in n unknowns is solvable It is a key result that is applied often in the sequel Theorem 58 Let A be an n X n matrixi If A is nonsingular then the system Ax b has a unique solution given by x A lb39 in particular the homogeneous system AX 0 has only the trivial solution x 0 If A is singular then the homoge neous system Ax 0 has in nitely many solutions and the nonhomogeneous system Ax b may have no solution or in nitely many solutions It is easy to show the rst part of the theorem when A is nonsingular using the machinery of matrix notationi If A is nonsingular then A 1 existsi Multiplying both sides of Ax b on the left by A 1 gives 147le Ailb I x Ailb x A lb which is the unique solution If A is singular one can appeal to a geometric ar gument in two dimensions That is if A is singular then det A 0 and the two lines represented by the two equations must be parallel can you show that Therefore they either coincide or they do not giving either in nitely many so lutions or no solution We remark that the method of nding and multiplying 172 5 Linear Systems by the inverse of the matrix A7 as above is not the most ef cient method for solving linear systems Row reduction methods introduced in high school al gebra and reviewed below provide an ef cient computational algorithm for solving large systems Example 59 Consider the homogeneous linear system 4 1 11 7 0 8 2 12 7 0 i The coef cient matrix has determinant zero so there will be in nitely many solutions The two equations represented by the system are 411 12 07 811 212 07 which are clearly not independent397 one is a multiple of the other Therefore we need only consider one of the equations say 411 12 0 With one equation in two unknowns we are free to pick a value for one of the variables and solve for the other Let 11 1 then 12 74 and we get a single solution x 174 More generally if we choose 11 a where a is any real parameter then 12 7401 Therefore all solutions are given by xltgtltiagtaltf4gt7 Thus all solutions are multiples of 17 74T and the solution set lies along the straight line through the origin de ned by this vector Geometrically the two equations represent two lines in the plane that coincidei Next we review the row reduction method for solving linear systems when n 3 Consider the algebraic system AX b or a1111 a1212 a1313 1 17 a2111 a2212 a2313 1 27 5 7 a3111 a3212 3313 113 At rst we assume the coef cient matrix A clj is nonsingular so that the system has a unique solution The basic idea is to transform the system into the simpler triangular form 51111 51212 51313 1 17 52212 523 1 27 aggzg 123 52 Matrices 173 This triangular system is easily solved by back substitution That is the third equation involves only one unknown and we can instantly nd 13 That value is substituted back into the second equation where we can then nd 12 and those two values are substituted back into the rst equation and we can nd 11 The process of transforming 57 into triangular form is carried out by three admissible operations that do not affect the solution structure 1 Any equation may be multiplied by a nonzero constant 2 Any two equations may be interchanged 3 Any equation may be replaced by that equation plus or minus a multiple of any other equation We observe that any equation in the system 517 is represented by its coef cients and the right side so we only need work with the numbers which saves writing We organize the numbers in an augmented array an L112 an b1 L121 L122 L123 b2 L131 L132 ass ha The admissible operations listed above translate into row operations on the augmented array any row may be multiplied by a nonzero constant any two rows may be interchanged and any row may be replaced by itself plus or minus any other rowi By performing these row operations we transform the augmented array into a triangular array with zeros in the lower left corner below the main diagonali The process is carried out one column at a time beginning from the left Example 510 Consider the system 11 12 13 07 211 7 213 2 11 7 12 13 6 The augmented array is 1 1 1 2 0 72 2 1 71 1 Begin working on the rst column to get zeros in the 21 and 31 positions by replacing the second and third rows by themselves plus multiples of the rst 174 5 Linear Systems rowi So we replace the second row by the second row minus twice the rst row and replace the third row by third row minus the rst rowi This gives 1 1 1 0 0 72 74 2 0 72 0 6 Next work on the second column to get a zero in the 32 position below the diagonal entryi Speci cally replace the third row by the third row minus the second row 1 1 1 0 0 72 74 2 0 0 4 4 This is triangular as desired To make the arithmetic easier multiply the third row by 14 and the second row by 712 to get 1 1 1 0 0 1 2 71 7 0 0 1 1 with ones on the diagonal This triangular augmented array represents the system 11 12 13 07 12 213 17 13 11 Using back substitution zg 17 12 73 and 11 2 which is the unique solution representing a point 277371 in R31 If the coef cient matrix A is singular we can end up with different types of triangular forms for example 1 1 1 1 1 1 1 1 1 1 1 1 0 1 7 0 0 6 7or 0 0 0 7 0 0 0 1 0 0 0 1 0 0 0 1 where the 96 denotes an entry These augmented arrays can be translated back into equationsi Depending upon the values of those entries we will get no solu tion the equations are inconsistent or in nitely many solutions As examples suppose there are three systems with triangular forms at the end of the process given by 1 1 3 0 1 0 3 3 1 2 0 1 0 1 2 5 7 0 0 1 1 7or 0 0 0 0 0 0 0 7 0 0 0 0 0 0 0 0 52 Matrices 175 There would be no solution for the rst system the last row states 0 7 and in nitely many solutions for the second and third systemsi Speci cally the second system would have solution 13 1 and 11 0 with 12 a which is arbitrary Therefore the solution to the second system could be written 11 0 0 0 12 a a 1 0 7 13 1 0 1 with a an arbitrary constant This represents a line in Rgi A line is a one dimensional geometrical object described in terms of one parameter The third system above reduced to 11 212 1 So we may pick 13 and 12 arbitrarily say 12 a and 13 b and then 11 17 2a The solution to the third system can then be written 11 17 2a 72 0 1 12 a a 1 12 0 0 7 I3 b 0 1 0 which is a plane in Rgi A plane is a twodimensional object in R3 requiring two parameters for its description The set of all solutions to a homogeneous system Ax 0 is called the nullspace of Al The nullspace may consist of a single point X 0 when A is nonsingular or it may be a line or plane passing through the origin in the case where A is singular Finally we introduce the notion of independence of column vectors A set of vectors is said to be a linearly independent set if any one of them cannot be written as a combination of some of the others We can express this statement mathematically as follows A set 17 of them of nvectors V17V27M7Vp is a linearly independent set if the equation1 C1V1C2V2 CPVP0 forces all the constants to be zero 7 that is cl Cg up 0 If all the constants are not forced to be zero then we say the set of vectors is linearly dependent In this case there would be at least one of the constants say CT which is not zero at which point we could solve for VT in terms of the remaining vectorsi Notice that two vectors are independent if one is not a multiple of the other 1 A sum of constant multiples of a set of Vectors is called a linear combination of those vectors 176 5 Linear Systems In the sequel we also need the notion of linear independence for vector functions A vector function in two dimensions has the form of a 2 vector whose entries are functions of time t397 for example 9W F lt W gt7 where t belongs to some interval I of time The vector function rt is the po sition vector and its arrowhead traces out a curve in the plane given by the parametric equations 1 1t y yt t 6 Il As observed in Section 51 solutions to twodimensional systems of differential equations are vector func tions Linear independence of a set of n vector functions r105 r2 t7 m rp t on an interval I means that if a linear combination of those vectors is set equal to zero for all t E I then the set of constants is forced to be zero In symbols 01r1t02r2tcprpt 07 t6 implies 0107 5207m75p0l Finally if a matrix has entries that are functions of t ie A At cl1396 then we de ne the derivative of the matrix as the matrix of derivatives or Amt aw Example 511 The two vector functions gt7 lt form a linearly independent set on the real line because one is not a multiple of the other Looked at differently if we set a linear combination of them equal to the zero vector ilel 51r1t 62136 0 and take t 0 then cl 502 07 701 07 which forces cl Cg 0 Because the linear combination is zero for all t we may take t 0 Example 512 The three vector functions r1lttgtlt2tgt7 r2lttgtlt gt7 r3tlt3sngt7 52 Matrices 177 form a linearly independent set on R because none can be written as a combi nation of the others That is if we take a linear combination and set it equal to zero 7 that is 51r1t 02r1t 53r1t 07 for all t E R then we are forced into cl Cg Cg 0 see Exercise 15 EXERCISES 13 710 2 H2 4gt7 Hg 7gt7 75 Find A B B 7 4AABBAA2BxABxA 1detBB3AI and detA 7 AI where A is a parameter 1 Let 2 With A given in Exercise 1 and b 21T solve the system Ax b using A ll Then solve the system by row reduction 3 Let 0 2 71 1 71 0 A 7 1 6 72 B 7 2 1 4 2 0 3 71 71 l 2 x 0 71 7 Find A 7 B B 7 4A BA A2 Bx det A A1 A 7 31 and detB 7 I Find all values of the parameter A that satisfy the equation det AiAI 07 where A is given in Exercise 1 2 71 A l lt 74 2 l Compute det Al Does A 1 exist Find all solutions to Ax 0 and plot the solution set in the plane F 5 Let 6 Use the row reduction method to determine all values m for which the algebraic system 2z3ym7 76179y5 has no solution a unique solution or in nitely many solutions 7 Use row reduction to determine the values of m for which the following system has in nitely many solutions I y 07 21 y 07 31 2y m2 0 178 5 Linear Systems H H H H 4 HH i If a square matrix A has all zeros either below its main diagonal or above its main diagonal show that detA equals the product of the elements on the main diagonal i Construct simple homogeneous systems Ax 0 of three equations in three unknowns that have a a unique solution b an in nitude of solutions lying on a line in R3 and c an in nitude of solutions lying on a plane in R3 Is there a case when there is no solution Let 0 2 71 A l 6 72 2 0 3 a Find det A by the expansion by minors formula using the rst column the second column and the third row Is A invertible Is A singular b Find the inverse of A and use it to solve Ax b where b l 0 4Ti c Solve Ax b in part b using row reduction Find all solutions to the homogeneous system Ax 0 if 72 0 2 A 2 74 0 0 4 72 i Use the de nition of linear independence to show that the 2 vectors 2 73 and 74 8T are linearly independent i Use the de nition to show that the 3vectors 010T 120T and 014T are linearly independent i Use the de nition to show that the 3vectors 101T 5710T and 7712T are linearly dependent i Verify the claim in Example 512 by taking two special values of ti Plot each of the following vector functions in the my plane where 700 lt tlt00i gt r2tltgtt r3tltt1gtequot r1t Show that these vector functions form a linearly independent set by setting clrl t 02r1 t 03r1t 0 and then choosing special values oft to force the constants to be zero 3cost 2sint 53 TwoDimensional Systems 179 17 Show that a 3 X 3 matrix A is invertible if and only if its three columns form an independent set of 3vectorsi 18 Find A t if cos t t2 0 At 262 l sin2t 0 T51 53 TwoDimensional Systems 531 Solutions and Linear Orbits A twodimensional linear system of differential equations m amby y cm dy Where a b c and d are constants can be Written compactly using vectors and went Alt gt7 the system can be written 1E gtlt Zgtlt3Efigtv Xt Z Zgtxti We often Write this simply as matricesi Denoting x Ax 58 Where we have suppressed the understood dependence of x on t We brie y reiterate the ideas introduced in the introduction Section Sill A solution to the system 58 on an interval is a vector function Xt mtytT that satis es the system on the required intervali We can graph mt and yt vsi t Which gives the state space representation or time series plots of the solu tion Alternatively a solution can be graphed as a parametric curve or vector function in the my plane We call the my plane the phase plane and we call a solution curve plotted in the my plane an orbit Observe that a solution is a vector function xt With components mt and In the phase plane the 180 5 Linear Systems orbit is represented in parametric form and is traced out as time proceedsi Thus time is not explicitly displayed in the phase plane representation but it is a parameter along the orbit An orbit is traced out in a speci c direction as time increases and we usually denote that direction by an arrow along the curve Furthermore time can always be shifted along a solution curvei That is if xt is a solution then Xt 7 c is a solution for any real number 5 and it represents the same solution curvei Our main objective is to nd the phase portrait or a plot of key orbits of the given systemi We are particularly interested in the equilibrium solutions of 58 These are the constant vector solutions x for which Ax 0 An equilibrium solution is represented in the phase plane as a point The vector eld vanishes at an equilibrium point The time series representation of an equilibrium solution is two constant functions If detA 0 then X 0 is the only equilibrium of 58 and it is represented by the origin 00 in the phase plane We say in this case that the origin is an isolated equilibriumi lf det A 07 then there will be an entire line of equilibrium solutions through the origin397 each point on the line represents an equilibrium solution and the equilibria are not isolated Equilibrium solutions are important because the interesting behavior of the orbits occurs near these solutions Equilibrium solutions are also called critical points by some authors Example 513 Consider the system I 72ziy y 2175y7 x i 7271 7 2 75 The coef cient determinant is nonzero so the only equilibrium solution is rep which we write as resented by the origin zt 0 yt 0 By substitution it is straightforward to check that x1lttgtlt 83 gt lt 3 gt lt 1 w is a solution Also 6741 l x2lttgtlt2ei4zgtlt2gte 53 TwoDimensional Systems 181 is a solution Each of these solutions has the form of a constant vector times a scalar exponential function of time t Why should we expect exponential solutions The two equations involve both I and y and their derivatives397 a solution must make everything cancel out and so each term must basically have the same formi Exponential functions and their derivatives both have the same form and therefore exponential functions for both I and y are likely candidates for solutions We graph these two independent solutions X1t and x2 t in the phase planer See gure 55 Each solution or orbit plots as a ray traced from in nity as time t approaches 700 into the origin as t approaches 00 The slopes of these raylike solutions are de ned by the constant vectors preceding the scalar exponential factor the latter of which has the effect of stretching or shrinking the vector Note that these two orbits approach the origin as time gets large but they never actually reach it Another way to look 4 20 Figure 55 X1t and XQt are shown as linear orbits rays entering the origin in the rst quadranti The re ection of those rays in the third quadrant are the solutions 7x1t and 7X2ti Note that all four of these linear orbits approach the origin as t e 00 because of the decaying exponential factor in the solution As t 4 foo backward in time all four of these linear orbits go to in nity at it is this If we eliminate the parameter t in the parametric representation 1 6 y 26 of X2 t7 say then y 21 which is a straight line in the my plane This orbit is on one ray of this straight line lying in the rst quadranti Solutions of 58 the form xt VeM where A is a real constant and V is 182 5 Linear Systems a constant real vector are called linear orbits because they plot as rays in the zy phase planet We are ready to make some observations about the structure of the solution set to the twodimensional linear system 58 All of these properties can be extended to three or even it dimensional systems 1 Superposition lfx1t and X2 t are any solutions and cl and 52 are any constants then the linear combination 51x1 t 52x t is a solution 2 General Solution If x1t and xt are two linear independent solu tions ie one is not a multiple of the other then all solutions are given by xt 01x1t CQXQ t where cl and Cg are arbitrary constantsi This combination is called the general solution of 5 8 3 ExistenceUniqueness The initial value problem X AX7 X050 X07 where x0 is a xed vector has a unique solution valid for all 700 lt t lt 00 The existenceuniqueness property actually guarantees that there are two independent solutions to a twodimensional systemi Let x1 be the unique so lution to the initial value problem X 1 AX17 x10 10T and X2 be the unique solution to the initial value problem x AXQ7 X007 1Ti These must be independent Otherwise they would be proportional and we would have KN 7090 for all t where k is a nonzero constantl But if we take t 0 we would have lt170gtT M071 which is a contradiction The question is how to determine two independent solutions so that we can obtain the general solution This is a central issue we address in the sequel One method to solve a twodimensional linear system is to eliminate one of the variables and reduce the problem to a single secondorder equationl Example 514 Method of Elimination Consider z 4x 7 By M 6zi7yl 53 TwoDimensional Systems 183 Differentiate the rst and then use the second to get I 41 7 3y 441 7 3y 7 361 7 7y 7219y72197z z 731 1017 which is a secondorder equation The characteristic equation is A23A710 0 with roots A 7572 Thus zt claim 02623 Then 1 4 2 yt 7gz g1 3616751 gage2 We can write the solution in vector form as xlttgtltsgt7lttgtczlt2zgtv In this form we can see that two independent vector solutions are X10 7 lt 3 7 W 7 lt gt and the general solution is a linear combination of these xt 01x1 t 52X2ti However simple this strategy appears in two dimensions it does not work as easily in higher dimensions nor does it expose methods that are easily adaptable to higherdimensional systemsi Therefore we do not often use the elimination methodi But we point out features of the phase planer Notice that X1 graphs as a linear orbit in the rst quadrant of the my phase plane along the ray de ned by of the vector 13 It enters the origin as t 7 0 because of the decaying exponential factor The other solution X2 also represents a linear orbit along the direction de ned by the vector 1723 This solution because of the increasing exponential factor 62 tends to in nity as t 7 00 Figure 56 shows the linear orbitsi Figure 57 shows several orbits on the phase diagram obtained by taking different values of the arbitrary constants in the general solution The structure of the orbital system near the origin where curves veer away and approach the linear orbits as time goes forward and backward is called a saddle point structure The linear orbits are sometimes called separatrices because they separate different types of orbitsi All orbits approach the separatrices as time gets large either negatively or positivelyi 184 5 Linear Systems 4 10 Figure 56 Linear orbits in Example 514 representing the solutions corre sponding to x1t and X2t and the companion orbits 7X10 and 7x2tl These linear orbits are called separatricesl Figure 57 Phase portrait for the system showing a saddle point at the origin 53 TwoDimensional Systems 185 532 The Eigenvalue Problem Now we introduce some general methods for the twodimensional system x Axi 59 We assume that det A y 0 so that the only equilibrium solution of the system is at the origin As examples have shown we should expect an exponentialtype solutioni Therefore we attempt to nd a solution of the form x vek 510 where A is a constant and V is a nonzero constant vector both to be determined Substituting x V6 and XAVeM into the 59 gives AveM AVe 7 AV AVi 511 Therefore if a A and V can be found that satisfy 511 then we will have determined a solution of the form 510 The vector equation 511 represents a wellknown problem in mathematics called the algebraic eigenvalue prob lemi The eigenvalue problem is to determine values of A for which 511 has a nontrivial solution V A value of A for which there is a nontrivial solution V is called an eigenvalue and a corresponding V associated with that eigenvalue is called an eigenvector The pair A7 V is called an eigenpairi Geometrically we think of the eigenvalue problem like this A represents a transformation that maps vectors in the plane to vectors in the plane397 a vector x gets transformed to a vector Axi An eigenvector of A is a special vector that is mapped to a multiple A of itself397 that is Az Azi In summary we have reduced the prob lem of nding solutions to a system of differential equations to the problem of nding solutions of an algebra problem every eigenpair gives a solution Geometrically if A is real the linear orbit representing this solution lies along a ray emanating from the origin along the direction de ned by the vector Vi If A lt 0 the solution approaches the origin along the ray and if A gt 0 the solution goes to in nity along the ray The situation is similar to that shown in gure 56 When there is a solution graphing as a linear orbit then there is automatically a second opposite linear orbit along the ray 7V This is because M is a solution then so is 7x 7V6 if x Ve To solve the eigenvalue problem we rewrite 511 as a homogeneous linear system A7A1v 0 512 186 5 Linear Systems By Theorem 58 this system will have the desired nontrivial solutions if the determinant of the coef cient matrix is zero or detA7AI or 513 Written out explicitly this system 512 has the form a 7 A 12 v1 7 0 c d 7 A 112 7 0 7 where the coef cient matrix A7AI is the matrix A with A subtracted from the diagonal elements Equation 513 is explicitly detlt a de gt a7Ad7A7cb0 or equivalently A2 7 a bA ad 7125 0 This last equation can be memorized easily if it is written A2 7 trAA detA 0 514 where trA ad is called the trace of A de ned to be the sum of the diagonal elements of Al Equation 514 is called the Characteristic equation associ ated with A and it is a quadratic equation in Al lts roots found by factoring or using the quadratic formula are the two eigenvaluesl The eigenvalues may be real and unequal real and equal or complex conjugatesl Once the eigenvalues are computed we can substitute them in turn into the system 512 to determine corresponding eigenvectors Vi Note that any multiple of an eigenvector is again an eigenvector for that same eigenvalue397 this follows from the calculation Acv cAv CAv Acvl Thus an eigenvector corresponding to a given eigenvalue is not unique397 we may multiply them by constantsl This is expected from Theorem 58 Some calculators display normalized eigenvectors of length one found by dividing by their length As noted the eigenvalues may be real and unequal real and equal or com plex numbers We now discuss these different cases 53 TwoDimensional Systems 187 533 Real Unequal Eigenvalues If the two eigenvalues are real and unequal say A1 and A27 then correspond ing eigenvectors v1 and V2 are independent and we obtain two independent solutions Vleklt and VgeM i The general solution of the system is then a linear combination of these two independent solutions Alt Agt CQVQE 7 xt clvle where cl and Cg are arbitrary constants Each of the independent solutions represents linear orbits in the phase plane which helps in plotting the phase diagrami All solutions orbits xt are linear combinations of the two inde pendent solutions with each speci c solution obtained by xing values of the arbitrary constants Example 515 Consider the linear system NIH 3 5 x lt 1 71 gt x 5115 The characteristic equation 5114 is 2 5 A 5A4 1 01 By the quadratic formula the eigenvalues are l A 75 721 Now we take each eigenvalue successively and substitute it into 5112 to obtain corresponding eigenvectorsi First for A 7 we get 71 v1 7 0 1 7 v2 0 7 which has a solution v1 v2T 12T1 Notice that any multiple of this eigen vector is again an eigenvector but all we need is one Therefore an eigenpair is 1 1 2 2 Now take A 721 The system 5112 becomes lt11gtltgtlt2gt7 188 5 Linear Systems which has solution v17 v2T 711 Thus another eigenpair is 473 The two eigenpairs give two independent solutions x1t lt gt6 and x2tlt 1 1 gt523 516 Each one plots along with its negative counterparts as a linear orbit in the phase plane entering the origin as time increases The general solution of the system 515 is l 7 fl XtC1lt2gt t2C2lt1 gt623 This is a twoparameter family of solution curves and the totality of all these solution curves or orbits represents the phase diagram in the my plane These orbits are shown in gure 58 Because both terms in the general solution decay as time increases all orbits enter the origin as t a ltxgti And as t gets large the term with 642 dominates the term with e Q i Therefore all orbits approach the origin along the direction 12 As 75 4 700 the orbits go to in nity397 for large negative times the term 6 2 dominates the term 6quot and the orbits become parallel to the direction 711 Each of the two basic solutions 516 represents linear orbits along rays in the directions of the eigenvectors When both eigenvalues are negative as in this case all orbits approach the origin in the direction of one ofthe eigenvectorsi When we obtain this type of phase plane structure we call the origin an asymptotically stable nodei When both eigenvalues are positive then the time direction along the orbits is reversed and we call the origin an unstable model The meaning of the term stable is discussed subsequentlyi An initial condition picks out one of the many orbits by xing values for the two arbitrary constantsi For example ifx0 14T or we want an orbit passing through the point 14 then Clltgtc2ltflgtltigtv giving cl 53 and 02 23 Therefore the unique solution to the initial value ltigtw2 ltigt 5 3 gait2 7 672 10642 672 problem is 53 TwoDimensional Systems 189 Figure 58 A node All orbits approach the origin tangent to the direction 17 2 as t 4gt oo Backwards in time as t a 700 the orbits become parallel to the direction 717 71 Notice the linear orbits Example 516 If a system has eigenpairs 27lt2gt7 24 with real eigenvalues of opposite sign then the general solution is xt cllt 3 gte2 52lt 1 gt633 In this case one of the eigenvalues is positive and one is negative Now there are two sets of opposite linear orbits one pair corresponding to 72 approaching the origin from the directions i37 2T and one pair corresponding to A 3 approaching in nity along the directions i7l5T The orbital structure is that of a saddle point refer to gure 57 and we anticipate saddle point structure when the eigenvalues are real and have opposite sign 534 Complex Eigenvalues If the eigenvalues of the matrix A are complex they must appear as complex conjugates or A a i bi The eigenvectors will be V W i iz Therefore 190 5 Linear Systems taking the eigenpair a bi7 W iz7 we obtain the complex solution W ize biti Recalling that the real and imaginary parts of a complex solution are real solutions we expand this complex solution using Euler s formula to get W ize quotem e quotW iz cos bt isin bt eatW cos bt 7 2 sin bt ieatW sin bt 2 cos ht Therefore two real independent solutions are x1t e quot W cos bt 7 2 sin bt7 X2 t e quotW sin bt 2 cos bt7 and the general solution is a combination of these Xt cleatW cos bt 7 Z sin bt 026 quot W sin bt 2 cos bit 517 In the case of complex eigenvalues we need not consider both eigenpairs397 each eigenpair leads to the same two independent solutions For complex eigenvalues there are no linear orbitsi The terms involving the trigonometric functions are periodic functions with period 27rb7 and they de ne orbits that rotate around the origin The factor 6 quot acts as an amplitude factor causing the rotating orbits to expand if a gt 0 and we obtain spiral orbits going away from the origin If a lt 0 the amplitude decays and the spiral orbits go into the origin In the complex eigenvalue case we say the origin is an asymptotically stable spiral point when a lt 0 and an unstable spiral point when a gt 0 If the eigenvalues of A are purely imaginary A ibi then the amplitude factor 6 quot in 517 is absent and the solutions are periodic of period QT given by xt cl W cos bt 7 2 sin bt 52 W sin bt 2 cos bit The orbits are closed cycles and plot as either concentric ellipses or concentric circles In this case we say the origin is a neutrally stable centeri Example 517 Let x 7 72 73 x 3 72 i The matrix A has eigenvalues A 72 i 32quot An eigenvector corresponding to A 72 3239 is V1 71 i Therefore a complex solution is N ltf1gte 23it 61 gt ilt 1 e 2tcos3tisin3t z 7 76 2 cos 3t 2 76 2 sin 3t 7 76 2 sin 3t 76 2 cos 3t 53 TwoDimensional Systems Therefore two linearly independent solutions are 72t 72t 6 cos 3t 6 s1n 3t x t x t i 1 lt 7672 sin3t gt7 2 lt fei2tcos3t gt The general solution is a linear combination of these two solutions xt clxlt 02x2 In the phase plane the orbits are spirals that approach the origin as t e 00 because the real part 72 of the eigenvalues is negative See gure 59 At the point 11 the tangent vector direction eld is 751 so the spirals are counterclockwise Figure 59 A stable spiral from Example 517 535 Real Repeated Eigenvalues One case remains when A has a repeated real eigenvalue A with a single eigen vector Vi Then x1 V6 is one solution representing a linear orbit and we need another independent solution We try a second solution of the form X2 e tv W7 where W is to be determined A more intuitive guess based on our experience with secondorder equations in Chapter 3 would have been eMtV7 but that does not work try it Substituting x2 into the system we get x e ve tvw7 Ax2 eA Atvwi 192 5 Linear Systems Therefore we obtain an algebraic system for W A 7 ADW v This system will always have a solution W and therefore we will have deter mined a second linearly independent solution In fact this system always has in nitely many solutions and all we have to do is nd one solution The vector w is called a generalized eigenvector Therefore the general solution to the linear system X Ax in the repeated eigenvalue case is xt elven CgeMOV w If the eigenvalue is negative the orbits enter the origin as t e 00 and they go to in nity as t4 fool If the eigenvalue is positive the orbits reverse direction in time In the case where the eigenvalues are equal the origin has a nodallike structure as in Example 513 When there is a single eigenvector associated with the repeated eigenvalue we often call the origin a degenerate node It may occur in a special case that a repeated eigenvalue A has two independent eigenvectors vectors V1 and V2 associated with it When this occurs the general solution is just xt clvleh CQVQEAti It happens when the two equations in the system are decoupled and the matrix is diagonal with equal elements on the diagonal In this exceptional case all of the orbits are linear orbits entering A lt 0 or leaving A gt 0 the origin397 we refer to the origin in this case as a starlike nodei Example 518 Consider the system x 7 21x 7714 The eigenvalues are A 3 3 and a corresponding eigenvector is V 11Ti x1t lt 1 gt633 Notice that this solution plots as a linear orbit coming out of the origin and Therefore one solution is approaching in nity along the direction 17 lTi There is automatically an oppo site orbit coming out of the origin and approaching in nity along the direction 7171 A second independent solution will have the form X2 egttv W where W satis es 71 l 1 A731Wlt711gtWlt1gti 53 TwoDimensional Systems 193 This equation has many solutions and so we choose W 7 0 7 1 Therefore a second solution has the form x2t e3 tvw e3 1 gttlt 1 lt Deg The general solution of the system is the linear combination X05 01X1t C2X2t If we append an initial condition for example xlt0gtlt3gt then we can determine the two constants cl and 02 We have 1 0 1 X0C1X10C2X20 cllt 1 gt52lt 1 gt lt 0 cl 17 Cg 71 Hence Therefore the solution to the initial value problem is given by W lt1gt lt i gte lt1gtlt tine gt lt 2696 gt As time goes forward t a 0 the orbits go to in nity and as time goes backward t a 700 the orbits enter the origin The origin is an unstable nodei How to draw a phase diagrami In general to draw a rough phase dia gram for a liner system all you need to know are the eigenvalues and eigenvec torsi If the eigenvalues are real then draw straight lines through the origin in the direction of the associated eigenvectorsi Label each ray of the line with an arrow that points inward toward the origin ifthe eigenvalue is negative and out ward if the eigenvalue is positive Then ll in the regions between these linear orbits with consistent solution curves paying attention to which eigendirec tion dominates as t 4 gt0 and t a 700 Real eigenvalues with the same sign give nodes and real eigenvalues of opposite signs give saddlesilf the eigenvalues are purely imaginary then the orbits are closed loops around the origin and if they are complex the orbits are spiralsi They spiral in if the eigenvalues have 194 5 Linear Systems negative real part and they spiral out if the eigenvalues have positive real part The direction clockwise or counterclockwise of the cycles or spirals can be de termined directly from the direction eld often by just plotting one vector in the vector eld Another helpful device to get an improved phase diagram is to plot the set of points where the vector eld is vertical the orbits have a vertical tangent and where the vector eld is horizontal the orbits have a horizontal tangentl These sets of points are found by setting zprime az by 0 and yprime or dy 0 respectively These straight lines are called the z and ynullclinesl Example 519 The system x 7 25x 7 720 has eigenvalues l i 32quot The orbits spiral outward because the real part of the eigenvalues l is positive They are clockwise because the second equation in the system is y 721 and so y decreases y lt 0 when I gt 0 Observe that the orbits are vertical as they cross the nullcline 21 5y 0 and they are horizontal as they cross the nullcline z 0 With this information the reader should be able to draw a rough phase diagraml 536 Stability We mentioned the word stability in the last section Now we extend the discus sion For the linear system x Ax an equilibrium solution is a constant vector solution xt x representing a point in the phase planet The zero vector x 0 the origin is always an equilibrium solution to a linear systeml Other equilibria will satisfy Ax 0 and thus the only time we get a nontrivial equilibrium solution is when detA 039 in this case there are in nitely many equilibrial lf detA y 0 then xquot 0 is the only equilibrium and it is called an isolated equilibriuml For the discussion in the remainder of this section we assume det A y 0 Suppose the system is in its zero equilibrium state lntuitively the equilib rium is stable if a small perturbation or disturbance does not cause the system to deviate too far from the equilibrium39 the equilibrium is unstable if a small disturbance causes the system to deviate far from its original equilibrium state We have seen in twodimensional systems that if the eigenvalues of the matrix A are both negative or have negative real parts then all orbits approach the origin as t 4gt ltxgtl In these cases we say that the origin is asymptotically 53 TwoDimensional Systems 195 stable node including degenerate and starlike nodes or an asymptotically stable spiral point If the eigenvalues are both positive have positive real parts or are real of opposite sign then some or all orbits that begin near the origin do not stay near the origin as t 7 00 and we say the origin an un stable nodeindexnode unstable including degenerate and starlike nodes an unstable spiral point and a saddle respectively If the eigenvalues are purely imaginary we obtain periodic solutions or closed cycles and the origin is a center In this case a small perturbation from the origin puts us on one of the elliptical orbits and we cycle near the origin39 we say a center is neu trally stable or just stable but not asymptotically stable Asymptotically stable equilibria are also called attractors or sinks and unstable equilibria are called repellers or sources Also we often refer to asymptotically sta ble spirals and nodes as just stable spirals and nodes7the word asymptotic is understood Summary We make some important summarizing observations that should be remembered for future use Chapter 6 For twodimensional systems it is easy to check stability of the origin and sometimes this is all we want to do The eigenvalues are roots of the characteristic equation A2 7 trA detA 0 By the quadratic formula A gum i trA2 7 4detAi One can easily checks the following facts H lf det A lt 0 then the eigenvalues are real and have opposite sign and the origin is a saddle 0 i If det A gt 0 then the eigenvalues are real with the same sign nodes or complex conjugates centers and spirals Nodes have trA2 7 4detA gt 0 and spirals have trA2 7 4det A lt 0 If trA2 7 4det A 0 then we obtain degenerate and starlike nodesi lf trA lt 0 then the nodes and spirals are stable and if trA gt 0 they are unstable lf trA 0 we obtain centersi 9 lf det A 0 then at least one of the eigenvalues is zero and there is a line of equilibriai An important result is that the origin is asymptotically stable and only if trA lt 0 and detA gt 0 196 5 Linear Systems EXERCISES 1 Find the eigenvalues and eigenvectors of the following matrices 71 4 2 3 2 78 s H4 a H1 72gt 2 Write the general solution of the linear system x Ax if A has eigenpairs 27 15T and 73 274 Sketch the linear orbits in the phase plane corresponding to these eigenpairsi Find the solution curve that satis es the initial condition x0 01T and plot it in the phase planer Do the same for the initial condition X0 76712 CA3 i Answer the questions in Exercise 2 for a system Whose eigenpairs are 76 172T and 71 175 4 For each system nd the general solution and sketch the phase portraiti Indicate the linear orbits if any and the direction of the solution curvesi aXlt1 bXlt 3 43gtx gtx dXlt 25 30gt we fXlti igtxi we 14gtx hgtxlt01 5i Solve the initial value problem Hi agtxy mac 53 TwoDimensional Systems 197 H H H 6 5 0 to Consider the system x 7 1 72 x 7 72 4 i a Find the equilibrium solutions and plot them in the phase planer b Find the eigenvalues and determine if there are linear orbitsi c Find the general solution and plot the phase portraiti i Determine the behavior of solutions near the origin for the system Xi 3a 711 for different values of the parameter a x For the systems in Exercise 4 characterize the origin as to type node de generate node starlike node center spiral saddle and stability unstable neutrally stable asymptotically stable i Consider the system 1 731ay y b172yi Are there values of a and b Where the solutions are closed cycles periodic orbits In an individual let I be the excess glucose concentration in the blood and y be the excess insulin concentration positive I and y values measure the concentrations above normal levels and negative values measure concen trations below normal levels These quantities are measured in mg per ml and insulin units per ml respectively and time is given in hours One simple model of glucosePinsulin dynamics is I 7am 7 by y or 7 dy Where 7am is the rate glucose is absorbed in the liver and 7by is the rate it is used in the muscle The rate or is the rate insulin is produced by the pancreas and idy is the rate degraded by the liver A set of values for the constants is a 3 b 4 3 5 0 2 and d 0i8i1f101 and y0 0 nd the glucose and insulin concentrations and graph time series plots over a 4 hour period i Find a twodimensional linear system Whose matrix has eigenvalues A 72 andA73i 198 5 Linear Systems 12 Rewrite the damped springmass equation mm 01 k1 0 as a sys tem of two rstorder equations for z and y Iquot Find the characteristic equation of the matrix for the system and show that it coincides with the characteristic equation associated with the secondorder DEi H 9 Consider an RCL circuit governed by LOU RC1 v 07 where v is the voltage on the capacitor Rewrite the equation as a twodimensional linear system and determine conditions on the constants R L and C for which the origin is an asymptotically stable spirali To what electrical response vt does this case correspond H g i What are the possible behaviors depending on 7 of the solutions to the linear system I VI 7 y7 y ziwv H 01 i Show that A 1 exists if and only if zero is not an eigenvalue of Al 4 on i For a 2 X 2 matrix show that the product of the two eigenvalues equals its determinant and the sum of the two eigenvalues equals its tracer H 7 For a 2 gtlt 2 matrix A of a linear system let p equal its trace and 4 equal its determinanti Sketch the set of points in the pqplane where the system has an asymptotically stable spiral at the origin Sketch the region where it has a saddle points 54 Nonhomogeneous Systems CUL quot to a t J39 39 1 linear system X Ax we now examine the nonhomogeneous system x Axft7 518 w lt 338 gt is a given vector function We think of this function as the driving force in the where systemi To ease the notation in writing the solution we de ne a fundamental matrix t as a 2 X 2 matrix whose columns are two independent solutions to the associated homogeneous system X Axi So the fundamental matrix is a 54 Nonhomogeneous Systems 199 square array that holds both vector solutions It is straightforward to show that t satis es the matrix equation t 14 t7 and that the general solution to the homogeneous equation x AX can therefore be written in the form xht tc7 T is an arbitrary constant vector The reader should do where C 5102 Exercise 1 presently which requires verifying these relations The variation of constants method introduced in Chapter 2 is applicable to a rstorder linear systemi Therefore we assume a solution to 518 of the form xt tct7 519 where we have varied the constant vector Ci Then using the product rule for differentiation which works for matrices x t tc t tct tc t A tct Ax ft A tct fti Comparison gives tc t ft or c t trl t we can invert the fundamental matrix because its determinant is nonzero a fact that follows from the independence of its columnsi Integrating the last equation from 0 to t then gives Ct Ot s lfsdsk7 where k is a arbitrary constant vectori Note that the integral of a vector func tions is de ned to be the vector consisting of the integrals of the components Substituting into 519 shows that the general solution to the nonhomogeneous equation 518 is xt tk t Z germane 520 As for a single rstorder linear DE this formula gives the general solution of 518 as a sum of the general solution to the homogeneous equation rst term and a particular solution to the nonhomogeneous equation second term Equation 520 is called the variation of parameters formula for systems It is equally valid for systems of any dimension with appropriate size increase in the vectors and matricesi It is sometimes a formidable task to calculate the solution 520 even in the twodimensional case It involves nding the two independent solutions to the homogeneous equation forming the fundamental matrix inverting the fundamental matrix and then integratingi 200 5 Linear Systems Example 520 Consider the nonhomogeneous system Hi we It is a straightforward exercise to nd the solution to the homogeneous system x lt 4 3 gtx 71 0 i The eigenpairs are 1 171T and 3 731 Therefore two independent solutions are e 736 all l L 7363 t e 0 6 St 1 1 e 363 1 63 363 1 6quot Sequot tdet439gt e e 72 e e 75 e gt 3 A fundamental matrix is and its inverse is By the variation of parameters formula 5120 the general solution is e 7363 t 1 6 5 36 s 0 tkltiez eSt gt0 lt 735 6735 gt lt S gtd5 1 et 736 1 336 s tk lt 76 63 gt0 lt 3673 gtd5 1 et 736 3 f se sds 43031quot 5 lt 76 63 gt lt fol 864st 1 et 736 gt lt 37 3t1e t gt 43 t k 7 2lt7et eSt 7673t kle 7 31626 t 7k1 tk2 31 7amp7 1f the nonhomogeneous term ft is relatively simple we can use the method of undetermined coef cients judicious guessing introduced for secondorder equations in Chapter 3 to nd the particular solution In this case we guess a particular solution depending upon the form of fti For example if both xt 54 Nonhomogeneous Systems 201 components are polynomials then we guess a particular solution with both components being polynomials that have the highest degree that appears If folt 2 then a guess for the particular solution would be a1t2 blt cl mmlt gt a2t2 bgt Cg Substitution into the nonhomogeneous system then determines the six con stants Generally if a term appears in one component of ft then the guess must have that term appear in all its components The method is successful on forcing terms with sines cosines polynomials exponentials and products and sums of those The rules are the same as for single equations But the calculations are tedious and a computer algebra system is often preferred Example 521 We use the method of undetermined coef cients to nd a particular solution to the equation in Example 520 The forcing function is 37 and therefore we guess a particular solution of the form X 7 at b p 7 ct d Substituting into the original system yields CH 3gtltgtlt3gtA Simplifying leads to the two equations a 4a Sct 412 3d 5 7b l 7 ati Comparing coef cients gives 4 l3 al7 127c 7 diiyi 202 5 Linear Systems EXERCISES 1 0 9 g Let 13 x2lt 3gt be independent solutions to the homogeneous equation x AX and let amt MO gt amt 1M0 be a fundamental matrixi Show by direct calculation and comparison of entries that t A ti Show that the general solution of the homoge neous system can be written equivalently as m lt 51x1 52x tc7 where C 5102T is an arbitrary constant vector 1 Two lakes of volume V1 and V2 initially have no contaminationi A toxic chemical flows into lake 1 at q T gallons per minute with a concentration 5 grams per gallon From lake 1 the mixed solution flows into lake 2 at 4 gallons per minute while it simultaneously flows out into a drainage ditch at 7 gallons per minute In lake 2 the the chemical mixture flows out at 4 gallons per minute If I and y denote the concentrations of the chemical in lake 1 and lake 2 respectively set up an initial value problem whose solu tion would give these two concentrations draw a compartmental diagrami What are the equilibrium concentrations in the lakes if any Find zt and Now change the problem by assuming the initial concentration in lake 1 is 10 and fresh water flows in Write down the initial value problem and qualitatively without solving describe the dynamics of this problem using eigenvalues Solve the initial value problem 1 7 3 71 z 1 y 7 1 1 y 2 7 1 t as 7 lt t is a fundamental matrix 10 ND gt7 lt gtew Hi i Solve the problem in Exercise 3 using undetermined coef cients to nd a particular solution o1 Nonhomogeneous Systems U1 7 5 0 i Consider the nonhomogeneous equation it 0 75 2 x 7 3 x e 7 710 Find the fundamental matrix and its inverse Find a particular solution to the system and the general solution In pharmaceutical studies it is important to model and track concentrations of chemicals and drugs in the blood and in the body tissues Let I and y denote the amounts in milligrams of a certain drug in the blood and in the tissues7 respectivelyi Assume that the drug in the blood is taken up by the tissues at rate T11 and is returned to the blood from the tissues at rate Tgyi At the same time the drug amount in the blood is continuously degraded by the liver at rate ngi Argue that the model equations which govern the drug amounts in the blood and tissues are I irlz 7 T31 my 2 mm 7 Tgyi Find the eigenvalues of the matrix and determine the response of the system to an initial dosage of 10 Io given intravenously with y0 0i Hint show both eigenvalues are negative In the preceding problem assume that the drug is administered intra venously and continuously at a constant rate Di What are the governing equations in this case What is the amount of the drug in the tissues after a long time An animal species of population P Pt has a per capita mortality rate mi The animals lay eggs at a rate of 12 eggs per day per animal The eggs hatch at a rate proportional to the number of eggs E Et397 each hatched egg gives rise to a single new animali a Write down model equations that govern P and E and carefully de scribe the dynamics of the system in the two cases I gt m and b lt mi U V Modify the model equations if at the same time an eggieating preda tor consumes the eggs at a constant rate of 7 eggs per day Solve the model equations in part b when I gt m and discuss the dynamics How would the model change if each hatched egg were multiyolked and gave rise to y animals 204 5 Linear Systems 55 ThreeDimensional Systems In this section we give some examples of solving three linear differential equa tions in three unknowns The method is the same as for twodimensional sys tems but now the matrix A for the system is 3 X 3 and there are three eigenval ues and so on We assume det A y 0 Eigenvalues A are found from the charac teristic equation detA 7 AI 07 which when written out is a cubic equation in A For each eigenvalue A we solve the homogeneous system A 7 ADV 0 to determine the associated eigenvectorsi We will have to worry about real complex and equal eigenvalues as in the twodimensional case Each eigenpair A7V gives a solution V6 which if A is real is a linear orbit lying on a ray in R3 in the direction de ned by the eigenvector Vi We need three independent solutions X1t7X2 tX3t to form the general solution which is the linear combination xt clxlt C2X2t 03X3t of those If all the eigenvalues are real and unequal then the eigenvectors will be independent and we will have three independent solutions397 this is the easy casei Other cases such as repeated roots and complex roots are discussed in the examples and in the exercises If all the eigenvalues are negative or have negative real part then all so lution curves approach 000 and the origin is an asymptotically stable equi libriumi If there is a positive eigenvalue or complex eigenvalues with positive real part then the origin is unstable because there is at least one orbit reced ing from the origin Threedimensional orbits can be drawn using computer software but the plots are often dif cult to visualize Examples illustrate the key ideas and we suggest the reader work through the missing detailsi Example 522 Consider the system I1 11 12 13 12 211 12 7 13 13 7811 7 512 7 313 with matrix 1 1 1 A 2 1 71 55 ThreeDimensional Systems 205 Eigenpairs of A are given by 73 74 0 71 4 72 5 2 1 2 7 71 73 74 x1 4 equot7 X2 5 e722 X3 em 2 7 71 Each represents a linear orbit The general solution is a linear combination of these three397 that is xt clxlt C2X2lttgt 53X3t1 The origin is unstable because of the positive eigenvalue Example 523 Consider N l l 1 0 2 0 3 0 xi 2 0 1 The eigenvalues found from det A 7 AI 07 are A 717 37 31 An eigenvector corresponding to A 71 is 17 07 71T and so is one solution To nd eigenvectors corresponding to the other eigenvalue a double root we form A 7 31V 07 or 72 0 2 v1 0 0 0 0 v2 0 2 0 72 v3 0 This system leads to the single equation v1 7 v3 07 With v2 arbitrary Letting v2 and v1 a we can Write the solution as as v1 1 0 112 a 0 8 1 7 v3 1 0 206 5 Linear Systems where a and B are arbitrary Therefore there are two independent eigenvectors associated with A 3 This gives two independent solutions Therefore the general solution is a linear combination of the three independent solutions we determined X05 Cixi C2X2 Csxs We remark that a given eigenvalue with multiplicity two may not yield two independent eigenvectors as was the case in the last example Then we must proceed differently to nd another independent solution such as the method given in Section 53 see Exercise 2c below Example 524 If the matrix for a threedimensional system X AX has one real eigenvalue A and two complex conjugate eigenvalues a i ib7 with associated eigenvectors V and W i iz7 respectively then the general solution is as is expected from Section 5327 xt elven 026 quot W cos bt 7 2 sin bt 036 quot W sin bt 2 cos bt EXERCISES 1 Find the eigenvalues and eigenvectors of the following matrices 2 3 2 3 4 1 0 A 0 6 2 397 B 2 0 2 397 C 0 1 0 0 0 71 4 2 3 2 Find the general solution of the following threedimensional systems 3 l 3 a x 75 73 73 x Hint A 4 is one eigenvalue 6 6 4 702 0 02 b X 02 704 0 x Hint A 71 is one eigenvalue 0 04 702 55 ThreeDimensional Systems 207 2 1 72 c x 71 0 0 xi Hint see Section 5313 0 2 72 1 0 1 d x 0 1 0 x 1 71 1 3 Find the general solution of the system 1 pr 7 y7 y I P97 2 722 Where p is a constant 4 Consider the system 0 1 2 x 1 0 2 xi 71 72 73 a Show that the eigenvalues are A 71 71 71 10 Find an eigenvector V1 associated With A 71 and obtain a solution to the systemi 0 V Show that a second independent solution has the form V2 tV1e t and nd V21 Show that a third independent solution has the form V3 W t2vlequot and nd V31 39D V Find the general solution and then solve the initial value problem x AX7 x0 010 6 Nonlinear Systems Nonlinear dynamics is common in nature Unlike linear systems where we can always nd explicit formulas for the solution nonlinear systems can seldom be solved For some nonlinear systems we even have to give up on existence and uniqueness So nonlinear dynamics is much more complicated than linear dynamics and therefore we rely heavily on qualitative methods to determine their dynamical behavior As for linear systems equilibrium solutions and their stability play a fundamental role in the analysis 61 Nonlinear Models 611 Phase Plane Phenomena A twodimensional nonlinear autonomous system has the general form 1 NW 61 y 9I7y7 62 where f and g are given functions of z and y that are assumed to have contin uous rst partial derivatives in some open region in the plane This regularity assumption on the rst partial derivatives guarantees that the initial value prob lem associated with 616 2 will have a unique solution through any point in the region Nonlinear systems arise naturally in mechanics circuit theory 210 6 Nonlinear Systems compartmental analysis reaction kinetics mathematical biology economics and other areas In fact in applications most systems are nonlineari Example 61 We have repeatedly noted that a secondorder equation can be reformulated as a rstorder systemi As a reminder consider Newton s second law of motion for a particle of mass m moving in one dimension m1 Fzz 7 where F is a force depending upon the position and the velocity Introducing the velocity y 1 as another state variable we obtain the equivalent rstorder system 1 7 y 7 mFI7yl Consequently we can study mechanical systems in an zy phase space rather than the traditional position7time space In this chapter we are less reliant on vector notation than for linear systems where vectors and matrices provide a natural language We review some general terminology of Chapter 5 A solution I zt y yt to 6162 can be represented graphically in two different ways see gure 51 in Chapter 5 We can plot I vs t and y vs t to obtain the time series plots showing how the states I and y vary with time t Or we can plot the parametric equations 1 1t y yt in the my phase plane A solution curve in the my plane is called an orbit On a solution curve in the phase plane time is a parameter and it may be shifted at will397 that is if z 1t y yt is a solution then I zt 7 c y yt 7 0 represents the same solution and same orbit for any constant cl This follows because the system is autonomous The initial value problem lVP consists of the solving the system 6il76i2 subject to the initial conditions 1050 a007 9050 yo Geometrically this means nding the orbit that goes through the point 10 yo at time to If the functions f and g are continuous and have continuous rst partial derivatives on R2 then the IVP has a unique solution Therefore two different solution curves cannot cross in the phase planet We always assume conditions that guarantee existence and uniqueness As is true for their linear counterparts there is an important geometric interpretation for nonlinear systems in terms of vector elds For a solution 61 Nonlinear Models 211 z m y ya we have mmm mmylttgtgt7gltzlttgt7ylttgtgtgtl Therefore at each point 17y in the plane the functions f and 9 de ne a vector fz ygzy that is the tangent vector to the orbit which passes through that point Thus the system 6162 generates a vector eld A different way to think about it is this The totality of all orbits form the flow of the vector eld lntuitively we think of the ow as uid particle paths with the vector eld representing the velocity of the particles at the various points A plot of several representative or key orbits in the zy plane is called a phase diagram of the systeml It is important that f and y do not depend explicitly upon timer Otherwise the vector eld would not be stationary and would change giving a different vector eld at each instant of time This would spoil a simple geometric approach to nonlinear systems Nonautonomous systems are much harder to deal with than autonomous onesl Among the most important solutions to 60762 are the constant so lutions or equilibrium solutions These are solutions zt ze yt ye where 18 and ye are constant Thus equilibrium solutions are found as solutions of the algebraic simultaneous system of equations ay 07 9179 0 The time series plots of an equilibrium solution are just constant solutions hor izontal lines in time In the phase plane an equilibrium solution is represented by a single point 187 We often refer to these as equilibrium pointsi Non linear systems may have several equilibrium points If an equilibrium point in the phase plane has the property that there is a small neighborhood about the point where there are no other equilibria then we say the equilibrium point is isolated Example 62 If a particle of mass m 1 moves on an zaxis under the in uence of a force Fz 312 7 17 then the equations of motion in the phase plane take the form I 1 97 y 312717 where the position I and the velocity y are functions of time t Here we can obtain the actual equation for the orbits in the zy phase plane in terms of z and y Dividing the two equations1 gives dydt 7 dy 7 31271 dzdt 7 dz y 1 Along an orbit x Mt y yt we also have y as a function of x or y Then the Chain rule dictates dgli fig 7 212 6 Nonlinear Systems Separating variables and integrating yields ydy3127ldx 1 if 13 7erE7 63 where we have chosen the letter E to denote the arbitrary constant of integra tion as we soon observe E stands for total energy This equation represents a family of orbits in the phase plane giving a relationship between position and velocity By dividing the equations as we did time dependence is lost on these orbitsl Equation 63 has an important physical meaning that is worth review ing The term y2 represents the kinetic energy onehalf the mass times the velocitysquared Secondly we recall that the potential energy Vz associated with a conservative force Fz is Vz 7 fFzdz or Fz 7dVdzi In the present case Vz 713 z where we have taken V 0 at z 0 The orbits 63 can be written 1 y2 713 1 7 E which states that the kinetic energy plus the potential energy is constant Therefore the orbits 63 represent constant energy curves The total en ergy E can be written in terms of the initial position and velocity as E yQ 0 7103 For each value of E we can plot the locus of points de ned by equation 63 To carry this out practically we may solve for y and y IS7IE7 y7 137IE Then we can plot the curves using a calculator or computer algebra systeml For write values of I that make the expression under the radical negative the curve will not be de ned Figure 61 shows several orbitsi Let us discuss their features There are two points z y 0 and z 7 y 07 where 1 y 0 These are two equilibrium solutions where the velocity is zero and the force is zero so the particle cannot be in motion The equilibrium solution I 7 E y 0 has the structure of a center and for initial values close to this equilibrium the system will oscillatei The other equilibrium 1 y 0 has the structure of an unstable saddle pointi Because 1 y for y gt 0 we have 1 gt 0 and the orbits are directed to the right in the upper halfplaner For y lt 0 we have 1 lt 0 and the orbits are directed to the left in the lower halfplaner For large initial energies the system does not oscillate but rather goes to z 00 y 00 that is the mass moves farther and farther to the right with faster speed 61 Nonlinear Models 213 Figure 61 Plots of the constant energy curves y2 7 13 z E in the zyphase planer These curves represent the orbits of the system and show how position and velocity relatel Time dependence is lost in this representation of the orbits Because 1 y the orbits are moving to the right I is increasing in the upper halfplane y gt 0 and to the left I is decreasing in the lower halfplane y lt 0 Example 63 Consider the simple nonlinear system If ya 64gt y ix 65gt 3 l l Clearly the origin 1 0 y 0 is the only equilibrium solution In this case we can divide the two equations and separate variables to get 2 2 ii IIl yy 3 integrating with respect to t gives y2ydt7 zzdtc where C is an arbitrary constant Changing variables in each integral y yt in the left integral and z zt in the right we obtain y2dy7 zdzc 214 6 Nonlinear Systems ys 712Ci Rearranging y lt07 12gt Consequently we have obtained the orbits for system 646 5 in terms of z and y These are easily plotted eg on a calculator for different values of C and they are shown in Figure 62 This technique illustrates a general 71 08 708 704 02 Figure 62 Phase diagram for z y27 y i zi Because 1 gt 0 all the orbits are moving to the right as time increases method for nding the equation of the orbits for simple equations in terms of the state variables alone divide the differential equations and integrate as far as possible With this technique however we lose information about how the states depend on time or how time varies along the orbits To nd solution curves in terms of time t we can write 64 as 1 y2 C 7 13237 which is a single differential equation for z We can separate variables but the result is not very satisfying because we get a complicated integrali This shows that time series solutions are not easily obtained for nonlinear problemsi Usually the qualitative behavior shown in the phase diagram is all we want If we do need time series plots we can obtain them using a numerical method which we discuss later 61 Nonlinear Models 215 We point out an important feature of the phase diagram shown in gure 62 The origin does not have the typical type of structure encountered in Chapter 5 for linear systems There we were able to completely characterize all equilibrium solutions as saddles spirals centers or nodes The origin for the nonlinear system 646 5 is not one of those Therefore nonlinear systems can have an unusual orbital structure near equilibria Why are the equilibrium solutions so important First much of the ac tion in the phase plane takes place near the equilibrium points so analysis of the ow near those points is insightful Second physical systems often seek out and migrate toward equilibria397 so equilibrium states can represent persistent states We think of z and y as representing two competing animal populationsl If a system is in an equilibrium state the two populations coexistl Those pop ulations will remain in the equilibrium states unless the system is perturbed This means that some event eg a bonanza or catastrophe would either add or subtract individuals from the populations without changing the under lying processes that govern the population dynamics If the in icted population changes are small the populations would be bumped to new values near the equilibriuml This brings up the stability issue Do the populations return to the coexistent state or do they change to another state If the populations return to the equilibrium then it is a persistent state and asymptotically stable If the populations move further away from the equilibrium then it is not persis tent and unstable If the populations remain close to the equilibrium but do not actually approach it then the equilibrium is neutrally stablel For each model it is important to discover the locally stable equilibrium states or persis tent states in order to understand the dynamics of the modell ln Example 62 the saddle point is unstable and the center is neutrally stable gure 61 and in Example 63 the equilibrium is unstable gure 62 For an unstable equi librium orbits that begin near the equilibrium do not remain nearl Examples of different types of stability are discussed in the sequel The emphasis in the preceding paragraph is on the word local That is what happens if small changes occur near an equilibrium not large changes Of course we really want to know what happens if an equilibrium is disturbed by all possible changes including an arbitrarily large change Often the adjec tives local and global are appended to stability statements to indicate what types of perturbations small or arbitrary are under investigation However we cannot usually solve a nonlinear system and so we cannot get an explicit resolution of global behaviorl Therefore we are content with analyzing local stability properties and not global stability properties As it turns out local stability can be determined because we can approximate the nonlinear system by a tractable linear system near equilibria Section 63 216 6 Nonlinear Systems EXERCISES 1 Consider the uncoupled nonlinear system I 12 y 7y a Find a relation between I and y that describes the orbits Are all the orbits contained in this relation for different values of the arbitrary constant b Sketch the vector eld at several points near the origin c Draw the phase diagram ls the equilibrium stable or unstable d Find the solutions I 1t y yt and plot typical time series Pick a single time series plot and draw the corresponding orbit in the phase plane 2 Consider the system 1 7 y 21 a Are there any equilibrium solutions b Find a relationship between I and y that must hold on any orbit and plot several orbits in the phase planer c From the orbits sketch the vector eld d Do any orbits touch the zaxis 3 Consider the nonlinear system 1 12 y2 7 4 y y 7 21 F 01 1533 a Find the two equilibria and plot them in the phase planer b On the plot in part a sketch the set of points where the vector eld is vertical up or down and the set of points where the vector eld is horizontal left or right Do parts a and b of the previous problems for the nonlinear system 1 y1 y y12i i A nonlinear model of the form 1 y 7 z 512 y y my has been proposed to describe cell differentiation Find all equilibrium so lutions Find all equilibria for the system 1 sin y y 21 i Consider the nonlinear system 1 y y 7z7y3i Show that the function Vzy 12 y2 decreases along any orbit iiei Vzt lt 0 and state why this proves that every orbit approaches the origin as t 7 00 61 Nonlinear Models 217 8 Consider the nonlinear system 1 12 7 y2 y z 7 y a Find and plot the equilibrium points in the phase planer Are they isolated b Show that on an orbit z y 1 Cey where C is some constant and plot several of these curves c Sketch the vector eld d Describe the fate of the orbit that begins at i 0 at t 0 as t 7 00 and as t7gt 700 e Draw a phase plane diagram being sure to indicate the directions of the orbits 612 The Lotka7V01terI39a Model Nonlinear equations play an central role in modeling population dynamics in ecologyi We formulate and study a model involving predator7prey dynamicsi Let I zt be the prey population and y yt be the predator population We can think of rabbits and foxes food sh and sharks or any consumer resource interaction including herbivores and plants If there is no predator we assume the prey dynamics is I T1 or exponential growth where 7 is the per capita growth rate In the absence of prey we assume that the predator dies via y 7 7mg where m is the per capita mortality rate When there are interactions we must include terms that decrease the prey population and increase the predator population To determine the form of the predation term we assume that the rate of predation or the the number of prey consumed per unit of time per predator is proportional to the number of prey That is the rate of predation is azi Thus if there are y predators then the rate that prey is decreased is azyi Note that the interaction term is proportional to my the product of the number of predators and the number of prey For example if there were 20 prey and 10 predators there would be 200 possible interactions Only a fraction of them a are assumed to result in a kill The parameter a depends upon the fraction of encounters and the success of the encounters The prey consumed cause a rate of increase in predators of Suzy where 5 is the conversion ef ciency of the predator population one prey consumed does not mean one predator borni Therefore we obtain the simplest model of predator7 prey interaction called the Lotka7Volterra model I TI 7 azy y 7mybzy7 218 6 Nonlinear Systems where b Ear To analyze the Lotka7Volterra model we factor the right sides of the equa tions to obtain 1 17quot 7 ay y y7m 121 66 Now it is simple to locate the equilibrial Setting the right sides equal to zero gives two solutions I 0 y 0 and z mb y Tai Thus in the phase plane the points 00 and mbJa represent equilibrial The origin rep resents extinction of both species and the nonzero equilibrium represents a coexistent state To determine properties of the orbits we usually plot curves in the my plane where the vector eld is vertical where z 0 and curves where the vector eld is horizontal y 0 These curves are called the null clinesl They are not usually orbits but rather the curves where the orbits cross vertically or horizontally The znullclines for 66 where z 07 are I 0 and y Tal Thus the orbits cross these two lines vertically The y nullclines where y 07 are y 0 and z mbl The orbits cross these lines horizontallyi Notice that the equilibrium solutions are the intersections of the z and y nullclinesi The nullclines partition the plane into regions where z and y have various signs and therefore we get a picture of the direction of the ow pattern See gure 63 Next along each nullcline we can nd the direction of y ynullcline x lt 0 x lt 0 ylt 0 lt ygt 0 ra I xnullcllne x gt 0 gt x gt 0 y lt 0 y gt 0 07 0 quot11 X Figure 63 Nullclines dashed and vector eld in regions between nullclinesi The I and y axes are nullclines as well as orbitsl the vector eld For example on the ray to the right of the equilibrium we have I gt mb y Tal We know the vector eld is vertical so we need only check the sign of yl We have y y7m bx Maxim 121 gt 0 so the vector eld points upwardi Similarly we can determine the directions along the 61 Nonlinear Models 219 other three raysi These are shown in the accompanying gure 63 Note that y 0 and z 0 both nullclines are also orbitsi For example when I 0 we have M imy or yt Ce m 7 when there are no prey the foxes die out Similarly when y 0 we have zt Ce so the rabbits increase in number Finally we can determine the direction of the vector eld in the regions between the nullclines either by selecting an arbitrary point in that region and calculating z and y or by just noting the sign of z and y in that region from information obtained from the systemi For example in the quadrant above and to the right of the nonzero equilibrium it is easy to see that z lt 0 and y gt 0397 so the vector eld points upward and to the left We can complete this task for each region and obtain the directions shown in gure 63 Having the direction of the vector eld along the nullclines and in the regions bounded by the nullclines tells us the directions of the solution curves or orbitsi Near 00 the orbits appear to veer away and the equilibrium has a saddle point structure The equilibrium 07 0 is unstable It appears that orbits circle around the nonzero equilibrium in a counterclockwise fashioni But at this time it is not clear if they form closed paths or spirals so more work is needed We attempt to obtain the equation of the orbits by dividing the two equa tions in 66 We get y 7 dy 7 ymbz gigi 17 704y Rearranging and integrating gives dyl z mdzo y z Carrying out the integration gives 39rlnyiaybzimlnzC7 which is the algebraic equation for the orbits It is obscure what these curves are because it is not possible to solve for either of the variables So cleverness is required If we exponentiate we get yreiay eCebxrim Now consider the y nullcline where z is xed at a value mb7 and x a positive C value ie x an orbit The right side of the last equation is a positive num ber A and so yr Aeayi If we plot both sides of this equation do thisliplot a power function and a growing exponential we observe that there can be at most two intersections397 therefore this equation can have at most two solutions for y Hence along the vertical line I mb there can be at most two cross ings 7 this means an orbit cannot spiral into or out from the equilibrium point because that would mean many values of y would be possible We conclude 220 6 Nonlinear Systems x xrxy yeryoxy Figure 64 Closed counterclockwise periodic orbits of the Lotka7Volterra predatoriprey model 1 z 7 my y 7y my The zaxis is an orbit leaving the origin and the y axis is an orbit entering the origin that the equilibrium is a center with closed periodic orbits encircling it A phase diagram is shown in gure 64 Time series plots of the prey and preda tor populations are shown in gure 65 When the prey population is high the predators have a high food source and their numbers start to increase thereby eventually driving down the prey populationl Then the prey population gets low ultimately reducing the number of predators because of lack of food Then the process repeats giving cycles The Lotka7Volterra model developed by Al Lotka and Vi Volterra in 1925 is the simplest model in ecology showing how populations can cycle and it was one of the rst strategic models to explain qualitative observations in natural systems Note that the nonzero equilibrium is neutrally stablel A small perturbation from equilibrium puts the populations on a periodic orbit that stays near the equilibriuml But the system does not return to that equilibriuml So the nonzero equilibrium is stable but not asymptotically stable The other equilibrium the origin which corresponds to extinction of both species is an unstable saddle point with the two coordinate axes as separatricesl

### BOOM! Enjoy Your Free Notes!

We've added these Notes to your profile, click here to view them now.

### You're already Subscribed!

Looks like you've already subscribed to StudySoup, you won't need to purchase another subscription to get this material. To access this material simply click 'View Full Document'

## Why people love StudySoup

#### "There's no way I would have passed my Organic Chemistry class this semester without the notes and study guides I got from StudySoup."

#### "I bought an awesome study guide, which helped me get an A in my Math 34B class this quarter!"

#### "Knowing I can count on the Elite Notetaker in my class allows me to focus on what the professor is saying instead of just scribbling notes the whole time and falling behind."

#### "It's a great way for students to improve their educational experience and it seemed like a product that everybody wants, so all the people participating are winning."

### Refund Policy

#### STUDYSOUP CANCELLATION POLICY

All subscriptions to StudySoup are paid in full at the time of subscribing. To change your credit card information or to cancel your subscription, go to "Edit Settings". All credit card information will be available there. If you should decide to cancel your subscription, it will continue to be valid until the next payment period, as all payments for the current period were made in advance. For special circumstances, please email support@studysoup.com

#### STUDYSOUP REFUND POLICY

StudySoup has more than 1 million course-specific study resources to help students study smarter. If you’re having trouble finding what you’re looking for, our customer support team can help you find what you need! Feel free to contact them here: support@studysoup.com

Recurring Subscriptions: If you have canceled your recurring subscription on the day of renewal and have not downloaded any documents, you may request a refund by submitting an email to support@studysoup.com

Satisfaction Guarantee: If you’re not satisfied with your subscription, you can contact us for further help. Contact must be made within 3 business days of your subscription purchase and your refund request will be subject for review.

Please Note: Refunds can never be provided more than 30 days after the initial purchase date regardless of your activity on the site.