Heath - Scientific Computing (523150), страница 82
Текст из файла (страница 82)
We will solve this ODEnumerically using a stepsize of h = 0.5.(a) Is this ODE stable?(b) Is Euler’s method stable for this ODE using this stepsize?(c) Compute the numerical value for the approximate solution at t = 0.5 given by Euler’smethod.(d ) Is the backward Euler method stable forthis ODE using this stepsize?(e) Compute the numerical value for the approximate solution at t = 0.5 given by thebackward Euler method.9.5 With an initial value of y0 = 1 at t0 = 0and a time step of h = 1, compute the approximate solution value y1 at time t1 = 1 for theODE y 0 = −y using each of the following twonumerical methods. (Your answers should benumbers, not formulas.)(a) Euler’s method(b) Backward Euler method9.6 For the ODE, initial value, and stepsizegiven in Example 9.10, prove that fixed-pointiteration for solving the implicit equation fory1 is in fact convergent.
What is the convergence rate?9.7 Consider the initial value problemy 00 = yfor t ≥ 0, with initial values y(0) = 1 andy 0 (0) = 2.(a) Express this second-order ODE as anequivalent system of two first-order ODEs.(b) What are the corresponding initial conditions for the system of ODEs in part a?(c) Is this a stable system of ODEs?(d ) Perform one step of Euler’s method for thisODE system using a stepsize of h = 0.5.(e) Is Euler’s method stable for this problemusing this stepsize?(f ) Is the backward Euler method stable forthis problem using this stepsize?9.8 Consider the initial value problem forthe ODE y 0 = −y 2 with the initial conditiony(0) = 1. We will use the backward Eulermethod to compute the approximate value ofthe solution y1 at time t1 = 0.1 (i.e., takeone step using the backward Euler methodwith stepsize h = 0.1 starting from y0 = 1at t0 = 0).
Since the backward Euler methodis implicit, and the ODE is nonlinear, we willneed to solve a nonlinear algebraic equation fory1 .304CHAPTER 9. INITIAL VALUE PROBLEMS FOR ODES(a) Write out that nonlinear algebraic equation for y1 .(b) Write out the Newton iteration for solvingthe nonlinear algebraic equation.(c) Obtain a starting guess for the Newton iteration by using one step of Euler’s method forthe ODE.(d ) Finally, compute an approximate value forthe solution y1 by using one iteration of Newton’s method for the nonlinear algebraic equation.(c) Self-starting9.9 For solving an ODE y 0 = f (t, y) numerically, each of the following methods,(1)9.11 The centered difference approximationyk+11= yk + [f (tk , yk )2+f (tk+1 , yk + f (tk , yk )h)]h,(2)yk+11= yk + [3f (tk , yk )2−f (tk−1 , yk−1 )]h,(3)yk+11= yk + [f (tk , yk )2+f (tk+1 , yk+1 )]his of second order, but the methods alsohave important differences.
For each propertylisted, state which of the three methods has orhave the given property.(a) Single-step method(b) Implicit method(d ) Unconditionally stable(e) Runge-Kutta type method(f ) Good for solving a stiff ODE9.10 Use the linear ODE y 0 = λy to analyzethe accuracy and stability of Heun’s method(see Section 9.6.2). In particular, verify thatthis method is second-order accurate, and describe or plot its stability region in the complexplane.y0 ≈yk+1 − yk−12hleads to the two-step leapfrog methodyk+1 = yk−1 + f (tk , yk )2hfor solving the ODE y 0 = f (t, y).
Determinethe order of accuracy and the stability regionof this method.9.12 Let A be an n × n matrix. Compare andcontrast the behavior of the linear differenceequationxk+1 = Axkwith that of the linear differential equationx0 = Ax.What is the general solution in each case?In each case, what property of the matrixA would imply that the solution remainsbounded for any starting vector x0 ? You mayassume that the matrix A is diagonalizable.Computer Problems9.1 The populations of two species, a preydenoted by y1 and predator denoted by y2 , canbe modeled by a system of ODEsy10y20= by1 − cy1 y2 ,= −dy2 + cy1 y2due to Lotka and Volterra.
The parameters band d govern the birth rate of prey and deathrate of predators, respectively, and the param-eter c governs the interaction of the two populations. With the parameter values b = 1,d = 10, and c = 1, and initial conditionsy1 (0) = 0.5 and y2 (0) = 1 (the populationsare normalized, and we treat them as continuous variables), use a library routine to solvethis system numerically, integrating to t = 10.Plot each of the two populations as a functionof time, and on a separate graph plot the tra-COMPUTER PROBLEMSjectory of the point (y1 (t), y2 (t)) in the planeas a function of time.
The latter is sometimescalled a “phase portrait.” Give a physical interpretation of the behavior you observe. Canyou find nonzero initial populations such thateither of the populations eventually becomesextinct?9.2 The Kermack-McKendrick model for thecourse of an epidemic in a population is givenby the system of ODEsy10y20y30= −cy1 y2 ,= cy1 y2 − dy2 ,= dy2 ,where y1 represents susceptibles, y2 representsinfectives in circulation, and y3 represents infectives removed by isolation, death, or recovery and immunity.
The parameters c and drepresent the infection rate and removal rate,respectively. Use a library routine to solve thissystem numerically, with the parameter valuesc = 1 and d = 5, and initial values y1 (0) = 95,y2 (0) = 5, y3 (0) = 0. Integrate from t = 0 tot = 1. Plot each solution component on thesame graph as a function of t. As expectedwith an epidemic, you should see the number of infectives grow at first, then diminish tozero. Experiment with other values for the parameters and initial conditions.
Can you findvalues for which the epidemic does not grow, orfor which the entire population is wiped out?9.3 Suppose that we have three chemicalspecies whose concentrations are denoted byy1 , y2 , and y3 . If the rate of the reactiony1 → y2 is proportional to y1 , and the rateof the reaction y2 → y3 is proportional to y2 ,then the concentrations are governed by thesystem of ODEsy10y20y30= −k1 y1 ,= k1 y1 − k2 y2 ,= k2 y2 ,where k1 and k2 are the rate constants for thetwo reactions.(a) What is the Jacobian matrix for this ODEsystem, and what are its eigenvalues? If therate constants are positive, is this system stable? Under what conditions will the system bestiff?305(b) Solve the ODE system numerically, assuming initial concentrations y1 (0) = y2 (0) =y3 (0) = 1.
Take k1 = 1 and experiment withvalues of k2 of varying magnitude, specifically,k2 = 10, 100, and 1000. For each value of k2 ,solve the system using a Runge-Kutta method,an Adams method, and a method designed forstiff systems, such as a backward differentiation formula. You may use library routinesfor this purpose, or you may wish to developyour own routines, perhaps using the classicalfourth-order Runge-Kutta method, the fourthorder Adams-Bashforth predictor and AdamsMoulton corrector, and the BDF formula givenin Section 9.6.
If you develop your own codes,a fixed stepsize will suffice for this exercise. Ifyou use library routines, compare the different methods with respect to their efficiency, asmeasured by function evaluations or executiontime, for a given accuracy. If you develop youown codes, compare the different methods withrespect to accuracy and stability for a givenstepsize. In each instance, integrate the ODEsystem from t = 0 until the solution is approximately in steady state, or until the method isclearly unstable or grossly inefficient.9.4 Experiment with several different libraryroutines having automatic stepsize selection tosolve the ODEy 0 = −200ty 2numerically.
Consider two different initial conditions, y(0) = 1 and y(−3) = 1/901, and ineach case compute the solution until t = 1.Monitor the stepsize used by the routines anddiscuss how and why it changes as the solution progresses. Explain the difference in behavior for the two different initial conditions.Compare the different routines with respect toefficiency for a given accuracy requirement.Rb9.5 A definite integral a f (t) dt can be evaluated by solving the equivalent ODE y 0 (t) =f (t), a ≤ t ≤ b, with initial condition y(a) = 0.The value of the integral is then simply y(b).Use a library ODE solver to evaluate eachdefinite integral in the first several ComputerProblems for Chapter 8, and compare its efficiency with that of an adaptive quadratureroutine for the same accuracy.306CHAPTER 9.
INITIAL VALUE PROBLEMS FOR ODES9.6 Homotopy methods for solving systemsof nonlinear algebraic equations parameterizethe solution space x(t) and then follow a trajectory from an initial guess to the final solution. As one example of this approach, forsolving a system of nonlinear equations f (x) =0, where f : Rn → Rn , with initial guess x0 , thefollowing ODE initial value problem is a continuous analogue of Newton’s method:x0 = −Jf−1 (x)f (x),x(0) = x0 ,where Jf is the Jacobian matrix of f , and ofcourse the inverse need not be computed explicitly.
Use this method to solve the nonlinear system given in Computer Problem 5.13.Starting from the given initial guess, integratethe resulting system of ODEs from t = 0 untila steady state is reached. Compare the resulting solution with that obtained by a conventional nonlinear system solver. Plot the trajectory of the components of x(t) from t = 0to the final solution. You may also want to trythis technique on some of the other ComputerProblems from Chapter 5.9.7 An important problem in classical mechanics is the motion of two bodies under mutual gravitational attraction.
Suppose that abody of mass m is orbiting a second body ofmuch larger mass M , such as the earth orbitingthe sun. From Newton’s laws of motion andgravitation, the orbital trajectory (x(t), y(t))is described by the system of second-orderODEsx00y 00= −GM x/r3 ,= −GM y/r3 ,where G is the gravitational constant and r =(x2 +y 2 )1/2 is the distance of the orbiting bodyfrom the center of mass of the two bodies.For this exercise, we choose units such thatGM = 1.(a) Use a library routine to solve this systemof ODEs with the initial conditionsx(0) = 1 − e,x0 (0) = 0y(0) = 0,1/21+e,y 0 (0) =1−ewhere e is the eccentricity of the resulting elliptical orbit, which has period 2π.
Try thevalues e = 0 (which should give a circular orbit), e = 0.5, and e = 0.9. For each case, solvethe ODE for at least one period and obtainoutput at enough intermediate points to drawa smooth plot of the orbital trajectory. Makeseparate plots of x versus t, y versus t, andy versus x. Experiment with different errortolerances to see how they affect the cost ofthe integration and how close the orbit comesto being closed. If you trace the trajectorythrough several periods, does the orbit tend towander or remain steady?(b) Check your numerical solutions in part a tosee how well they conserve the following quantities, which should remain constant:Conservation of energy:1(x0 )2 + (y 0 )2−2rConservation of angular momentum:x y 0 − y x09.8 Consider a restricted form of the threebody problem in which a body of small massorbits two other bodies with much largermasses, such as an Apollo spacecraft orbitingthe earth-moon system. We will use a twodimensional coordinate system in the plane determined by the three bodies, with the originat the center of mass of the two larger bodies,and the coordinate system rotating so that thetwo larger bodies appear fixed.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
