WWII

WM

W

cm THE NUMERICAL SOLUTION OF LINEAR
AND MONLSNEAR SYHEM MflBELS

 

The-sis for the Dogs“ of Ph. D.
‘MiCHlGAN STATE UMVERSSTY
Dinkar Shankar Rana
19.62

This is to certify that the

thesis entitled

ON THE NUIJHSRICAL SOLUTION OF LINEAR AND NOI‘EINEAR SYSTEM

MODELS
presented by

Dinkar S. Rane

has been accepted towards fulfillment
of the requirements for

Mdegree inmal Engineering

    

 

/7 ‘7’ .2
,fl -( ,..--"'. ‘ .
fl, , ‘ r . fi/ 4/:-
V’vvzf- 4"”, I" 7 r .71 \‘31 W,
/ Major professor ,
t.»

Date November 9, 1962

0-169

LIBRARY

Michigan State
University

 

ABSTRACT

ON THE NUMERICAL SOLUTION OF LINEAR AND NONLINEAR

SYSTEMtMODELS

Dinkar Shankar Rane

ON THE NUMERICAL SOLUTION OF

LINEAR AND NONLINEAR SYSTEM MODELS

Dinkar Shankar Rene

AN ABSTRACT OF A THESIS

Submitted to
Michigan State University
in partial fulfillment of the requirements

for the degree of

DOCTOR OF PHILOSOPHY

Department of Electrical Engineering

1962

ON THE NUKERICAL SOLUTION OF LINEAR AND NONLINEAR SYSTEM

MODELS

Abstract of Thesis

by

Dinkar S. Rane

The general availability of digital computers has lead
to the development of various numerical techniques for solving
linear and nonlinear systems. The relative merit of each such
technique varies from problem to problem. However, in all cases
one of the primary considerations is computation time and numerical
accuracy. The consideration of computation time has lead to
the develOpment of this thesis. The numerical accuracy has been
increased in most cases or at least kept to the same degree as in
the existing methods.

The mathematical model of the system used is assumed
to be given in normal or standard form. Whenever components
with nonlinear algebraic terminal eguations are prestnt, the mathemat-

ical model of the system may appear in a semi-normal form. The

general mathematical model is then taken as

 

 

 

X [F1(X,Y,z)
Y = F2(X,Y,Z)
I o F3(X,Y,Z)

Abstract Dinbar S. Rane

JLC.

For linea system, i.e., when

a recursive formula is developed for the numerical solution. This

.‘

formula can meet any prescribed degree of accuracy wither; increasing

\

J- V

the computation time per step. Theoretically th; computation time

is one-fourth that of the standard Range—Kutta method.

Some mixed, linear and nonlinear forms are also considered.
The Runge-Kutta method has been modified to increase the damputational
speed. The numerical accuracy has been increased in some cases while
the original degree of accuracy has been unchanged in others. The

n

overall e Tect is a faster numerical solution with a high degree

0

or accuracy. The results are compared with the results by other
existing methods such as the Runge-Kutta method, predictor-corrector

methods and others.

0

Finally the semi-normal form of the mathematical model 18

considered end two prc~ .ures are described for obtaining a numerical
solution. 'hese procedures also yield a shorter computation time with-
ou. loss in accure- .

‘,:erical examples a‘. iven for the purpose of comparing the
cetual ,lasive wafttat3UK "fine and numerical accuracy and indeed

Support the t. .' 'lcal u;velopment.

ON THE NUMERICAL SOLUTION OF

LINEAR AND NONLINEAR SYSTEM MODELS

By

Dinkar Shankar Rene

'A THESIS

Submitted to
Michigan State University
in partial fulfillment of the requirements

for the degree of

DOCTOR OF PHILOSOPHY

Department of Electrical Engineering

1962

ACKI‘IONLEDGEMEI‘IT

The auth r is grateful to Dr. H.E. Koenig, his thesis
adviser and major professor, for His constant guidance during the
preparation of this thesis and during the pe*iod of investigation
precseding its completion and to Dr. R.J. Reid and Dr. G.P. Weeg
for their helpful suggestions. He also wisnes to thank Dr. L.N. VonTersch
and Professor J.W. D nnell for their constant encouragement.

~: 1- A ->:« 9:. -:<- -x-~>e -x.x- -x- -x- ->.L as

Section

\Jl

TABLE OF CONTENTS

Page

INTRODUCTION . . . . . . . . . . . . . . . . . . . . . . l

LINEAR SYSTEMS OF DIFFERENTIAL EQUATIONS.

KI.)

MIXED LINEAR AND NONLINEAR DIFFERENTIAL EQUATIONS. . . . . 18
NONLINEAR ALGEERAIC EQUATIONS. . . . . . . . . . . . . . . 3C
mammEsmummm. ... .... ... .... ... .. M+

CONCLUSION. . . . . . . . . . . . . . . . . . . . . . . .. 53

iii

LIST OF TABLES

Table Page
5.1 SOLUTION OF LINEAR SYSTEM OF EQUATIONS ................... A6

5.2 SOLUTION OF MIXED LINEAR AND NONLINEAR EQUATIONS

WITH CLOSED LINE.R SET ................................... us
5.3 SOLUTION OF A CLASS OF LINEAR AND NONLINEAR

SYSTEMS .................................................. 52
5.A SOLUTION OF ANOTHER CLASS OF LINEAR AND NONLINEAR

SYSTEMS .................................................. 55

5.5 SOLUTION OF SYSTEM MODEL WITH NONLINEAR

ALGEBRAIC EQUATIONS ...................................... 57

Appendix

B.

L ST OF APPEI‘IDICES

RECURSSIVE FORI-IULA FOR VARIOUS DEGREES

OF ACCURACY .

RIHIGE-KUTTA I-‘EFHOD AND A PREDICTOR

CONNECTOR I'IJI‘NOD. . . .

I. INTRODUCTION

A prerequisite for any system design procedure, is some
kind of mathenatical model of the system. The model may be obtained
Ly various transform techniques such as the Laplace transform, or in
the time domain as a set of differential equations. The transform
techniques have been applied very effectively in the area of linear
systems but do not apply to nonlinear systems. The time-domain analysis,
however, can be extended to include nonlinear systems in general.

A time-domain mathematical model of physical systems
particularly suited to computer solutions 115 in the normal or standard
form. General procedures 'er formulating mathematical models of
physical systems in this form are given by Nirth(l).

The objective of this thesis is to develop high speed
numerical methods for obtaining solutions to the type of mathematical

model found in the analysis of physical syssems. The general form of

the seminorma] form

p.
H.
U)

the mathematical model considere

 

 

 

 

- C "' ‘ 1 —h L]
X F1(A,Y,Z) Llrl(t)
d r
E? Y = r2(X,Y,A) + B2E9(t) (1 l)
o F.(X.Y.Z B.E-(t)
_j 3 5).)

 

 

2, F3 are vector

functions and the last equation is an algebraic equation nonlinear in 2.

where X, Y, El(t) and E?(t) are vectors; Fl, F
J

IN
...;-

The problem, then is to develOp a method for obtaining a
numerical solution of nonlinear differential and algebraic equations.

Various methods have been suggested for obtaining the ntmerical

The most CONJOD of these metnous

solutions to such problems.
are Runge-Kutta, Num rov's and the various predictor-corrector schemes,

such as those of Nystrom, Milne, Hartree etc. All these methods

are compared with the procedures of this thesis.

In this th sis Specia fO‘ms of the general mathematical
model are considered. In section 2, linear differential systems are
considered and a recurssive formula is develOped for the solution. The
solution at nth interval of time in terms of initial conditions is

shown to have a form resemtling that of the nth state vector in

Markov Prccesscs.

In section 3, combinations of linear and nonlinear differ-
ential systems are cons'dered. A procedure applicable to a special

class of nonlinearity is considered first with sussequent extention to

a more general form.

In section A the procedure is extended to certain classes
) L
of mathematical models involving both nonlinear algebraic and differ-

ential equations.

II.

 

 

Let the mathematical
P — r
X a
1 ll
X2 {121
,1 .1-
\u. L) g o
x a‘
L k El
.1»

 

 

LIIHLKR

model of a linear

H
R)

{D

J— ,

 

 

l) for the vectors

 

DIFFERENTIAL EQUATICNS

nL;

 

n '

of the form

‘(h), X(2h) ...

 

(2

 

IL—

X(nh) . . .

a’faylor's series ex ansion, and then this solution

or in vector notation
3— [x(t)]
at
”no solution
is obtained first as
is put in a corcise

2.1 Taylors Series Exyansion

 

form to obtain a recurssion

formula.

An analytic function f(t+h), where h is a constant, can be

GXpanded in terms of

to t as

.o( .1-
. I

its

SLICC

t) +

RHFJ

essive derivatives with res

.‘ J-
per

J.

b) +
r‘ '5
i: .3

 

