Просмотреть файл в отдельном окне: 4f31388e9b6bd4e9c18c71c6a99ec679.pdf

6

Ordinary Differential

Equations — Initial and

Boundary Value Problems

6.1 INTRODUCTION

An example of historical interest in solving an unknown function which is governed

by an ordinary differential equation and an initial condition is the case of finding

y(x) from:

dy

=y

dx

y( x = 0 ) = 1

and

(1)

As we all know, y(x) = ex. In fact, the exponential function ex is defined by an

infinite series:

ex = 1 +

x1 x 2

+

+… =

1! 2!

?

?

i=0

xi

i!

(2)

To prove that Equation 2 indeed is the solution for y satisfying Equation 1, here

we apply an iterative procedure of successive integration using a counter k. First,

we integrate both sides of Equation 1 with respect to x:

?

x

0

dy

dx =

dx

x

? y dx

(3)

0

After substituting the initial condition of y(x = 0) = 1, we obtain:

?

x

y(x) = 1 + y dx

(4)

0

So we are expected to find an unknown y(x) which is to be obtained by integrating it? Numerically, we can do it by assuming a y(x) initially (k = 1) equal to

1, investigate how Equation 4 would help us to obtain the next (k = 2), guessed y(x),

and hope eventually the iterative process would lead us to a solution. The iterative

equation, therefore, is for k = 1,2,…

© 2001 by CRC Press LLC

?

x

y( k +1) (x) = 1 + y( k )dx

(5)

0

The results are y(1) = 1, y(2) = 1 + x, y(3) = 1 + x + (x2/2!), and eventually the

final answer is the infinite series given by Equation 2.

What really need to be discussed in this chapter is not to obtain an analytical

expression for y(x) by solving Equation 1 and rather to compute the numerical values

of y(x) when x is equally incremented. That is, for a selected x increment, x (or,

stepsize h), to find yiy(xi) for i = 1,2,… until x reaches a terminating value of xe

and xi = (i–1) x. A simplest method to find y2 is to approximate the derivative of

y(x) by using the forward difference at x1. That is, according to the notation used

in Chapter 4, we can have:

dy

dx

x = x1

=?

?y1 y 2 ? y1

=

= y1

?x

?x

(6)

Or, y2 = (1 + x)y1. In fact, this result can be extended to any xi to obtain yi + 1 =

1 + (x)yi. For the general case when the right-hand side of Equation 1 is equal to a

prescribed function f(x), we arrive at the Euler’s formula yi + 1 = yi + f(xi) x. Euler’s

formula is easy to apply but is inaccurate. For example, using a x = 0.1 in solving

Equation 1 with the Euler’s formula, it leads to y2 = (1 + 0.1)y1 = 1.1, y3 = (1 +

0.1)y2 = 1.21 when the exact values are y2 = e0.1 = 1.1052 and y3 = e0.2 = 1.2214,

respectively. The computational errors accumulate very rapidly.

In this chapter, we shall introduce the most commonly adopted method of RungeKutta for solution of the initial-value problems governed by ordinary differential

equation(s). For the fourth-order Runge-Kutta method, the error per each computational step is of order h5 where h is the stepsize. Converting the higher-order

ordinary differential equation(s) into the standardized form using the state variables

will be illustrated and computer programs will be developed for numerical solution

of the problem.

Engineering problems which are governed by ordinary differential equations and

also some associated conditions at certain boundaries will be also be discussed.

Numerical methods of solution based on the Runge-Kutta procedure and the finitedifference approximation will both be explained.

6.2 PROGRAM RUNGEKUT — APPLICATION

OF THE RUNGE-KUTTA METHOD

FOR SOLVING THE INITIAL-VALUE PROBLEMS

Program RungeKut is designed for solving the initial-value problems governed by

ordinary differential equations using the fourth-order Runge-Kutta method. There

are numerous physical problems which are mathematically governed by a set of

ordinary differential equations (ODE) involving many unknown functions. These

unknown functions are all dependent of a variable t. Supplementing to this set of

© 2001 by CRC Press LLC

FIGURE 1. The often cited vibration problem shown requires the changes of the elevation

x and velocity v to be calculated.

ordinary differential equations are the initial conditions of the dependent functions

when t is equal to zero. For example, the often cited vibration1 problem shown in

Figure 1 requires the changes of the elevation x and velocity v to be calculated using

the equations:

(

)

m d 2 x dt 2 + c(dx dt ) + kx = f (t )

(1)

dx dt = v

(2)

and

where m is the mass, c is the damping coefficient, k is the spring constant, t is the

time, and f(t) is a disturbing force applied to the mass. When the physical parameters m, c, and k, and the history of the applied force f(t) are specified, the complete

histories of the mass’ elevation x and velocity v can be calculated analytically, or,

numerically if the initial elevation x(t = 0) and v(t = 0) are known. If m, c, and k

remain unchanged throughout the period of investigation and f(t) is a commonly

encountered function, Equation 1 can be solved analytically.1 Otherwise, a numerical

method has to be applied to obtain approximate solution of Equation 1.

Many numerical methods are available for solving such initial-value problems

governed by ordinary differential equations. Most of the numerical methods require

that the governing differential equation be rearranged into a standard form of:

© 2001 by CRC Press LLC

dx1 dt = F1 (x1, x 2 ,…, x n , t; parameters)

dx 2 dt = F2 (x1, x 2 ,…, x n , t; parameters)

(3)

. . . . . . . . . .

dx n dt = Fn (x1, x 2 ,…, x n , t; parameters)

For example, the variables x and v in Equations 1 and 2 are to be renamed x1

and x2, respectively. Equation 1 is to be rewritten as:

m(dv dt ) + cv + kx = f (t )

and then as:

dv dt = [f (t ) ? cv ? kx] m

and finally as:

dx 2 dt = F2 (x1, x 2 , t; m, c, k )

Meanwhile. Equation 2 is rewritten as:

dx1 dt = F1 (x1, x 2 , t; m, c, k )

Or, more systematically the problem is described by the equations:

dx1 dt = F1 (x1, x 2 , t; m, c, k ) = x 2

dx 2 dt = F2 (x1, x 2 , t; m, c, k ) = [f (t ) ? kx1 ? cx 2 ] m

(4)

and having the initial conditions x1(t = 0) and x2(t = 0) prescribed.

Runge-Kutta method is a commonly used method for numerical solution of a

system of first-order ordinary differential equations expressed in the general form

of (3). It is to be introduced and illustrated with a number of practical applications.

RUNGE-KUTTA METHOD (FOURTH-ORDER)

Consider the problem of finding x and y values at t>0 when they are governed

by the equations:

(dx dt ) ? 4x + (dy dt ) ? 6y = ?2

© 2001 by CRC Press LLC

(5)

and

2(dx dt ) + x + 3(dy dt ) + 2 y = 0

(6)

when initially their values are x(t = 0) = 7 and y(t = 0) = –4. The analytical solutions

are obtainable2 and they are:

x = 5e t + 2

and

y = ?3e t ? 1

(7)

To solve the problem numerically, Equations 5 and 6 need to be decoupled and

expressed in the form of Equation 3. Cramer’s rule can be applied by treating dx/dt

and dy/dt as two unknowns and x and y as parameters, the converted standard form

after changing x to x1 and y to x2 is:

dx1 dt = F1 (x1, x 2 , t; constants), x1 (0) = 7

(8)

dx 2 dt = F2 (x1, x 2 , t; constants), x 2 (0) = ?4

(9)

F1 (x1, x 2 , t;constants) = ?6 + 13x1 + 20 x 2

(10)

F2 (x1, x 2 , t;constants) = 4 ? 9x1 ? 14 x 2

(11)

where:

and

Numerical solution of x1 and x2 for t>0 is to use a selected time increment t

(often referred to as the stepsize h for the independent variable t). Denote t0 as the

