CH-06 (523176), страница 4
Текст из файла (страница 4)
If we assume avalue for 0, say A, the Runge-Kutta process will be able to generate the followingtablexx 1=y (1)00hy1(1)2hy2(1)……jhyj(1)……NhyN(1)x 2=θ (1)Aθ1(1)θ2(1)…θj(1)…θN(1)© 2001 by CRC Press LLC(5)If yN(1) = 0, then the value A selected for (x = 0) is correct and the y and values listed in Equation 5 are the results for the selected stepsize. If yN(1) is notequal to zero, then the value incorrectly selected, we have to make a second try byletting (x = 0) = B to obtain a second table by application of the Runge-Kuttamethod. Let the second table be denoted as:xx 1=y (2)00hy1(2)2hy2(2)……jhyj(2)……NhyN(2)x 2=θ (2)Bθ1(2)θ2(2)…θj(2)…θN(2)(6)Again, if yN(2)= 0, then the value B selected for θ(x = 0) is correct and the y andθ values listed in Equation 6 are the results for the selected stepsize.
Otherwise, ifthe problem is linear, the solutions can be obtained by linearly combining the twotrial results as, for j = 1,2,…,N:Yj = αy(j1) + βy(j2 )andθ j = αθ(j1) + βθ(j2 )(7,8)where the weighting coefficients and are to be determined by solving the equations:α +β =1andαy(N1) + βy(N2 ) = 0(9,10)Equation 10 is derived from the boundary condition y(x = Nh = L) = 0 and basedon Equation 7. Equation 9 needs more explanation because it cannot be derived ify(x = 0) = 0. Let us assume that for the general case, y(x = 0) = . Then Equation7 gives + = which leads to Equation 9. Using Cramer’s rule, we can easilyobtain:α = y(N2 ) Dandβ = − y(N1) D(11,12)and(1)D = y(2N − yN(13)NUMERICAL EXAMPLESLet us consider the problem of a loaded beam as shown in Figure 8. Thecrosssection of the beam has a width of 1 cm and a height of 2 cm which results ina moment of inertia I = 2/3 cm4.
The reactions at the left and right supports can becomputed to be 5/3 N and 25/3 N, respectively. Based on these data, it can be shownthat the equations for the internal bending moments are:for 0 ≤ x ≤ 20 cm5x 3,M( x ) = 651 2for 20 < x ≤ 30 cm−200 + 3 x − 2 x ,© 2001 by CRC Press LLC(14)If the beam has a Young’s modulus of elasticity E = 2x107 N/cm2, we may decideon a stepsize of h = 1 cm and proceed to prepare a computer program using thefourth-order Runge-Kutta method to generate two trial solutions and then linearlycombining to arrive at the desired distributions of the deflected shape y(x) and slope(x).
The FORTRAN version of this program called OdeBvpRK to be presentedlater has produced the following display on screen:© 2001 by CRC Press LLCFORTRAN VERSION© 2001 by CRC Press LLCThe Subprogram FUNCTION F which defines the initial-value problem is codedin accordance with Equation 14. The two trial initial slopes are selected as equal to0.1 and 0.2. The trial results are kept in the three-dimensional variable XT, in whichthe deflection y(k)(j) for the kth try at station x = xj = jh is stored in XT(1,j,k) whereasthe slope there is stored in XT(2,j,k) for j = 1,2,…,30 and k = 1,2. Such a threesubscripts arrangement facilitates the calling of the subroutine RKN becauseXT(1,KS,NTRY) is transmitted as XIN(1) and automatically the next valueXT(2,KS, NTRY) as XIN(2), and the computed results XOUT(1) and XOUT(2) areto be stored as XT(1,KS + 1,NTRY) and XT(2,KS + 1,NTRY), respectively.
Noticethat there are only two dependent variables, NV = 2.After the weighting coefficients (ALPHA in the program) and (BETA) havebeen calculated, the final distributions of the deflection and slope are saved in firstand second rows of the two-dimensional variable X, respectively. It should beemphatically noted that the solutions obtained is only good for the selected stepsizeh = 1 cm. Whether it is accurate or not remains to be tested by using finer stepsizesand by repeated application of the Runge-Kutta methods.It can be shown that the maximum deflection of the loaded beam is equal to–2.019 cm and is obtained when the stepsize is continuously halved and two consecutively calculated values is different less than 0.0001 cm in magnitude.
Theneeded modification of the above listed program to include this change in the stepsizeand testing of the difference in the maximum deflection is left as homework for thereader.© 2001 by CRC Press LLCQUICKBASIC VERSIONMATLAB APPLICATIONSIn the program RungeKut, MATLAB is used for solving initial value problemsby application of its m function ode45 based on the fourth- and fifth-order RungeKutta method. Here, this function can be employed twice to solve a boundary-valueproblem governed by linear ordinary differential equations. To demonstrate theprocedure, the sample problem discussed in FORTRAN and QuickBASIC versionsof the program OdeBvpRK, with the aid of function BVPF.m listed in the subdirectory <mFiles>, can be solved by interactive MATLAB operations as follows:© 2001 by CRC Press LLCNotice that format compact enables the display to use fewer spacings; YT1 andYT2 keep the two trial solutions, and ode45 automatically determines the beststepsize which if used directly will result in a coarse plot as shown in Figure 9.
Theplot showing solid-line curve for the deflection and broken-line curve for the slopecan, however, be refined by linear interpolation of the data (X,Yanswer) and expanding X and the two columns of Yanswer into new data arrays of XSpline, YSpline,and YPSpline, respectively. Toward that end, the m function spline in MATLAB isto be applied as follows:© 2001 by CRC Press LLCFIGURE 9.The result of plotting the spline curves is shown in Figure 10.MATHEMATICA APPLICATIONSThe Runge-Kutta method, particularly the most popular fourth-order method,can be applied for solution of boundary-value problem governed by ordinary differential equation(s). Here, only the application of this method is elaborated; readersare therefore referred to program RungeKut to review the method itself and thedevelopment of related programs and subprograms.
The boundary-value problem isto be solved by continuously guessing the initial condition(s) which are not provideduntil all boundary conditions are satisfied if the problem is nonlinear. When theproblem is linear, then only a finite number of guesses are necessary. A system oftwo first-order ordinary differential equations which governed the loaded elasticbeam problem previously solved in the MATLAB application is here adapted todemonstrate the application of the Runge-Kutta method.© 2001 by CRC Press LLCFIGURE 10. The result of plotting the spline curves.In[1]: = EI = 4.*10^7/3; F[X_] = If[X>20,(–200 + 65*X/3X^2/2)/EI, 5*X/3/EI]In[2]: = Id1 = (NDSolve[{Y[X] = = YP[X], YP'[X] = = F[X], Y[0] = = 0,YP[0] = = 0.1}, {Y,YP}, {X, 0, 30}])In[3]: = Y30Trial1 = Y[30]/.
Id1Out[3] = {3.00052}EI value and F(X) are defined in In[1]. In[2] specifies the two first-order ordinarydifferential equations involving the deflection Y and slope YP, describes the correctinitial condition Y(X = 0) = 0, gives a guessed slope Y(X = 0) = 0.1, and decideson the limit of investigation from X = 0 to X = 30. In[3] interpolates the ending Yvalue by using the data obtained in Id1. A second trial is then to follow as:In[4]: = Id2 = (NDSolve[{Y[X] == YP[X], YP'[X] = = F[X],Y[0] == 0, YP[0] == 0.2},{Y,YP}, {X, 0, 30}])© 2001 by CRC Press LLCIn[5]: = Y30Trial2 = Y[30]/.
Id2Out[5] = {6.00052}Linear combination of the two trial solutions is now possible by calculating acorrect Y(X = 0) which should be equal to the value given in Out[8].In[6]: = d = Y30Trial2Y30Trial1; a = Y30Trial2/d; b = -Y30Trial1/d;In[7]: = Print[“Alpha = “,a,” Beta = “,b]Out[7] = Alpha = {2.00017} Beta = {–1.00017}In[8]: = YP0 = 0.1*a + 0.2*bOut[8] = {–0.0000174888}Finally, the actual deflection and slope can be obtained by providing the correctset of initial conditions and again applying the Runge-Kutta method.In[9]: = Id = (NDSolve[{Y[X] = = YP[X], YP'[X] = = F[X],Y[0] = = 0, YP[0] = = –0.0000174888},{Y,YP}, {X,0,30}])In[10]: = (Do[Print[“X = “, Xv, “ Y = “, Y[Xv]/.Id, “ DY/DX = “,YP[Xv]/.Id], {Xv,0,30}])Out[10] =X=X=X=X=X=X=X=X=X=X=X=X=X=X=X=© 2001 by CRC Press LLC01234567891011121314Y=Y=Y=Y=Y=Y=Y=Y=Y=Y=Y=Y=Y=Y=Y={0.}{–0.0000174645}{–0.0000348041}{–0.0000518936}{–0.0000686077}{–0.0000848222}{–0.000100412}{–0.000115251}{–0.000129215}{–0.00014218}{–0.000154019}{–0.000164608}{–0.000173788}{–0.000181454}{–0.000187493}DY/DX =DY/DX =DY/DX =DY/DX =DY/DX =DY/DX =DY/DX =DY/DX =DY/DX =DY/DX =DY/DX =DY/DX =DY/DX =DY/DX =DY/DX ={–0.0000174888}{–0.0000174262}{–0.0000172387}{–0.0000169262}{–0.0000164887}{–0.0000159262}{–0.0000152387}{–0.0000144262}{–0.0000134887}{–0.0000124262}{–0.0000112387}{–0.00000992619}{–0.00000848869}{–0.00000692619}{–0.00000523869}X=X=X=X=X=X=X=X=X=X=X=X=X=X=X=X=15161718192021222324252627282930Y=Y=Y=Y=Y=Y=Y=Y=Y=Y=Y=Y=Y=Y=Y=Y={–0.000191786}{–0.00019422}{–0.00019468}{–0.000193037}{–0.000189163}{–0.000182912}{–0.000174247}{–0.00016319}{–0.00014959}{–0.000133573}{–0.00011547}{–0.0000951574}{–0.0000725124}{–0.0000483635}{–0.0000233292}{0.0000023099}DY/DX =DY/DX =DY/DX =DY/DX =DY/DX =DY/DX =DY/DX =DY/DX =DY/DX =DY/DX =DY/DX =DY/DX =DY/DX =DY/DX =DY/DX =DY/DX ={–0.00000342619}{–0.00000148869}{0.000000573812}{0.00000276006}{0.00000506839}{0.0000074989}{0.0000100026}{0.0000125002}{0.0000149974}{0.0000174176}{0.0000196352}{0.0000216325}{0.0000233628}{0.0000246913}{0.0000255217}{0.0000258107}6.4 PROGRAM ODEBVPFD — APPLICATIONOF FINITE DIFFERENCE METHODFOR SOLVING BOUNDARY-VALUE PROBLEMSThe program OdeBvpFD is designed for numerically solving boundary-valueproblems governed by the ordinary differential equation which are to be replacedfinite-difference equations.
To illustrate the procedure involved, let us consider theproblem of an annular membrane which is tightened by a uniform tension T andrigidly mounted along its inner and outer boundaries, R = Ri and R = Ro, respectively.As shown in Figure 11, it is then inflated by application of a uniform pressure p. Thedeformation of the membrane, Z(R), when its magnitude is small enough not toaffect the tension T, can be determined by solving the ordinary differential equation6Z(R) satisfies Equation 1 is for Ri<R<Ro and the boundary conditions.d 2 Z 1 dZp+=−dR 2 R dRTZ( R i ) = 0andZ( R o ) = 0(1)(2)If the finite-difference approximation is to be applied for solving Equation 1,we will be seeking not for the expression Z(R) but for the numerical values at aselected stations of R in the interval Ri<R<Ro, say N. Let these stations be designatedas Rk for k = 1 to N and the lateral displacements of the membrane as ZkZ(Rk).Using the first-order and second-order central differences (see the program DiffTabl),© 2001 by CRC Press LLCFIGURE 11.
An annular membrane tightened by a uniform tension T and rigidly mountedalong its inner and outer boundaries, R = Ri and R = Ro, respectively, and then inflated byapplication of a uniform pressure p.the first and second derivatives of Z(R) appearing in Equation 1 at Rk can beapproximated as, respectively:Z − Z k −1dZat R k k +1dR2∆R(3)Z − 2 Z k + Z k +1d2Zat R k k −12dR( ∆R)2(4)andwhere R is the increment in R and is related to the decided number of station Nby the equation:∆R = ( R o − R i ) ( N + 1)© 2001 by CRC Press LLC(5)Substituting Equations 3 and 4 into Equation 1, we obtain for R = Rk that:Z k −1 − 2 Z k + Z k +1 1 Z k +1 − Z k −1p−=−RkT2 ∆R( ∆R)2Multiplying both sides by (R)2 and collecting terms, we can have for k =1,2,…,N:2∆R ∆R − p( ∆R )1 + Z k −1 − 2 Z k + 1 − Z k +1 =T 2R k 2R k (6)The two boundaries are ZoZ(R = Ri) and ZN + 1Z(R = Ro), at which Z = 0.The above equation thus lead to the matrix equation:[C]{Z} = {R}(7)where:{Z} = [Z1 Z 2 … Z N ](8){R} = − p( ∆R)2 [1 1…1]T T(9)Tand if the elements for the coefficient matrix [C] are denoted as ck,j, based onEquation 6 they are to be calculated using the formulas:c k ,k = −2,for k = 1, 2,… Nc k ,k −1 = 1 − ( ∆R 2 R k ), for k = 2, 3,…, Nc k ,k +1 = 1 + ( ∆R 2 R k ), for k = 1, 2,…, N − 1c k , j = 0,(10)elsewhereandR k = k∆R(11)After having calculated [C] and {R}, Equation 7 can be solved by calling thesubroutine GAUSS for {Z}.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
