Heath - Scientific Computing (523150), страница 74
Текст из файла (страница 74)
Successively refine a uniform mesh until a givenerror tolerance is met. Estimate the error ateach stage by comparing the values obtainedfor consecutive mesh sizes. What kind of datastructure is needed for reusing previously computed function values?(b) Write an adaptive quadrature routine usingthe composite Simpson rule. Successively refine only those subintervals that have not yetmet an error tolerance. What kind of datastructure is needed for keeping track of whichsubintervals have converged?After debugging, test your routines using someof the integrals in the previous problems andcompare the results with those previously obtained. How does the efficiency of your adaptive routine compare with that of your nonadaptive routine?8.14 Select an automatic adaptive quadratureroutine and try to devise an integrand functionfor which it gives an answer that is completelywrong.
(Hint: This problem may require atleast one round of trial and error.) Can youdevise a smooth function for which the adaptive routine is seriously in error?8.15 (a) Solve the integral equationZ01(s2 + t2 )1/2 u(t) dt =(s2 + 1)3/2 − s33on the interval [0, 1] by discretizing the integral using the composite Simpson quadraturerule with n equally spaced points tj , and alsousing the same n points for the si . Solve theresulting linear system Ax = y using a libraryroutine for Gaussian elimination with partialpivoting. Experiment with various values for nin the range from 3 to 15, comparing your results with the known unique solution, u(t) = t.Which value of n gives the best results? Canyou explain why?(a) The unit square, i.e., 0 ≤ x ≤ 1, 0 ≤ y ≤ 1.(b) For each value of n in part a, compute thecondition number of the matrix A. How doesit behave as a function of n?(b) The quarter of the unit disc lying in thefirst quadrant, i.e., x2 + y 2 ≤ 1, x ≥ 0, y ≥ 0.(c) Repeat part a, this time solving the linearsystem using the singular value decomposition,over each of the following regions:274CHAPTER 8.
NUMERICAL INTEGRATION AND DIFFERENTIATIONbut omit any “small” singular values. Try various thresholds for truncating the singular values, and again compare your results with theknown true solution.(d ) Repeat part a, this time using the methodof regularization. Experiment with variousvalues for the regularization parameter µ todetermine which value yields the best resultsfor a given value of n. For each value of µ, plota point on a two-dimensional graph whose axesare the norm of the solution and the norm ofthe residual.
What is the shape of the curvetraced out as µ varies? Does this shape suggestan optimal value for µ?(e) Repeat part a, this time using an optimization routine to minimize ky − Axk22 subject tothe constraint that the components of the solution must be nonnegative. Again, compareyour results with the known true solution.(f ) Repeat part e, this time imposing the additional constraint that the solution be monotonically increasing, i.e., x1 ≥ 0 and xi −xi−1 ≥ 0,i = 2, .
. . , n. How much difference does thismake in approximating the true solution?8.16 In this exercise we will experimentwith numerical differentiation using data fromComputer Problem 3.1:ty0.01.01.02.72.05.83.06.64.07.55.09.9For each of the following methods for estimating the derivative, compute the derivative ofthe original data and also experiment with randomly perturbing the y values to determinethe sensitivity of the resulting derivative estimates. For each method, comment on boththe reasonableness of the derivative estimatesand their sensitivity to perturbations.
Notethat the data are monotonically increasing, soone might expect the derivative always to bepositive.(a) For n = 0, 1, . . . , 5, fit a polynomial of degree n by least squares to the data, then differentiate the resulting polynomial and evaluatethe derivative at each of the given t values.(b) Interpolate the data with a cubic spline,differentiate the resulting piecewise cubic polynomial, and evaluate the derivative at each ofthe given t values (some spline routines provide the derivative automatically, but it canbe done manually if necessary).(c) Repeat part b, this time using a smoothingspline routine. Experiment with various levelsof smoothing, using whatever mechanism forcontrolling the degree of smoothing that theroutine provides.(d ) Interpolate the data with a monotonic Hermite cubic, differentiate the resulting piecewisecubic polynomial, and evaluate the derivativeat each of the given t values.Chapter 9Initial Value Problems for OrdinaryDifferential Equations9.1Ordinary Differential EquationsWe now turn to the study of differential equations, that is, equations involving derivativesof the unknown solution.
We have previously considered only algebraic equations, for whichthe unknown solution is a discrete vector in a finite-dimensional space. For a differentialequation, on the other hand, the unknown solution is a continuous function in an infinitedimensional space. Our approach to solving differential equations numerically will be basedon finite-dimensional approximations, a process called discretization. We will replace differential equations with algebraic equations whose solutions approximate those of the givendifferential equations.First, we establish some notation and definitions.
A system of ordinary differentialequations (ODEs) has the general formy 0 (t) = f (t, y(t)),where t is a real variable, y: R → Rn is a vector-valued function of t, f : Rn+1 → Rn , andy 0 (t) = dy(t)/dt denotes the derivative with respect to t, i.e., y 0 (t) dy1 (t)/dt dy2 (t)/dt ...=.1 y20 (t) ...0yn (t)dyn (t)/dtThus, we have a system of coupled differential equations in which we are given the functionf and we wish to determine the unknown function y. An important special case, which wewill often consider for simplicity, is n = 1, i.e., a single scalar ODE.2752769.1.1CHAPTER 9. INITIAL VALUE PROBLEMS FOR ODESInitial Value ProblemsAn ordinary differential equation y 0 = f (t, y) by itself does not determine a unique solutionfunction because the equation merely specifies the slopes of the solution components y 0 (t)at each point but not the actual solution value y(t) at any point.
Thus, in general, there isan infinite family of functions that satisfy the differential equation, provided f is sufficientlysmooth.To single out a particular solution, we must specify the value, usually denoted by y0 , ofthe solution function at some point, usually denoted by t0 . Thus, part of the given problemdata is the requirement thaty(t0 ) = y0 .This additional requirement determines a unique solution to the ODE, provided that f iscontinuously differentiable. Because the independent variable t usually represents time, wethink of t0 as the initial time and y0 as the initial value. Hence, this is termed an initialvalue problem. The ODE governs the dynamic evolution of the system in time from itsinitial state y0 at time t0 onward, and we seek a function y(t) that describes the state ofthe system as a function of time.Example 9.1 Initial Value Problem.
Consider the scalar ordinary differential equationy 0 = y. This is an ODE of the form y 0 = f (t, y), where in this example f (t, y) = y. Thefamily of solutions for this equation is given by y(t) = cet , where c is any real constant.If we impose an initial condition, such as requiring that y(t0 ) = y0 , then this will singleout the unique particular solution that satisfies the initial condition. For this example, ift0 = 0, then we get c = y0 , which means that the solution is y(t) = y0 et . Some membersof the family of solutions for this equation are sketched in Fig. 9.1, including the particularsolution that satisfies the given initial condition.y.................................................................0............................................................................................................................................................................................................................................................................................................................................
........................................................................................................................................................................................................................................................................................................................0 ..........................................................................................................................................................................................................................................................................................................................................................................y =yy •t0tFigure 9.1: The family of solution curves for the ODE y 0 = y.9.1.2Higher-Order ODEsIf the first derivative is the highest-order derivative of the solution function appearing inthe equation, an ODE is said to be of first order.
Equations with higher-order derivatives9.1. ORDINARY DIFFERENTIAL EQUATIONS277occur frequently in practice but can be transformed into an equivalent first-order system asfollows. For example, given an nth order scalar equationu(n) = f (t, u, u0 , . . . , u(n−1) ),define the n new unknowns y1 (t) = u, y2 (t) = u0 , . .
. , yn (t) = u(n−1) , so that the originalequation becomes the first-order system of n equations 0 y1y20 y y3 2 .. ...= . . 0 yn−1 ynf (t, y1 , y2 , . . . , yn )yn0Thus, in general, a scalar ODE of order n is equivalent to a system of n first-order ODEs.If a system of ODEs contains equations having higher-order derivatives, then each suchcomponent equation can be transformed into an equivalent first-order system in this samemanner. For example, a system of two second-order equations would yield an equivalentsystem of four first-order equations.
For this reason, most ODE software is designed to solveonly first-order equations, and we will also restrict our attention to first-order equations indiscussing numerical solution methods.Example 9.2 Newton’s Second Law. To illustrate the transformation of a higherorder ODE into an equivalent system of first-order ODEs, consider Newton’s Second Lawof Motion, F = ma, in one dimension.
This is a second-order ODE, since the acceleration ais the second derivative of the position coordinate, which we denote by x. Thus, the ODEhas the formx00 = F/m,where F and m represent force and mass, respectively. To transform this second-order,scalar ODE into a system of two first-order ODEs, we define two new functions y1 = x andy2 = x0 . This step gives us the system of two first-order equations 0 y1y2=.y20F/mWe can now use a method for first-order equations to solve this system. When we do so,the first component of the solution y1 will be the same as the solution x to the originalsecond-order equation.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