initial instant t = 0 and tj + 1 as the instant after j increments of time, that is, tj + 1 =

(j + 1)h. If the values for x1 and x2 at tj, denoted as x1,j and x2,j respectively, are

already known, the fourth-order Runge-Kutta method is to use the following formulas to calculate x1 and x2 at tj + 1, denoted as x1,j + 1 and x2,j + 1:

(

)

x i, j+1 = x i, j + p i,1 + 2 p i,2 + 2 p i,3 + p i,4 6

(12)

for i = 1,2. The p’s in Equation 12 are the Runge-Kutta parameters to be calculated

using the functions F1 and F2 by adjusting the values of the variables x1 and x2 at

tj. The formulas for calculating these p’s are, for i = 1,2

(

p i,1 = hFi t j , x1, j , x 2, j

© 2001 by CRC Press LLC

)

(13)

[

(

)

(

)]

(14)

[

(

)

(

)]

(15)

p i,2 = hFi t j + ( h 2), x1, j + p1,1 2 , x 2, j + p2,1 2

p i,3 = hFi t j + ( h 2), x1, j + p1,2 2 , x 2, j + p2,2 2

(

p i,4 = hFi t j + h, x1, j + p1,3 , x 2, j + p2,3

)

(16)

Equations 12 to 16 are to be used to generate x1 and x2 values at tj for j = 1,2,3,…

which can be tabulated as:

t

x1

x2

0

7

–4

h

?

?

2h

?

?

3h

?

?

.

.

.

.

.

.

.

.

.

jh

x1,j

x2,j

.

.

.

.

.

.

.

.

.

te

?

?

where te is the ending value of t at which the computation is to be terminated. The

first pair of values to be filled into the above table is for x1 and x2 at t = h (j = 1).

Based first on Equations 13 to 16 and then Equation 12, the actual computations for

h = 0.1 and at t1 go as follows:

p1,1 = hF1(t0,x1,0,x2,0) = 0.1F1(0,7,–4) = 0.1(–6 + 91-80) = 0.5

p2,1 = hF2(t0,x1,0,x2,0) = 0.1F2(0,7,–4) = 0.1(4-63 + 56) = –0.3

p1,2 = hF1(t0 + 0.05,x1,0 + 0.25,x2,0–0.15) = 0.1F1(.05,7.25,–4.15)

= 0.1(–6 + 13x7.25-20x4.15) = 0.525

p2,2 = hF2(t0 + 0.05,x1,0 + 0.25,x2,0–0.15) = 0.1F2(.05,7.25,–4.15)

= 0.1(4-9x7.25 + 14x4.15) = –0.315

p1,3 = hF1(t0 + 0.05,x1,0 + 0.2625,x2,0–0.1575)

= .1F1(.05,7.2625,–4.1575)

= 0.1(–6 + 13x7.2625-20x4.1575) = .52625

p2,3 = hF2(t0 + 0.05,x1,0 + 0.2625,x2,0–0.1575)

= .1F2(.05,7.2625,–4.1575) = 0.1(4–9x7.2625 + 14x4.1575)

= –0.31575

p1,4 = hF1(t0 + 0.1,x1,0 + 0.52625,x2,0–0.31575)

= 0.1F1(0.1,7.52625,–4.31575)

= 0.1(–6 + 13x7.52625–20x4.31575) = 0.552625

p2,4 = hF2(t0 + 0.1,x1,0 + 0.52625,x2,0–0.31575)

= .1F2(.1,7.52625,–4.31575) = .1(4–9x7.52625 + 14x4.31575) = –0.331575

x1,1 = 7 + (0.5 + 2x0.525 + 2x0.52625 + 0.552625)/6

= 7 + (3.155125)/6 = 7.5258541

(17)

x2,1 = –4 + (–0.3–2x0.315–2x0.31575–0.331575)/6

= –4 + (–1.893075)/6 = –4.3155125

(18)

The exact solution calculated by using Equation 7 are:

x1,1 = 7.5258355

© 2001 by CRC Press LLC

and

x2,1 = –4.3155013

(19)

The errors are 0.000247% and 0.000260% for x1 and x2, respectively. Per-step

error for the fourth-order Runge-Kutta method is difficult to estimate because the

method is derived by matching terms in Equation 12 with Taylor-series expansions

of x1 and x2 about ti through and including the h4 terms. But approximately, the perstep error is of order h5. For better accuracy, the fifth-order Runge-Kutta method

should be applied. For general use, the classic fourth-order Runge-Kutta method is,

however, easier to develop a computer program which is to be discussed next.

SUBROUTINE RKN

A subroutine called RKN has been written for applying the fourth-order RungeKutta method to solve the initial-value problems governed by a set of first-order

ordinary differential equations. It has been coded according to the procedure

described in the preceding section. That is, the equations must be in the form of

Equation 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 seven

arguments for this subroutine, the first four are input arguments where the last is an

output argument. The fifth argument P keeps the Runge-Kutta parameters generated

in this subroutine. The sixth argument XT is needed for adjusting the input argument

XIN. These two arguments, P and XT, are included for handling the general case

of N variables. Listing them as arguments makes possible to specify them as matrix

and vector of adjustable sizes.

FORTRAN VERSION

© 2001 by CRC Press LLC

QUICKBASIC VERSION

PROGRAM RUNGEKUT

Program RungeKut which calls the subroutine RKN is to be run interactively

by specifying the inputs through the keyboard. Displayed messages on screen instruct

user 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 following

inputs 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 only

9 functions can be handled by the program RungeKut, but the user should

be able to expand the program to accommodate any N value which is

greater 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 of

the independent variable, denoted as t0 and te (T0 and TEND in the

program RungeKut), respectively. It is over this range, the values of the

N dependent variables are to be calculated.

(6) Enter the stepsize, h (DT in the program RungeKut), with which the

independent variable is to be incremented.

(7) Enter the N initial values.

© 2001 by CRC Press LLC

QUICKBASIC VERSION

FORTRAN VERSION

© 2001 by CRC Press LLC

FUNCTION F

According to Equations 12 to 15, the Runge-Kutta parameters pi,j (matrix P in

the program RungeKut) are calculated using two FOR-NEXT loops — an I loop

covering 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 defined

by Equations 8 to 11. We create a supporting function program F as follows:

QuickBASIC Version

© 2001 by CRC Press LLC

FORTRAN Version

The 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, and

the 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 by

Equations 8 to 11 is manifested by a copy of the screen display:

© 2001 by CRC Press LLC

SAMPLE APPLICATIONS

OF THE PROGRAM

RUNGEKUT

As a first example, consider the problem of a beam shown in Figure 2 which is

built into the wall so that it is not allowed to displace or rotate at the left end. For

consideration of the general case when the beam is loaded by (1) a uniformly

distributed 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 the

free end of the beam, x = 40 cm. It is of concern to the structural engineers to know

how the beam will be deformed. The equation for finding the deflected curve of the

beam, 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 along

the beam. It can be shown that the moment distribution for the loading shown in

Figure 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 LLC

for

for

for

for

0 x 10 cm

10 x 20 cm

20 < x 30 cm

30 < x 40 cm

(22)

FIGURE 2. The problem of a beam, which is built into the wall so that it is not allowed to

displace or rotate at the left end.

To convert the problem into the form of Equation 3, we replace y, dy/dx, and x

by x1, x2, and t, respectively. Knowing Equation 21 and that the initial conditions

are y = 0 and dy/dx = 0 at x = 0, we can obtain from Equation 21 the following

system of first-order ordinary differential equations:

dx1 dx = F1 (x, x1, x 2 ) = x 2 ,

x1 (x = 0) = 0

(23)

and

dx 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 Equation

21 whereas the beam stiffness EI is to be set equal to 2x105 N/cm2 in numerical

