Heath - Scientific Computing (523150), страница 84
Текст из файла (страница 84)
Using a stepsize of h = 0.5, we first step fromt0 = 0 to t1 = 0.5. The classical fourth-order Runge-Kutta method gives the approximatesolution value at t11y (1) = y (0) + (k1 + 2k2 + 2k3 + k4 )6 10.62500.50.500.68750.875=++2+2+=.0.750.75001.5001.75010.06Next we step from t1 = 0.5 to t2 = 1, obtaining 10.6252.02.00.8751.251.4375(2)y =++2+2+=,1.7501.5002.252.25003.04.06so we have hit the value y(1) = 2 instead of the desired value y(1) = 1. We try again, thistime with an initial slope of y 0 (0) = −1, obtaining 10−0.5−0.50−0.3125−1.25−0.375(1)y =++2+2+=−10.00.750.75001.50−0.2506andy(2) 1−0.375−0.1250.250.43751.00.0=++2+=,+2−0.2502.25003.02.01.5002.256so we have hit the value y(1) = 0 instead of the desired value y(1) = 1.
We now havethe initial slope bracketed between −1 and 1. We omit the further iterations necessary toidentify the correct initial slope, which turns out to be y 0 (0) = 0: 100.00.000.18750.3750.125(1)y =++2+=+200.75001.5000.7500.00.756andy(2) 10.1250.3750.750.93751.51.0=++2+=,+23.00.7501.5002.252.25003.06312CHAPTER 10.
BOUNDARY VALUE PROBLEMS FOR ODES2.01.51.00.50.0 •−0.5• ← first attempt.................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................• ← target••• ← second attempt•0.51.0Figure 10.2: Shooting method for two-point boundary value problem in Example 10.1.so we have indeed hit the target solution value, y(1) = 1.
These results are illustrated inFig. 10.2.A potential difficulty with the shooting method is that the initial value problem may beill-conditioned, perhaps owing to diverging solution curves over part of the domain, whichmay make it difficult to hit the desired target. A remedy is provided by multiple shooting,in which the interval [a, b] is divided into subintervals and shooting is carried out overeach subinterval.
Requiring continuity at the internal mesh points provides the boundaryconditions for the individual subproblems. Therefore, multiple shooting results in a systemof nonlinear equations to solve rather than just a single scalar equation.10.3Superposition MethodAnother way of replacing boundary value problems with initial value problems is the superposition method . Consider the homogeneous, linear, second-order ODEy 00 = p(t)y + q(t)y 0 ,a ≤ t ≤ b,with boundary conditionsy(a) = α,y(b) = β.The solution can be expressed as a superposition (i.e., linear combination) of two independent solutions, which can be obtained numerically by solving the equation with each of thetwo sets of initial conditionsy(a) = 1,y 0 (a) = 0and y(a) = 0,y 0 (a) = 1.This method becomes somewhat more complicated if the equation is inhomogeneous, andmuch more complicated if the equation is nonlinear.
Moreover, it may be necessary to useorthogonalization to maintain independence of the solutions computed for the initial valueproblems.10.4Finite Difference MethodBoth the shooting and superposition methods convert boundary value problems into initialvalue problems. We now consider methods that approximate boundary value problems10.4. FINITE DIFFERENCE METHOD313directly by systems of algebraic equations. Finite difference methods convert boundaryvalue problems into systems of algebraic equations by replacing any derivatives that appearwith finite difference approximations. For example, to solve the two-point boundary valueproblemy 00 = f (t, y, y 0 ), a ≤ t ≤ b,with boundary conditionsy(a) = α,y(b) = β,we first divide the interval [a, b] into n equally spaced subintervals.
Let ti = a + ih, i =0, 1, . . . , n, where h = (b − a)/n. We seek an approximation yi ≈ y(ti ) at each of the meshpoints ti . We already have y0 = α and yn = β. We next replace the derivatives with finitedifference approximations (see Section 8.7.1), such asy 0 (ti ) ≈yi+1 − yi−12hand y 00 (ti ) ≈yi+1 − 2yi + yi−1,h2choosing the finite difference formulas so that they have the same order truncation error,in this case O(h2 ), since the accuracy will be limited by the least accurate formula. Thisreplacement yields a system of algebraic equationsyi+1 − yi−1yi+1 − 2yi + yi−1 − h2 f ti , yi ,=02hto be solved for the unknowns yi , i = 1, .
. . , n − 1. This system of equations may be linearor nonlinear, depending on whether f is linear or nonlinear in y and y 0 . In this example,each equation in the system involves only three adjacent unknowns, which means that thematrix of the linear system—or the Jacobian matrix in the nonlinear case—is tridiagonal,thereby saving on both work and storage compared with a general system of equations. Suchsavings are generally true of finite difference methods: they yield sparse systems becauseeach equation involves only a few variables.Example 10.2 Finite Difference Method. We demonstrate the finite difference methodby using it to solve the two-point boundary value problem of Example 10.1,y 00 = 6t,0 ≤ t ≤ 1,with boundary conditionsy(0) = 0and y(1) = 1.To illustrate the concepts involved, yet keep computation to a minimum, we will computean approximate solution at a single mesh point in the interval [0, 1], namely t = 0.5. Thus,including the boundary points, we have three mesh points, t0 = 0, t1 = 0.5, and t2 = 1.From the boundary conditions, we know that y0 = y(t0 ) = 0 and y2 = y(t2 ) = 1, andwe seek an approximate solution y1 ≈ y(t1 ).
Approximating the second derivative by astandard finite difference quotient at the point t1 gives the equationy2 − 2y1 + y0= f (t1 , y1 , y10 ).h2314CHAPTER 10. BOUNDARY VALUE PROBLEMS FOR ODESSubstituting the boundary data, mesh size, and right-hand side for this example, we obtain1 − 2y1 + 0= 6t1 ,(0.5)2or4 − 8y1 = 6(0.5) = 3,so that1= 0.125,8which agrees with the approximate solution at t = 0.5 that we computed by the shootingmethod in Example 10.1.y(0.5) ≈ y1 =In a practical problem, a much smaller stepsize and many more mesh points wouldbe required to achieve acceptable accuracy, and we would therefore obtain a system ofequations to solve for the approximate solution values at the mesh points, rather than asingle equation as in this example. Nevertheless, this system would still be easy to solvebecause it would be tridiagonal.10.5Finite Element MethodAnother approach to reducing boundary value problems to algebraic systems is the finiteelement method.
Finite element methods approximate the solution to a boundary valueproblem by a linear combination of a finite collection of basis functions φi , typically piecewise polynomials, which for historical reasons are called “elements.” The approximationtherefore has the formnXy(t) ≈ u(t) =xi φi (t).i=1The coefficients xi are determined by imposing one of several possible requirements on theresidual, which is defined, as usual, to be the difference between the left and right sidesof the differential equation. For this reason, these methods are also known as weightedresidual methods.
Each of the three most commonly used criteria leads to a different classof methods:• Collocation: The residual is zero (i.e., the differential equation is satisfied exactly) at ndiscrete points.• Galerkin: The residual is orthogonal to the space spanned by the basis functions.• Rayleigh-Ritz : The residual is minimized in a weighted least squares sense.The latter two criteria are often equivalent. It may be helpful in understanding them torecall Fig. 3.2: the true solution to the differential equation does not in general lie in thespace spanned by the basis functions, so we seek an approximate solution (i.e., a linearcombination of basis functions) such that the residual is minimized, or is orthogonal to thespace spanned by the basis functions. Because they are based on an inner product on afunction space, these two methods involve the computation of integrals, either analyticallyor by some quadrature rule.10.5.
FINITE ELEMENT METHOD315Each of these three criteria leads to a system of equations to be solved for the coefficientsxi . The system of equations may be linear or nonlinear, depending on whether f is linearor nonlinear. The system will be sparse, and hence require much less work and storage,if the elements are “local,” which means that each basis function is zero throughout mostof the domain of the problem and that they have little overlap. Typical examples in onedimension are B-splines, such as the piecewise linear “hat” functions (see Section 7.3.4).The resulting sparse matrix of the system, called the stiffness matrix , is assembled elementby element and is a sum of contributions from each element.A related family of methods, called spectral methods, uses eigenfunctions of the differential operator as basis functions for expanding the approximate solution (e.g., trigonometricseries for the second derivative operator).
Similar use of other basis functions, such asLegendre or Chebyshev polynomials, leads to a pseudospectral method .Example 10.3 Collocation Method. We first illustrate the finite element method byusing collocation to solve the two-point boundary value problem of Example 10.1,y 00 = 6t,0 ≤ t ≤ 1,with boundary conditionsy(0) = 0and y(1) = 1.With the finite element method, we approximate the solution of the ODE by a functionrather than by a table of approximate values. Specifically, using the collocation method,we seek a function u(t) that satisfies the boundary conditions and also satisfies the ODEexactly at a discrete set of mesh points in the interval.
Again, for simplicity, we will useonly one interior mesh point, namely t = 0.5. For illustrative purposes, the function wechoose is a quadratic polynomial represented in the monomial basis, sou(t) = x0 + x1 t + x2 t2 .Note thatu0 (t) = x1 + 2x2 t,and u00 (t) = 2x2 .In determining the coefficients xi , we will enforce the boundary conditions at the endpointsof the interval, and the ODE at the point t = 0.5. For a general second-order two-pointboundary value problemy 00 = f (t, y, y 0 ), a ≤ t ≤ b,with boundary conditionsy(a) = αand y(b) = β,these requirements give the three equations,x0 + x1 a + x2 a2 = α,x0 + x1 b + x2 b2 = β,and u00 (t) = f (t, u(t), u0 (t)).Substituting the data and functions for this example, we obtain the system of three equationsx0 = 0,x1 + x2 = 1,and 2x2 = 6(0.5) = 3,316CHAPTER 10.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