fl
‘1

_u_

cascuuently using vector notation, (2.3) can be OXpanded as

X(t+.

A h2 ,2 ‘n ,n
‘ _ ‘ u L a L h u +
: 2((L)+LE'E'X(LJ)+7—7X(o)+...+F—EX(V)+.

dt' ° dt

(2.h)

where the indicated derivatives of the vector X(t) as calculated from

(

o o) n
L- n 1... L15

1"

 

 

"\

A X(t) + BE(t) 2.5.1)

A“[3x(t) + rs(t)1 + ArE(t) + PE"(t)

A“ [xx(c) + EE(L)] + 3‘ ss'(t) + ABE"(t) + BE"'(t)

ALL X(t) + A3 EE(t) + AI DE'(t) + A133"“) + EE"'(t)

(23.5.lt)

In equations (2.5.1) through (2.5.h) the derivatives of E(t) can

. . , . _ h H .
be obtained approx1mately in terms of E(t), E(t + 2) etc., by using
the forward differences. The third and higher derivatives involve

h i , . . fi ‘_
E(t + 3- ): E(t + 2h) etc. However, an apprOXimation 15 made for this,

. . . , fir h
so that all the derivatives are eXpressed in terms of L(t),:4t +-3),

E(t + h).

 

E's) = 8-;- m) 2’ 5
us __ h
2
P l
= h [E(t + 3 h) - E(t)] (2°C)
h" ‘1 It 4- ~‘\ 1 1: 1 -
L(U):7-‘EL(U)—E[E (ti‘E)-E(t)]
I). , l r‘
z ...2 [mt +11) - 2 E(t + 51:) + E(t)] (2.7)
h
n. _ 2 ' n ,_ h '71"
E (t)=*ELE(L+-2- '(t)]
9 3p 11
.... —: [E(L + -,—.—) - 3 E(t + h) + 3 E(t + 2,) - E(t)]

dh . . . .
Using quadratic extrapolation E(t + f—) is approximated 1n
2

1.
terms of E (t + h), E(t + %) and E(t), so that
4 3h lg 5 _ L L h l A
E(b +-§—) —'§ E(t + n) - 2 E(s + 2) + 5 E(t) (2 9)

01‘

Where

er—J
l
m
:_J
Q J
Lo

. l l .. l .2 {i l : 3
+n[-,—U+—,:nz\+—-g.11\ 4-21-11AJIBm)
‘3 J la c_‘*'|r

1 1. C a 1 3 3
l[ 7: [U] + —]1[A] + —— 1 [A + 7,— 1] [A] ] [sefi 1

h [ — [u] + 7 [A] J [a] (2.1:)

If X(nh) is denoted as X(n) and X(nh + h) as X(n + 1),

then a the discrete points, t = nh, (2.12) becomes

X(n+1) = [t] X(n) + s E(n) + 30 E(n + ) + a E(n+1) (2.17)

mflhd

1 3

Eauation (2.17) is a recursion formula for a solution to the system of

.5

equations (2.2) i.e.

x(1) = [a] X(o) + ,31 E(O) +3.15%) + 33 12(1) (2.1.3.1)
x(2) = [a] X(l) + s1 E(l) + :0 E(%) + ;3 E(2) (2 1 a)

and so on.

1,

Note that the matrix [a] is of order n x n and [31], [£2

and [31] are of order n X m. These matrices are calculated once for all
3
at the bevinning.

Vt)

Inrmugrpractical problems the input signals are all of the

P 81(6) Fol-q
EE(t) = = c(t) = [E] e(t) (2 1))
311(t) ”1.
L _ _.‘;

 

 

 

 

-Q—

[3 ] E(n) = [3 ] [E] e(n)

If we let

P
R)

"C a
U
L—a
r—\
I?!
L—J
II
P”
U.)
A
[‘J
I" D
C.)
LA.)
V

then the recursion relation in (2.17) reduces to

X(r+l) = [a] X(n) + X e(n) + x2 e(n + %) + X e(n+1) (2.21)

1 3

where X x2, x are column matrices.

1’ 3

2.2 Comparison with the“ Solution Procedures

 

There are many numerical procedures available for a solution
to a set of linear differential equations. Most common and widely used

of these is the Runge-Kutta method. The predicotr-corrector methods

! W

of Milne, Hartree and NyStrom and otners, are also used for special

'1 .1.

forms of linear difierential equations. However, the procedure given

1. ..' ,,.. .. 1,. ' -. ,s -' ,1 4.1.1 ...," p .. 4 .L -' (:9
asove gives a much lassel OJIUblOU and tag accuracy Can be maintained

above any "esired derree. The dearee of accuracy is a function only

Q

of the eXpressions for [a], [Q ], [pm], [99] as shown in Appendix A.

l c. J

The recurssive formula 2.12 remains unchan ed. Thus retardless
1.2 J L)

of the degree of accur-cy the computation time is the same. Actual

‘l 'v

computer programming indeed proves this point.

0.1.?

InSpection or bfle widely used Runge-Kutta method shows

0
fl . o . o -\ o o o L‘
that ior a system with n variables and m ariVing functions, (An +

O

hnm + 5n) multiplicatiCIs and (An + lln) additions are required.
. - 2

On the other hand, in the procedure of equation (2.17), (n + jnm)

multiplications and (n" + jnm) additions are required. This gives

('3 q
. q 4.". , i , . . o j V \. L .
a saving of (;n + nm + pn) multiplications and (3n + nm + lln)

7 f 0 _ .' , _- ,‘ - - _ ,. .21 ' ..I. ' . ,0 s- -, . - ' i .1.
ac itions. inis redresenss a relative time rector of appr Ximately

L

(x

l G ' . .."fi. " - . —
4 to l ior re atively shall m. Even ior the largest value of m,

namely m = n, the tihe iaCter is better than 2 to 1. In practice,

‘ ' r '1‘ ‘» . ‘V ‘7' ‘ J- o -. ‘ '1 ' I14 . ‘I‘ ‘ J‘ 1 » ’ I'I an J‘ > ' . I
however, this ahvansaje is reduced sli ntly «ecause oi th time

)

required to print out the intermediate results.

The bookkeeping required in the procedure given above is
also considerably less as compared to that in the Runge-Kutta method.
In Runge-Kutta method the solution at the nth interval must be retained.

At the same time, the values of the matrices K as defined

Km K K

1., .4) 3)
in Appendix C are also to be stored. The evaluation of K1’ K2, Ki, KIL
3 e.

.1.

is not symmetric, hence this requires additional bookeeping. On

the other hand in the solution given above it i.

(J
:3
(D
O
(D
U)
(n
f.)
’1

J- V
1 y a J U11 I J fj ‘j c o I -3
only the solution as n interval and he coeiiiCient matrices [a] [91]

J [S ] once calculated are used throughout the solution of a given

A. '7
‘4 .3

[5

problem.

The comparison of the procedure given above with the virious
predictor correcter methods also shows an advantage. The formulas for
a typical pr-dictor-corrector method surfesteu ov Hamming are given in
Appendix B. A careful study of these form: las shows hat the first and
second derivatives of the variables have to be evaluate d twice for each
step, t.;resy in re as in" the computation time considerably over that
required in the above procedure. A further study indicates tlat the
va ue of X(t ), X'(t) and X"(t) a‘ the points (n-1) and n are required
to evaluate the solution at (n+1). This rerre"erts a substantial amount

bookkeeping he nee puts a limit on the size of problem that can be
solved. Another disadzantage in the method of Lemming as compaied to
a ..4. (3H

J. 1 ‘ 0 ~ ‘r‘ ‘ .0 ‘ 1’" ’ _‘ r. J. 1 .. r! O n
that given asove is the us or sue oCCODd( 3 Hi atives while the s.ssem

,

U
(‘0 Jw . h -. ~ fl ‘1 w 4'1“ a

01 equations (2.a) involves nly the fiist de'l ivatives.

n 1

¢:.3 State Vector and State Matrix

 

ution at any

J
J
A
P.)
r_J
~J
V
(D
*3
w
(.4
C)
C)
(D
U}
C -
L)
C
('1
F4

The recurssive rOlJJHu
int rval, in terms of the solution at the preceeain; interval and the
driving function at all the previous intervals and mid—inis ervals. Lf
recursive substitution, for :.u1a (2.17) can be used to develop an
expression for the solution at any interval in terms of only he initie 1
conditions on the dependent variables azid the 'riving function at all
the previous intervals and mi'-intervals. Such an expression is very
conver ient to use in many cas es. Moreover, in the special case o
:onstant input, this xiiessicn has a very interesting fo 1*m. Even in
the general case of time varying driving functions the rcuif* d expression

4-1

has a striking similarity wit} exoress ion in the Value- Deteriination 0 er-

ation of Markov Process.

The recurssive formula 2.17 can he rewritten in detail

. the intervals h, 2h, ... nh as:

X(l) = 2 X(o) + 91 E(o) + 22 E(%) + 23 2(1) (2.22)
x(.:) = c X(l) +931 13(1) + ,2 E(l-,l;) + :33 E(l)

) + 20 2(2) 1 2.23

. _ l
+elms+p ms)+2non mam

N
A

:5
v

H

01"

If we let

then

a X(n-l) + 31 E(n-l) + 92

-12-

'3 E(J) +2 E(J +

aura

2 h l 2 2
l E(l) + sq e(LZ) + )3 E(a) J
c. a. J
\ ,— 7-: 1- W
, e(n—.) + sf L(lg) + ad E(n-l) J
l L 2 a
1
E(n-l) + r\ E(nefi) + 33 E(n) ]