calculation of the deflection of the beam using program RungeKut.

To apply program RungeKut for solving the deflection equation, y(x) which has

been renamed as x1(t), we need to define in the QuickBASIC program RungeKut

two 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 QuickBASIC

programming. In view of Equation 21, F2 in Equation 24 has to be defined by

modifying the subprogram function, here we illustrate it with a FORTRAN version:

© 2001 by CRC Press LLC

FUNCTION BM(X)

IF ((X.LE.0.).AND.(X.LT.10.)) BM = –1121./3 + 24*X-X**2/2

IF ((X.GE.10.).AND.(X.LT.20.)) BM = –871./3 + 4*T +

T**2T**3/30

IF ((X.GE.20.).AND.(X.LT.30.)) BM = 4*X–157

IF ((X.GE.30.).AND.(X.LT.40.)) BM = 4*X–160

RETURN

END

In fact, each problem will require such arrangement because these function

statements and subprogram function describe the particular features of the problem

being 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 errors

in % of the analytical solution ( = –0.68917) are tabulated below:

Stepsize h, cm

y at x = 40 cm,

5

2

1

.5

.25

.1

.05

Error, %

–0.79961

–0.73499

–0.71236

–0.70083

–0.69501

–0.69151

–0.69033

16.0

6.65

3.36

1.69

.847

.340

.168

This 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 loads

directly. This approach saves the reader from deriving the expressions for the bending

moment. Readers interested in such an approach should solve Problem 4.

A NONLINEAR OSCILLATION PROBLEM SOLVED

BY

RUNGEKUT

The numerical solution using the Runge-Kutta method can be further demonstrated by solving a nonlinear problem of two connected masses m1 and m2 as shown

in 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, their

instantaneous 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 in

term of y as:

[

? = tan ?1 (y + h ) b

]

and

[

L = (y + h) + b2

2

]

0.5

(25,26)

If we denote the cable tension as T, then the Newton’s second law applied to

the two masses leads to m1g–2Tsin? = m1d2y/dt2 and T-m2g = m2d2z/dt2. By eliminating T, we obtain:

© 2001 by CRC Press LLC

FIGURE 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.

? 2m2

?

d 2y 2m2

d 2z

+

sin ? 2 = g?1 ?

sin ??

2

dt

m1

dt

m

?

?

1

(27)

The displacements y and z are restricted by the condition that the cable length

must remain unchanged. That is z = 2{[(y + h)2 + b2]1/2[h2 + b2]1/2}. This relationship

can be differentiated with respect to t to obtain another equation relating d2y/dt2 and

d2z/dt2 which is:

2

d 2 z ?2

dLdy 2 ?? dy ?

d2y ?

= 2 (y + h)

+ ?

+ (y + h) 2 ?

2

dt ?

dt

L

dtdt L ?? dt ?

(28)

Equation 28 can be substituted into Equation 27 to obtain:

?? 2 2 m 2

d2y

sin ?

=

?gL ?

2

dt

L(2 y + h + L ) ??

m1

2

? 2

dLdy

? dy ? ? ??

gL

y

h

L

?

+

+

2

2

(

)

?

? dt ? ?? ??

dtdt

?

?

(29)

By letting x1 = y, x2 = dy/dt, x3 = z, and x4 = dz/dt, then according to Equation

3 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. It can be shown

that the final form of the system of four first-order differential equations are:

dx1

= x2 ,

dt

dx3

= x4 ,

dt

© 2001 by CRC Press LLC

2

2

dx 2 gL (1 ? 2 rm sin ?) + 4 rm v sin ?

,

=

2

dt

L 1 + 4 r sin ?

(

dx 4 ?2 v

dx

=

+ 2 2 sin ?

dt

L

dt

2

and

)

(30)

FIGURE 4. A numerical case where b = 3.6 m, h = .15 m, m2/m1 = 0.8, and initial conditions

y = z = dy/dt = dz/dt = 0 has been investigated and the results for -y, z, -dy/dt, and dz/dt have

been plotted for 0?t?12 seconds.

where:

rm =

m2

m1

and

dL

v2 = 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 conditions

y = z = dy/dt = dz/dt = 0 has been investigated and the results for -y, z, -dy/dt, and

dz/dt have been plotted in Figure 4 for 0?t?12 seconds. The reason that negative

values 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 varies

between 0 and 3.8358 m. The oscillation has a period approximately equal to

2x6.6831 = 13.3662 seconds, so the frequency is about 0.47 rad/sec. The masses

reach 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 the

oscillation for the studied period are listed below:

© 2001 by CRC Press LLC

MATLAB APPLICATION

MATLAB has a file called ode45.m which implement the fourth- and fifthorder Runge-Kutta integration. Here, we demonstrate how the sample problem used

in the FORTRAN and QuickBASIC versions can also be solved by use of the m

file. The forcing functions given in Equations 10 and 11 are first prepared as follows:

© 2001 by CRC Press LLC

If 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 of

investigation 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 which

the forcing functions are defined. The second and third argument the initial and final

values of time, respectively. The fourth argument is a vector containing the initial

values of the dependent variables. The resulting display, after rearranging T and X

side-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 created

to 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), and

spring 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 LLC

FIGURE 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 by

entering 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 2

in the program NewRaphG) of about 0.28 cm at approximately t = 2.5 seconds and

finally settles to a static deflection of 0.25 cm, and that the maximum ascending

velocity is about 0.18 cm/sec.

As another example of dynamic analysis in the field of fluid mechanics, Figure 8

shows the flow of a fluid between two connected tanks. The valve settings control

the amount of flows, qi for i = 1,2,3. The levels of the tanks h1 and h2 change in

time depending on these settings and also on the supply rates Q1 and Q2 and the

discharge 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 the

formulas:

A1

dh1

1

= Q1 ?

(h ? h3 )

dt

R1 1

© 2001 by CRC Press LLC

and

A2

dh 2

1

= Q2 ?

(h ? h3 )

dt

R2 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 commands

shown below.

where A’s are the cross-sectional areas of the tanks, and h3 is the pressure head at

the junction indicated in Figure 7 and is related to the discharge rate q3 by the

equation:

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 differential

equations in h1 and h2 to be:

© 2001 by CRC Press LLC

FIGURE 7. The flow of a fluid between two connected tanks.

A1

dh1

= Q1 ? a1h1 + a 3h 2

dt

and

A2

dh 2

= Q 2 ? a 2 h 2 + a 3h1

dt

a1 = ( R 2 + R 3 ) ?

and

a 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 both

tanks.4

MATHEMATICA APPLICATIONS

Mathematica 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 LLC

Out[1] = {{X1 -> InterpolatingFunction[{0., 25.}, ],

X2 -> InterpolatingFunction[{0., 25.}, ]}}

In[1] shows that NDSolve has three arguments: the first argument defines the

two ordinary differential equations and the initial conditions, the second argument

lists the dependent variables, and the third argument specifies the independent

variable and the range of investigation. Id is a name selected for easy later reference

of the results generated. Out[1] indicates that interpolation functions have been

generated for X1 and X2 for t in the range from 0 to 25. To print the values of X1

and 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 = 0

t=

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 LLC

X1 = {0.}

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

X1 =

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 MATLAB

application. In[2] shows that the replacement operator/. is employed in X1[tv]/.Id

which requires all t appearing in the resulting interpolating function for X1(t) created

in the Id statement to be substituted by the value of tv. The looping DO statement

instructs that the tv values be changed from a minimum of 0 and a maximum of 25

using an increment of 1. To have a closer look of the overshoot region, In[2] can

be modified to yield

In[3]: = (Do[Print[“t = “,tv,” X1 = “,X1[tv]/. Id,” X2 = “,X2[tv]/. Id],

{tv, 2, 3, 0.1}])

Out[3] =

t=2

t=

t=

t=

t=

t=

t=

t=

