Pao Y.C. Engineering Analysis: Interactive Methods and Programs with FORTRAN, QuickBASIC, MATLAB, and Mathematica

Подождите немного. Документ загружается.

X = 15 Y = {–0.000191786} DY/DX = {–0.00000342619}

X = 16 Y = {–0.00019422} DY/DX = {–0.00000148869}

X = 17 Y = {–0.00019468} DY/DX = {0.000000573812}

X = 18 Y = {–0.000193037} DY/DX = {0.00000276006}

X = 19 Y = {–0.000189163} DY/DX = {0.00000506839}

X = 20 Y = {–0.000182912} DY/DX = {0.0000074989}

X = 21 Y = {–0.000174247} DY/DX = {0.0000100026}

X = 22 Y = {–0.00016319} DY/DX = {0.0000125002}

X = 23 Y = {–0.00014959} DY/DX = {0.0000149974}

X = 24 Y = {–0.000133573} DY/DX = {0.0000174176}

X = 25 Y = {–0.00011547} DY/DX = {0.0000196352}

X = 26 Y = {–0.0000951574} DY/DX = {0.0000216325}

X = 27 Y = {–0.0000725124} DY/DX = {0.0000233628}

X = 28 Y = {–0.0000483635} DY/DX = {0.0000246913}

X = 29 Y = {–0.0000233292} DY/DX = {0.0000255217}

X = 30 Y = {0.0000023099} DY/DX = {0.0000258107}

6.4 PROGRAM ODEBVPFD — APPLICATION

OF FINITE DIFFERENCE METHOD

FOR SOLVING BOUNDARY-VALUE PROBLEMS

The program OdeBvpFD is designed for numerically solving boundary-value

problems governed by the ordinary differential equation which are to be replaced

ﬁnite-difference equations. To illustrate the procedure involved, let us consider the

problem of an annular membrane which is tightened by a uniform tension T and

rigidly mounted along its inner and outer boundaries, R = R

and R = R

, respectively.

As shown in Figure 11, it is then inﬂated by application of a uniform pressure p. The

deformation of the membrane, Z(R), when its magnitude is small enough not to

affect the tension T, can be determined by solving the ordinary differential equation

Z(R) satisﬁes Equation 1 is for R

<R<R

and the boundary conditions.

(1)

(2)

If the ﬁnite-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 a

selected stations of R in the interval R

<R<R

, say N. Let these stations be designated

as R

