CH-06 (523176), страница 2
Текст из файла (страница 2)
It has been coded according to the proceduredescribed in the preceding section. That is, the equations must be in the form ofEquation 3 by having the first derivatives of the dependent variables (x1 through xN)all on the left sides of the equations and the right sides be called F1 through FN.These functions are to be defined in a Function subprogram F.The FORTRAN version of Subroutine RKN is listed below.
There are sevenarguments for this subroutine, the first four are input arguments where the last is anoutput argument. The fifth argument P keeps the Runge-Kutta parameters generatedin this subroutine. The sixth argument XT is needed for adjusting the input argumentXIN. These two arguments, P and XT, are included for handling the general caseof N variables. Listing them as arguments makes possible to specify them as matrixand vector of adjustable sizes.FORTRAN VERSION© 2001 by CRC Press LLCQUICKBASIC VERSIONPROGRAM RUNGEKUTProgram RungeKut which calls the subroutine RKN is to be run interactivelyby specifying the inputs through the keyboard.
Displayed messages on screen instructuser how to input the necessary data and describe the problem in proper sequence.From the provided listing of the program RungeKut, user will find that the followinginputs and editing need to be executed in the sequence specified:(1) Input the number of variables, N, involved.(2) Define the N functions on the right sides of Equation 3, F1 through FN,by editing the DEF statements starting from statement 161. Presently only9 functions can be handled by the program RungeKut, but the user shouldbe able to expand the program to accommodate any N value which isgreater than 9 by renumbering the program and adding more DEF statements.(3) Type RUN 161 to run the program.(4) Reenter the N value.(5) Enter the beginning (not necessary equal to zero!) and ending values ofthe independent variable, denoted as t0 and te (T0 and TEND in theprogram RungeKut), respectively.
It is over this range, the values of theN dependent variables are to be calculated.(6) Enter the stepsize, h (DT in the program RungeKut), with which theindependent variable is to be incremented.(7) Enter the N initial values.© 2001 by CRC Press LLCQUICKBASIC VERSIONFORTRAN VERSION© 2001 by CRC Press LLCFUNCTION FAccording to Equations 12 to 15, the Runge-Kutta parameters pi,j (matrix P inthe program RungeKut) are calculated using two FOR-NEXT loops — an I loopcovering N sets of variables and a J loop covering the four parameters in each set.As an illustrative example, let us apply program RungeKut for the problem definedby Equations 8 to 11. We create a supporting function program F as follows:QuickBASIC Version© 2001 by CRC Press LLCFORTRAN VersionThe computation can then commence by entering through keyboard the beginning and ending values of t, t0 = 0 and te = 3, respectively, the stepsize h = 0.1, andthe initial values x1,0 = 7 and x2,0 = –4.
The complete sequence of question-andanswer steps in running the program RungeKut for the problem described byEquations 8 to 11 is manifested by a copy of the screen display:© 2001 by CRC Press LLCSAMPLE APPLICATIONSOF THE PROGRAMRUNGEKUTAs a first example, consider the problem of a beam shown in Figure 2 which isbuilt into the wall so that it is not allowed to displace or rotate at the left end. Forconsideration of the general case when the beam is loaded by (1) a uniformlydistributed load of 1 N/cm over the leftmost quarter-length, x between 0 and 10 cm,(2) a linearly varying distributed load of 0 at x = 10 cm and 2 N/cm at x = 20 cm,(3) a moment of 3 N-cm at x = 30 cm, and (4) a concentrated load of 4 N at thefree end of the beam, x = 40 cm. It is of concern to the structural engineers to knowhow the beam will be deformed.
The equation for finding the deflected curve of thebeam, usually denoted as y(x), is:3()EI d 2 y dx 2 = M(x)(21)where EI is the beam stiffness and M(x) is the variation of bending moment alongthe beam. It can be shown that the moment distribution for the loading shown inFigure 4 can be described by the equations:−1121 3 + 24 x − x 2 2,−871 3 + 4 x + x 2 − x3 30,M=4 x − 157,4 x − 160,© 2001 by CRC Press LLCforforforfor0 x 10 cm10 x 20 cm20 < x 30 cm30 < x 40 cm(22)FIGURE 2.
The problem of a beam, which is built into the wall so that it is not allowed todisplace or rotate at the left end.To convert the problem into the form of Equation 3, we replace y, dy/dx, and xby x1, x2, and t, respectively. Knowing Equation 21 and that the initial conditionsare y = 0 and dy/dx = 0 at x = 0, we can obtain from Equation 21 the followingsystem of first-order ordinary differential equations:dx1 dx = F1 (x, x1, x 2 ) = x 2 ,x1 (x = 0) = 0(23)anddx 2 dx = F2 (x, x1, x 2 ) = M EI,x 2 (x = 0) = 0(24)For Equation 24, the moment distribution has already been described by Equation21 whereas the beam stiffness EI is to be set equal to 2x105 N/cm2 in numericalcalculation of the deflection of the beam using program RungeKut.To apply program RungeKut for solving the deflection equation, y(x) which hasbeen renamed as x1(t), we need to define in the QuickBASIC program RungeKuttwo functions:DEF FNF1(X) = X(2)DEF FNF2(X) = BM(X)/(2*10^5)Notice that the bending moment M is represented by BM in QuickBASICprogramming.
In view of Equation 21, F2 in Equation 24 has to be defined bymodifying the subprogram function, here we illustrate it with a FORTRAN version:© 2001 by CRC Press LLCFUNCTION BM(X)IF ((X.LE.0.).AND.(X.LT.10.)) BM = –1121./3 + 24*X-X**2/2IF ((X.GE.10.).AND.(X.LT.20.)) BM = –871./3 + 4*T +T**2T**3/30IF ((X.GE.20.).AND.(X.LT.30.)) BM = 4*X–157IF ((X.GE.30.).AND.(X.LT.40.)) BM = 4*X–160RETURNENDIn fact, each problem will require such arrangement because these functionstatements and subprogram function describe the particular features of the problembeing analyzed.
The computed deflection at the free end of the beam, y(x = 40 cm),by application of the program RungeKut using different stepsizes, and the errorsin % of the analytical solution ( = –0.68917) are tabulated below:Stepsize h, cmy at x = 40 cm,521.5.25.1.05Error, %–0.79961–0.73499–0.71236–0.70083–0.69501–0.69151–0.6903316.06.653.361.69.847.340.168This problem can also be arranged into a set of four first-order ordinary differential equations and can be solved by using the expressions for the distributed loadsdirectly. This approach saves the reader from deriving the expressions for the bendingmoment. Readers interested in such an approach should solve Problem 4.A NONLINEAR OSCILLATION PROBLEM SOLVEDBYRUNGEKUTThe numerical solution using the Runge-Kutta method can be further demonstrated by solving a nonlinear problem of two connected masses m1 and m2 as shownin Figure 3.
A cable of constant length passing frictionless rings is used. Initially,both masses are held at rest at positions shown. When they are released, theirinstantaneous positions can be denoted as y(t) and z(t), respectively. The instantaneous angle and the length of the inclined portion of the cable can be expressed interm of y as:[θ = tan −1 (y + h ) b]and[L = (y + h) + b22]0.5(25,26)If we denote the cable tension as T, then the Newton’s second law applied tothe two masses leads to m1g–2Tsinθ = m1d2y/dt2 and T-m2g = m2d2z/dt2.
By eliminating T, we obtain:© 2001 by CRC Press LLCFIGURE 3. The numerical solution using the Runge-Kutta method can be further demonstrated by solving a nonlinear problem of two connected masses m1 and m2. 2m2d 2y 2m2d 2z+sin θ 2 = g1 −sin θ2dtm1dtm1(27)The displacements y and z are restricted by the condition that the cable lengthmust remain unchanged. That is z = 2{[(y + h)2 + b2]1/2[h2 + b2]1/2}. This relationshipcan be differentiated with respect to t to obtain another equation relating d2y/dt2 andd2z/dt2 which is:2d 2 z −2dLdy 2 dy d2y = 2 (y + h)+ + (y + h) 2 2dt dtLdtdt L dt (28)Equation 28 can be substituted into Equation 27 to obtain: 2 2 m 2d2ysin θ=gL −2dtL(2 y + h + L ) m12 2dLdy dy gLyhL−++22() dt dtdt(29)By letting x1 = y, x2 = dy/dt, x3 = z, and x4 = dz/dt, then according to Equation3 we have F1 = x2, F2 to be constructed using the right-hand side of Equation 29,F3 = x4, and F4 to be constructed using the right-side of Equation 28.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