t=

t=

t=

X1 = {0.267432}

2.1

2.2

2.3

2.4

2.5

2.6

2.7

2.8

2.9

3.

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 overshoot

occurs at t = 2.6 and X1 has a maximum value equal to 0.28502.

6.3 PROGRAM ODEBVPRK — APPLICATION OF RUNGE-KUTTA

METHOD FOR SOLVING BOUNDARY-VALUE PROBLEMS

The program OdeBvpRK is designed for numerically solving the linear boundaryvalue problems governed by the ordinary differential equation by superposition of

two solutions obtained by application of the Runge-Kutta fourth-order method. To

explain the procedure involved, consider the problem of a loaded beam shown in

Figure 8. Mathematically, the deflection y(x) satisfies the well-known flexural equation.5

d2y M

=

dx 2 EI

(1)

where M is the internal moment distribution, E is the Young’s modulus and I is the

moment of inertia of the cross section of the beam. For the general case, M, E and

I can be function of x. The boundary conditions of this problem are:

© 2001 by CRC Press LLC

FIGURE 8. The problem of a loaded beam.

y(x = 0) = 0

and

y(x = L ) = 0

(2)

where L is the length of the beam. The problem is to determine the resulting

deflection y(x). Knowing y, the moment and shearing force, V, distributions can

subsequently be determined based on Equation 1 and V = dM/dx. The final objective

is 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 if

in additional to the initial condition given in (2), y(x = 0) = 0, we also know the

slope, , at x = 0. But, we can always make a guess and hope that by making better

and better guesses the trial process will eventually lead to one which satisfies the

other boundary condition, namely y(x = L) = 0 given in Equation 3. In fact, if the

problem is linear, all we need to do is making two guesses and linearly combine

these two trial solutions to obtain the solution y(x).

Let us first convert the governing differential Equation 1 into two first-order

equations as:

? dx1

? dx = x 2 ,

? dx 2 M

= ,

?

? dx EI

x1 (x = 0) = y(x = 0) = 0

(3,4)

x 2 (x = 0) = ?(x = 0) = ?0

To apply the fourth-order Runge-Kutta method, we have to first decide on a

stepsize h, for example h = L/N we then plan to calculate the deflections at N + 1

locations, x + jh for j = 1,2,…N since j = 0 is the initial location. If we assume a

value for 0, say A, the Runge-Kutta process will be able to generate the following

table

x

x 1=y (1)

0

0

h

y1(1)

2h

y2(1)

…

…

jh

yj(1)

…

…

Nh

