Shampine, Allen, Pruess - Fundamentals of Numerical Computing (523185), страница 47
Текст из файла (страница 47)
Output parameters are X, Y, the integration was advanced to X and Y isthe solution there; H, the step size suitable for the next step; FLAG, a flag reportingwhat the code did; STEP, the actual length of the step taken (output X minus inputX); and YCOEFF, an array of coefficient values for quintic Hermite interpolation to beused by the routine Yvalue.Some of the variables in the call list require more explanation. The initial step sizeH informs the code of the scale of the problem.
It should be small enough to capturechanges in the solution near the initial point. Also, the sign of H indicates the directionof integration because the code will try to step to X + H. After the first call, the codeprovides a suitable H for the next call.The variable FIRST enables the code to initialize itself. The start of a new problemis indicated by input of FIRST = .TRUE. in FORTRAN and first = 1 in C or C++.
Afterthe first call, the code sets FIRST = .FALSE. in FORTRAN and first = 0 in C and C++for subsequent calls. The error parameters TOL and THRES (or threshold in C andC++) tell the code how accurately the solution is to be computed. The vector THRESmust have dimension at least NEQ in the calling program. All components of THRESmust be nonnegative. The relative error tolerance TOL must satisfy10u < TOL < 0.016.6 THE CODE RKE237where u is the unit roundoff of the machine. The tolerances are used by the code in alocal error test at each step that requires roughly that|local error| < TOL max(|Y(I)|, THRES(I))for component I of the vector Y.
Setting THRES(I ) = 0 results in a pure relative errortest on the component. On the first call to the code, if some Y(I ) = 0, the correspondingTHRES(I) must be strictly positive. The size of the solution component is carefullydefined so that vanishing of the true solution at the current step is very unlikely tocause trouble. Any such trouble can be avoided by a positive value of THRES(I ) .The code will not attempt to compute a solution at an accuracy unreasonable forthe computer being used. It will report if this situation arises. To continue the integration after such a report, TOL and/or THRES must be increased. Note that Rke is anefficient code for moderate relative accuracies.
For more than, say, six-digit accuracy,other methods are likely to be more efficient.The true (global) error is the difference between the true solution of the initialvalue problem and the computed one. Nearly all initial value codes, including this one,control only the local error at each step, not the global error. Moreover, controllingthe local error in a relative sense does not necessarily result in the global error beingcontrolled in a relative sense. Roughly speaking, the codes produce a solution Y(x)that satisfies the differential equation with a residual R(x),that is usually bounded in norm by the error tolerances.
Usually the true accuracyof the computed Y is comparable to the error tolerances. This code will usually, butnot always, deliver a more accurate solution when the problem is solved again withsmaller tolerances. By comparing two such solutions a fairly reliable idea of the trueerror in the solution at the larger tolerances can be obtained.The principal task of the code is to integrate one step from X toward X + H. Routine Rke is organized so that subsequent calls to continue the integration involve little(if any) additional effort. The status of the integration is reported by the value of theFLAG parameter.
After a successful step the routine Yvalue is used to approximatethe solution at points within the span of the step. A typical call isCALL YVALUE(NEQ, X, STEP, YCOEFF, POINT, U)in the FORTRAN version andYvalue(neq, x, step, ycoeff, point, u);in the C and C++ versions. Input parameters are NEQ, X, STEP, YCOEFF (as returnedfrom Rke) and POINT, the point at which a solution is desired. The output is U(*), thevector of solution components at P. Routine Yvalue can be used only after a successfulstep by Rke and should be used only to interpolate the solution values on the interval[ X - STEP,X].238CHARTER 6ORDINARY DIFFERENTIAL EQUATIONSExample 6.8.To illustrate Rke, we solve the initial value problemY´1Y´2= Y1,= -Y2,Y l(0)Y 2(0)==11on the interval [0,1] and print the solution at x = 1.
The problem has the solutionYl(x) = ex, Y2(x) = e-x. As the solution component Y1 is increasing, a pure relativeerror test is appropriate and we set THRES(1) = 0. Y2, on the other hand, is decreasingand we choose THRES(2) = 10-5, which results in an absolute error test for small |Y2|.1.00000000000000XOUT =The numerical solution is 2.7182755628071 3.6787784616084E-012.71828182845903.6787944117144E-01The true solution isExample 6.9.
This example i llustrates the use of Yvalue in conjunction with Rke.The initial value problemis integrated over the interval [0,10] and the solution tabulated at x = 0, 1, 2,. . . , 10.Note that Rke must be called before Yvalue. The output is as follows.XOUTXOUTXOUTXOUTXOUTXOUTXOUTXOUTXOUTXOUTXOUT= 0.00= 1.00= 2.00= 3.00= 4.00= 5.00= 6.00= 7.00= 8.00= 9.00=10.00Y(1)Y(1)Y(1)Y(1)Y(1)Y(1)Y(1)Y(1)Y(1)Y(1)Y(1)= 1.000000= 1.298484= 0.421178=-1.634813=-1.743960=-0.878664= 1.187072= 1.933030= 1.245572=-0.329599=-2.008260Y(2)Y(2)Y(2)Y(2)Y(2)Y(2)Y(2)Y(2)Y(2)Y(2)Y(2)= 1.000000=-0.367034=-1.488951=-1.485472= 0.568922= 1.258102= 2.521700=-0.406833=-0.963183=-2.467240=-0.034172EXERCISESUnless otherwise indicated, use 10-5 for tolerancevalues and 10-7 for the threshold values for the computing exercises.in Exercise 6.4.(a) y´ = -½y3, y(0) = 1; b = 3(b) y´ = -2xy2, y(0) = 1; b = 1-66.16 Use Rke with a tolerance of 10 and threshold valuesof 1 to calculate y(b) in the following cases.
Check thecomputed output with the exact answers that are given(c) y´ = ¼y(1 - y/20), y(0) = 1; b = 5(d) y´ = 100(sinx - y), y(0) = 0; b = 16.6 THE CODE RKE239(e) y´ = (15 cos 10x)/y, y(0) = 2; b = π/4.6.19 Use Rke to approximate the solution to the initialvalue problemAre the true (global) errors within the toleranceon the local errors? Which problems needed the mostP´(t) = 0.001P(t)[1000(1 - 0.3 cos πt/6) -P(t)],steps? Why do you think they are more difficult thanthe others?with P(0) = 250; sketch its graph for 0 < t < 36.
(If6.17 An important equation of nonlinear mechanics is vanyou have worked Exercise 5.10, compare results.)der Pol’s equation:6.20 The response of a motor controlled by a governor canx´´(t) + ε( x2(t) - 1)x´(t) + x(t) = 0be modeled byfor ε > 0. Regardless of the initial conditions, all solutions of this equation converge to a unique periodicsolution (a stable limit cycle).
For ε = 1, choose someinitial conditions t0, x(t0), and x´(t0) and integrate theequation numerically until you have apparently converged to the limit cycle. A convenient way to viewthis is to plot x´(t) against x(t), a phase plane plot. Inthe phase plane, a periodic solution corresponds to aclosed curve.6.18 Deriving the equations for the quintic interpolatingpolynomial used in Rke/Yvalue is not difficult.
Write2345U(p) = a + bz + cz + dz + ez + fzfor p in [x - ∆, x], whereThe motor should approach a constant (steady-state)Assume s(0) = s´(0) = u´(0) =speed asθ (0) = 0, u(0) = 50, v(0) = w(0) = x(0) = 75.(a) Evaluate v(t) for t = 0, 25, 50, 75, 100, 150, 200,250, 300, 400, 500.(b) What doesappear to be? You cancheck this by working out the steady-state solution(the constant solution of the differential equation).6.21 Consider the initial value problem(a) Apply the interpolation conditionsy´´´ - y´´ sinx - 2y´ cosx + ysinx = lnx,y(1) = A1,y´(1) = A2,y´´(1) = A3.Show that the solution y(x) satisfies the first integralrelationy´´(x) - y´(x) sinx - y(x) cosx = c2 + xl nx - xto generate six linear equations in the six unknowns a,b, c, d, e, and f.and the second integral relation(b) Solve the linear system to get(In the code, α =and γ = UR-UL .)What are c1, c2 in terms of A1, A2, A3? Choose valuesfor A1, A2, and A3 and integrate this problem numerically.
Monitor the accuracy of your solution by seeinghow well it satisfies the integral relations. Argue thatif the integral relations are nearly satisfied, then thenumerical solution may or may not be accurate, butthat if they are not satisfied, the numerical solutionmust be inaccurate.240CHAPTER 6ORDINARY DIFFERENTIAL EQUATIONS6.22 The Jacobian elliptic functions sn(x), cn(x), and dn(x)satisfy the initial value problemy´1y´2y´3===and the following table:y 2 y 3 , y1(0) = 0-y 1 y 3 , y2(0) = 1-k 2 y 1 y 2 , y 3 (0) = 1where k2 is a parameter between 0 and 1 and y1(x) = 6.23 A simple model of the human heartbeat givessn(x), y2(x) = cn(x), and y3(x) = dn(x).ε x´ = -(x3 - Ax + c)Evaluate these functions numerically.
Check youraccuracy by monitoring the relationsc´ = x,sn2(x) + cn2(x) = 1dn2(x) + k2sn2(x) = 1dn2(x) - k2cn2(x) = 1 - k2.where x(t) is the displacement from equilibrium of themuscle fiber, c(t) is the concentration of a chemicalcontrol, and ε and A are positive constants. Solutionsare expected to be periodic. This can be seen by plotting the solution in the phase plane (x along the horizontal axis, c on the vertical), which should produce aclosed curve. Assume that ε = 1.0 and A = 3.Argue that if these relations are well satisfied numerically, you cannot conclude that the computed functions are accurate, rather that their errors are correlated. If the relations are not satisfied, the functions(a) Calculate x(t) and c(t) for 0 < t < 12 starting withmust be inaccurate. Thus, this test is a necessary testx(0) = 0.1, c(0) = 0.1.
Sketch the output in the phasefor accuracy but it is not sufficient.plane. What does the period appear to be?The Jacobian elliptic functions are periodic. You2(b) Repeat (a) with x(0) = 0.87, c(0) = 2.1.can get the true solutions for k = 0.51 from the factthat the period is 4K, where K = 1.86264 08023 6.24 Devise a step size strategy for Euler’s method with a32738 55203 02812 20579 ···. If tj = jK, j =local error estimator based on Heun’s method.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
