Heath - Scientific Computing (523150), страница 75
Текст из файла (страница 75)
In addition, we will also get the second component y2 , which is thesame as the velocity x0 . In three dimensions, Newton’s Law would comprise three secondorder equations, one for each spatial coordinate, and would yield an equivalent system ofsix first-order equations.9.1.3Stable and Unstable ODEsRoughly speaking, if the members of the solution family for an ODE move away fromeach other with time, then the equation is said to be unstable; but if the members of278CHAPTER 9. INITIAL VALUE PROBLEMS FOR ODESthe solution family move closer to each other with time, then the equation is said to bestable.
If the solution curves are neither converging nor diverging (i.e., they remain nearbybut do not actually come together), then the equation is said to be neutrally stable. Thisdefinition of stability for ODEs is consistent with the general concept of stability discussedin Section 1.2.7 in that it reflects the sensitivity of a solution of the ODE to perturbations. Asmall perturbation to a solution of a stable equation will be damped out with time becausethe solution curves are converging, whereas for an unstable equation the perturbation willgrow with time because the solution curves are diverging.The stability of a cone provides a helpful geometric analogy. If a cone resting on itscircular base is slightly perturbed, it will return to its original position; the position isstable.
If a cone is balanced on its point, any slight perturbation will cause it to fall; theposition is unstable. If a cone is resting on its slanted side, then a slight perturbation willmove the cone to a new position nearby; the position is neutrally stable.Note that the concept of stability of an ODE depends on the entire family of solutions,not just on some particular solution. Moreover, both stable and unstable behavior canoccur in different portions of the domain of interest for the same equation.This qualitative concept of stability for an ODE y 0 = f (t, y) can be made more precisequantitatively by considering the Jacobian matrix Jf (t, y) with entries{Jf (t, y)}ij = ∂fi (t, y)/∂yj .If any of the eigenvalues of this matrix have positive real parts, then the equation is unstable.If all of the eigenvalues have negative real parts, then the equation is stable.
If one or moreeigenvalues have zero real parts, and all of the remainder have negative real parts, then theequation is neutrally stable. Since the entries of Jf are functions of t and y, its eigenvaluesmay vary with time, and hence the stability of the equation may vary from region to region.For scalar ODEs, which we will focus on for simplicity, stability of an ODE is determinedby the sign of its Jacobian, which is scalar valued in that case.Example 9.3 Unstable ODE. In Example 9.1, we considered the scalar ODE y 0 = y andsketched its family of solution curves y(t) = cet in Fig. 9.1. From the exponential growthof the solutions, we know that the solution curves for this equation move away from eachother as time increases, as we see in Fig. 9.1.
We can therefore conclude that the equationis unstable. More rigorously, we note that the Jacobian of f (i.e., ∂f /∂y) is positive (infact, it is the constant 1), so the equation is unstable.Example 9.4 Stable ODE. Let us now consider a different scalar equation, namely,y 0 = −y. The family of solutions for this equation is given by y(t) = ce−t , where c is anyreal constant. For this equation we see that the Jacobian of f is negative (∂f /∂y = −1),so the equation is stable. We also can see this from the exponential decay of the solutions,as shown in Fig. 9.2, in which some members of the solution family for this equation aredrawn.Example 9.5 Neutrally Stable ODE.
Consider the scalar ODE y 0 = a, for a given9.1. ORDINARY DIFFERENTIAL EQUATIONS279y.................................. ......................0............................................................................. .......................................................................................................................................................................................................................... .................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................y = −yt0tFigure 9.2: The family of solution curves for the ODE y 0 = −y.constant a.
The family of solutions is given by y(t) = at + c, where c is any real constant.Thus, the solution curves, as illustrated for a = 12 in Fig. 9.3, are parallel straight lines thatare neither converging nor diverging, and hence the equation is neutrally stable. Note that∂f /∂y = 0 for this equation, consistent with its neutral stability. Note also that the issuethat determines stability is not whether the solution curves are increasing or decreasing(either case can apply for this equation, depending on whether a is positive or negative)but rather the relationship of the solution curves to each other.y......................................................................10................................................2...................................................................................................................................................................................................................................................................................................................
.............................................................................................................................................................................................. ................................................................................................................................................................ ....................................................................................................................... .........................................................................
.............................................................................................................................................................................................................................y =t0tFigure 9.3: The family of solution curves for the ODE y 0 = 12 .Example 9.6 Linear System of ODEs. A linear, homogeneous system of ODEs withconstant coefficients has the formy 0 = Ay,where A is an n × n matrix.
Suppose we have the initial condition y(0) = y0 . Let theeigenvalues of A be denoted by λi , and the corresponding eigenvectors by ui , i = 1, . . . , n.For simplicity, assume that the eigenvectors are linearly independent, so that we can express280CHAPTER 9. INITIAL VALUE PROBLEMS FOR ODESy0 as a linear combinationy0 =nXαi ui .i=1Then it is easily confirmed thaty(t) =nXαi ui eλi ti=1is a solution to the ODE that satisfies the initial condition. We see that eigenvalues of Awith positive real parts yield exponentially growing solution components, eigenvalues withnegative real parts yield exponentially decaying solution components, and pure imaginaryeigenvalues with zero real parts yield oscillatory solution components.
These are consistent with our definitions of instability, stability, and neutral stability, respectively, as theJacobian J = A for this problem.9.2Numerical Solution of ODEsAn analytical solution of an ODE is a closed-form formula for computing the value of thesolution function at any point t. In contrast, a numerical solution of an ODE is a table ofapproximate values of the solution function at a discrete set of points. Such a numericalsolution is obtained by simulating the behavior of the system governed by the differentialequation. Approximate solution values are generated step by step in discrete incrementsmoving across the interval in which the solution is sought.
For this reason, numericalmethods for solving ODEs are sometimes called discrete variable methods.In stepping from one discrete point to the next, we will in general incur some error, whichmeans that our new approximate solution value will lie on a different member of the familyof solution curves for the ODE from the one on which we started. The stability or instabilityof the equation determines in part whether such errors are magnified or diminished withtime.9.2.1Euler’s MethodA numerical solution of an ODE is generated by simulating the behavior of the systemgoverned by the ODE.
Starting at t0 with the given initial value, we wish to track thetrajectory dictated by the ODE. Evaluating f (t0 , y0 ) tells us the slope of the trajectory atthat point. We use this information to predict the value y1 of the solution at some futuretime t1 = t0 + h for some suitably chosen increment h.The simplest example of this approach is Euler’s method . Consider the Taylor seriesy 00 (t) 2h + ···2y 00 (t) 2= y(t) + f (t, y(t))h +h + ···.2y(t + h) = y(t) + y 0 (t)h +9.2. NUMERICAL SOLUTION OF ODES281Euler’s method is derived by dropping terms of second and higher order to obtain theapproximate solution valueyk+1 = yk + f (tk , yk )hk ,which allows us to step from time tk to time tk+1 = tk + hk . Equivalently, if we replacethe derivative in the differential equation y 0 = f (t, y) with a finite difference quotient, weobtain an approximating algebraic equationyk+1 − yk= f (tk , yk ),hkwhich gives Euler’s method when solved for yk+1 .