yN(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 not

equal to zero, then the value incorrectly selected, we have to make a second try by

letting (x = 0) = B to obtain a second table by application of the Runge-Kutta

method. Let the second table be denoted as:

x

x 1=y (2)

0

0

h

y1(2)

2h

y2(2)

…

…

jh

yj(2)

…

…

Nh

yN(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, if

the problem is linear, the solutions can be obtained by linearly combining the two

trial 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:

? +? =1

and

?y(N1) + ?y(N2 ) = 0

(9,10)

Equation 10 is derived from the boundary condition y(x = Nh = L) = 0 and based

on Equation 7. Equation 9 needs more explanation because it cannot be derived if

y(x = 0) = 0. Let us assume that for the general case, y(x = 0) = . Then Equation

7 gives + = which leads to Equation 9. Using Cramer’s rule, we can easily

obtain:

? = y(N2 ) D

and

? = ? y(N1) D

(11,12)

and

(1)

D = y(2

N ? yN

(13)

NUMERICAL EXAMPLES

Let us consider the problem of a loaded beam as shown in Figure 8. The

crosssection of the beam has a width of 1 cm and a height of 2 cm which results in

a moment of inertia I = 2/3 cm4. The reactions at the left and right supports can be

computed to be 5/3 N and 25/3 N, respectively. Based on these data, it can be shown

that the equations for the internal bending moments are:

for 0 ? x ? 20 cm

?5x 3,

?

M( x ) = ?

65

1 2

for 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 decide

on a stepsize of h = 1 cm and proceed to prepare a computer program using the

fourth-order Runge-Kutta method to generate two trial solutions and then linearly

combining 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 presented

later has produced the following display on screen:

© 2001 by CRC Press LLC

FORTRAN VERSION

© 2001 by CRC Press LLC

The Subprogram FUNCTION F which defines the initial-value problem is coded

in accordance with Equation 14. The two trial initial slopes are selected as equal to

0.1 and 0.2. The trial results are kept in the three-dimensional variable XT, in which

the deflection y(k)(j) for the kth try at station x = xj = jh is stored in XT(1,j,k) whereas

the 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 because

XT(1,KS,NTRY) is transmitted as XIN(1) and automatically the next value

XT(2,KS, NTRY) as XIN(2), and the computed results XOUT(1) and XOUT(2) are

to be stored as XT(1,KS + 1,NTRY) and XT(2,KS + 1,NTRY), respectively. Notice

that there are only two dependent variables, NV = 2.

After the weighting coefficients (ALPHA in the program) and (BETA) have

been calculated, the final distributions of the deflection and slope are saved in first

and second rows of the two-dimensional variable X, respectively. It should be

emphatically noted that the solutions obtained is only good for the selected stepsize

h = 1 cm. Whether it is accurate or not remains to be tested by using finer stepsizes

and 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. The

needed modification of the above listed program to include this change in the stepsize

and testing of the difference in the maximum deflection is left as homework for the

reader.

© 2001 by CRC Press LLC

QUICKBASIC VERSION

MATLAB APPLICATIONS

In the program RungeKut, MATLAB is used for solving initial value problems

by 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-value

problem governed by linear ordinary differential equations. To demonstrate the

procedure, the sample problem discussed in FORTRAN and QuickBASIC versions

of the program OdeBvpRK, with the aid of function BVPF.m listed in the subdirectory , can be solved by interactive MATLAB operations as follows:

© 2001 by CRC Press LLC

Notice that format compact enables the display to use fewer spacings; YT1 and

YT2 keep the two trial solutions, and ode45 automatically determines the best

stepsize which if used directly will result in a coarse plot as shown in Figure 9. The

plot showing solid-line curve for the deflection and broken-line curve for the slope

can, 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 is

to be applied as follows:

© 2001 by CRC Press LLC

FIGURE 9.

The result of plotting the spline curves is shown in Figure 10.

MATHEMATICA APPLICATIONS

The 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; readers

are therefore referred to program RungeKut to review the method itself and the

development of related programs and subprograms. The boundary-value problem is

to be solved by continuously guessing the initial condition(s) which are not provided

until all boundary conditions are satisfied if the problem is nonlinear. When the

problem is linear, then only a finite number of guesses are necessary. A system of

two first-order ordinary differential equations which governed the loaded elastic

beam problem previously solved in the MATLAB application is here adapted to

demonstrate the application of the Runge-Kutta method.

© 2001 by CRC Press LLC

FIGURE 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]/. Id1

Out[3] = {3.00052}

EI value and F(X) are defined in In[1]. In[2] specifies the two first-order ordinary

differential equations involving the deflection Y and slope YP, describes the correct

initial condition Y(X = 0) = 0, gives a guessed slope Y (X = 0) = 0.1, and decides

on the limit of investigation from X = 0 to X = 30. In[3] interpolates the ending Y

value 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 LLC

In[5]: = Y30Trial2 = Y[30]/. Id2

Out[5] = {6.00052}

Linear combination of the two trial solutions is now possible by calculating a

correct 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*b

Out[8] = {–0.0000174888}

Finally, the actual deflection and slope can be obtained by providing the correct

set 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 LLC

0

1

2

3

4

5

6

7

8

9

10

11

12

13

14

Y=

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=

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

Y=

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 — 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

finite-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 = Ri and R = Ro, respectively.

As shown in Figure 11, it is then inflated 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 equation6

Z(R) satisfies Equation 1 is for Ri

Ordinary Differential

Equations — Initial and

Boundary Value Problems

6.1 INTRODUCTION

An example of historical interest in solving an unknown function which is governed

by an ordinary differential equation and an initial condition is the case of finding

y(x) from:

dy

=y

dx

y( x = 0 ) = 1

and

(1)

As we all know, y(x) = ex. In fact, the exponential function ex is defined by an

infinite series:

ex = 1 +

x1 x 2

+

+… =

1! 2!

?

?

i=0

xi

i!

(2)

To prove that Equation 2 indeed is the solution for y satisfying Equation 1, here

we apply an iterative procedure of successive integration using a counter k. First,

we integrate both sides of Equation 1 with respect to x:

?

x

0

dy

dx =

dx

x

? y dx

(3)

0

After substituting the initial condition of y(x = 0) = 1, we obtain:

?

x

y(x) = 1 + y dx

(4)

0

So we are expected to find an unknown y(x) which is to be obtained by integrating it? Numerically, we can do it by assuming a y(x) initially (k = 1) equal to

1, investigate how Equation 4 would help us to obtain the next (k = 2), guessed y(x),

and hope eventually the iterative process would lead us to a solution. The iterative

equation, therefore, is for k = 1,2,…

© 2001 by CRC Press LLC

?

x

y( k +1) (x) = 1 + y( k )dx

(5)

0

The results are y(1) = 1, y(2) = 1 + x, y(3) = 1 + x + (x2/2!), and eventually the

final answer is the infinite series given by Equation 2.

What really need to be discussed in this chapter is not to obtain an analytical

expression for y(x) by solving Equation 1 and rather to compute the numerical values

of y(x) when x is equally incremented. That is, for a selected x increment, x (or,

stepsize h), to find yiy(xi) for i = 1,2,… until x reaches a terminating value of xe

and xi = (i–1) x. A simplest method to find y2 is to approximate the derivative of

y(x) by using the forward difference at x1. That is, according to the notation used

in Chapter 4, we can have:

dy

dx

x = x1

=?

?y1 y 2 ? y1

=

= y1

?x

?x

(6)

Or, y2 = (1 + x)y1. In fact, this result can be extended to any xi to obtain yi + 1 =

1 + (x)yi. For the general case when the right-hand side of Equation 1 is equal to a

prescribed function f(x), we arrive at the Euler’s formula yi + 1 = yi + f(xi) x. Euler’s

formula is easy to apply but is inaccurate. For example, using a x = 0.1 in solving

Equation 1 with the Euler’s formula, it leads to y2 = (1 + 0.1)y1 = 1.1, y3 = (1 +

0.1)y2 = 1.21 when the exact values are y2 = e0.1 = 1.1052 and y3 = e0.2 = 1.2214,

respectively. The computational errors accumulate very rapidly.

In this chapter, we shall introduce the most commonly adopted method of RungeKutta for solution of the initial-value problems governed by ordinary differential

equation(s). For the fourth-order Runge-Kutta method, the error per each computational step is of order h5 where h is the stepsize. Converting the higher-order

ordinary differential equation(s) into the standardized form using the state variables

will be illustrated and computer programs will be developed for numerical solution

of the problem.

Engineering problems which are governed by ordinary differential equations and

also some associated conditions at certain boundaries will be also be discussed.

Numerical methods of solution based on the Runge-Kutta procedure and the finitedifference approximation will both be explained.

6.2 PROGRAM RUNGEKUT — APPLICATION

OF THE RUNGE-KUTTA METHOD

FOR SOLVING THE INITIAL-VALUE PROBLEMS

Program RungeKut is designed for solving the initial-value problems governed by

ordinary differential equations using the fourth-order Runge-Kutta method. There

are numerous physical problems which are mathematically governed by a set of

ordinary differential equations (ODE) involving many unknown functions. These

unknown functions are all dependent of a variable t. Supplementing to this set of

© 2001 by CRC Press LLC

FIGURE 1. The often cited vibration problem shown requires the changes of the elevation

x and velocity v to be calculated.

ordinary differential equations are the initial conditions of the dependent functions

when t is equal to zero. For example, the often cited vibration1 problem shown in

Figure 1 requires the changes of the elevation x and velocity v to be calculated using

the equations:

(

)

m d 2 x dt 2 + c(dx dt ) + kx = f (t )

(1)

dx dt = v

(2)

and

where m is the mass, c is the damping coefficient, k is the spring constant, t is the

time, and f(t) is a disturbing force applied to the mass. When the physical parameters m, c, and k, and the history of the applied force f(t) are specified, the complete

histories of the mass’ elevation x and velocity v can be calculated analytically, or,

numerically if the initial elevation x(t = 0) and v(t = 0) are known. If m, c, and k

remain unchanged throughout the period of investigation and f(t) is a commonly

encountered function, Equation 1 can be solved analytically.1 Otherwise, a numerical

method has to be applied to obtain approximate solution of Equation 1.

Many numerical methods are available for solving such initial-value problems

governed by ordinary differential equations. Most of the numerical methods require

that the governing differential equation be rearranged into a standard form of:

© 2001 by CRC Press LLC

dx1 dt = F1 (x1, x 2 ,…, x n , t; parameters)

dx 2 dt = F2 (x1, x 2 ,…, x n , t; parameters)

(3)

. . . . . . . . . .

dx n dt = Fn (x1, x 2 ,…, x n , t; parameters)

For example, the variables x and v in Equations 1 and 2 are to be renamed x1

and x2, respectively. Equation 1 is to be rewritten as:

m(dv dt ) + cv + kx = f (t )

and then as:

dv dt = [f (t ) ? cv ? kx] m

and finally as:

dx 2 dt = F2 (x1, x 2 , t; m, c, k )

Meanwhile. Equation 2 is rewritten as:

dx1 dt = F1 (x1, x 2 , t; m, c, k )

Or, more systematically the problem is described by the equations:

dx1 dt = F1 (x1, x 2 , t; m, c, k ) = x 2

dx 2 dt = F2 (x1, x 2 , t; m, c, k ) = [f (t ) ? kx1 ? cx 2 ] m

(4)

and having the initial conditions x1(t = 0) and x2(t = 0) prescribed.

Runge-Kutta method is a commonly used method for numerical solution of a

system of first-order ordinary differential equations expressed in the general form

of (3). It is to be introduced and illustrated with a number of practical applications.

RUNGE-KUTTA METHOD (FOURTH-ORDER)

Consider the problem of finding x and y values at t>0 when they are governed

by the equations:

(dx dt ) ? 4x + (dy dt ) ? 6y = ?2

© 2001 by CRC Press LLC

(5)

and

2(dx dt ) + x + 3(dy dt ) + 2 y = 0

(6)

when initially their values are x(t = 0) = 7 and y(t = 0) = –4. The analytical solutions

are obtainable2 and they are:

x = 5e t + 2

and

y = ?3e t ? 1

(7)

To solve the problem numerically, Equations 5 and 6 need to be decoupled and

expressed in the form of Equation 3. Cramer’s rule can be applied by treating dx/dt

and dy/dt as two unknowns and x and y as parameters, the converted standard form

after changing x to x1 and y to x2 is:

dx1 dt = F1 (x1, x 2 , t; constants), x1 (0) = 7

(8)

dx 2 dt = F2 (x1, x 2 , t; constants), x 2 (0) = ?4

(9)

F1 (x1, x 2 , t;constants) = ?6 + 13x1 + 20 x 2

(10)

F2 (x1, x 2 , t;constants) = 4 ? 9x1 ? 14 x 2

(11)

where:

and

Numerical solution of x1 and x2 for t>0 is to use a selected time increment t

(often referred to as the stepsize h for the independent variable t). Denote t0 as the

initial instant t = 0 and tj + 1 as the instant after j increments of time, that is, tj + 1 =

(j + 1)h. If the values for x1 and x2 at tj, denoted as x1,j and x2,j respectively, are

already known, the fourth-order Runge-Kutta method is to use the following formulas to calculate x1 and x2 at tj + 1, denoted as x1,j + 1 and x2,j + 1:

(

)

x i, j+1 = x i, j + p i,1 + 2 p i,2 + 2 p i,3 + p i,4 6

(12)

for i = 1,2. The p’s in Equation 12 are the Runge-Kutta parameters to be calculated

using the functions F1 and F2 by adjusting the values of the variables x1 and x2 at

tj. The formulas for calculating these p’s are, for i = 1,2

(

p i,1 = hFi t j , x1, j , x 2, j

© 2001 by CRC Press LLC

)

(13)

[

(

)

(

)]

(14)

[

(

)

(

)]

(15)

p i,2 = hFi t j + ( h 2), x1, j + p1,1 2 , x 2, j + p2,1 2

p i,3 = hFi t j + ( h 2), x1, j + p1,2 2 , x 2, j + p2,2 2

(

p i,4 = hFi t j + h, x1, j + p1,3 , x 2, j + p2,3

)

(16)

Equations 12 to 16 are to be used to generate x1 and x2 values at tj for j = 1,2,3,…

which can be tabulated as:

t

x1

x2

0

7

–4

h

?

?

2h

?

?

3h

?

?

.

.

.

.

.

.

.

.

.

jh

x1,j

x2,j

.

.

.

.

.

.

.

.

.

te

?

?

where te is the ending value of t at which the computation is to be terminated. The

first pair of values to be filled into the above table is for x1 and x2 at t = h (j = 1).

Based first on Equations 13 to 16 and then Equation 12, the actual computations for

h = 0.1 and at t1 go as follows:

p1,1 = hF1(t0,x1,0,x2,0) = 0.1F1(0,7,–4) = 0.1(–6 + 91-80) = 0.5

p2,1 = hF2(t0,x1,0,x2,0) = 0.1F2(0,7,–4) = 0.1(4-63 + 56) = –0.3

p1,2 = hF1(t0 + 0.05,x1,0 + 0.25,x2,0–0.15) = 0.1F1(.05,7.25,–4.15)

= 0.1(–6 + 13x7.25-20x4.15) = 0.525

p2,2 = hF2(t0 + 0.05,x1,0 + 0.25,x2,0–0.15) = 0.1F2(.05,7.25,–4.15)

= 0.1(4-9x7.25 + 14x4.15) = –0.315

p1,3 = hF1(t0 + 0.05,x1,0 + 0.2625,x2,0–0.1575)

= .1F1(.05,7.2625,–4.1575)

= 0.1(–6 + 13x7.2625-20x4.1575) = .52625

p2,3 = hF2(t0 + 0.05,x1,0 + 0.2625,x2,0–0.1575)

= .1F2(.05,7.2625,–4.1575) = 0.1(4–9x7.2625 + 14x4.1575)

= –0.31575

p1,4 = hF1(t0 + 0.1,x1,0 + 0.52625,x2,0–0.31575)

= 0.1F1(0.1,7.52625,–4.31575)

= 0.1(–6 + 13x7.52625–20x4.31575) = 0.552625

p2,4 = hF2(t0 + 0.1,x1,0 + 0.52625,x2,0–0.31575)

= .1F2(.1,7.52625,–4.31575) = .1(4–9x7.52625 + 14x4.31575) = –0.331575

x1,1 = 7 + (0.5 + 2x0.525 + 2x0.52625 + 0.552625)/6

= 7 + (3.155125)/6 = 7.5258541

(17)

x2,1 = –4 + (–0.3–2x0.315–2x0.31575–0.331575)/6

= –4 + (–1.893075)/6 = –4.3155125

(18)

The exact solution calculated by using Equation 7 are:

x1,1 = 7.5258355

© 2001 by CRC Press LLC

and

x2,1 = –4.3155013

(19)

The errors are 0.000247% and 0.000260% for x1 and x2, respectively. Per-step

error for the fourth-order Runge-Kutta method is difficult to estimate because the

method is derived by matching terms in Equation 12 with Taylor-series expansions

of x1 and x2 about ti through and including the h4 terms. But approximately, the perstep error is of order h5. For better accuracy, the fifth-order Runge-Kutta method

should be applied. For general use, the classic fourth-order Runge-Kutta method is,

however, easier to develop a computer program which is to be discussed next.

SUBROUTINE RKN

A subroutine called RKN has been written for applying the fourth-order RungeKutta method to solve the initial-value problems governed by a set of first-order

ordinary differential equations. It has been coded according to the procedure

described in the preceding section. That is, the equations must be in the form of

Equation 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 seven

arguments for this subroutine, the first four are input arguments where the last is an

output argument. The fifth argument P keeps the Runge-Kutta parameters generated

in this subroutine. The sixth argument XT is needed for adjusting the input argument

XIN. These two arguments, P and XT, are included for handling the general case

of N variables. Listing them as arguments makes possible to specify them as matrix

and vector of adjustable sizes.

FORTRAN VERSION

© 2001 by CRC Press LLC

QUICKBASIC VERSION

PROGRAM RUNGEKUT

Program RungeKut which calls the subroutine RKN is to be run interactively

by specifying the inputs through the keyboard. Displayed messages on screen instruct

user 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 following

inputs 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 only

9 functions can be handled by the program RungeKut, but the user should

be able to expand the program to accommodate any N value which is

greater 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 of

the independent variable, denoted as t0 and te (T0 and TEND in the

program RungeKut), respectively. It is over this range, the values of the

N dependent variables are to be calculated.

(6) Enter the stepsize, h (DT in the program RungeKut), with which the

independent variable is to be incremented.

(7) Enter the N initial values.

© 2001 by CRC Press LLC

QUICKBASIC VERSION

FORTRAN VERSION

© 2001 by CRC Press LLC

FUNCTION F

According to Equations 12 to 15, the Runge-Kutta parameters pi,j (matrix P in

the program RungeKut) are calculated using two FOR-NEXT loops — an I loop

covering 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 defined

by Equations 8 to 11. We create a supporting function program F as follows:

QuickBASIC Version

© 2001 by CRC Press LLC

FORTRAN Version

The 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, and

the 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 by

Equations 8 to 11 is manifested by a copy of the screen display:

© 2001 by CRC Press LLC

SAMPLE APPLICATIONS

OF THE PROGRAM

RUNGEKUT

As a first example, consider the problem of a beam shown in Figure 2 which is

built into the wall so that it is not allowed to displace or rotate at the left end. For

consideration of the general case when the beam is loaded by (1) a uniformly

distributed 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 the

free end of the beam, x = 40 cm. It is of concern to the structural engineers to know

how the beam will be deformed. The equation for finding the deflected curve of the

beam, 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 along

the beam. It can be shown that the moment distribution for the loading shown in

Figure 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 LLC

for

for

for

for

0 x 10 cm

10 x 20 cm

20 < x 30 cm

30 < x 40 cm

(22)

FIGURE 2. The problem of a beam, which is built into the wall so that it is not allowed to

displace or rotate at the left end.

To convert the problem into the form of Equation 3, we replace y, dy/dx, and x

by x1, x2, and t, respectively. Knowing Equation 21 and that the initial conditions

are y = 0 and dy/dx = 0 at x = 0, we can obtain from Equation 21 the following

system of first-order ordinary differential equations:

dx1 dx = F1 (x, x1, x 2 ) = x 2 ,

x1 (x = 0) = 0

(23)

and

dx 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 Equation

21 whereas the beam stiffness EI is to be set equal to 2x105 N/cm2 in numerical

calculation of the deflection of the beam using program RungeKut.

To apply program RungeKut for solving the deflection equation, y(x) which has

been renamed as x1(t), we need to define in the QuickBASIC program RungeKut

two 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 QuickBASIC

programming. In view of Equation 21, F2 in Equation 24 has to be defined by

modifying the subprogram function, here we illustrate it with a FORTRAN version:

© 2001 by CRC Press LLC

FUNCTION BM(X)

IF ((X.LE.0.).AND.(X.LT.10.)) BM = –1121./3 + 24*X-X**2/2

IF ((X.GE.10.).AND.(X.LT.20.)) BM = –871./3 + 4*T +

T**2T**3/30

IF ((X.GE.20.).AND.(X.LT.30.)) BM = 4*X–157

IF ((X.GE.30.).AND.(X.LT.40.)) BM = 4*X–160

RETURN

END

In fact, each problem will require such arrangement because these function

statements and subprogram function describe the particular features of the problem

being 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 errors

in % of the analytical solution ( = –0.68917) are tabulated below:

Stepsize h, cm

y at x = 40 cm,

5

2

1

.5

.25

.1

.05

Error, %

–0.79961

–0.73499

–0.71236

–0.70083

–0.69501

–0.69151

–0.69033

16.0

6.65

3.36

1.69

.847

.340

.168

This 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 loads

directly. This approach saves the reader from deriving the expressions for the bending

moment. Readers interested in such an approach should solve Problem 4.

A NONLINEAR OSCILLATION PROBLEM SOLVED

BY

RUNGEKUT

The numerical solution using the Runge-Kutta method can be further demonstrated by solving a nonlinear problem of two connected masses m1 and m2 as shown

in 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, their

instantaneous 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 in

term of y as:

[

? = tan ?1 (y + h ) b

]

and

[

L = (y + h) + b2

2

]

0.5

(25,26)

If we denote the cable tension as T, then the Newton’s second law applied to

the two masses leads to m1g–2Tsin? = m1d2y/dt2 and T-m2g = m2d2z/dt2. By eliminating T, we obtain:

© 2001 by CRC Press LLC

FIGURE 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.

? 2m2

?

d 2y 2m2

d 2z

+

sin ? 2 = g?1 ?

sin ??

2

dt

m1

dt

m

?

?

1

(27)

The displacements y and z are restricted by the condition that the cable length

must remain unchanged. That is z = 2{[(y + h)2 + b2]1/2[h2 + b2]1/2}. This relationship

can be differentiated with respect to t to obtain another equation relating d2y/dt2 and

d2z/dt2 which is:

2

d 2 z ?2

dLdy 2 ?? dy ?

d2y ?

= 2 (y + h)

+ ?

+ (y + h) 2 ?

2

dt ?

dt

L

dtdt L ?? dt ?

(28)

Equation 28 can be substituted into Equation 27 to obtain:

?? 2 2 m 2

d2y

sin ?

=

?gL ?

2

dt

L(2 y + h + L ) ??

m1

2

? 2

dLdy

? dy ? ? ??

gL

y

h

L

?

+

+

2

2

(

)

?

? dt ? ?? ??

dtdt

?

?

(29)

By letting x1 = y, x2 = dy/dt, x3 = z, and x4 = dz/dt, then according to Equation

3 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. It can be shown

that the final form of the system of four first-order differential equations are:

dx1

= x2 ,

dt

dx3

= x4 ,

dt

© 2001 by CRC Press LLC

2

2

dx 2 gL (1 ? 2 rm sin ?) + 4 rm v sin ?

,

=

2

dt

L 1 + 4 r sin ?

(

dx 4 ?2 v

dx

=

+ 2 2 sin ?

dt

L

dt

2

and

)

(30)

FIGURE 4. A numerical case where b = 3.6 m, h = .15 m, m2/m1 = 0.8, and initial conditions

y = z = dy/dt = dz/dt = 0 has been investigated and the results for -y, z, -dy/dt, and dz/dt have

been plotted for 0?t?12 seconds.

where:

rm =

m2

m1

and

dL

v2 = 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 conditions

y = z = dy/dt = dz/dt = 0 has been investigated and the results for -y, z, -dy/dt, and

dz/dt have been plotted in Figure 4 for 0?t?12 seconds. The reason that negative

values 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 varies

between 0 and 3.8358 m. The oscillation has a period approximately equal to

2x6.6831 = 13.3662 seconds, so the frequency is about 0.47 rad/sec. The masses

reach 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 the

oscillation for the studied period are listed below:

© 2001 by CRC Press LLC

MATLAB APPLICATION

MATLAB has a file called ode45.m which implement the fourth- and fifthorder Runge-Kutta integration. Here, we demonstrate how the sample problem used

in the FORTRAN and QuickBASIC versions can also be solved by use of the m

file. The forcing functions given in Equations 10 and 11 are first prepared as follows:

© 2001 by CRC Press LLC

If 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 of

investigation 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 which

the forcing functions are defined. The second and third argument the initial and final

values of time, respectively. The fourth argument is a vector containing the initial

values of the dependent variables. The resulting display, after rearranging T and X

side-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 created

to 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), and

spring 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 LLC

FIGURE 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 by

entering 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 2

in the program NewRaphG) of about 0.28 cm at approximately t = 2.5 seconds and