for k = 1 to N and the lateral displacements of the membrane as ZkZ(R

Using the ﬁrst-order and second-order central differences (see the program DiffTabl),

dR R

+=−

ZR and ZR

()

=00

the ﬁrst and second derivatives of Z(R) appearing in Equation 1 at R

can be

approximated as, respectively:

(3)

and

(4)

where R is the increment in R and is related to the decided number of station N

by the equation:

(5)

FIGURE 11. An annular membrane tightened by a uniform tension T and rigidly mounted

along its inner and outer boundaries, R = R

and R = R

, respectively, and then inﬂated by

application of a uniform pressure p.

at R



+−

−

2∆

at R

ZZZ

kkk



−+

()

∆

∆RRRN

=−

()

Substituting Equations 3 and 4 into Equation 1, we obtain for R = R

that:

Multiplying both sides by (R)

and collecting terms, we can have for k =

1,2,…,N:

(6)

The two boundaries are Z

Z(R = R

) and Z

N + 1

Z(R = R

), at which Z = 0.

The above equation thus lead to the matrix equation:

(7)

where:

(8)

(9)

and if the elements for the coefﬁcient matrix [C] are denoted as c

k,j

, based on

Equation 6 they are to be calculated using the formulas:

(10)

and

(11)

After having calculated [C] and {R}, Equation 7 can be solved by calling the

subroutine GAUSS for {Z}. This solution is only for the selected value of N. N

should be continuously increased to test if the maximum Z value would be affected.

The program OdeBvpFD has been developed in both QuickBASIC and FOR-

TRAN versions based on the above described procedure for solving the boundary

problems governed by ordinary differential equations. It is designed for the general

case where the coefﬁcient matrix [C] and the right-hand side vector {R} of the

matrix equation derived from the ﬁnite-difference approximation, [C]{Z} = {R}, are

ZZZ

kkk

kk−++−

−+

()

−

=−

∆













−+−













−

()

−+

∆∆∆R

CZ R

[]

{}

ZZZZ

{}

=…

[]

RpR T

{}

=−

()

…

[]

∆

11 1

kk,

, , ,=− = …212for k

cRR N

kk k

,,,,

−

=−

()

=…

()

=…−

12 23

12 12 1

∆

for k

elsewhere

RkR

=∆

to be deﬁned by the user. In addition, the user needs to specify the number of station,

N, the stepsize, R, and the location of the left boundary, R

. Here, Z and R are

referred as the dependent and independent variables, respectively using the mem-

brane problem only for the convenience of explanation; user could have other

notations as dependent and independent variables.

As a second example, we consider the simply supported beam of length 3L

which has a solid section at its middle portion and hollow sections at both end

portions as shown in Figure 12. It is of interest to know how the beam will deﬂect

under its own weight. This is the case of a beam subjected to uniformly distributed

loads. For simplicity, let us assume that the distributed loads to be w

and w

, in

N/cm

, for the middle and end portions, respectively. To determine the deﬂection

y(x) for 0<x<3L, we need the following relevant equations which are available in a

standard textbook on mechanics of materials

(12–15)

where , M, and V are the slope, bending moment, and shearing force, respectively.

E is the Young’s modulus and I is the moment of inertia of the beam. The distributed

load w is considered as positive when it is applied upward in the direction of the

positive y-axis and the moment is considered as positive if it causes the beam to

bend concave up. The shearing forces are considered as positive if they are related

to M according to Equation 14.

By successive substitutions, Equations 12–15 can yield an equation which relates

the deﬂection y directly to the distributed load w as:

(16)

FIGURE 12. As a second example, we consider the simply supported beam of length 3L which

has a solid section at its middle portion and hollow sections at both end portions as shown

w=== = =θ

, , ,













This is a fourth-order ordinary differential equation and needs fourth supple-

mentary initial, boundary, or mixed conditions in order to completely solve for the

deﬂection y(x). For the simply supported beam subjected to its own weight, the four

boundary conditions are:

(17)

If EI is not a function of x, the central-difference approximation for Equation

16 at x = x

is:

(18)

If N in-between stations are selected for determination of the deﬂections there,

then the two boundaries are at k = 0 and k = N + 1. In view of the subscripts k–2

and k + 2 in Equation 18, k can thus only take the values 2 through N–1. For solving

through y

, we hence need two more equations which are the two moment

conditions in Equation 17. At x

= 0, the second-order, forward-difference formula

can be applied to give

Since y = 0 at the left support x = x

= 0, the above equation yields

(19)

Similarly, at the right support x = x

N + 1

= 3L, the second-order, backward-

difference formula can be applied to give

Since y = 0 at the right support x = x

N + 1

= 3L, the above equation yields

(20)

Meanwhile, the boundary condition y = 0 at the left support x = x

= 0 can be

substituted into Equation 18 for k = 2 to obtain:

(21)

y and M EI

at x and x L=====000

at x

yyyyy

at x

kkkkk

21 12

464



−− ++

−+−+

()

∆

MEI

at x

yyy

EI at x=

−+

()

210

0 

∆

−+=20

MEI

at x

yyy

EI at x

NnN

−+

()

+−

0 

∆

NN−

−=

−+ − +=

()

46 4

1234

yy yy x

at x∆

Similarly, the boundary condition y = 0 at the right support x = x

N + 1

= 3L can

be substituted into Equation 18 for k = N–1 to obtain:

(22)

Equation 18 now should have the modiﬁed form of, for k = 3,4,…,N–2:

(23)

In matrix form, Equations 19–23 can be written as [C]{Y} = {R} where {Y} =

• • • y

]

and the elements of the coefﬁcient matrix [C] denoted as c

i,j

are to

be calculated using the formulas:

(24)

For the right-side vector {R} in the matrix equation [C]{Y} = {R}, its elements

denoted as r

are to be calculated using the formula:

(25)

where w(x

) means the distributed load at x = x