.32. (J+l)]
J

Comparing this to the value-determination formula in the

Markov nroeess, namelv

We see that they are identical in form.

If for time varying driving functions we let

. _. (.2 . 1 .2 . 2
E(J) = :21 3(a) + 22 E(J + 2) + '23 M.) + l) (2.20)

r‘. ‘

Equations (2.2;) and (2.31) indicate that the linear
mathematical model of a control system given by (2.2) can be put

in the same form as the mathematical models of Markov Processes.

Thus, the ontimisation nroccdures using dynamic programming given by

( C) . .

Howard can U3 aeplied dir-ctly to deterministic systems. The possibil—

ity of cxploitiwg these techniiues through this type of recurssion

formula represents an area worthy of further investigation.

-111-

In many control systems applications the only information

required is the outgut signals

{.3

r! “\ 4r, \ -'- '. , , - \ v. .‘ J— ‘V . 3w V
s the iii: eldn oi LAA? in.ut blnjfll.b.
" ‘2‘ \ ‘ ‘ “\ .. -. I ‘j "w ""' \ '. 1* ’. ‘ ". V' ‘

In otnei woids, only the iesyonse oi the system is of interest.

This is also the case when numerical solution to a higher order

differential eqaation is require . Let the differential equation

’1 < D I} 1
L 4‘. (L f‘ {1
i——-+ Cn-l ————— x + Cn—2 ————— x + ... + C-

 

x + COX + f(t) = 0

W4- 1 - w -* *4,
A f’v I
(.L p \.ki} Gt k4. LI

 

dX . C K d 2 n
at‘ = [C , m = .L) ... n — zhi-l (2;.3;)
(AU (— _n"-L

Equation (2.32) when written in normal form then becomes

rt; _ th. -]
i

. .
4’;

X If

3
d __
at : :

 

 

 

 

 

 

 

 

 

 

[—5 1 o ... o o "I T — F‘s——
O O l . . . 3 C x C
J i
O C 3 . . . 3 3 56 L
= ; : : : : I + I 2.33!)
O O O . . . O l xn-2 O
c0 -c1 ~c2 . . . -cn-2 -cn-l xn-l -f(t)
__ _i __ __J __ J

In such a case only the solution for x is of interest.

In control systems, the response can be obtained conveniently from

Q

the earlier results an' particularly from equation (2.31). Let

this set of equations oe arran ed so that the output variables

Xo(t) are the first-h elements of the state vector X(t) and let

E(n) denote the first h-rows of [a] . Then from equation 2.31 the

n~l
Xo(n) = E(n) X(o) + :E: R (n—J—l) E(j) (2-35

1
U

~16-

r1~ ,. .. - .. ..°. — - 12, .+ . ' ‘2 _ ».
ine choression for the out us rariasles

ran
kuo

can to written
7 _

Xo(l) ..(i)

E(o)
Xo(2) R(l)

R(2)

E(nmi)

.R(n) E(n-l)

I—

lic)(:i)

 

 

 

 

 

In the sfiecial case where there is only one in

J. h ’.
(‘1 ’3‘) 1 (‘v
s) ukud .Ls)

-. . .,.W , 1‘ .. "'7‘. .1.-. ,. .
.Litnyie, \hhfll the txritrol s;

7-

J'}0 -2 2|. '121‘2
Ll4fSU 45

rows of 9L8

— ~-.-. 1r —_ J‘jw'
l\\h) g'CHi‘BbCIloo u.;k3

'fiprj ”.3 ’3,’
k.» ska n- u JL

) seco es

>.
A
r
V
ii
to
V
C)
J
a)

O)

XCfl +

 

 

 

 

 

“ J- J- — 1. ’7
[(0 (it; v '— L1, (.11)

E(ii-Z)

.. '
.-._3'~r/.3
\‘ 4-1.‘(LA’1

C)

C)

'I)“‘;" “- 1‘; .
.1. .2 \lei :1 ‘-‘ :1.

1“
x [a]‘, i.e.

C)

) R(O)

 

 

 

E(o)
Eu)
E(qa)
IE(II~]-)
— _l

..r
(2.;e)

“.+ p-nsml 22
2.2/LL La 01L); Lg». ) kwi)

Qital computer,

—1I ‘_ “'1'
o ei(o)
o 01(0)
0 '3 (2)

C0

 

 

 

«I ~ ' C‘ . l 1' ~ J--_. .<'
131‘: s.) 18 a CO 11.321 1..C1ol'l...

[X] = [R] [X(oH + is) ['37.] (9.3:)

O
F.

lflieitz [S] :is IJJUOl"CTl£Rl:tUJlr i1] tin: su.gxxtricn.s t) . — c), l' , .

Eiiation (2.3 ) is a yarticularly convenient model to use in the analysis
tems such as sampled data control systens.

In such a vstem the vectors of discrete Maanitudes c, ’ ' = l 2 ...
I ) t.) l ) ’ )

I \ “v1 ‘ ’ 1‘ r‘ '. V 4' VJ' --~.- "'— ' V -' . - ' 1 7 . r- . 4.- ‘nrs J.“ ‘I‘.
regiesent the discrete ouslus levels of Lao ululeal comjucei.

III. MIXED LINEAR AND M NLINEAR DIFFERENTIAL EQUATIONS
The vast majority of jhysica_ systems of practical in
in engineering contain hoth linear and nonlinear components, i.e.,
the system seldom contains all nonlinear components. When such is the
case it is sometimes possible to select the variables (i.e. a

formulation tree in th svstem :raoh so that the normal form model
./ C ..

rw-nvvnn c‘ ' t‘n ~ w
appears in he form

I I
A
LU

...—J
v

1’" 1-1 ' ' ,. _, . ,2 a. K , 2"" 4.3, :- ij .r» m.
wncre g is a lineai facupl luncsion oi tnc rein

and F is nonlinear vector function. For such a system a substantial
savina in computer tfme can be realized ey combining the procedur
c'ven in section 2 with other known procedures.

3.1 Modified nungc-Kutta Solution for Mixed Linear and Nonlinear Syst—ms

4- (2

With Closed Linear See

 

Solution to the second set of (3.1) can be obtained
independent of the first set. Let it be obtained by the procedure
of section 2.

The modification of the Runge-Kutta method required to use
this procedure is explained by first examining the standard Runse—Kutta

method .

18

-19-

The standard Range-Kutta procedure for obtaining a solution

><
I...)
A
:3
+
H
v
u
><
A
s
+
Cflh‘
W
.1.
h)
71
+

1 - 1'? 2K3 + K2] (3.3)

_ _:_L_ r ’3 r) 1. 3 ‘l

X2(n+l) — X2(n) + 6 LLl + 2L2 + 2L3 + Lh‘ (3.4)
K = F [A (n), Xn(n), nhj 2 i

1 cf ‘1‘?)
L1 = . Jr [A (n), nn] 3 C)
K _ u v W (n) + i K r ) _ L a + l )1] 2 C)

2 — L (1 2 l, {\2 I1 + , I11; 2 A 3
L = i r [x (.) + l L nh + s h] (3 is)

2 2 s a 1’ 2

A
U)
—J
V

L3 = h F2 [X2(n) +-% L2, nh +‘% h] (3 11)

K = h F [Xl(n) + K3, X2(n) + L2, nh + h] (3.;)

— h F [Xq(n) + L3, nh + h] (3.12)

.1.

It can be seen that in order to obtain a solution to the

first set, the quantities Ll’ LO, L L3L must l“e obtained from the
L— t

3)
second set. However, if the recurrsive formula of section 2 is used
to obtain a solution to the second set, the quantities L1’ 0, 2, a

