IN,. 0
TECHNICAL REPORT T7866
U.S. ARMY
TRAPEZOIDAL CONVOLUTION REVISITED: RK CONVOLUTION OR THE DIGITAL SIMULATION
MISSILEOF CONTINUOUS
RESEARCH
DEVELOPMTE(IT COTIMiA fD C3
SYSTEMS VIA ZTRANSFORMS
Richard E. DicksonDictre
Technology Laboratory
15 June 1979
LLj
Gm~
Approved for Public Release,
d stone Arsenal, Alabama 35809
s.79 MVIFOM
000. 1 APA ?
Distribution Unlimited.
10 12 012
DISPOSITION INSTRUCTIONS DESTROY THIS REPORT WHEN IT IS NO LONGER NEEDED. DO NOT RETURN IT TO THE ORIGINATOR.
DISCLAIMER THE FINDINGS IN THIS REPORT ARE NOT TO BE CONSTRUED AS AN OFFICIAL DEPARTMENT OF THE ARMY POSITION UNLESS SO DESIGNATED BY OTHER AUTHORIZED DOCUMENTS.
TRADE NAMES USE OF TRADE NAMES OR MANUFACTURERS IN THIS REPORT DOES NOT CONSTITUTE AN OFFICIAL INDORSEMENT OR APPROVAL OF THE USE OF SUCH COMMERCIAL HARDWARE OR SQFTWARE.
. . ...•
, .. ..
......
.
..
.
.
.
.
.
.
..
.
.....
..
.
.
,
.......
i
Unclassified SECURITY CLASSIFICATION OF TNIS PAGE ffffe DM.
Mtovor__________________
REZAD INSINUCflONS
REPORT DOCMENTATION PAGE I.
BEFORE COMPLETING FORM GOVT ACCESSION NO. S. RECIPICNT'S CATALOG MUR
REPOT HUNG.
TYPE OF REPORT a PERIOD COVERED
eveied
4. TITLE fand Subeft.I
I Technical Repeto
RK Convolution or the Digital Simulation *~~of Continuous Systems Via ZTransf arms I .3
1.
(
a CO
AUTOR(.)S
TRAC
OTA0
Richard E. Dickson1. S. PERFORMING ORGANIZATION NAME AND ADDRESS
PORAM ELEMENT.PROJECT. TASK AREA A WORK UNIT NUMBERS
I*.
$
Commande" US Army Missile Research and Development Coimmand ATTN: DRDMITD Redstone Arsenal, Alabama 35809
/ 49A4 ar um*7
16
I I. CONTROLLING OFFICE NAME AND ADDRESS
Commander US Army Missile Research and Development Comman('
7
Z1
DRDMITI"UDROPAE
ATTN:
73_____________
RedstoneArsenLlabama_35809 AGENCY NAME S1AODRESS(II dElgorenI Dim CoritrolIng Office) Tr. MONI1TORINGO
IS. SECURITY CLASS. (of Iftl. MOON)
Unclassified ISM.
16.
DECL JUSIFICATION/ DOWNGRADING SC,*EOGULE
DISTRIBUTION STATEMENT (of this RoE)
Approved for Public Release; Distribution Unlimited.
IT. DISTRIBUTION STATEMENT (of the abstract eitomd In Sleek 20.II fflevt
1S.
SUPPLEMENTARY NOTES
it.
KCEY
WORDS (Cmntinu
on reverse aide It noe...inv and Idomit?~ by mock inbot)
Digital Simulation Integration Initial Conditions
.rZtransforms RagazziniZadeb Identity Convolution A _5
ACT
(MNI
heom Report)
M e.,Ov IND N 0=908
=9 0000Pt
loek monkc)
A proof of the RagazziniZadeb identity which uses the properties of the Dirac delta is outlined. Then the mean value theorem of the integral calculus and various numerical integrators are used to develop convolution approximations of use in the digital simulation of continuous systems. The accuracy of the approximations is examined for the single integrator; also, there are sample problems showing the incorporation of initial conditions.
IDD I ism~ VS
3
EIIO
FMV
~~,
IOSL
Unclassified S~5CUMTV CLASSFICATION OF roIS PANE (111M
=ore .
1.4.
Section
CO TN4Page
1. Introduction............................................
5
2. Some Preliminaries .....................................
5
3. Mean Value Convolution.................................
9
4. RK Convolution........................................
14
5. Analysis of Errors......................................
23
6. Some Sample Problems ..................................
34
A. B. C. D.
Single Integration .................................. Double Integration .................................. Single Pole Filter ................................... A Forced Damped Oscillator..........................
7. Results and Conclusions ................................
34 38 41 47 5
Appendix A. Transform Table for Selected Functions and Distributions .......................................
51
Appendix B. Dirac Delta Distribution......................
5
Appendix C. Numerical Differentiation ....................
57
Appendix D. RK(N) Convolution..........................
69
Reference .............................................
72
TABLES Page
Table 1. Single Integrators for Various Order Holds .................... 2.
Convolution Approximations for Z[g(s)f(s)] .....................
16 24
3. Checks of the Convolution Approximations for f(s) = 1/s and g(t) a Constant, Ramp or Parabola .............. 24 4. Ratios of Approximate Solution to the Exact for a Constant, Ramp or Parabolic Input to a Single Integrator .....
25
5. Accuracy for an Exponential Input to an Integrator ..........
30
6. Accuracy for a Sine Wave Input to an Integrator .............. 31 7. Accuracy for a Sine Wave Input to a Tunable Integrator ....... 32 8.
Accuracy at the Shannon Limit (wT =.7r)
.........................
33
9. Convolution Approximations for Single Integration ............ 35 10. Recurrences for Single Integration ............................. 38 11. Convolution Approximations for Double Integration ........... 40 12. Recurrences for Double Integration ............................ 41 13. Convolution Approximations for a Single Pole Filter ........... 43
2
TABLES (CONCLUDED) 14.
Recurrences for a Single Pole Filter ............................ 43
15.
Recurrences for a LeadLag Filter
(g±b 16.
................................................
4
Recurrences for a LeadLag Filter
+ ./a
.................................................... 45
Recurrence f r( S
a) ........................................ 46
18.
Startup Step(s) for Oscillator .................................. 48
19.
General Recurrences for Oscillator .......................... 48
CI
Integrators Via Log Approximation ......................... 60
C2
Tunings Independent of "wT"............................... 63
C3
Differentiation of a Cosine Wave, "Amplitude/Phase" ........ 64
C4 C5
Differentiation of a Cosine Wave ........................... 65 Differentiation of a Cosine Wave,wT<
C6
Differentiation of a Cosine Wave, wT<«< .................... 67
C7
Differentiation of a Cosine Wave, wT «r .................... 68
3
fte
ACKNOWLEDGMENT
The author wishes to acknowledge the debt he owes Professor Charles A. Halijak for his briefings on "Numerical Transforms and Digital Simulation of Dynamical Systems" and "Advanced ZTransforms" held at Redstone Arsenal during the summers of 1976 and 1977, respectively. His introduction to the power of trapezoidal convolution and instilling of a questioning attitude toward the Ztransform literature (including Halijak)
have stood the author in good stead.
4
1. INTRODUCTION A decade ago some facts and fallacies in digital simulation were discussed in an interesting article [ I]. Some of the facts appear to be fallacies as that author had warned might be the case. The fallacies arose, in part, because of premature use of the general recurrence; that is, improper incorporation of the initial conditions. It is, of course, unfair to "sharp shoot" after a decade, but some of the "facts" (fallacies) are still present in today's ztransform literature. These and other questions were addressed in an earlier report (21 documenting the results of an initial investigation. Several suggestions were made as to avenues for further possible development. They were as follows: a)
Place ztransforms on a firm foundation using distribution theory,
b) Determine when modified ztransforms and/ or tunable convolution are advantageous and in what combination. This would include the use of higher order holds. c)
Analyze the effects of tuning for other inputs and other transfer functions.
Some progress has been made, and, though in no way complete at this time, it was felt that the results to date were worth reporting if only to stimulate further interest. There is the story of a physicist and a mathematician in an airplane which flew over a flock of sheep, all white save one, who was black. The physicist proceeded to theorize about the number of black sheep in the universe. The mathematician knew there was one flock with a sheep who was black ON TOP. An engineer would probably choose to ignore the color of the sheep since it would have no effect on the mutton. This author is a physicist.
2. SOME PRELIMINARIES In this section certain details as to the notation used and some conjectures concerning the use of distributions in ztransforms will be presented. 5
That the Laplace transform of the Dirac delta* is a constant and the Laplace transform of the nth derivative of the Dirac delta Iss"(see Appendixes A and B) would certainly seem to indicate an intimate connection between Laplace transforms and distribution theory.
For reasons which will become apparent, the Laplace transform will be defined as (31, L(t(t))
 f f(t) u(t) est dt
(1)
where u(t) is the Heaviside unit step (whose derivative is the Dirac delta). Of course Equation (1) readily reduces to 00
f(s) = f f(t) eStdt
(2)
Note that the introduction of the unit step, Equation (1), remov i the need for defining the behavior of "f(t)" prior to "t" equal zero. The lower lim. of integration in Equation (2) is contained in the unit step in Equation (1). For example, consider convolution: O.D
g(t) * f(t)= f g(t)
u(t) f(it) u(1t) dt,
(3)
f(t)
(4)
Co
=
f
g(t)
u(t) u(rt) dt.
The term, "u(t) u(rt)," will be recognized as the unit pulse from 0 to r. Equation (4) becomes "1
g(t) * f(t)  f g(t) f(tt) dt. 0
(5) *The )irac delta distribution is not a function. Dirac to distinguish it from the Kroneker delta.
6
It should be noted that the Dirac delta is the unity element in convolution,
go.(6)
g(t) 6(1t) dt 
For a discussion of the convolution of distributions, see Kecs and Teodorescu (4]. The Dirac delta may also be thought of as a sampler
fg(t)
6(nTt) dt
=
(7)
g(nT)
If one has the sequence 151
f Wt
(8)
f(t) 5(nT't)
n0
its Laplace transform is
n0
Taking
z
(10)
eS
one has
n0
4
7
The ztransform may be thought of as a discrete Laplace transform, that is. a transformation to the sccqucncy domain. The notation
Z(f(s)l
(12)
L[f (t)1
will be used in the following. To be consistent, let
Zff(s)]
=
ffT(nT)
(13)
ti(nT) z
and
Zjf(t)) =
f(t)
3(n'l't)
(14)
(t)] W
(15)
Since g(s) f(s) = L[g(t)
it follows that Zig(s) f(z)j
Szil u(nT) f [g(nTt) u(nTt) 'n=cx>
f(kT)
u(t) 6(kt)] dt
kff0
COo
(16)
From the properties of the Dirac delta, [6,71 and Appendix B, 8
a
(17)
g (iFrkT) f (kr)
.z
*
gtnrkT) z
(17
(18)
t(kT) z
nO) ,<(0
which is the Cauchy product of power series
= t g(nT)z n.=o ti II k=U (kT )"
=
g(z)
(19)
,
(20)
t(z),
the RagazziniZadeh identity. Ztg(s) f (S)J
(z)f(z)
but unfortunately ([gts) f(s)
# g(z)
f(z)
that is g(z)
# y(z),
and therein lies the difficulty in applying ztransforms to continuous systems, that is, the digital simulation of continuous systems.
3.
MEAN VALUE CONVOLUTION
A very useful relationship would be the solution of Zig(s)
z n u(nT)
f(s)]  I n
f g(t) u(t)

9
f(nTt) u(nTt) dt
(21)
where both "g(t)" and "A~t)" are functions. There are no Dirac deltas to "remove" the integral as was the case with the RagazziniZadeh identity. The mean value theorem of the integral calculus [5] guarantees there is some "k such that (k+1)T
f
g(t) f(nTt) dt
=
T g(kT+ 6 kT) f(nTkT6 T)
(22)
kT
and Equation (21) may be rewritten ni T g(kT+6 I.T)f(nT.kT 6 kT)
z n=O
Since
"&k
(23)
k0O
is
most likely different for each "k" it would be difficult to proceed.
Assuming
6(24)k=
that
is.
"45k"
is
the same for all "W', one has
00 n1 T g(kT6T) f(nTkT6T) xznI k=O n0O
(25)
and as before 1k
OD
Ty
nk f(nTkT6T) z"
g(kT+6T) z
(26)
n0 k=O
One problem, the sum is to "n

I" not "n" and to obtain the form of the Cauchy
product a term must be added and subtracted, 10
H
g~k+6T) '~
f (,vTSr) zf
(6T
(27)
and, finally, from the definition of the modified ztransform
ZIg( .) f()l
T
g(z,)If(z.,6)
 f
(28)
T)i
which will be referred to as Mean Value Convolution, MVC for short.
Since the assumption. Equation (24). is suspect, some checks are warranted. Consider the case where .(29)
s
and (30)

From Equation (28) and Appendix A
(31)
1z
Checking against Appendix A for I /s' Equation (31) is found to be exact. For g(s) the same and
111
f(s).
__
(32) Equation (28) and Appendix A yield
S
(33)
which simplifies to
T z (tt
) nzj
Z) 3
(34)
For a small time step. one would suspect that
6. (35) would be reasonable and since
(36)
one would have
2 2(1 z)
3
(37)
which for I/s' is exact.
12
Interchanging the roles of g(s) and f(s) in Equation (28), one has
+(~ Z G
T (T
]
[
(38)
and simplifying
T'2z (nj+(l :au) 3
(1 z)
)
(39)
In this case n = 6= 2(40)
and Equation (39) becomes
T2z (1+z)
(41)
2(1  z)
which for I/s'is again exact. That the above checks are exact is not surprising since the first would correspond to integration of a constant, and the second to integration of a ramp or the double integration of a constant. For these cases
k
62
(42)
=
is what one would expect. So far, so good, but I/ s4 presents difficulties whose discussion will be deferred until a solution to the difficulty is available in a later section. 13
_II__
_
_
_
_
_
_
_
I__ _ _
_
_
_
_.__
_
_
_
_
.. _....
_
_
4.
RK CONVOLUTION
In the previous section the mean value theorem was invoked to "remove" the integral. In this section the approach will be to "numerically" integrate. Since the approach of this report is that of ztransforms a brief summary of integrators with holds is in order. Consult Jury 161. Smith [51 or Rosko 171 for details.
Z[K(s) f(s)
(4:3)
 y(z) f(z)
may be approximated by
.(l(s)\ Z [g (S)t" €)] "
g(z)
"(44) Z
1
that is
Y60 = S(Z) F(
1+
(45)
where "n"is the order of the hold. For an integrator, "f(s) = I/s", and a zero order hold, "n
(Z) W)
T z v.(z)
14
0", one has
(46)
Eular integration or the Left Riemann Sum. For a first order hold, "n = I",
(z)
g2(z I~) 1
(Tz)
(47)
+zg(z) 22(1  z)
trapezoidal integration. For a second order hold, "n = 2",
Z(4) g(Z) /
Z()
(T3 (1+
4z + z2 T i + 4z + z
\60I_
(T 2 Z(1;4)
gZ
g(Z)
3 G1+ z) (1z)
Simpson's rule. Smith's "tunable" integrators [51 may be obtained by using the modified ztransform. The "tunable" integrators for zero and first order holds will be found in Table I. One may proceed to higher order holds (Table 1), but these integrators do not appear practical.
A digression: Of course, holds may be used with plants other than integrators; for example, a single pole filter, f(s)
(49)
s__.
15
+
J+
0 0
4
+
+
0
~
CN
2~
1 N
N
~
+
I +
E0
(NI
4
16

+
+ +
+ N.~

4
l , 4
4.4
7
A single pole filter and a zero order hold, "n
g(z) Z(i)
= 0",
yield
(50)
( T z g(z) a(lz e aT)
If one chooses to implement Figure 1,
Hol
g (S)
T
T
Figure 1. An Integrator with feedback. an integrator with feedback for a zero order hold, one would have (s"()
\
(51)
T z g(z)
"1z(1aT)
which would amount to using a (I/0) Pade' approximation for e T,
aT
e
1aT,
(52)
in Equation (50). Dividing Equation (51) by (50) one has the ratio*
(aT
1_aT aT aT )1ze)
*Subtract one and multiply by one hundred for percnt emr.
17
Applying the final value theorem,
.oc. the ratio is one, as desired, but
1 (n
applying the initial value theorem, / 0 (n0). one has the ratio
Since "a" isgiven, "T" must be chosen so that Equation (52) is as near one as desired, or use Equation (50) instead. For a single pole filter and a first order hold one has

LY,,,' (z)
aT 1)aT
aT
7(11) + (,
,1
aT

aT
j'.
)
(53)
)
It is common practice in digital simulation to follow the analog computer practice of reducing problems to interconnected single integrator'. This is not necessary nor may it be desirable, but enough of this digression. In the following, several integrators (Table 1) will be used to find a,)proximate solutions to Equation (21). There is some repetition of earlier material 121 which is included for completeness. Using Eular integration, Equation (46), in Equation (21) one has* (k+) T g(r) f(n'ri) i z" Y" f k0 kT 1n)O A"
n1
, it=
n
)
T' g(kT) T
f(nTkT)
(54)
n=O k=O
As before, adding and substracting...
=
' in Tz 11=0
k0
k() g
f(nTkT)g(nT)f 0
*How recurrences are developed istreated inthe sample problem sction.
18
(55)
The Cauchy product...,
=
T
W) zk
f(nT) zj T fO
(56)
g(nT) zn
and from the definition of the ztransform,
Z [(s) f(s)]= T g(z) (f(z)f
,
(57)
Eular convolution. If the modified ztransform is used (2], for a zero order hold one has oo I
n1 (n g(nT+T) f(nTkTT) + (1l) g(nT) f(nTkT)1
z T
n=O
(58)
k=O
and then, as above,
co n
T zn I g(kT) f(nTkT)n=O kf[O
Znfrio(1)
T
Zn g(nT
f(a) + (l2.) f
go
nO
n=O
(59)
and finally,
Z[g(s) f(s)]
T g(z) f(z)T[Tg
o
19
f(z) + (in) g(z) f]
(60)
Eular tunable convolution. Proceeding on with the trapezoidal integrator,
Z
Tn1
2
n.0
(61)
[g(kT+T) f(nTkTt) + g(kT) f(nTkT)] k0O
and after manipulating indicies
x T
g(kT) f(nTkT)
g0f(nT) + g(nT) f 1(62)

and rearranging 00 nk
nk
g(kT) z TY n0 k=O
T~ [o iz
f(nTkT)z
n g WT) + goiz
f(nT
(63)
and finally.
T f (z)
Z [g(s) f (s)]
g (z)
g (Z) + g
[f
f (Z)1
(64)
trapezoidal convolution [7,8,91. Equation (64) may also be written
Z[g ()f(s)]
1
( (g (z)
2

g)
f (z) +g (Z)
0
(f (z)

f0 )]
(65)
For the modified ztransform and first order hold (Table 1) one has
jz n=()
+
~'
I g(kT+2TI)
f (n1k1*2T')
k=0) I
(1 + 2ri
2rK2) g (kT+T)

f (nTkTT)
+
( 12 r + i2 ) g(kt) f(nTkT]
(66) Changing summing indices
z i
11
2
nj,,2'+]
g~kT)f(nTkT)
+
~nT
(1+1j2y )Y+kT
k='
ki
ni +(121j+rj 2 )j g(kT) f (Tk) kO
(67)
and adding and subtracting the necessary terms, nn12 ,
TY z' Y g(kT) f(nTkT) n=U (k=U(1Og

(
+ 2n + q2g
2r
1
+
2
f(nT)
[gn
+T
~
f(T)
T
g(T)
f(nTT)1I
(68)
21
From the definition of the ztrangform (and the Cauchy product), finally,
Zlgls) f(s)l
T 9(7) t (.)
+ T2 [g z.I

'2(0
f()
+
2
g(T)
ri

11")go
f z,1)
f W)
.
(69)
Trapezoidal tunable convolution. Simpson's rule (48),
3(1 +z) (1  z)
presents indexing problems since the recurrence would be
X~ = X
+ :Ik
+ 4ki
+
k(70)
Also, Simpson's rule is known to be sensitive to noise [6,10]. Another form, used in third order RungeKutta integration, is
T
+4 (1/2z+
(71)
whose recurrence is
Xn
=X
+ ni
+ 4 in1/2+
n1
22
(72)
Equation (71) may be factored into 3 k1z
3 \1 z7/(3
that is,
2n .n 1
3 i 2
i. )J+ 23 T 6 .n + :.n1
' n1/2
(74)
Note that Equations (73) and (74) are a linear combination of trapezoidal integration, and Mean Value integration with 6 = i 2. Therefore RK convolution isa linear combination of trapeioidal convolution and Mean Value Convolution,
Z [g(s) r(s)]
g(Z)  go) f(z)
+ 4 g(z,1/2) If(z,1/2) f(T/2)] + g(z) (f(z) 
f)]
(75)
S. ANALYSIS OF ERRORS Earlier the ratio between the zero order hold integrator with feedback and the zero order hold single pole filter was taken to compare these different simulation approaches, Equation (51). The input was not specified. Now the ratio of the convolution approximation to the exact solution will be taken for given inputs. The accuracy is highly dependent upon the input. Checks for the various convolution approximations, Table 2, similar to those done for MVC, Equations (29) through (42), will be found in Table 3. The ratios to the exact ztransform will be found in Table 4.
23
TABLE 2. CONVOLUTION APPROXIMATIONS FOR Z[g (a)I s)
M.V. C
T 9(Z,6) [f (z,6)
E.C.
T g (z) If (Z)  f 0 ]
E.T.C.
1
T.
[g(z)

2 T)
1 [(g(,)

g.)
T.T.C.

f (6T)]I
f(Z) + g(z)
f(z) + g(z)
(gz
(f(z)


2 (1ni)
f(z) +
4
g(T)
g(z,1/2)
f)
fo)]
g0] f (Z) + 9(z) [f (Z)
(1+2r1n')
+ ri2[g(z,1) f(T)
RKC.
(f (z)

(12rri)]
f(z,1)]I
fz,/)(T1]+()(fZ)f 0)
TABLE 3. CHECKS OF THE CONVOLUTION APPROXIMATIONS FOR 1(s) g(t) A CONSTANT, RAMP OR PARABOLA
2 m~~.Tz
M..(1Z)
EC E..(1z)
T2z(l+z)
2
zT2Z2 2
Tz (1Z) 2
T.C.
RK.Ti
RK.(1Z)
2
24
=AND
3 q 2 T z(1+6z+z )
(iz) 4
2(1z) 3
8
(1i) 3
Tz(1+z) 2 (1i)4
T2z(1+z) 2 (1z) 3
T3z(1+2z+z2 ) 4 (10)4
T2Z(1+Z) 2 (1Z) 3
T3z(1+4z+z2) 6 (1z) 4
TABLE 4. RATIOS OF APPROXIMATE SOLUTION TO THE EXACT FOR A CONSTANT. RAMP AND PARABOLIC INPUT TO A SINGLE INTEGRATOR
ts)
N.VC
T.
I
)7
C.)
ki
I+
~
kit
I 1wc uto
4z +
I

(kit
\%otuld he teio ot
lukld
tie
Iumotmol
miltnM
me.tm
vcta lon.T o z
0 the ratio
tahe a411111) and paawaholac inputs. 1I (ihe **t in the nuivneraloa is
in~terprtetd na it sahift, the ratious wo~uld he two and three,~ reaptetively. 'rhe approxiantion for the ramp waaa dettermned~~ with
~] ~
(76)
Initerchantging gtaa) and fAs) in Equation (57) one haui


and the ratio would be
25i
(77)
But recalling that
(79)
T
is Rectangular integration, one may conclude the difference between Eular integration and rectangular integration, the Right Riemann Sum, is a time step. This difficulty arises because of all the approximations in Table 2, only Eular Convolution lacks symmetry, and convolution should be symmetric.
One could proceed to higher powers of time for inputs, but b ing only able to determine the ratio at the initial time is not very enlightening. A time step by time step approach is possible 121 but does not seem practical. Fortunately some transcendental functions work very well as inputs for ratio analysis. For example, consider an exponential input to an inte, rator. That is,
f (s)
g(t)
(80)
=a e at(
(81)
and, it follows that
(82)
g~
taking the ztransform of Equations (80) and (82),
(83)
f(z) 1z 26
K(Z)
=(84)
1
aT
v
and substituting into Equation (57) yields 1)'
(1
~z~ea)(6 _8) z___ _aT

z
\z
for Eular Convolution. Taking the ztransform of )
(I
aT
_____
\s~~aJ(Z) (I
)
(86)
zT)
yields the exact solution. Taking the ratio of Equation (85) and (86), aT 1 at
(87)
Note that no ..z's" remain. When this occurs, the assumption,
(24)
= 6
appears quite reasonable since the ratio does not depend upon the time step 121. Factoring out an
"eaT/,," in
the denominator, Equation (87) becomes 27
a
aT e T/ (eaT/2
2

aT/2)
(88)
and, since, aT /2 e
aT/2
= sinh
(9
+ cosh aT/2
(89)
and
sinh
aT/2
)
eaT/2 T2(90)
(eaT/2 1/2 12
after a few manipulations one has
(aT /2) [ctnh
(aT/2)
+ 1]
,
(91)
the ratio for an exponential input to an Eular integrator. If "a" is imaginary, that is,
a  iW
(92)
Equation (91) becomes
(wT/2) [ctn (wT/2)
+ i]
(93)
28
Since this "ratio" is complex, one must determine its amplitude and phase. Multiplying Equation (93) by its complex conjugate and taking the square root, one has
(wAT/2)
[ctn2(wT/2)+
(94)
111/2
which readily simplifies to
(95)
(wT/2) sin(wT/2)
for the amplitude ratio of a sine wave input to an Eular integrator. Dividing the imaginary part by the real in Equation (93), and taking the arc tangent yields
(wT/2)
(96)
for the phase error of a sine wave input to a Eular integrator. One would proceed in a similar manner for the other integrators of interest; see Reference 2 (6,and 10) for details. The results are summarized in Tables 5 and 6. For a low frequency sine wave input the Eular integrator has half the error of the Trapezoidal integrator, but unfortunately it has a linear phase error. For 6 = I/2, MVC does not have the phase error and since there would be no advantage to choosing a "B"other than one half, it is worthy of consideration. Also,the ratio of MVC error to TC error is 1:2, and the errors are of opposite sign which serves to explain their mixture in RKC. The linear phase error of EC and/or the fact that MVC samples at mid interval for 6 = 1/2 may explain the observation that in some simulations the difference between the simulation and exact answer is a shift in time. 29
TABLE S.
ACCURACY FOR AN EXPONENTIAL INPUT TO AN INTEGRATOR
TMANSFORI4
RA It)
KII
I ZI z eAT)
N.y.C.
TI
j.
aT
1.C.
az(TI (IZ)(I ze AT
E.T.C.
AT kh(I 
(aTIZ
Simpson
KK.
)e )AT
I(4. +
T)
_k
siiih
(d'/
2
)(Vtnh (aT/2 I + 2(,
6
+
)
cndT211I
6fl
4
al
+
(1
(aT/2') ctith ( aT/2 I+
t.
_____
AT(I + 4e.T2+e& )z 6(12)(1ge 
j2.I)~~ dT/ .I
(AT) 6
cosha)0
2
, osh AT sinh (aT/2)
30
AA)
4
a T1.
(aT)?
I, +
12
TABLE 6.
ACCURACY FOR A SINE WAVE INPUT TO AN INTEGRATOR
PHASE
RATlO
0
Exact
1 + lT)
(toT/2) sin ( wT/2 )24
M.nV.C.
TC.
Simpson
RC
(twT/2)
ctn(wT/2 )
/w•
2
wT/2
1 + (T) 24
(toT/2) sin (wT2)
E.C.
(1/2  6) wT
2
1
o(w)
cg
\
+ cos T tT) 3 (2sin
1+
6T2
+ sin cos (wT/2)+
4
(wT/2)
31
45 45
Saiw)0 720
0
N
N
I3
3
0
0
N 
3
0 I4
N
'0
~ 4
'
I~, N
N
N I
4
4
0
ma
+1

lal
0
0

*
0 '0
4
z
'4
o
4
'4, 3
4

+ 4
4
4
4
ma 4
ma
N N
N..
a
4


+
Nfr~
2 4
0 IL
2
+

N
N
~
I3 N N
4.1
U
L~.l
o4
N
4
+
N
4
N
U
.
.
N
N
3
C
£~~J
'a ~
N
3
I. 3 N
ma
'0
S 4 I.
N.
+1
+1
N
N N..
*
U
U
*
*
*4
a 'C
32
*
C.)
4
*
f.4
'C
Though it is unwise to sample at the Shannon limit of two samples per cycle, it is instructive to examine the ratios and phases at that limit, Table 8. The noisy qualities of the Simpson integrator are apparent. Since RKC takes twice as many samples as Simpson. it has no difficulty at the limit; its ratio is two thirds that of MVC since TC's is zero. Trapezoidal integration appears particularly bad at the limit since it might stabilize an otherwise unstable system. The last entry was tuned [2] for unity ratio at the limit but it has a ninety degree phase error. The implications in determining gain and phase margins are obvious. TABLE 8.
ACCURACY AT THE SHANNON LIMIT (wT
RATIO
1
Exact
= 1.57...
M.V.C.
=
7r)
PHASE
(0 dB)
0
(2 dB)
(1/2

6)
TTT
E.C.
2 = 1.57.,.
(2 dB)
2
T.C.
0
( 'dB)
0
Simpson
W
dB)
0
RKC.=
1.05...
E.T.C. E.T.C.
6
6E.T.C.
(0.2 dB)
2
( )1/2
0
+ 7T
1.28
(1 dB)
+ 
1
(0 dB)
+
2
6
33
6.
SOME SAMPLE PROBLEMS
Now to the meat of the subject, the numerical solution of differential equations. To incorporate nonzero initial conditions, one should proceed from the differential equation using n1
nIs X(s)
L[X

I

IX
ni
CM
(97)
To attempt to proceed directly from the transfer function is difficult because to obtain the transfer function, the ratio of output to input, it was assul.led that the initial conditions were all zero. A. Single Integration For "n = 1" Equation (97) becomes
X(s)

s X(s)  X
(98)
0
and solving for "X(s)" one has
X(s)  *(s) + S
S
(99)
0
Taking the ztransform of Equation (99),
X(z)

Z~()+
x0 Z(I)
(100)
There is no difficulty in taking the ztransform of the initial conditions since they are constants in the frequency domain; they are multiplied by Dirac deltas in the time domain [23. 34
The input, "X(s)", is another matter and requires the application of the convolution approximations, Table 2. Now to develop the recurrence for MVC, Table 9,
(z)
TABLE 9.
X(z)
= T z X(z,1/2) + X0
(101)
CONVOLUTION APPROXIMATIONS FOR SINGLE INTEGRATION
(X (s
M.V.C.
T X(z,
E.C.
(1 z) [
E.T.C.
TIn +
1/2)
" 1
or T k(z)
(1r)z] X(z)

Tr
1z
T.C.
RKC.
A+ 2 (LtZ)j Tz
o 1 z

T2 1z X
6(1z)[.*(z) + 4 (i(z,

1/2) X(T/2))4 (X(z)
X)]
and substituting the definition of the (modified) ztransform, one has
SX(nT) zn  X(nT) z n0
(nT+T) zn+l + X
T nffO
35
2o
(102)
Now the important step, equating coefficients of like powers of 'T', one has =O.
Xn
(
Xn
(103)
+ T
,n1/
n>O
(104)
It is important to note that the general recurrence, Equation (104), does not apply at nO'0. For EC there are two possibilities because of the lack of symmetry.
(1z) X(z)
=T
(X(z)k I 0
+
X
(105)
0
becomes
00 zn 1 X (nT)z n0
X(nT)
zr1~+1 00z+X(16 z =TJ k(nT) n=10
0~+X(~
and equating coefficients of like powers of "z", X0
Xn
X0
(107)
(108)
Xn1 + T knn>0
Rectangular Integration! For
(I z)
X(z)
 T z(z)
+X
,
36
(109)
(II0
n+1
) no X(nT) and, ...
,
z

X(nl)
z
=
T
n+ . X(nT) z n=O
l
+ Xo
(110)
finally, x= x
(111)
+ T Xn1,
Xn = Xn
(112)
n>O
'\Eular Integration.
Proceeding in a similar manner would produce Table 10. It might appear at first that MVC and RKC are not symmetrical since, if the rolls of "g(s)" and "f(s)" are interchanged, one would have
[X(z
T

1/2) 
(113)
X(T/2)]
1z
and T 6(1z) [(k(z)  90) + 4z X(z,1/2) + z X(z)]
,
(114)
respectively, but these approximations lead to the same recurrences found in Table '9. Of course, these recurrences, Table 10, are those used to develop the convolution
approximations in an earlier section, and are not unexpected. At least they demonstrate consistency.
37
RECURRENCES FOR SINGLE INTEGRATION
TABLE 10.
Xo = Xo
M.V.C.
X n = Xn_ 1 + T knl/2
E.C.
X n = Xn_ 1 + T Xn
E.T.C.
X n = Xn 1 + T[rj X
T.C.
X
T.T.C
Xn
=Xn
RKC.
X
= X
B.
=X
n
n1
i
+ 1(X 2
+ T
n
+ (irn)
+ k
I2ri)
+ 6 [X
X 11
n1
X+
(l +2r)
Xn1+ r]2(X
_
2k
+ X
k
+ 4X1/2 + X
Double Integration
Double integration is interesting because it illustrates some of the problems associated with initial conditions and start up. Equation (97) for "n = 2" is
X(s) = s2X(s)

sX
0
 x
(115)
0
and solving for "X(s)", one has txo
X(s) =
+ X0 + s0
(116)
38
Taking the ztransform of Equation (116),
X(Z) 
(
Z(
+'(,).(117)
Before proceeding, consider the case where R(t) 0,that is.
X(W

X
X(18
+
The exact ztransform is X(z)
TzX(
. Tz
+
X
(119)
(Z) 1
(t or, rearranging,
(
Z) X(z)
=
T z X + (
(120)
z)X
Substituting the definition of the ztransform and equating coefficients of like powers of "z", X0 . X0 ,
(121)
X1 X 0 +T
(122)
Xn
2 XnI
(123)
2' n> 1.
It is important to note that the general recurrence, Equation (123), may be first applied at "n = 2" and DOES NOT apply at "n = 1". As the order of the system 39
increases the number of steps before the general recurrence may be applied increases! For further discussion on this point see Reference 2, in particular, Appendix D. With a little mathematical induction one may readily demonstrate that
S=0 + n
'r
10 0
(124)
Mathematical induction serves to further illustrate the point made above. In addition to showing a relationship for "n" is valid for "n + I'"one must also find some "n" for which it is valid. For equation (124) the smallest such "n" is"n I but, for Equation (123) it is "n = 2". Table 11 contains the required approximations convolution and Table 12 the recurrences. Note how the input enters the startup step, "n=I", in particular for TC. TABLE 11.
CONVOLUTION APPROXIMATIONS FOR DOUBLE INTEGRATION
(T/2)1
( z) .'I M.V.C..
E.C.E.C
T.C.
12)
[X(z
2 T0_z z { 2
_ R0
(z)
Tz) 2 [I(z)

R 1/21
G Z)
RKC.
T2 (1+z)[(z,1I2) 2
40

X(T/2)] +z(X(z)Xo)
TABLE 12.
RECURRENCES FOR DOUBLE INTEGRATION x () = x
0
1...X = vX X x 0+
M.V.c.
X E.C.
Xn 2
1/2 R11 T2 11  Xn2 r x + T2 r
X1  X + T X = 2 X n n 1 X
T.C.
=X + T 1 0
X
RKC.
+1
o x
2n
n/2 + Rn] 3
2 ""
+T n2
n 1
,
n>1
> >1
+ 112 T 2 0
0
+ Tn
2 X
X1  X0
xn
T2 2
o
k0+T
1
2 k
/
+ k 0] /6
2  [Xn
n
+ 2
n1 +
n
, n>l 2
C. Single Pole Filter Now for a bone with a little more meat on it, the single pole filter whose
differential equation is
k (t)
= a
X(t) + g(t).
(125)
Applying Equation (97), one has
s X(s)x
Wa X(S) + g(s),
(126)
41
and solving for "X(s)", g(s) + X
X(s)
(127)
i S + a
For X. = 0, one would have the transfer function
X(t)
g(s)
18 (128)
_I s+
taking the ztransform of Equation (127),
(129)
a +
+ X
X(z) =/s))
Proceeding as before, one has Tables 13 and 14. Like the single integrator there is no startup step but note how past values are "decayed" based on how "old" they are. A leadlag, in Equation (130), is sometimes required. (130)
rIsbA g(s)
x(S)
To incorporate initial conditions convert Equation (130) back to the time domain,
i(t)
+ a X(t)

(t)
(131)
+ b g(t),
and then apply Equation (97),
a X(s)
X
X(s)

a g(0)
g 0 + b g(s)
42
(132)
4
TABLE 13.
CONVOLUTION APPROXIMATIONS FOR A SINGLE POLE FILTER
z(o') (g(z, 1/2)
aT 1z e
M.V.C.
Tze
E.C.
(T/2))

aT g(z)
aT )
(1 
aT )
T( 1
a(T+ze) fg(z)  go ]
T.C. (1z e
aT RKC.
T aT 1(1+z 6(1 ze )
TABLE 14.

e aT)g(z)

g
+ 4e
2(g(x, 1/2) g(1/2T))]
RECURRENCES FOR A SINGLE POLE FILTER Xo  Xo
*aT/2
aT
Xn1+ T e
X..C n e
Xn . e  a T Xn_ + T g
E.C. C.T
e aT eaT
. n
RKC.
nini
l
gnI/
x
=aTx
n =e
[g +
n
T 6
n
43
+4 e aT/2 T
gn 1]
2 nn]
1
+e aT • n1J _
and solving Equation (132) for "X(s)" (s+b) g(s)

g0+
(133
X(S)
(133)_________
taking the ztransform of Equation (133),
X (z)

r/ +
(
+ (X0 80 ) z+.(14
Noting that s +b 1 s+b
(ab) S+ a
,(135)
Equation (134) may be written
X(Z)= g (Z) Z [( jb~gs]+(

g)
z(
(136)
The approximation required in Equation (136) is the same as the single pole filter with the coefficient "(ab)". Table 15 contains the recurrences. If the leadlag is of the form, I +s/b . a Is + bj 1+Ts/a b k ;"+a
(137)
then the recurrences in Table 16 would apply. If one had the form
Y1+/a
aI
(138)
s +a)
44
I RECURRENCES FOR A LEADLAG FILTER,
TABLE IS.
aT T (ab) eg
M.V.C
Xn =
E.C.
Xn ' gn +
T.C.
Xn = (1T(ab)/2)

aT[X_
+
aT
n2
Ini gn
(I+T(ab))gn1]
1
(+T(ab)/2)gn]
eaT[xn
g+
aT/2
RKC.
+e

aT
gn1/2
(ab)e
gn3
. n=6(ab))
(ab))
(I+
I+
TABLE 16. RECURRENCES FOR A LEADLAG FILTER,
M4.V. C.
x
E.C.
X
T (a

a
(
n
n
b

b) eaTI 2
b
g
1
I nl
(1T (ab)) gn +e n
)+
gn1/2


b
Yn = "b a n ++ e aT[x n1  b( T a1 1n+ T(
X
T.C.
(
RC
RKC.
Xn
a
(I
(ab)
+e aT [X
a( b
=
nI
2aT (ab) e
)
gn
3

+
)
(ab) 6(4
45
n1
g nlj1 b1n g
(ab)) gn + eaTnI  A
 A
eaT Exn~ a
gn1
(I+T (ab)) g aT / 2
n1/2
s/a
~
1
gn1
in Table 16 lot
(ab)

(139)
a ,
(140)
a
Unfortunately Table 17 does not lend itself to differentiation since the requirement that "aT<
TABLE17. RECURRENCE FOR
(as)
T aT/2 M.V.C.
Xn = a(g n  aT e
E.C.
Xn
T.C.
n
RKC.
aT
a (1aT)gn+ eaT (X
a (1
Xn = a(1
+eaT (Xn 1

g n
1
an
a(1+aT/2)
2a2 TeaT/2
3
a (1 + aT/6) gn1
46

ag)
aT/2 ) gn+ e  aT (Xn
+aT/6)
(Xn
gn 1 / 2)+e
)
gn 1 )
D.
A Forced Damped Oscillator
This problem was discussed in some detail in Reference 2, and only some results will be presented here. The transfer function
21 s
+ 2
1>>
2 WS + i 0
(141)
° 0
is represented as 1 (s+a) 2+W2
(142)
where
a
=
w
(143)
and W
(1C2)1/2
(144)
.
the startup step(s) will be found in Table 18 and the general recurrence(s) in Table 19. Figure2 serves to illustrate the difference between the innerconnected integrator approach, Figure2(a), and that suggested herein, Figure2(b). The velocity, "X(s)", would be found using the single pole filter recurrences, Table 14. Unless "I(s)"and "X(s)" are needed for some other computation or output, they need not be computed.
47
•
.4
TABLE 18. STARTUP STEP(S) FOR OSCILLATOR X X 0
0
aT+Isin wT) X
eaT (cos
0
1
M.V.C.
+ TeaT/
2
sin wT/2
+
e aT sin wT w
0....
g
g1/2
E.C.
+0
T.C.
+ T T
aT/2
+ j e
RKC
TABLE 19.
sin wT
aT
E.C.
+ T •
+Te
g
+ e  aT/2
cos wT/2 go
GENERAL RECURRENCES FOR OSCILLATOR, n>l
=2e aT cos wT Xn1
M.V.C.
wT/2 [2
sin
aT 2
 aT a
/
2a T
sioCu12)
%2
['n1/2
.... +
e'a aT gn3/]
sin wT 8nl U)
T.C. T.C
+ TeT sin WT
+ T + 2
eaT/2
gn1
sin(fL/2
(&T/1,2 48L1+
48
S cos WT12
aT &a _1,.
]
Figure 2a. Interconnected single Integrators.
Y
Figure 2. A fored damped oscillator. 49
Its)
7.
RESULTS AND CONCLUSIONS
The Heaviside unit step and the Dirac delta were shown to be useful in outlining a proof of the RagazziniZadeh identity. The mean value theorem of the integral calculus and various numerical integrators were used to develop approximations of use in the digital simulation of continuous systems. The accuracy of these approximations was then studied for various inputs to a single integrator. The inputs considered were a constant, ramp, parabolic, exponential and sine wave. Recurrences were then developed for single and double integration to illustrate incorporation of initial conditions. Recurrences for a leadlag transfer function and a forced damped oscillator were also presented. The sample problems were simple to illustrate the approach. The approach has been applied to the real time digital portion of the hybrid simulation of a high spin rate air defense missile [11]. The simulation is presently in operation. Another application would be subroutines of often countered transfer functions for use in simulation models such as Reference 12. Further effort on the recommendations[2] is required. It is hoped that in trying to determine the facts that additional fallacies have not been introduced. Comments are welcomed.
50
APPENDIX A. TRANSFORM TABLE FOR SELECTED FUNCTIONS' AND DISTRIBUTIONS
*Healy. M.. Tables of Lu4place. Fourier and Ztransforms. London: W. and R. Chambers Ltd.. 1967. Cadzo*, l.A., DisreteTime SYstems. Englewood Cliffs. New Jersey,: Prentice Hall, 1973.
51
+
to + C1
I N
. c' I
+
'4.4
t
N
4
"4 %I
CN + r"
N

41
*
4)
1
41
N
+
4
N
N
N
Ij1
41taI
N
+
4t
N
441
0 sw
%52
*
N4440
100
IQ00
0 #
"
+'61
+ ('4
+
04
('
#4
u *a
42
'a 0
++
+ 0
0
0 SN
+N
+
i+
V4.4
.
sr
94
*
0
.4
I*44
41
.*N
.4
ul
aa 3. "Uag4 N4 +
P4
3
4
+
'N~. U

1
N
a

o
0 
%

54
I 

do
APPENDIX B.
DIRAC DELTA DISTRIBUTION
55
APPENDIX B. DIRAC DELTA DISTRIBUTION*
J
f(t) 8(1:
(Bi1)
nT)dt  f(nT)

60
where f(t) is a welldefined funiction at t nT. 8(t) N(at)
(B2)
5 (t)
(1)(B3

Xlg;(nT) 1 6( t
Ig (t9)

g(nT) nT), 9(T
t 8(1:)  0 f(t)b(t
f
(t

0
(B4)
(B5)
nT)

f(nT)B(t
r)b(t

nT)dt



nT)

8(.r

(B6) T)
(B7)
b
8(t)
f (t)dt (1),
f(o)(B8
*Dirac, P. A. M., The Principles of Quantum Mechanics, London: Oxford University Press, 1968, Messiah, A., Quantum Mechanics, New York: John Wiley and Sons, 1964.
56
APPENDIX C.
NUMERICAL DIFFERENTIATION
57
Differentiation is to be avoided on analog and digital computers whenever possible, but there are occasion ... Equation (97) for "n
I" is
and taking the ztransform
UZS X(s))
k(z)
X(O)
.(C2)
Since e St
z
(C3)
it follows that S
=
Qn(1/z)
(C4)
and Equation (C2) becomes
k(z) =
(zQn(1/z)
X(s))
X(O)
.(C5)
Applying the RagazziniZadeh identity, Equation (C5) becomes
i(z)
=1Qn(i/z)
X(z)

X(o)
.(C6)
58
One of the possible series for the logarithm is
'n)
(1
=
z)

n=1
(C7)
n
and substituting into Equation (C6) and equating coefficients of like powers of z: 
(C.8)
0
0 SX
l
1
I X
 X
T
'(C9)
x 2 X T
1
=
X=
x 3 X
2
T
+ [(X3 
(x 2 X
)(x 2T
1
+ +
x2 )
(X
3

(X 2
 X

2
1
x)
(C10)
)(xXl)
2T
X 1 )]
2

x1 )

(X
X
(CIl)
3T
This series was chosen because each term is a finite difference, and there is no attempt to use any term unless the whole term can be used. That is, for a given "X.,"
n1z)
n
Z)
(C12)
ii Since all past values of "X," are used, a general recurrence never occurs!
59
One might truncate the series at one, two or three terms, i.e., X x n XT
'
2T 3X
n
(C13)
, n>l1

(C14)
,
2
+9
ni
11 Xn n
l*ll2
18X
3IX
,
gn2 O
X.
T

n >0
Xn3
Xn2 6T
'
n
(C15)
2
respectively. Of course, integrators may be developed in a similar fashion (Table Ci).
TABLE C1
INTEGRATORS VIA LOG APPROXIMATION
n(i/zf
z
Rectangular
T
lz
lz
11 +A
A1
2 k 
2i +z
T_ +
CLASSICAL NAME
Eular
T z iz
z
1

1 = T s
fn(l/z)
1
.1z +
"
.
3T 8(
]
(1 + z)3  z)
60
Trapezoidal
Simpson's 3/8
Since some checks are warranted, consider the differentiation of a cosine wave. Equation (C13) would become
cos ntoT

COS (n
C
Qi____
which may be reduced to
(
in WT/2
(C17
/sin(nwT wT/2
WT/2
Of course the answer is "w sin(nwT)". The mean value theorem of the differential calculus states
n nlnT
xni

0O
(
(n1)T
If in Equation (C18) Ti =
(C19)
112
then in Equation (C17)
sn
wT/2 wT/2
(C20)
may be interpreted as an amplitude error.
61
Since there is tunable integration [5), why not tunable differentiation?
YOi Z )+ Tz
(1

y)(1 T

z)
= y
2*y)z Tz
I

(1

y)z 2(C21)
Invoking the mean value theorem, for a cosine wave
w sin (nwT

nwTf)
cos (n + I)wT + (1  2y) cos nwT
=y

w
[(i
wT ) sin nwT
1y)cos (n
+ (2y
1(

oWT)
I~T(C22)
cos nw]j
(C23)
and finally
tan rTwT(2y 1)
tan
(C24)
wT/2
Therefore

I
tan
121 ta
ta
WT/ T2(C25)
62
or tan r :T"
,v]/(1
t1
1/21
)
1 /
(C26)
In small angles, "i' T.  i T

12 1
,
(C27)
or 1/2 + rY
For any
"L,
(C28)
T,"
TABLE C2 TUNINGS INDEPENDENT OF "wT" nY 
1/2
0
0
1/2 1
112
When the other differentiators of interest are checked in a manner similar to Equation (C13), Tables (C3) and (C4) would resi It. The only difference between the tableL is one of point of view. When small angle approximations are warranted, Tables (C5) and (C6) are appropriate. One may conclude that for the first two differentiators, "y = 0" and "y = I," a phase shifted interpretation is more appropriate, while for the last two, "n = 2" and "n = 3," it is not. At the Shannon limit, Table (C7), the first two are 36%low while the last two have suffered a serious "phase" error. Differentiation is to be avoided. 63
IL
CA
ILH
33 z
c+ +
+
;A0
3
0' 3
Ji
4
k
C
3
CH
A. P
II
m
3
(N3
H
zw
I
(N
E30
V44
4
3
H64
0
3
03 H 3
';A
H
4
w
u3
()
U)
3
W
IL
HH
4C,
CH
4
+
I
iii
(
3
w
cnI
IA.
3
4
'
C'4
H
H
3
4C.
04'
650
4
4.L"4r
4
H
F
HH
En
I
F
N
F
( in
1
1
I
I
C,
40 MI
4
I N
o z6
iz
m H
4
n
H 1r

H
TABLE C4
DIFFERENTIATION OF A COSINE WAVE, wT'
i
r)r 0

[l
o' in ,ncsT js2)s
(  j t)


2
7'(n .,j,, OT) WT 1)
(()4
67
=o

3+n w
TABLE C7
1
T z
DIFFERENTIATION OF A COSINE WAVE, wT  7
tn~l)w
sin nwT
Tz (y = 1/)
0W2n i~nTr2
w (2/i)
y~/
z12
sin nwT +T/
T
1X(z
I)n (1
(n =2)
w (4/a)
cos nwT
(n
w(20/31)
cos nwT
=3)
68
.......
APPENDIX D. RK(N) CONVOLUTION
69
Simpson's 3/8 integration rule is
X(z)
=k
3T (1 + Z)O (z) + X0 z
(D1)
and the ratio for a sine wave input is
(D2)
[Cos 3wT/2 + 3 coswTd/2
4\2
sin 3wT/2J
Note the rather poor behavior at wT =27T / 3; BOOM! One of the forms used in fourth order RungeKutta integration is X(O)
1 1+z1/3 )3.
(D3)
+ 
X()=X(z)
and the ratio for a sine wave input is
I ()T\coswT/2 + 3 cos w/
At the Shannon limit, "wT
7
"the
D4
ratio is 1.0202621
70
...
Equation (D3) may be written
X(Z)
=~
[z)
(T
I.
[32
+
3 [Tz2/3]
()+XO
As with RK(3), Equation (D5) is a linear combination of trapezoidal integration and mean value integration but with Yi = 1/3, 2/3 in this case. Therefore
Z18 (S)
S
8 [
Z
(f (z)
+ 3 g(z, 2/3)
+ f(Z) (g(z)
ff(0)) + 3 g (z1 /3)
(f (z,.2/3)
g(O)]
(f (z ./3)  f (T/3))*
f(2T/3)) f
(D6)
RK(4) Convolution.
Higher order convolutions, RK(N)C, could be developed in this fashion using trapezoidol integration/convolution and mean va'ue integration/convolution.
71
REFERENCES Pennsylvania: AddisonWesley Publishing Company, 1972.
1. Gaskill, R.A., "Fact and Fallacy Digital Simulation," in Simulation edited by J. McLeod, New York: McGraw Hill, 1968.
8.
using Trapeoidal Convotion, ITM
2. Dickson. R.E.. Tunable Integration and Tunable Trapezoidal Convoltion A Potpourri,Technical Report TD
Systems Division, 26 August 1960.
7712, US Army Missile Research and Development Command, Redstone Arsenal, Alabama, 5 May 1977.
9r Atens 9. Halijak C.A.. t SanpledIata S and Genratmg Fun'tions. NASA CR61349. March 1971.
3. Sneddon, I.N., The Use of Integral Transforms, New York: McGraw Hill. 1972. 4.
Halijak, C.A., DigitalApproxiination of Solutions of 1)1ferential Equations
10.
Hamming, R.W.. Numerical Methods for Scientist and Engineers. New York: McGraw Hill Book Company. 1973.
P.P.. W. andof Teodorcscu. Kecs. the Theory of Application Ditribtions in ehanics D is tribu tio n s in M e c h a n ic s .
II.
Tunbridge Wells. Kent: Abrams Press, 1974.
Grimes. V.S.. Jr., "Real Time Digital o lq h ol n i r m , Model f ep Rolling .4ir Framey Technical Report T7875, US Army Missile Research and Development Command. Redstone Alabama. 31 July 1978.
5. Smith, J.M., Mathematical Modeling and Digital Simulation for Engineers and Scientists, New York: John Wiley & Sons, 1977.
12.
Arsenal.
Holmes. W.M.. Kesting. J. and Mitchell, E.E.L., Modular Aissile Simulation Model Development Using Advanced Continuous Simulation Language, Technical Report TD7714. US Army Missile Research and Development Command. Redstone Arsenal, Alabama. 3 June 1977.
6. Jury. E.L., Theory and Application of the ZTransform Method, New York: John Wiley & Sons, 1966.
7. Rosko, J.S.. Digital Simulation of Physical Systems, Reading, 72
DISTRIBUTION
Defense Documentation Center
No
No
Copies
Copies Commander
12
Armament Development Teet Center (AFSC) ATTN: ADA (L Parks) Eglin Air Force Bae", Florida 32542
Cameron Station Alexandria, Virginia 23144 11T Research Institute ATTN: GACIAC 10 West 35th Street Chicago, Illinois 60616 Professor C A Halijak Department of Electrical Engineering University of Alabama in Huntsville P. 0 Box 1247 Huntsville, Alabama 35807
1
J. M. Smith & Associates P. 0. Box 324 Greenbelt. Marytand 20770
2
Bolt Beranek and Newman, Inc. ATTN: DrL S. Baron Dr. J1.E. Berliner 50 Moulton Street Cambridge, Massachusetts 02138
DRSMILP. Mr. Voigt DROMIT. Dr, Kobler DRDMI.TE1 DRDMITO. D. Cifiax DRDMITD1 DRDMITDD1 DRDMITDK1 DROM(TOS ORDMITDR1 DRDMtTDF. V. Grimes1 DROMITOD. R. Dickson DRDMITL DRDMITKI DRDMITB DOMITBD DROMITI (Reference Copy)I DROMITI (Record Copy)
I I
73
I
100 I 5
3 I