finally settles to a static deflection of 0.25 cm, and that the maximum ascending

velocity is about 0.18 cm/sec.

As another example of dynamic analysis in the field of fluid mechanics, Figure 8

shows the flow of a fluid between two connected tanks. The valve settings control

the amount of flows, qi for i = 1,2,3. The levels of the tanks h1 and h2 change in

time depending on these settings and also on the supply rates Q1 and Q2 and the

discharge 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 the

formulas:

A1

dh1

1

= Q1 ?

(h ? h3 )

dt

R1 1

© 2001 by CRC Press LLC

and

A2

dh 2

1

= Q2 ?

(h ? h3 )

dt

R2 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 commands

shown below.

where A’s are the cross-sectional areas of the tanks, and h3 is the pressure head at

the junction indicated in Figure 7 and is related to the discharge rate q3 by the

equation:

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 differential

equations in h1 and h2 to be:

© 2001 by CRC Press LLC

FIGURE 7. The flow of a fluid between two connected tanks.

A1

dh1

= Q1 ? a1h1 + a 3h 2

dt

and

A2

dh 2

= Q 2 ? a 2 h 2 + a 3h1

dt

a1 = ( R 2 + R 3 ) ?

and

a 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 both

tanks.4

MATHEMATICA APPLICATIONS

