CH-06 (523176), страница 3
Текст из файла (страница 3)
It can be shownthat the final form of the system of four first-order differential equations are:dx1= x2 ,dtdx3= x4 ,dt© 2001 by CRC Press LLC22dx 2 gL (1 − 2 rm sin θ) + 4 rm v sin θ,=2dtL 1 + 4 r sin θ(dx 4 −2 vdx=+ 2 2 sin θdtLdt2and)(30)FIGURE 4. A numerical case where b = 3.6 m, h = .15 m, m2/m1 = 0.8, and initial conditionsy = z = dy/dt = dz/dt = 0 has been investigated and the results for -y, z, -dy/dt, and dz/dt havebeen plotted for 0≤t≤12 seconds.where:rm =m2m1anddLv2 = x 2 sin θ − x 2 dt(31,32)A numerical case where b = 3.6 m, h = .15 m, m2/m1 = 0.8, and initial conditionsy = z = dy/dt = dz/dt = 0 has been investigated and the results for -y, z, -dy/dt, anddz/dt have been plotted in Figure 4 for 0≤t≤12 seconds. The reason that negativevalues of y and dy/dt are used is that the mass m1 is moving downward as positive.It can be observed from Figure 4 that y varies between 0 and 3.0704 m, and z variesbetween 0 and 3.8358 m.
The oscillation has a period approximately equal to2x6.6831 = 13.3662 seconds, so the frequency is about 0.47 rad/sec. The massesreach their maximum speeds, dy/dtmax = 1.4811 and dz/dtmax = 1.8336 in m/sec,when y = .5ymax = 1.5352 and z = .5zmax = 1.9179 m, respectively. Details for theoscillation for the studied period are listed below:© 2001 by CRC Press LLCMATLAB APPLICATIONMATLAB has a file called ode45.m which implement the fourth- and fifthorder Runge-Kutta integration. Here, we demonstrate how the sample problem usedin the FORTRAN and QuickBASIC versions can also be solved by use of the mfile. The forcing functions given in Equations 10 and 11 are first prepared as follows:© 2001 by CRC Press LLCIf this file FunF.m is stored on a disk which has been inserted in disk drive A,ode45.m is to be applied with appropriate initial conditions and a time interval ofinvestigation as follows:>> T0 = 0; Tend = 3; X0 = [7;–4]; [T,X] = ode45(‘A:Funf’,T0,Tend,X0); plot(T,X)Notice that ode45 has four arguments.
The first argument is the m file in whichthe forcing functions are defined. The second and third argument the initial and finalvalues of time, respectively. The fourth argument is a vector containing the initialvalues of the dependent variables. The resulting display, after rearranging T and Xside-by-side for saving space instead of one after the other, is:The plots of X(1) and X(2) vs. T using the solid and broken lines, respectively,are shown in Figure 5.As another example, MATLAB is applied to obtain the displacement and velocity histories of a vibration system, Figure 1. First, a m file FunMCK.m is createdto describe this system as:Notice that the first and second variables X(1) and X(2) are displacement (x)and velocity (v = dx/dt), respectively, and the mass (m), damping constant (c), andspring constant (k) are taken as 2 N-sec2/cm, 3 N-sec/cm, and 4 N/cm, respectively.For a system which is initially at rest and disturbed by a constant force of F(t) = 1© 2001 by CRC Press LLCFIGURE 5.
The plots of X(1) and X(2) vs. T using the solid and broken lines, respectively.N, the MATLAB solution for 0≤t≤25 seconds shown in Figure 6 is obtained byentering the commands:>> T0 = 0; Tend = 25; X0 = [0;0]; X = [0;0]; [T,X] = ode45(‘a:FunMCK’,T0,Tend,X0)We can observe from Figure 6 that the mass has a overshoot (referring to Figure 2in the program NewRaphG) of about 0.28 cm at approximately t = 2.5 seconds andfinally settles to a static deflection of 0.25 cm, and that the maximum ascendingvelocity is about 0.18 cm/sec.As another example of dynamic analysis in the field of fluid mechanics, Figure 8shows the flow of a fluid between two connected tanks. The valve settings controlthe amount of flows, qi for i = 1,2,3.
The levels of the tanks h1 and h2 change intime depending on these settings and also on the supply rates Q1 and Q2 and thedischarge rate q3. Expressing the valve settings in terms of the resistances Ri for i =1,2,3, the conservation of masses requires that the flow rates be computed with theformulas:A1dh11= Q1 −(h − h3 )dtR1 1© 2001 by CRC Press LLCandA2dh 21= Q2 −(h − h3 )dtR2 2(33,34)FIGURE 6. For a system initially at rest and disturbed by a constant force of F(t) = 1 N,the MATLAB solution for 0≤t≤25 seconds shown here is obtained by entering the commandsshown below.where A’s are the cross-sectional areas of the tanks, and h3 is the pressure head atthe junction indicated in Figure 7 and is related to the discharge rate q3 by theequation:q 3 = q1 + q 2 = h 3 R 3 = ( h1 − h 3 ) R1 + ( h 2 − h 3 ) R 2(35)Or, h3 can be written in terms of R’s and h1 and h2 as:h 3 = R 3 ( R 2 h1 + R1h 2 ) ( R1R 2 + R1R 3 + R 2 R 3 )(36)By eliminating h3 terms from Equations 33 and 34, we obtain two differentialequations in h1 and h2 to be:© 2001 by CRC Press LLCFIGURE 7.
The flow of a fluid between two connected tanks.A1dh1= Q1 − a1h1 + a 3h 2dtandA2dh 2= Q 2 − a 2 h 2 + a 3h1dta1 = ( R 2 + R 3 ) ∆anda 2 = ( R1 + R 3 ) ∆(37,38)where:(39,40)and∆ = R1R 2 + R 2 R 3 + R 3R1(41)By assigning values for the parameters involved in the above problem, RungeKutta method can again be applied effectively for computing the fluid levels in bothtanks.4MATHEMATICA APPLICATIONSMathematica solves a set of ordinary differential equations based on the RungeKutta method by a function called NDSolve. The following run illustrates its interactive application using the same example in the MATLAB presentation of Figure 7:In[1]: = Id = (NDSolve[{X1’[t] = = X2[t], X2’[t] = = .5*(1–3*X2[t]–4*X1[t]),X1[0] = = 0, X2[0] = = 0},{X1,X2}, {t,0,25}])© 2001 by CRC Press LLCOut[1] = {{X1 -> InterpolatingFunction[{0., 25.}, <>],X2 -> InterpolatingFunction[{0., 25.}, <>]}}In[1] shows that NDSolve has three arguments: the first argument defines thetwo ordinary differential equations and the initial conditions, the second argumentlists the dependent variables, and the third argument specifies the independentvariable and the range of investigation.
Id is a name selected for easy later referenceof the results generated. Out[1] indicates that interpolation functions have beengenerated for X1 and X2 for t in the range from 0 to 25. To print the values of X1and X2, Id can be referred to interpolate the needs as follows:In[2]: = (Do[Print["t =", tv," X1 =", X1[tv]/. Id, " X2 =", X2[tv]/. Id],{tv, 0, 25, 1}])Out[2] = t = 0t=t=t=t=t=t=t=t=t=t=t=t=t=t=t=t=t=t=t=t=t=t=t=t=t=© 2001 by CRC Press LLCX1 = {0.}12345678910111213141516171819202122232425X1 =X1 =X1 =X1 =X1 =X1 =X1 =X1 =X1 =X1 =X1 =X1 =X1 =X1 =X1 =X1 =X1 =X1 =X1 =X1 =X1 =X1 =X1 =X1 =X1 =X2 = {0.}{0.13827}{0.267432}{0.280915}{0.256722}{0.245409}{0.246925}{0.249969}{0.250675}{0.250238}{0.249931}{0.249923}{0.24999}{0.250014}{0.250006}{0.249999}{0.249998}{0.25}{0.25}{0.25}{0.25}{0.25}{0.25}{0.25}{0.25}{0.25}X2 =X2 =X2 =X2 =X2 =X2 =X2 =X2 =X2 =X2 =X2 =X2 =X2 =X2 =X2 =X2 =X2 =X2 =X2 =X2 =X2 =X2 =X2 =X2 =X2 ={0.183528}{0.0629972}{–0.0193286}{–0.0206924}{–0.00278874}{0.00365961}{0.00187841}{–0.000172054}{–0.000477574}{–0.000125015}{0.0000639464}{0.0000489215}{4.83211 10–7}{–0.0000104109}{–3.5685 10–6}{1.38473 10–6}{1.35708 10–6}{1.39173 10–7}{–4.40766 10–7}{–4.61875 10–7}{–1.93084 10–8}{2.47763 10–8}{9.81364 10–8}{6.7364 10–8}{3.57288 10–8}These results are in agreement with the plotted curves shown in the MATLABapplication.
In[2] shows that the replacement operator/. is employed in X1[tv]/.Idwhich requires all t appearing in the resulting interpolating function for X1(t) createdin the Id statement to be substituted by the value of tv. The looping DO statementinstructs that the tv values be changed from a minimum of 0 and a maximum of 25using an increment of 1. To have a closer look of the overshoot region, In[2] canbe modified to yieldIn[3]: = (Do[Print[“t = “,tv,” X1 = “,X1[tv]/. Id,” X2 = “,X2[tv]/. Id],{tv, 2, 3, 0.1}])Out[3] =t=2t=t=t=t=t=t=t=t=t=t=X1 = {0.267432}2.12.22.32.42.52.62.72.82.93.X1 =X1 =X1 =X1 =X1 =X1 =X1 =X1 =X1 =X1 ={0.273097}{0.277544}{0.280862}{0.283145}{0.284496}{0.285019}{0.284819}{0.284002}{0.282669}{0.280915}X2 = {0.0629972}X2 =X2 =X2 =X2 =X2 =X2 =X2 =X2 =X2 =X2 ={0.0504262}{0.0386711}{0.0278364}{0.0179948}{0.00919032}{0.00144176}{–0.00525428}{–0.0109202}{–0.0155945}{–0.0193286}This detailed study using a time increment of 0.1 reaffirms that the overshootoccurs at t = 2.6 and X1 has a maximum value equal to 0.28502.6.3 PROGRAM ODEBVPRK — APPLICATION OF RUNGE-KUTTAMETHOD FOR SOLVING BOUNDARY-VALUE PROBLEMSThe program OdeBvpRK is designed for numerically solving the linear boundaryvalue problems governed by the ordinary differential equation by superposition oftwo solutions obtained by application of the Runge-Kutta fourth-order method.
Toexplain the procedure involved, consider the problem of a loaded beam shown inFigure 8. Mathematically, the deflection y(x) satisfies the well-known flexural equation.5d2y M=dx 2 EI(1)where M is the internal moment distribution, E is the Young’s modulus and I is themoment of inertia of the cross section of the beam. For the general case, M, E andI can be function of x. The boundary conditions of this problem are:© 2001 by CRC Press LLCFIGURE 8.
The problem of a loaded beam.y(x = 0) = 0andy(x = L ) = 0(2)where L is the length of the beam. The problem is to determine the resultingdeflection y(x). Knowing y, the moment and shearing force, V, distributions cansubsequently be determined based on Equation 1 and V = dM/dx. The final objectiveis to calculate the stress distributions in the loaded beam using the M and V results.The Runge-Kutta method for solving an initial problem can be applied here ifin additional to the initial condition given in (2), y(x = 0) = 0, we also know theslope, , at x = 0.
But, we can always make a guess and hope that by making betterand better guesses the trial process will eventually lead to one which satisfies theother boundary condition, namely y(x = L) = 0 given in Equation 3. In fact, if theproblem is linear, all we need to do is making two guesses and linearly combinethese two trial solutions to obtain the solution y(x).Let us first convert the governing differential Equation 1 into two first-orderequations as: dx1 dx = x 2 , dx 2 M= , dx EIx1 (x = 0) = y(x = 0) = 0(3,4)x 2 (x = 0) = θ(x = 0) = θ0To apply the fourth-order Runge-Kutta method, we have to first decide on astepsize h, for example h = L/N we then plan to calculate the deflections at N + 1locations, x + jh for j = 1,2,…N since j = 0 is the initial location.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