. Consider the beam shown in

Figure 12 with dead loads w

and w

. That is:

(26)

Notice that the average value of w

and w

at the junctions of hollow and solid

portions of the beam. Expressions for the EI product can be given similar to those

for w. In a sample application of the QuickBASIC version of the program OdeB-

vpFD, we will give two subprogram functions which are prepared for the beam

problem shown in Figure 12 based on Equations 25 and 26.

An alternative approach for solving the simply supported beam by ﬁnite-differ-

ence approximation of Equation 13 and using the two boundary conditions y = 0 at

yyyyx

at x

NNNN N−−− −

−+−=

()

321

464∆

yyyyy x

at x

kkkkk k−− ++

−+−+=

()

21 12

464 ∆

cc cc

c for i N

c c for i N

NN NN

ii ii

11 12 1

6231

4231

1342

,, , ,

, ,,,

==− = =

==…−

==− =…−

== =…−

−

−+

elsewhere

r r and r x

at x for i N

0231== =

()

=…− , , ,∆

w x w for x L and L x L

w x w for L x L

w x w w for x L and x L

()

=≤<<≤

()

=<<

()

′

023

x = 0 and x = 3L is given as a homework problem. This approach requires the

moment equation M(x) be derived prior to the solution of the matrix equation

[C]{Y} = {R}.

FORTRAN VERSION

Sample Application

Consider the membrane problem shown in Figure 1. Let the inner radius R

equal to 3 in, the stepsize or radial increment R be equal to 9.5 in, and the number

of stations N between the two boundaries be selected as equal to 11 (that is, the

other boundary is at R = R

+ (N + 1) R = 3 + 12x0.5 = R

= 9 in.) If the tension

T is equal to 100 lbs/in and the pressure p is equal to 5 lbs/in

, a function subprogram

RI can then be accordingly prepared as listed in the program OdeBvpFD. An

interactive run of this program has resulted in a display on screen as follows:

QUICKBASIC VERSION

Sample Applications

Consider the same membrane problem as for the sample application using the

FORTRAN version of the program OdeBvpFD. The QuickBASIC version has a

COMMON SHARED statement allowing N to be shared in the subprograms. Also,

this version has been expanded to solving up to 119 stations. It gives an interactive

run as follows:

The simply supported beam under its own weight is considered as a second

example. For w

= w

= –2 and a uniform EI value equal to 1 (that is when the beam

has a uniform cross section without the hollow end portions), L = 100, x = 30,

and N = 9 stations between the supports, the subprograms CIJ and RI are prepared

according to Equations 24 to 26 as follows:

It should be noted that the arguments X and DX for FUNCTION CIJ are not

actually involved in any of the statements therein. They are kept so that the statements

in the program OdeBvpFD involving CIJ can remain as general as possible to

accommodate other problems which may need these arguments as linkages. The

interactive application of the above two FUNCTION subprograms is demonstrated

below.

The analytical solution of this beam problem can be easily obtained. From any

textbook on mechanics of materials (see footnote *), the maximum deﬂection for a

beam of length 3L and subjected to a uniformly distributed load w is y

max

5w(3L)

/384EI. Since w = –2, L = 100, and EI = 1 are used in the above illustrative

run, y

max

is equal to 5x(–2)x(300)

/384 = 2.109x10

. The computed result of y at x

(for N = 9, the 5th station is the mid-length of the beam) by the program OdeBvpFD

is equal to –1.21500x10

which has an error about 42%. A homework problem is

given for readers to exercise different N and x values, which will show that for N

equal to 19, 29, 59, 99, and 119, the respective y

max

values are –1.63285x10

–1.78535x10

, — 1.93665x10

, –1.92011x10

, and –1.83x10

. Notice that the accu-

racy continue to improve until the roundoff errors begin to affect the solution when

N becomes large, for such cases the double precision should be used in solving the

matrix equation [C]{Y} = {R} by the Gaussian elimination method.

MATLAB APPLICATIONS

The two sample problems discussed in FORTRAN and QuickBASIC version

of the program OdeBvpFD can be executed interactively by MATLAB with the

commands entered from keyboard as follows:

Membrane Problem