(— 3 w-

are not available.

To show how the Runge-Kutta method can be conveniently

.l.1

modified so as to avoid this difficulty, sujpose a solution to the

se en? set has been oLtained by a recurrsive fornula of section 2

'3

i
\

o ~ ~ 1 .21 5
at the pOints t = h, an, 3;, ..., nn, (n+l)h and t = Eh, 3n, En,...,

 

2n+l .
( ) n

O
c.

Consider now the recursion formula

{3 Xl(n+l) = X + 2 K; + 2 K5 + K']

J

C‘ x .
.VV-aAC ‘53

y
(

f) 1
K' = h Fl [x (n) +-—g , x2 (an), nn + Q; (3.16)

h

>—..a

1 F [x (n) + K}, x (n+l), nh + h] (3.17)
1 1 3 2

It can be noted that equations (3.13) through (3.17) are
similar to equations (3.h) through (3.3) except for the following

changes.

1
L is renlaeed bv X (n+-
- u IN I)

. (.

L.

\hJPJHJ
{—1

Iv.)

. . l
L_ is replaced sy Xq(n+;)
L ,

C.

ft

X0(n) + Lq is re laced by X, (n+l)

The following discussion shows that with this modification there is

no loss of accuracy yet that there is a saving in comrutation time.

Geometrieally, the Runge-Kutta formula actually calculates
0 +0 ’1 J. o a l a L.
the derivatives or X(n) at p01nts nh, nn +-§ h, nn + h. The solution
at point (n+l) is then calculated by adding to X(n), the weighted
average of the derivatives multiplied by the time interval h.

The modified procedure given above is exactly equivalent to this

geometric interpretation. Matheaatieally, F [Xl(t), X2(t), t] can be

1

rewritten as F[Xl(t), t] since X is also a time varying function

2

and is known in terms 01 t. Hence

 

Ki = h F] [Xl(n), X2(n), nh] = Kl (3.1-3)
K. 1 h
K' = 11 F1 [Xl(n) + 2, X2(n+ E)’ nh E]
K1
1 h
— h Fl [Xl(n) + 2 , nh + 2]
2: 1%? (3.19)
1
, _ [ K2 h
K3 - h Fl Xl(n) + —E, X2(n+ 3), nh +'§]
= h F [X (n) +-2 ' nh + 1—1]
1 1 2 2’ 2
5 K
3 (3.20)
K11: h Fl [Xl(n) + K'3, X0 (n+l), nh +11]
= ‘ ' gI .2
n Fl [Xl(*1) + K3, nh + h] KLL (3 1)

If the number of equations in the nonlinear and linear
sets are 3 and k, canectively, the computation time by the standard
RungC-Kutta method is h (3+k) units while by the procedure described,

it is (hj+2k) units.

An alternate way of modifying the Range—Kutta method to
realize even a larger time saving is to obtain a numerical solution
to the second set (linear set) by the recurssive formula of section
2, at time h, 2h, 3h ... nh (n+l) h. Solution to the first set is

then obtainer as follows.

X2(n+l), Xm(n) are known, calculate the average

i = % [X (n) + Xq(n+l)l (3.22)

Consider the recursion formula

1 l H H r" 3’ H

1 : 4‘( x ‘ "7- 2 LL: Q.2’D

Xl(1+l) l(*1) + 0 [K1 + K2 + K3 + KM] (a a)
idieina
K: = h Fl [Xl(n), X2(n), th (3.21)
Kg = h Fl [Xl(n) + 2' i, X”, nn + g] (3.23
n _ l — . h 2 A

K3 ~11 F1[Xl(n) + 5 Kg, X2, nn + E] (3.26
Kg = h Fl [Xl(n) + g, X2(n+l), nh + h] (3 27)

{‘0

Equations (3.23 through (3.27) are si ilar to equations

(3.h) through (3. ) exceyt for the following changes

><
+

PJH4
L"
H
0

replaced be2

N
+
rah;
t“
H.

0)

replaced by X2

X0 + L3 is replaced Ly XQ(n+l)
L. J C.
The following discussion shows that the procedure described
ts at at he sum as ice f a 13a 1 as .e sta-dar7 ”e-tztta
111w qJ- leafl’ t'ho’) "I“ fiS“’ ) O" CCLT Cur tlfifi flI.) CL RUIIU FL 1"
method and that there is [gain a saving in computation time. The

rain based on the reometrical inter-

u x.)

irst part of the statement is a

pretation of the Range-Kutta method. The value K as calculated by

l...—I

equation (3.5) is an approximation of the first derivative of X at nh

‘ . 1 _v, 0 _ 0 " ‘I ’ i I. r“* v 0 ‘3 ‘ 0 f w ‘
when K0! K1 are auls1yl1ea by n as ~.nomi 1.. equation (3.c) and (3.7)
(— J

-u w 0 0 - I '3 ~ 0 o l w

anu the result is an approximation to the .crivatives of X at h + g n.
(—
however, the solution is not known at nh + :h but is approximated by
l;

l . - . . . . .
) + — Ll. In the QuOTC procedure this approx1mation lS taken

(

><

PO

nNeJ S

, C‘
8. o

2
[Xq(n) + Xh(n+l)]. This is, in fact a m eh Letter approximation
C. L

J ..L

than the one in the Standard Range-Kusta method.

Mathematically we have

N
n
DJH4

[Xq(n) + X2(n+l)]

-2h-

but from the recursion formula (2.1?)

PJHA

X2(n+l) = a X2(n) + file(n) + 32 e(n + ) + 93 e(n + 1)

>0
n
IbIH

[X2(n) + X2(n) + (a-u) X2(n) + 31 e(n) + [32(n + -:Li)

+ s e(n+1)]

’2)
J

u
>
A
s
V
+
FJH4
rune
\.
O
A
s
V
"L

(«..-uj Mn) +

Therefore

1 . l u
(i—U) X2(n) [U+hA + 2 n A + C L A + 2H h A -U] X2(n)

2 l 2 2 +
= h[A + % h A + E h A + é: h3 A J Xq(n)

Q
L“

<~

l J‘Jl"
t is also to be noted that
(a-U) X201) a L2 (3.31)

Similarly it can be shown that

Xq(n+l) e X2(n) + L? (3.32)

f the numbers of equations in the nonlinear and linear
set are again 3 and h, reSpectively, the computation time by standard
Runge-Kutta method is h(j+k) units. While by the second modified
procedure given above is (hj+k) units, or a saving of 3k units.

3.2 Modified Rance-Kutta Procedure for M'Xed Linear and Nonlinear Svstems.

1.) o

In this section a more general set of mixed linear and
nonlinear equations in the general form is considered. Let the

mathematical model of the system be given as

)

—

 

 

 

 

F Xl Fl (X1, X2, c)
d
— = (3.33)
dt
X2 F2 (x1, X2, t)
_ __J _ _
Where
Fl (X1, X2 t) is the nonlinear vector function and
F2 (X1, X2, t) is a linear vector function of the form
F2 (x1, x2) = A1 x1 + A2 x2 (3.311)

The mathematical models of the vast majority of physical systems are
of this form since only rarely are all components in the system all
linear or all nonlinear. This form of mathematical model might also
arise from nonlinear differential equations of the third or higher order.

This is shown by a simple example -

 

 

 

1n n-l ,n-2

(L C + f r }' (-L ( ‘ "L ] _ r)
‘ un-l ’ 1 n-2 ’ ’ at ’ ’ ‘

Ct at at

7 2 7n-l

1L4". (A. 41 (A. X

rL‘henif --=X,-—-=X ,... = X

dt 1 2 2 n-l n-l

at dt

The higher order equation reduces to the following set of first

order equations

 

 

 

 

r “ . "
X [‘X1
X1 ‘72
9. _
dt ’
Xn-2 Xn-l
Xn_l -r (X, X1, , Xn_0, Xn_l, t)
_ J ._ ._

(3.36)

Using the standard Runge-Kutta procedure given in Appendix C,

we have for the linear equations

d
dt X2 — Al Xl + A2 X

the recursion formula
X2(nh+h) = X2(n) + AX

2

where

IH

AX2 = ”I[Ll + 2 L2 + 2 L3 + Lu]

(\

(3.37)

-27-

and .
L1 = h Al Xl(n) + A2 X2(n) (3.3a
K1 L1
1.2 = h Al [x1(n) + —2- + A2 [ x2(n) + 2—1 (3.39)
K”) L”)
= f ._e ‘ .2 1
L3 h Al [Xl(n) + 2 ] + A2 [X2(n) + 2 ] (3.10)
L1+ = h Al [Xl(n) + K3 1 + A2 [X2(n) + L3] (3.u1)
Substituting (3.38) through (3.hl) into (3.37) gives
_ l l 2 2 1 3 3
AX2 _ h [ U + 2 h A2 + 6 h A£,+'§E h A2] [Ale + A2X2]
= % K2U +-% h A2 +-% h2 A:) A1 Kl + (U +-% h A2) A1 K2
+ A K ] (3-19)
where l 3
X1 = h F [Xl(n), X201), nh] (3.113)
K L
h
K2 = h F [Xl(n) + 2; , X2(n) +-§§ , nh +-§ ] (3.uu)
2 L2 h
K3 = h F [Xl(n) + 5- , X2(n) + 5- , nh + 5 1 (3-45)
Ku = h F [Xl(n) + K3 , X2(n) + L3 , nh + h J (3.u6)

-28-

nspcctflgu of equations (3.h3) through (3.hC) shows

H

that K1’ K0, K are of the same order; and at least for the present,
C

3

will be assumed to be equal for this special case, i.e., K3 *5 K2“! K1 “—"
h F [Xl(n), X2(n), nh] (3.h7). The error at each step can be calculated

in terms of the parameters and can be corrected later if it becomes

too large. With this approximation (3.h2) can be rewritten as

AXE = h a [1“.le(11) + A2X2(n)] + h; AlKl (3.11!)
where
_. _ 1 2 2 1. 3 3
a — [U + h A2 + 6 h A2 + 2h h A2 ] (3.h9)
B:%[hU+hA2+%h2A§] (3-50)

The recurssion formula therefore becomes

X2(n+1) = X2(n) + h a [Ale(n) + A2X2(n)] + h s AlKl (3.51)

The error in equation (3.51) is expressed in terms of K

1}
K2 where KJ’ j = 1,2,3,h are defined as
Kl = h F [Xl(n), X2(n), nh] (3.52)
K 1 h
K2 = h F [Xl(n) + 2— , X2(n) + 2 AXE’ nh + 2 J (3.53)

K0
K3 = h F [Xl(n) + E‘: , X2(n) + E AXE, nh + —] (33-51")
K2+ — h F [Xl(n) + K3, X2(n+l), nh + h 1 (3-55)

and
l n [7.2
Xl(n+l) _ Xl(n) + 6 [K1 + 2K2 + CK3 + K)L J (3.,L)
The error in (3.56) is
1 l . f ,, -
e(n+1) = E h [a U + 7 h AC] [Al] .Kl - KC] (3.57)

If the error ezas given by (3.57) is too large, the correction may

be made by subtracting;e(n+l) from X (n+l) as calculated by 3.56

2
and the calculations in equations (3.52) through (3.57) repeated.
Only one such repetition is normally required and would not even be
necessary for a carefully chosen step size.

The computation time for this procedure is substantually
less than the time for the Runge-Kutta (standard) method. This
eSpeeially is true if the mathematical model has a small number
of nonlinear equations and a large number of linear equations.

Let these numbers be 3 and k respectively. Then the computation
time is (uj+k) units as compared to h\j+k) units by the Standard

Runge-Kutta Method.

I
Lt)
C)
I

3.3 More on the Solution of Mixed Linear and Nonlinear Systems

 

 

 

Another special, but very common class of nonlinear
differential equations encountered in the study of physical systems

is the form

d _ p Q

EX _ .(X) + G(X) (3-5v)
where

F(X) = AX (3 59)
and

G(X) is nonlinear in X containin~ quadraticcn’higher
degree terms in X. The procedure to be develOped for solving the
equations involves first a linear approximation to (3.53) followed
by a change in the parameters to take into account the nonlinearity,
i.e., the nonlinear function G(X) is neglected as a starting point
only. Let the linear approximation to (3.53) be taken as

g? X(t) = AX(t) (3.60)

From the results of section two the solution to (3.60) is

X(t+h) = aX(t) (3.61)

-31-

where

1
a = U + hA + % h2A2 +‘% 3A3 + %E huA* + ... (3.62)

The form of the solution to (3.5;) is then taken as

k»)
C\.
J7

v

X(t+h) = as t+h) (

Equation (3.6M) is obtained by Taylor's series expansion of X(t+h)
which gives:

t+h

X(t+h) = a X(t) + a [‘S G [X(t)] dtJ (3-65)

t

Equation (3.65) then gives the recurssive formula

X(n+l) = a[x(n) + M(n+l)] (3.66)
where

m(t) = ~(G[X(t)] dt (3.67)

Numerical integration of (3.67) is obtained by any of the standard
numerical methods, such as Runge-Kutta method. The scheme then

reduces to the following modified Runge-Kutta Procedure.

Procedure:
A numerical solution to a set of nonlinear differential

equations (3.53) is given by the recurssive formula

X(n+l) = a[x(n) + M(n+l)) (

LA)
c\
._

Where

and

Example:

Then if

Let

rflm

M(n+l) = % [Ml(n+l) + 2M2(n+l) + 2M3(n+l) + Mu(n+l)]

Ml(n+l) = h G [X(n)]

M2(n+l) = h G [X(n) +-— Ml(n+l)]
M3(n+l) = h G [X(n) +._ M2(n+l)]
Mu(n+l) = h G [X(n) + M3(n+l)]

The exact equation for motion of a pendulum is

 

 

Q.
r’.
F0

(3.70)

(3.71)

—4
rs
v

(3.

(3.73)

(3.7M)

(3.75)

(3.76)

-33..

and

 

._ 53
_ - M). -611) (3.77)

The system of first order equations therefore is

 

 

 

 

 

 

 

 

F“ F ' " ' ‘
x1 0 l xl ( O
51. ___ +

X2 _‘k 0 A - X24 _€X1j

or in the notations of this section
- X = AX + G(X) (3.79)

Solution by the above Runge-Kutta procedure is realized by the

recurssive formula

X(n+l) = a[x(n) + M(n+l)] (3.80)

where
M(n+l) = g- [Ml(n+l) + 2l«1,)(n+l) + 2143(n+l) + Mu(n+l)]
and o — (3.:31)
Ml(n+l) = h 3
9% Xi(n) (3004)
..J
G
Mg(n+l) = h (3.33)
gulp.) + 7:- 142(n+1)1i

 

 

-3h-

 

 

?" O '7
M (n+l) = h (3-34)
3 k 1 3
j[xl(n) + §M2(n+l)] J
0 7
r.1u(n+1) = h (3.35)

  

_[ )1 (n) + M3(I1+l)]1

1
Z 1

 

 

NONLINEAR ALGEBRAIC EQUATIONS

IV

contain components modeled in terms of non-
possible to develop

possible to

When systems

linear algebraic equations it is not always
may not be

from the system model.

tem model in normal form, i.e., it

a sys

model, in this case is referred to as being in
seminormal form and for the purposes of our discussion is represented

eliminate the nonlinear algebraic equations

A

The mathematical

 

 

 

 

in the form
X F (X,Y, t)
d :
dt (u.1)
o G (X,Y, t)
t _ L _
Where G(X,Y, t) is a vector function nonlinear in X and Y or at
least nonlinear in Y; so that Y cannot be expressed explicitely in
If G(X,Y, t) is linear in Y, the vector Y can
first equation and the model reduced to normal

3A)
143

terms of X and t.

be eliminated from t
Many attempts have been made to obtain a solution to the
None of the procedures, however, is general.
h The

form.

nonalgebraic equations.

Each method is applicable only to a certain class of proslems.
procedures to be used then depends on the nature of the nonlinearities

involved in a given problem.

-36-

1.1 Method 2; Differentiation

 

(1)

One method of solution related to recent work by Wirth
is to transform the nonlinear algebraic equations into differential

equations by partial differentiation. The procedure is as follows:

The algebraic equation in (h.l) viz.

G(X,Y, t) = o

is differented with respect to t, so that

d DG d< 3G dY BO
__ t :: —— —— + "'_ """ = °
dt G(X’Y’ ) -ax at -sr dt + ‘at -» O (A 2)
’BG 'bG . _ . .1
Where_;§- and 3?. IS the Jacobian of the vector G with respect to

the vectors Y and X reSpectively.

For example consider the set of equations

—— — — a

 

 

xl i fl (x1, x2, x3, yl, yE, t )
‘2 f2 (X1’ Ae’ X3’ y1’ ye’ t )
9'... x = f (x,x,,x,y,vn,t) (1+.3)
at 3 3 1 a 3 1 c
O g1 (x1, x2, x3, yl, y2, t )
L O \ [— 32 ( X1: X2; X3; yl, ye; t )

 

 

The form of equation (4.2) for this system of equations is

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

— T T— fl _
’3 X1 3X2 3 X1 (1‘0 tyl NE. 71. B t
dxn
-—-“_(': + + :: 0
(it
BX 9339) ‘34:“ LL11 'B y ‘DZY r: at 2 L
l (.1 :5 L L'. J

 

 

G . . .
If the matrixggg is nonSingular, then the solution (h.h) for the

u o EY I
derivative vector 5; is

 

El _ _ [Ifl [19 £35 + 3.5.3 1
3t BY 3X ch; 3t
_ tG-l gfi 3G
.. - [a] [5X am. t) + M 1
= P (X. Y. t) (11.5)
.. . . ?G -l . L . .
The conditions under which [-—J exists are established in

?Y

reference (1). The initial conditions on Y must satisfy the nonlinear

algebraic equation

iy considering (h.5) in conjunction with (h.l) the mixed system
of nonlinear algebraic and differentail equations are transformed

into the set of nonlinear differential equations in normal form

fo,r, t)

me

 

.—

 

l P(X,Y, t)
obtain

 

 

‘ections 2 and 3 can be used to

fa

The numerical procedures of

a solution.

 

 

 

 

 

Example:
Let the mathematical model of a system be given as
_ 1 —' '- P '1
X f(x, y, t) ax + by
i z : (1+.S)
dt
O g(x, y, t) x - sin y + six mt
t we have

 

Taking the time derivative of g with respect to

dg dx dy

'7— = —" - CO" ‘\ , + 5'

at dt ° J at ‘” 00°‘Dt

or
-E¥ = (ax + by + a>cos wt / cos v
LI
)
(4-9)

t)

_ P (X: 31')

L0

The nonlinear differential equations to be solved then are

 

 

 

 

[— _ '1
x W ax + by
i
dt .__ (1+ . 10)
y p( 1-1, 3'; t)
L .. _ _
with initial conditions satisfying the equation
f [x(o), y(o), O] = x(o) - sin y(o) + sin (0) = O (h.ll)

n . .
y o) = /C is such a value so that

Therefore, the initial conditions for (h.lO) are

x(o) = 0.5 and y(o) = n/E.

h.2 fin Alternate Procedure

 

In many a problem the technique as indicated in section
h.l may not be applicable for the reason that Jacobiandgg may not
exist or is singular. Even if the Jacobian is nonsingular the final
mathematical fo~m may be very complex and, as such, is difficult

handle numerically. The following numerical procedure overcomes this

difficulty.

-ho-

Again consider the mathematical model

— qr ..—

X F (X) Y) t) W

3— = (u.1i)
O G (X) Y: t)

 

 

 

 

The basic procedure involves an estimate of X(n+l) from values

of X(n) and X(n). Before considering the details of how this estimate
is to be made let it be assumed for the present that an estimated

value of X(n+l) is available. The problem then reduces to the solution

of the set of nonlinear algebraic equations

G [X(n+l), Y(n+l), nh+h] = o (3.12)

for Y(n+l). The Newton-Raphson can perhaps be applied to obtain
this solution using Y(n) as starting values. However, if the true
solution for the vector Y(n+l) is considerably different from

Y(n) this method may not converge, or it may converge very slowly,
thereby requiring an excessive number of iterations to arrive at a
solution. In an attempt to avoid this difficulty an estimate of
Y(n+l) is made using quadratic extrapolation. Using Y(n+l) to

represent this estimate we have

g<n+1> Y(n) + [Y(n) - Y(n-l)]+-% tr<n) - 2Y<n> + Y(n-2)1

g Y(n) - 2 Y(n-l) +-l Y(n-2) (3.13)

2

-111...

Then

G [X(n+l), Y(n+l), nh+h] = G [X(n+l),'f(n+l), nh+h]

r: [Y(n+l) - Y(mlHE-g- [X(n+i), Y(n+l), nh+h]

and as a first degree approximation

3
B

Q

Y(n+l) = Y(n+l) - [ (n+1)]‘l'c‘(n+1) (2.42.)

+<

This last expression is precisely the Newton-Rarhson formula. The
modification that Y(n+l) as evaluated by (b.13) is used instead of
Y(n). On the basis of a typical example it appears that only one
such calculation is necessary for convergence. Although it may first
appear from the form of (n.1u) that a first order approximation is
used, it is actually a third order approximation,

since Y(n+l) has been estimated by quadratic extrapolation.

The disadvantage of this method appears to be in the evaluation of

. ‘30 -l
the inverse [3Y1

at every step. However it may be possible to
evaluate this inverse analytically.
The estimate of X(n+l) required in the above procedure can be

established by any of the predictor type recursion formulas, but

Milne's formula seems to be more apprOpriate.

-hg-

Milne's formula is

£§,[2 F(n-2) — F(n-l) + 2 F(n)] (5-15)

X(n+l) = X(n-3) +

After the vector Y(n+l) is calculated by (b.1h) the esti-

mated veetor X(n+l) can be corrected by any of the integration
use estimated solutions. Simpson's formula is one

formulas that

such formula and gives

[X'(n+l) + M X'(n) + X'(n-l)]

bull—3

X(n-l) +

X(n+l)

f§(n+1) + u F(n) + F(n-l) 1 (3.16)

wlta

x(n-i) +

It can be note1 that the procedure developed in this section
requires solutions at the preceeding step . This is a handicap for

t can be overcome by using a smallar interval

the first few steps.
for the first few steps and by use of lower order methods that do
not require the solutions at the preceeding steps. For example let

F x-l Ff (x,y,t)T
(b.17)

..I

. v 1 4-
8 (159.15“)

 

 

 

 

-h3-

Tue first three steps x(l), x(2), x(3) are obtained from the first

~..l

. . ., . .l
dexree approximation Wltn a temporary (local) step size 01 :h, taus,
L;

Mh=xe>+gfimw.mw.m

 

(kit)
1 '1
l L3[X(:§)) 3(0)) '5']
;P(:) = Lf(C) - .3 l h
:3? £3[X(f_5): D/(O): '2‘]
mn=xd>+§amp.xp.9
(It-l9)

 

1.:
A
}_J
V
ll

:4
A
RNPJ
V
l

V. EXAMPLE SOLUTIONS

In this section the computer results for several simple
cases have been given for the purpose of comparing computation times

1

with some of the other methoas.

 

 

5.1 Numerical Solution of 3 Linear System in Tormal Form

—" —1 _ -l 0 q - V l — l 1

“l l

 

(i
C
1...:
I
DJ
C)

 

 

 

 

 

 

 

 

Using the procedure of section 2.1

2
Q
A
:3
v
'C)

xl(n+l) all ll 1

 

x2(n+l) 621 a x0(n) 801

PO
PO

 

 

 

 

 

 

 

vhe re
1 2 l 3 l h
f’ = l - h,t — A - h +-—r h
Vall 2 Z 24
a = 0
l2
’3. 2 2: 3 5 I;
an = h - i h + 4 h - 7 n
a1 2 2 O
'3 r) A
1 J L
022 — l - 2 h + 2 n - - h +‘E h

an

The numerical results for the above linear system as obtained by
three different methods are given in Table (5.1) along with the
computation time. Calculations were made with a time increment of

0.001 units, with results p‘inted out every.10th interval.

-16-

TABLE 5.1

SOLUTION OF LINEAR SYSTEM OF EQUATION

Recurrsive Method of Section 2.1

Computation time;

Time Units x1 x2
0.00 2.00000 3.00000
0.10 1.90u8u 2.63293
0.20 1.81873 2.32u21
0.30 1.7u082 2.06u0u
0.u0 1.67032 1.8uu31
0.50 1.60653 1.65835
0.60 1.5u881 1.50060
0.70 1.u9659 1.366u8
0.80 1.uu933 1.25217
0.90 1.u0657 1.15u32
1.00 1.36788 1 07088

Computation time: 2 minutes, #5 seconds
Predict-Correct Method of Appendix B

Time Units x1 x2
0.00 2.00000 3.00000
0.10 1.9ousu 2.63293
0.20 1.81873 2.32u21
0.30 1.7u082 2.06u0u
0.10 1.67032 1.8uu31
0.50 1.60653 1.65835
0.60 1.54881 1.50060
0.70 1.u9659 1.366u8
0.80 1.uu933 1.25217
0.90 1.u0657 1.15132
1.00 1.36788 1.07088

8 minutes 5 seconds

TABLE 5.1 Cont'd.

Runse-Kutta Method

Time Units x1 :2
0.00 2.00000 3.00000
0.10 1.9ou3u 2.63293
0.20 1.81873 2.32u21
0.30 1.7h082 2.06h0u
0.u0 1.67032 1.8uh31
0.50 1.60653 1.65836
0.60 1.5133 1.50061
0 {0 1.u9658 1 36650
0.30 1.uu931 1.25219
0.90 1.h0655 1.15u3u
1.00 1.36785 1.07091

Computation time: 8 minutes 35 seconds

5.2 Numerical solution to mixed linear and nonlinear system-closed
linear set.

The mathematical model under consideration is

{-x {-x x I ' 0 T

l 1 2

d
3:; X2 = X2 + l (5.3)
J- _ 1 L _

The numerical solution as obtained by the two procedures develOped in

 

 

 

 

 

L1.

section 3.1 and by the Rungc-Kutta method along with the computation time

for each is given in Table (5.2). The calculations were made with a time

k

1

increment of 0.001 units. The results were printed out at every 10th

step with every 10th printed result being given in Table (5.2).

Time Units

.00
.Ol
.02

00000000000
0
\D

-18-

TABLE 5.2

SOLUTION OF MIXED LINEAR AND NONLINEAR EQUATIONS

Procedure 1

+ + + + + + + + + + +

0\ F‘ F1 F4 F4 FJ F4 F4 F4 21 F1 F4

Computation time:

Time Units

0.00
0.01
0.02
0.03
0.0u
0.05
0.06
0.07
0.08
0.09
0.10

x1

.00000
.01015
.02061
.03139
.0u250
.0539u
.0657h
.07790
.090u3
.1033u
.11665

Procedure 2

+ + + + + + + + + + +

Computation time:

x1

1.00000
1.01015
1.02061
1.03139
1.04250
1.05395
1.0657u
1.07790
1.09013
1.10335
1.1666

WITH CLOSED LINEAR SET

+ + + + + + + + + + +

+ + + + + + + + + +

+

x2

.00000
.02010
.ououo
.06091
.08162
.1025u
.12367
.1u502
.16657
.18835
.2.031

F4 r4 F” F’ F4 F4 F” F’ F’ F’ F’

minutes, 50 seconds

x2

.00000
.02010
.0u040
.06091
.08162
.1025u
.12367
.0u502
1.16657
1.18835
1.2103u

F4 24 F4 F4 F4 F“ F4 F4

6 minutes, 5 seconds

C OUT; i UPC (1

T ,...-..\ ' i4. f 41‘” ‘_
Run e-hutta hetnoe

L. .
x.)

Time Units

3.0 +
0.01 +
0.32 +
0.33 +
0.05 +
8.35 +
0.0C F
0.87 +

Computation time:

‘1

1.00000
1.0101:
1.02061
1.03139
1.0h250
1.05395
1.0557h
1.07790
1.090u3
1.10335

1.11666

.rA
ALL-

1.33330

'\"‘, /
1- 013:-U_K/

1.12h67
1.1u502

1.16657

’\ . , -1
a m1nutes 15 seconds.

C. n
— H .-
/\)

 

_o [w] .‘ -._-‘ p- ’\ ..L o
Consiua: LUC sec or e;uacions
P— ? — -'

9— = (3.1+)

:\
{\D
l
H
I
k.)
D.)

 

 

 

 

4H, 1. J. .3 . .. f‘ C‘ ,. J. 3 , ’3 I‘
In a“; natations 01 Section 3.2

+ sin 100t (5.5)

 

thJ
b'
5.,
,2
+
('\It—’
a.
,_)
i
+
:3“
p.
c
I—

¥~—J
p“)
t")
1H
,2
k)
L—J

rake
A
\J‘n
_\]
v

T
(\h4

I
n
Cflh‘
I:
C}
+
3‘
rUD
+
#424
5’
3::

P.)

u
H
f'_|
p.
I
H4
5 J
+
£424
5...:
._1
L.—
A
\11
L)
v

‘ -4. c \ I ' . r‘ N r‘ J a -. -1*\/q ... ‘ ,1 I a.,>\ 7 -.
The numeiical scldtions as attaineu by the piocecuie

of Section 3.2 an‘ by the Range—Kutta scheme are given in Table 5.3.1

f\

along with computation times. The time increment used is 0.0001
units and results are printed out at every 10th step; with every

10t1 printed results given in Table 5.3.

SOLUTION OF A CLASS OF MIXED LINEAR AND NONLINEAR SYSTEM

-52-

MODIFIED RUNCE-KUTTA PROCEDURE

Time Units x1 42
0.01 0.000000 1.000000
0.01 0.0103u73 0.5;-999
0.02 0.0213 33 0.579991
0.03 0.032970u 0 969569
0.0a 0.0u50701 0.959930
0.05 0.0576121 0.9h9 67
0.06 0.07002u7 0. 397 0
0.07 0.0023371 0.929671
0.03 0.09h5375 0.9195u0
0.09 0.106uyo .0.;:;3;0
0.10 0.11.013 0.095221

Computation time: 10 minutes, 15 seconds

tandard Runge-Kutta Procedure

Time Units x1 x2
0.00 0.000000 1.000000
0.01 0 0103u73 0 981999
0.02 0 0213633 0 97,991
0.03 (1032970h 0 963970
0.0M 0.0u50701 0.950030
0.05 0.576122 0.9u9363
0.06 0.07002h7 0.939 32
0.07 0.0 23375 0.929672
0.08 0.09h5376 0.9195u2
0.09 0.1065 0 0.909392
0.10 0.113013 0. 99226

Computation time: 12 minutes, 55 seconds

5.1.

Mixed Linear and Nonlinear Mathematical Model Given in Equation (3.3)

 

 

The mathematical model under consideration is

 

 

 

 

 

 

 

 

 

 

 

 

I j _— 0 l — " _
x1 Al
d __
:— — (5.9)
at 5 3
XL —10 0 KA -: ‘
2 2 3 l
L _, L _ L J
Dy notations of section 3.3;
1"
0 N F0 T
A = and~G(x) = (5-10)
'2
-10 0 EX:
J .
L _ 1- _
Then
- 2 25 u
all — l - 511 + -—6 h
= 2 3 5 .5
C212 h 3 h + Z I1
= _ 2 3 2 5
0:21 10 (h 3 h + 6 h )
_ .2 2.2 1*
22? - l - 5h + 6 h

-5h-

The numerical solution to (5.9) by N dified Runge-Kutta
Method and Standa‘d Runge-Kutta Method is given in Table 5.A. The

time increment used in 0.001 time units and results are irinted out
T

at every 10th step; with every 10th printed result given in Table 5.A.

-55..

TABLE 5.4
SOLUTIONS OF ANOTHER CLASS OF LINEAR AND NONLINEAR SYSTEMS

Modified Range-Kutta Procedure of Section 3.3

Time Units x1 x2
0.00 + 0.500000 0.000000
0.10 + 0.h76216 - 0.u72181
0.20 + 0.h06958 - 0.902567
0.30 + 0.298456 - 1.250510
0.u0 + 0.160822 - 1.h79810
0.50 + 0.007335 - 1.56hu10
0.60 - 0.1u6878 - 1.h9u180
0.70 - 0.286602 - 1.277550
0.80 - 0.398322 - 0.93930h
0.90 - 0.u71588 - 0.51u959
1.00 - 0.u99789 - 0.04u926

Computation time: 7 minutes, 10 seconds

Standard Runge-Kutta Procedure

Time Units x1 x2
0.00 + 0.500000 0.000000
0.10 + 0.h76216 - 0.h72181
0.20 + 0.h06958 - 0.902567
0.30 + 0.298u56 - 1.250510
0.u0 + 0.160822 - 1.h79810
0.50 + 0.007335 - 1.56uh10
0.60 - 0.1u6878 - 1.191180
0.70 - 0.286602 - 1.277550
0.80 - 0.398322 - 0.939303
0.90 - 0.u71589 - 0.51u958
1.00 - 0.h99790 - 0.0uu925

Computation time: 10 minutes, 50 seconds

C‘
Q

.4 Viz-\"fi
‘sli_)_,|._L

7—1.
(J)

 

considered'

The mathematical model

 

 

 

 

 

 

 

x -x + cos y
d 2: .
73' 6'
U
L,0 L L - sin y
/
x(o) = 0.5, y(o) = 0.5236
According to the method of tifferentiation 5.11) reduces to
- 1 - 1
x -x + cos y
d
dt =
y (-x+cos y)/ cos y
(5.12)

 

“\r
J

-I{ + COS

 

 

v
1 - ("/60S y)
x(o) = 0.5. y(o) = 0.5235
The solutions of (5.11) by the two procedures of section A,
' The time inc ement
t printed out at ever;

.1th every 10th resal

1

w w ,-
.JCKJL

along with the computation time are given in Table 5.5.

o
3-8 L.

of 0.001 units
result and is given in Table 5.5.

10th

-57..

TABLE 5.5
SOLUTION 0F SYSTEM MODEL WITH NONLINEAR AIGEERAIC EQUATIONS

Method of Differentiation

Time x y

0.00 + 0.500000 + 0.523600
0.10 0.533832 0.563127
0.20 0.562574 0.597498
0.30 0.586874 0.627193
0.40 0.607325 0.652691
0.50 0.624470 0.674454
0.60 0.638793 0.692929
0.70 0.650721 0.708534
0.80 0.660628 0.721656
0.90 0.668837 0.732645
1.00 0.675627 0.741816

Computation time:

P.
S
(D

.00
.10
.20
.30
.40
.50
.60
.70
.80
.90
.00

i—‘OOOOOOOOOO

Computation Time:

Alternate Procedure

X

0.500000
0.533832
0.562575
0.586874
0.607325
0.624470
0.638793
0.650721
0.660628
0.668838
0.675627

12 minutes, 30 seconds

y

0.523600
0.563126
0-597597
0.627192
0.652689
0.674453
0.692928
0°708533
0.721655
0.732644
0.741815

7 minutes, 5 seconds

VI. C NSLUSION

Sevgral procedures for numerical solutions of
mathematical models are given in this thesis. In section two, a
rec1rssfvc iorxula is developed for linear syssems. The computation
time by this 101‘ .ual is one fourth of tl'lat 12y the Runjje-Kutta Izzethod.
Moreaver, any ~;iren dsgree of accuracy can 1e maintained without
increasing the com utntion time per step. In the case of linear
systems, the expression for the nth stage solution vector in terms
of the initial conditions vector and input veceors has a strik'ng
similarity to the mathematical model of :iscrete state Markov
Process. Thus the "Dynamic Programming" techniques eveloped by

Howard and others for Optimizing the systems also apply to the

control systems.

In section 3, mixed linear and nonlinear systems are
considered. The procedures level ped for the various tyres of linear
and nonlinear combinations are essentially: nadifi ed unge-Kutta
Procedures, in which the computation time has been reduced considerably.
This saving in time is achieved by applying the results of section 2
to the subsets in the system.

In section A, systems containing nonlinear algebraic
eguatior s are considered Two procedu‘es are given for the solution

q

of such a system. The Inc hod of differentiation as given in section 4.1

.J-‘ a .:_, .1.
0110 Eli. oCI‘H/fl CC

4..

gives better results, *ut is relatively slower as compared .0

r~ocedure of section 4.2. Bot1 the procedures, however, seem to

4'.

b4

converge rapidly.

t_.
.l.

\n
e L)

-59..

In section 5, severa1 examples of solutions by the
procedures developed in this thesis are given and compared with
I. I L) ..

the solutions ottaincd by several other methocs. These comparisons,

indeed, are consistent with the theoretical results.

F0

(:0

10.

My 'RENCES

 

WIRTH J.L., Time Domain Models of Physical Systems and Existence
of Solutions, Technical Report No. 1, NSF, G-20949, dichigan

State University, East Lansing, 1962.

EILNE W.E., Numerical Solution of Differential Equations, John

Wiley and Sons Inc., 1953.
HOUSEHOLDER, A.S., Principles of Numerical Analysis, McGraw Hill, 1953
HARTREE, D.R., Numerical Analysis, OxfoniUniversity Press, 1953.

IUUHING, R.W., Numerical Methods for Scientists and Engineers,

McGraw 2111, 1962.

HOWARD, R.A., Dynamic Programming and Markov Processes, The

Technology Press and John Wiley and Sons Inc. 1960.

BELLMAN R., Stability Theory of Differential Equations, McGraw Hill,
1953.

BUCKINGHAM, R.A., Numerical Methods, Pitman, 1957.

KUNZ, K.S., Numerical Analysis, McGraw Kill, 1957.

BENNETT A.A., Milne W.E and Bateman N., Numerical Integration of

Differential Equations, Dover Publications, 1956.

Recurssive formula

C‘. .1- ° q
030 Section ¢.l

m1 - 1 a .7 _
iuira aegree

X(n+l) = [6,

WEE l'C

f—l
t
5...;
u
*7
f‘fi
our)
c:

LolH

h[—}

7—!

.1)

El
II

J
ll
"fi‘
P
fl

F0
H
[C

T)
u
5’
PHP4
c:

c:
+

:3
;)
+

PDH4

Fl
r—w
>>

U+—E-m+-1-:1;A‘+;—,-
3 1~ t

Appendix A

 

rees of accuracy.

X(n) + 91 E(n) + Qc E(u + %) + 9‘ E(n+l)

i
l— J

2 1:3

L. _.L

(A.

C\
v

-a
V

 

 

-LQ-
A.3 Fifth degree accuracy.
X(n+l) = a X(n) + 91 E(n) + 32 E(n+ %) + a? E(n+l) (A.ll)
wgeie
r3 , ’1 :1 i ‘ r r
[V] = [U + LA + % nwAh + l hJA“ + %r nknk F E, h A ] (A.12)
C. 1 L“' 1(4'3
fix -— 2]... .3... ..L 7 1: 2 .1.— 13n-3 1 all. 2}- . '3
[,1] h [ 1* U + 10 nA “'35 L A + pi i n 4 1,3 1 A 1 (A.13)
r) 2‘
[:32] = 11[- “i? U + %:"5 11A + £3 11 A + EA 1173.3 ] (Anlh)
[ ] [ 3 U 2 hA l h2 2
___. ..;._ —_ —— r-
83 h 15 + 15 + 30 A 1 (A 1;)
D.h Sixth degree accuracy.
X(n+l) =13 X(n) + 51 E(n) + sq E(n + %) + a“ E(n+l) (A.1C
. 1 .2,2 1 .3 3 1 .2 h 1 S 1 .6 C
[T] = [ U + nA + E n n + E n A + 2n n A + 125 h5A + 725 n A ]
(A.17)
19 11 . 1 2 2 11 ,3 3 1 .h A 1 5 5
" = —_,. —} ‘n ~9— J —— 1 » —— 1 /
[31] h [3: U + 36 nA + 9 1 A + 3:0 n 1 + 1,0 1 A + 750 1 A ]
(A.13)
O 2 j ’1 \
[.32] = h [-J- U + ;1L_) ‘A + %‘Z all + 1—1—5 hit) + .373. JAM] (£1.19)
. r 2 2 “
[:33] : h[%%U+éZ-}1A+%—CIA 4-55.}1J1A3] (11.20)

APPENDIX B

 

D.l Runge Kutta Hath d.
d ‘ a
Let E? X(t) = F [X(t), t] (2.1)

Then X(n+l) is calculated from X(n) and 3.1 hy the recursion equation

-- .1; C“, ’3 A.
X(n+l) _ X(n) + 6 [K1 + 2K2 + 2K3 + KM ] (2.2)
the re

Kl = hF [X(n), nh] (B 3)

K = hF [X(n) +i K ) nh +-3 1 (B u)
2 2 2 ’ 2 °
. l h 2

K0 = “F [x(n) +-; KG), nh +-;] (n.5)

J l- L. (-
K1+ = hF [X(n) + K3), nh + h 1 (3.6)
J

For a linear system the numerical solution Ly the Runge-Kutta

.1

method is equivalent to the fourth order Taylor series expansion when
the function F is not a direct function of the independent variable t.

When, F is a direct function of t, the truncation error factor is

much higher than the exact expansion.

C3

-cu-

B.2 A Predictor - Corrector Method.

.1.

A very commonly used procedure for the numerical solution b0

linear differential equations is the predictor-corrector method.
There are many variations of this method, but the following one sug—

gested by Milne and modified by Hamming is used as a comparison with

the procedure of section 2.1, and so is included here.

The details of the general prector-corrector m thods can be found in

Reference 5.

Let the differential equation be

d
Ex(t) : £0913)

Then
x(m) = X(n) £23 [1,. (n+1) + X'(n)] + 131—; [X"(n) - p" (n+1) 1
(2.3)
where
P(n+1) = X(n-l) + 211 x: (n+1) + 112 [éi X"(n) + g x" (n-1)] (3.9)

W1 USE e311