Mathematica 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 LLC

Out[1] = {{X1 -> InterpolatingFunction[{0., 25.}, ],

X2 -> InterpolatingFunction[{0., 25.}, ]}}

In[1] shows that NDSolve has three arguments: the first argument defines the

two ordinary differential equations and the initial conditions, the second argument

lists the dependent variables, and the third argument specifies the independent

variable and the range of investigation. Id is a name selected for easy later reference

of the results generated. Out[1] indicates that interpolation functions have been

generated for X1 and X2 for t in the range from 0 to 25. To print the values of X1

and 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 = 0

t=

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 LLC

X1 = {0.}

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

X1 =

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 MATLAB

application. In[2] shows that the replacement operator/. is employed in X1[tv]/.Id

which requires all t appearing in the resulting interpolating function for X1(t) created

in the Id statement to be substituted by the value of tv. The looping DO statement

instructs that the tv values be changed from a minimum of 0 and a maximum of 25

using an increment of 1. To have a closer look of the overshoot region, In[2] can

be modified to yield

In[3]: = (Do[Print[“t = “,tv,” X1 = “,X1[tv]/. Id,” X2 = “,X2[tv]/. Id],

{tv, 2, 3, 0.1}])

Out[3] =

t=2

t=

t=

t=

t=

t=

t=

t=

t=

t=

t=

X1 = {0.267432}

2.1

2.2

2.3

2.4

2.5

2.6

2.7

2.8

2.9

3.

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 overshoot

occurs at t = 2.6 and X1 has a maximum value equal to 0.28502.

6.3 PROGRAM ODEBVPRK — APPLICATION OF RUNGE-KUTTA

METHOD FOR SOLVING BOUNDARY-VALUE PROBLEMS

The program OdeBvpRK is designed for numerically solving the linear boundaryvalue problems governed by the ordinary differential equation by superposition of

two solutions obtained by application of the Runge-Kutta fourth-order method. To

explain the procedure involved, consider the problem of a loaded beam shown in

Figure 8. Mathematically, the deflection y(x) satisfies the well-known flexural equation.5

d2y M

=

dx 2 EI

(1)

where M is the internal moment distribution, E is the Young’s modulus and I is the

moment of inertia of the cross section of the beam. For the general case, M, E and

I can be function of x. The boundary conditions of this problem are:

© 2001 by CRC Press LLC

FIGURE 8. The problem of a loaded beam.

y(x = 0) = 0

and

y(x = L ) = 0

(2)

where L is the length of the beam. The problem is to determine the resulting

deflection y(x). Knowing y, the moment and shearing force, V, distributions can

subsequently be determined based on Equation 1 and V = dM/dx. The final objective

is 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 if

in additional to the initial condition given in (2), y(x = 0) = 0, we also know the

slope, , at x = 0. But, we can always make a guess and hope that by making better

and better guesses the trial process will eventually lead to one which satisfies the

other boundary condition, namely y(x = L) = 0 given in Equation 3. In fact, if the

problem is linear, all we need to do is making two guesses and linearly combine

these two trial solutions to obtain the solution y(x).

Let us first convert the governing differential Equation 1 into two first-order

equations as:

? dx1

? dx = x 2 ,

? dx 2 M

= ,

?

? dx EI

x1 (x = 0) = y(x = 0) = 0

(3,4)

x 2 (x = 0) = ?(x = 0) = ?0

To apply the fourth-order Runge-Kutta method, we have to first decide on a

stepsize h, for example h = L/N we then plan to calculate the deflections at N + 1

locations, x + jh for j = 1,2,…N since j = 0 is the initial location. If we assume a

value for 0, say A, the Runge-Kutta process will be able to generate the following

table

x

x 1=y (1)

0

0

h

y1(1)

2h

y2(1)

…

…

jh

yj(1)

…

…

Nh

yN(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 not

equal to zero, then the value incorrectly selected, we have to make a second try by

letting (x = 0) = B to obtain a second table by application of the Runge-Kutta

method. Let the second table be denoted as:

x

x 1=y (2)

0

0

h

y1(2)

2h

y2(2)

…

…

jh

yj(2)

…

…

Nh

yN(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, if

the problem is linear, the solutions can be obtained by linearly combining the two

trial 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:

? +? =1

and

?y(N1) + ?y(N2 ) = 0

(9,10)

Equation 10 is derived from the boundary condition y(x = Nh = L) = 0 and based

on Equation 7. Equation 9 needs more explanation because it cannot be derived if

y(x = 0) = 0. Let us assume that for the general case, y(x = 0) = . Then Equation

7 gives + = which leads to Equation 9. Using Cramer’s rule, we can easily

obtain:

? = y(N2 ) D

and

? = ? y(N1) D

(11,12)

and

(1)

D = y(2

N ? yN

(13)

NUMERICAL EXAMPLES

Let us consider the problem of a loaded beam as shown in Figure 8. The

crosssection of the beam has a width of 1 cm and a height of 2 cm which results in

a moment of inertia I = 2/3 cm4. The reactions at the left and right supports can be

computed to be 5/3 N and 25/3 N, respectively. Based on these data, it can be shown

that the equations for the internal bending moments are:

for 0 ? x ? 20 cm

?5x 3,

?

M( x ) = ?

65

1 2

for 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 decide

on a stepsize of h = 1 cm and proceed to prepare a computer program using the

fourth-order Runge-Kutta method to generate two trial solutions and then linearly

combining 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 presented

later has produced the following display on screen:

© 2001 by CRC Press LLC

FORTRAN VERSION

© 2001 by CRC Press LLC

The Subprogram FUNCTION F which defines the initial-value problem is coded

in accordance with Equation 14. The two trial initial slopes are selected as equal to

0.1 and 0.2. The trial results are kept in the three-dimensional variable XT, in which

the deflection y(k)(j) for the kth try at station x = xj = jh is stored in XT(1,j,k) whereas

the 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 because

XT(1,KS,NTRY) is transmitted as XIN(1) and automatically the next value

XT(2,KS, NTRY) as XIN(2), and the computed results XOUT(1) and XOUT(2) are

to be stored as XT(1,KS + 1,NTRY) and XT(2,KS + 1,NTRY), respectively. Notice

that there are only two dependent variables, NV = 2.

After the weighting coefficients (ALPHA in the program) and (BETA) have

been calculated, the final distributions of the deflection and slope are saved in first

and second rows of the two-dimensional variable X, respectively. It should be

emphatically noted that the solutions obtained is only good for the selected stepsize

h = 1 cm. Whether it is accurate or not remains to be tested by using finer stepsizes

and 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. The

needed modification of the above listed program to include this change in the stepsize

and testing of the difference in the maximum deflection is left as homework for the

reader.

© 2001 by CRC Press LLC

QUICKBASIC VERSION

MATLAB APPLICATIONS

In the program RungeKut, MATLAB is used for solving initial value problems

by 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-value

problem governed by linear ordinary differential equations. To demonstrate the

procedure, the sample problem discussed in FORTRAN and QuickBASIC versions

of the program OdeBvpRK, with the aid of function BVPF.m listed in the subdirectory , can be solved by interactive MATLAB operations as follows:

© 2001 by CRC Press LLC

Notice that format compact enables the display to use fewer spacings; YT1 and

YT2 keep the two trial solutions, and ode45 automatically determines the best

stepsize which if used directly will result in a coarse plot as shown in Figure 9. The

plot showing solid-line curve for the deflection and broken-line curve for the slope

can, 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 is

to be applied as follows:

© 2001 by CRC Press LLC

FIGURE 9.

The result of plotting the spline curves is shown in Figure 10.

MATHEMATICA APPLICATIONS

The 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; readers

are therefore referred to program RungeKut to review the method itself and the

development of related programs and subprograms. The boundary-value problem is

to be solved by continuously guessing the initial condition(s) which are not provided

until all boundary conditions are satisfied if the problem is nonlinear. When the

problem is linear, then only a finite number of guesses are necessary. A system of

two first-order ordinary differential equations which governed the loaded elastic

beam problem previously solved in the MATLAB application is here adapted to

demonstrate the application of the Runge-Kutta method.

© 2001 by CRC Press LLC

FIGURE 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]/. Id1

Out[3] = {3.00052}

EI value and F(X) are defined in In[1]. In[2] specifies the two first-order ordinary

differential equations involving the deflection Y and slope YP, describes the correct

initial condition Y(X = 0) = 0, gives a guessed slope Y (X = 0) = 0.1, and decides

on the limit of investigation from X = 0 to X = 30. In[3] interpolates the ending Y

value 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 LLC

In[5]: = Y30Trial2 = Y[30]/. Id2

Out[5] = {6.00052}

Linear combination of the two trial solutions is now possible by calculating a

correct 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*b

Out[8] = {–0.0000174888}

Finally, the actual deflection and slope can be obtained by providing the correct

set 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 LLC

0

1

2

3

4

5

6

7

8

9

10

11

12

13

14

Y=

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=

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

Y=

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 — 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

finite-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 = Ri and R = Ro, respectively.

As shown in Figure 11, it is then inflated 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 equation6

Z(R) satisfies Equation 1 is for Ri