Trapezoidal Convolution Revisited - Defense Technical Information

Loading...

IN,. 0

TECHNICAL REPORT T-78-66

U.S. ARMY

TRAPEZOIDAL CONVOLUTION REVISITED: R-K CONVOLUTION OR THE DIGITAL SIMULATION

MISSILEOF CONTINUOUS

RESEARCH

DEVELOPMTE(IT COTIMiA fD C3

SYSTEMS VIA Z-TRANSFORMS

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

R-K Convolution or the Digital Simulation *~~of Continuous Systems Via Z-Transf 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: DRDMI-TD 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

DRDMI-TI"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

.rZ-transforms Ragazzini-Zadeb Identity Convolution A _5

ACT

(MNI

heom Report)

M e.,Ov IND N 0=908

=9 0000Pt

loek monkc)

A proof of the Ragazzini-Zadeb 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. R-K 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. R-K(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 Lead-Lag Filter

(g±b 16.

................................................

4

Recurrences for a Lead-Lag Filter

+ ./a

.................................................... 45

Recurrence f r( S

a) ........................................ 46

18.

Startup Step(s) for Oscillator .................................. 48

19.

General Recurrences for Oscillator .......................... 48

C-I

Integrators Via Log Approximation ......................... 60

C-2

Tunings Independent of "wT"............................... 63

C-3

Differentiation of a Cosine Wave, "Amplitude/Phase" ........ 64

C-4 C-5

Differentiation of a Cosine Wave ........................... 65 Differentiation of a Cosine Wave,wT<
C-6

Differentiation of a Cosine Wave, wT<«< .................... 67

C-7

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 Z-Transforms" 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 Z-transform 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 z-transform 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 z-transforms on a firm foundation using distribution theory,

b) Determine when modified z-transforms 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 z-transforms will be presented. 5

That the Laplace transform of the Dirac delta* is a constant and the Laplace transform of the n-th 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) e-st 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) e-Stdt

(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(i-t) u(1-t) dt,

(3)

f(-t)

(4)

Co

=

f

g(t)

u(t) u(r-t) dt.

The term, "u(t) u(r-t)," will be recognized as the unit pulse from 0 to r-. Equation (4) becomes "1

g(t) * f(t) - f g(t) f(t-t) 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(1-t) 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(nT-t) dt

=

(7)

g(nT)

If one has the sequence 151

f Wt

(8)

f(t) 5(nT'-t)

n-0

its Laplace transform is

n-0

Taking

z

(10)

e-S

one has

n-0

4

7

The z-transform 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(nT-t) u(nT-t) 'n=-cx>

f(kT)

u(t) 6(k-t)] dt

kff-0

COo

(16)

From the properties of the Dirac delta, [6,71 and Appendix B, 8

-a-

(17)

g (iFr-kT) f (kr)

.z

*

gtnr-kT) z

(17

(18)

t(kT) z

n-O) ,<(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 Ragazzini-Zadeh 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 z-transforms 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(nT-t) u(nT-t) 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 Ragazzini-Zadeh identity. The mean value theorem of the integral calculus [5] guarantees there is some "k such that (k+1)T

f

g(t) f(nT-t) dt

=

T g(kT+ 6 kT) f(nT-kT-6 T)

(22)

kT

and Equation (21) may be rewritten n-i 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 n-1 T g(kT-6T) f(nT-kT-6T) xznI k=O n0O

(25)

and as before -1k

OD

Ty

nk f(nT-kT-6T) z"

g(kT+6T) z

(26)

n-0 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 (,vT-Sr) z-f

(-6T

(27)

and, finally, from the definition of the modified z-transform

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.

R-K 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 z-transforms 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) (1-z)

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 +

E-0

---(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(l-z 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)

-"1-z(1-aT)

which would amount to using a (I/0) Pade' approximation for e T,

-aT

e

-1-aT,

(52)

in Equation (50). Dividing Equation (51) by (50) one has the ratio*

(aT

1_-aT aT aT )1-ze)

*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 (n-0). 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(1--1) +- (--,

,-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'r-i) i z" Y" f k-0 kT 1n)O A"

n-1

, it=

n-

)

T' g(kT) T

f(nT-kT)

(54)

n=O k=O

As before, adding and substracting...

=

' in Tz 11=0

k-0

k() g

f(nT-kT)-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 z-transform,

Z [(s) f(s)]= T g(z) (f(z)-f

,

(57)

Eular convolution. If the modified z-transform is used (2], for a zero order hold one has oo I

n-1 (n g(nT+T) f(nT-kT-T) + (1-l) g(nT) f(nT-kT)1

z T

n=O

(58)

k=O

and then, as above,

co n

T zn I g(kT) f(nT-kT)n=O kf[O

Znfrio(1)

T

Zn g(nT

f(a) + (l-2.) f

go

n-O

n=O

(59)

and finally,

Z[g(s) f(s)]

T g(z) f(z)-T[Tg

o

19

f(z) + (i-n) g(z) f]

(60)

Eular tunable convolution. Proceeding on with the trapezoidal integrator,

Z

Tn-1

2

n.0

(61)

[g(kT+T) f(nT-kT-t) + g(kT) f(nT-kT)] k0O

and after manipulating indicies

x T

g(kT) f(nT-kT)

g0f(nT) + g(nT) f 1(62)

-

and rearranging 00 nk

n-k

g(kT) z TY n-0 k=O

T~ [o iz

f(nT-kT)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 z-transform and first order hold (Table 1) one has

jz n=()

+

~'

I g(kT+2TI)

f (n1-k1*-2T')

k=0) I

(1 + 2ri

2rK2) g (kT+T)

-

f (nT-kT-T)

+

( 1-2 r + i2 ) g(kt) f(nT-kT]

(66) Changing summing indices

z i

11-

2

nj,,2'+]

g~kT)f(nT-kT)

+

~n-T

(1+1j2y )Y+kT

k=-'

k-i

n-i +(1-21j+rj 2 )j g(kT) f (Tk) k-O

(67)

and adding and subtracting the necessary terms, nn12 ,

TY z' Y g(kT) f(nT-kT) n=U (k=U-(1Og

-

(

+ 2n + q2g

2r

1

+

2

f(nT)

[gn

+T

~

f(-T)

T

-g(T)

f(nT-T)1I

(68)

21

From the definition of the z-trangform (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 k1---z

3 \1 -z7/(3

that is,

2n .n- 1

-3 i- 2

i. )J+ 23 T 6- .n + :.n-1

' n-1/2

(74)

Note that Equations (73) and (74) are a linear combination of trapezoidal integration, and Mean Value integration with 6 = i 2. Therefore R-K 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 z-transform 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 (1-ni)

f(z) +

4

g(T)

g(z,1/2)

f)

fo)]

g0] f (Z) + 9(z) [f (Z)

-(1+2r1-n')

+ ri2[g(z,1) f(-T)-

R-KC.

(f (z)

-

(1-2rri)]

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..(1-Z)

EC E..(1-z)

T2z(l+z)

2

zT2Z2 2

Tz (1-Z) 2

T.C.

RK.Ti

RK.(1-Z)

2

24

=AND

3 q 2 T z(1+6z+z )

(i-z) 4

2(1-z) 3

8

(1-i) 3

Tz(1+z) 2 (1-i)4

T2z(1+z) 2 (1-z) 3

T3z(1+2z+z2 ) 4 (1-0)4

T2Z(1+Z) 2 (1-Z) 3

T3z(1+4z+z2) 6 (1-z) 4

TABLE 4. RATIOS OF APPROXIMATE SOLUTION TO THE EXACT FOR A CONSTANT. RAMP AND PARABOLIC INPUT TO A SINGLE INTEGRATOR

ts)

N.V-C

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 z-transform of Equations (80) and (82),

(83)

f(z) 1-z 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 z-transform 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

"e-aT/,," 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

)

e-aT/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 R-KC. 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 (I-Z)(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(1-2)(1-ge -

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

R-C

(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 R-KC 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

R-KC.=

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 non-zero initial conditions, one should proceed from the differential equation using n-1

nIs X(s)

L[X-

-

I

--

IX

n--i

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 z-transform of Equation (99),

X(z)

-

Z~()+

x0 Z(I)

(100)

There is no difficulty in taking the z-transform 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)

(1-r)z] X(z)

-

Tr

1-z

T.C.

R-KC.

A+ 2 (LtZ)j T-z

o 1 -z

-

T2 1-z X

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

-

1/2)- X(-T/2))4 (X(z)

X)]

and substituting the definition of the (modified) z-transform, one has

SX(nT) zn - X(nT) z n-0

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

(1-z) X(z)

=T

(X(z)-k I 0

+

X

(105)

0

becomes

00 zn 1 X (nT)z n-0

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)

Xn-1 + 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 Xn-1,

Xn = Xn-

(112)

n>O

'\Eular Integration.

Proceeding in a similar manner would produce Table 10. It might appear at first that MVC and R-KC 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)]

1-z

and T 6(1-z) [(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 kn-l/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

R-KC.

X

= X

B.

=X

n

n-1

i

+ 1(X 2

+ T

n

+ (i-rn)

+ k

I-2ri)

+ -6 [X

X 11

n-1

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 z-transform 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 z-transform 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 z-transform and equating coefficients of like powers of "z", X0 . X0 ,

(121)

X1 -X 0 +T

(122)

Xn

2 Xn-I

(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 T-0_z z { 2

_ R0

-(z)

Tz) 2 [I(z)

-

R 1/21

G -Z)

R-KC.

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 1-1 - Xn-2 r x + T2 r

X1 - X + T X = 2 X n n- 1 X

T.C.

=X + T 1 0

X

R-KC.

+-1

-o x

2n

n-/2 + Rn-] 3

2 ""

+T n-2

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

n-1 +

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

W-a 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 z-transform 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 lead-lag, 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 1-z 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. (1-z e

aT R-KC.

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

Xn-1+ T e

X..C n -e

Xn . e - a T Xn_ + T g

E.C. C.T

e- aT e-aT

-. n

R-KC.

n-in-i

l

gn-I/

x

=aTx

n =e

[g +

n

T 6

n

43

+4 e -aT/2 T

gn- 1]

2 nn-]

1

+e -aT • n-1J _

and solving Equation (132) for "X(s)" (s+b) g(s)

-

g0+

(133

X(S)

(133)_________

taking the z-transform of Equation (133),

X (z)

-

r/ +

(

+ (X0 -80 ) z+.(14

Noting that s +b -1 s+b

(a-b) 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 "(a-b)". Table 15 contains the recurrences. If the lead-lag 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 LEAD-LAG FILTER,

TABLE IS.

aT T (a-b) eg

M.V.C

Xn =

E.C.

Xn ' gn +

T.C.

Xn = (1-T(a-b)/2)

-

aT[X_

+

aT

n--2

In-i- gn-

(I+T(a-b))gn-1]

1-

(+T(a-b)/2)gn]

eaT[xn-

g+

-aT/2

R-KC.

+e

-

aT

gn-1/2

(a-b)e

gn-3

. n=-6(a-b))

(a-b))

-(I+

I+

TABLE 16. RECURRENCES FOR A LEAD-LAG FILTER,

M4.V. C.

x

E.C.

X

T (a

-

a

(

n

n

b

-

b) e-aTI 2

b

g

1

I nl

(1-T (a-b)) gn +e n

)+

gn-1/2

-

-

b

Y-n = "b a n ++ e aT[x n-1 - b( T a1 1n+ T(

X

T.C.

(

R-C

R-KC.

Xn

a-

(I-

(a-b)

+e aT [X

a( b

=

n-I-

2aT (a-b) e

)

gn

3

-

+

)

(a-b) 6(-4

45

n-1

g n-lj1 -b1n g

(a-b)) gn + eaTn-I - A

- A

eaT Exn~ a

gn-1

(I+T (a-b)) g aT / 2

n-1/2

s/a

~

1

gn-1

in Table 16 lot

(a-b)

-

(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

R-KC.

-aT

a (1-aT)gn+ eaT (X

a (1-

Xn = a(1

+e-aT (Xn 1

-

g n-

1

an

-a(1+aT/2)

2a2 Te-aT/2

3

a (1 + aT/6) gn-1

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

(1-C2)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.

+ Te-aT/

2

sin wT/2

+

e -aT sin wT w

0....

g

g1/2

E.C.

+0

T.C.

+ T T

-aT/2

+ -j e

R-KC

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 Xn-1

M.V.C.

wT/2 [2

sin

aT 2

- aT a

/

2a T

sioCu12)

%-2

['n-1/2

.... +

e'a -aT gn-3/]

sin wT 8nl U)

T.C. T.C

+ Te-T sin WT

+ T + 2

e-aT/2

gn-1

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 Ragazzini-Zadeh 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 lead-lag 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 Z-transforms. 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:

(B-i1)

nT)dt - f(nT)

-

60

where f(t) is a well-defined funiction at t -nT. 8(t) N(at)

(B-2)

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

(B-4)

(B-5)

nT)

-

f(nT)B(t

r)b(t

-

nT)dt

-

-

-

nT)

-

8(.r

-

(B-6) T)

(B-7)

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 z-transform

UZS X(s))

k(z)

-X(O)

.(C-2)

Since e St

z

(C-3)

it follows that S

=

Qn(1/z)

(C-4)

and Equation (C-2) becomes

k(z) =

(z-Qn(1/z)

X(s))

X(O)

.(C-5)

Applying the Ragazzini-Zadeh identity, Equation (C-5) becomes

i(z)

=1Qn(i/z)

X(z)

-

X(o)

.(C-6)

58

One of the possible series for the logarithm is

'n)

(1

=

z)

-

n=1

(C-7)

n

and substituting into Equation (C-6) and equating coefficients of like powers of z: --

(C.8)

0

0 SX

l

1

I X

- X

T

'(C-9)

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)

(C-10)

)-(x-Xl)

2T

X 1 )]

2-

-

x1 )

-

(X

-X

(C-Il)

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.,"

n1-z)

n

Z)

(C-12)

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

(C-13)

, n>l1

-

(C-14)

,

-2

+9

n-i

11 Xn n

l*ll2

-18X

3IX

,

gn-2 O

X.

T

--

n >0

Xn-3

Xn-2 6T

'

n

(C-15)

2

respectively. Of course, integrators may be developed in a similar fashion (Table C-i).

TABLE C-1

INTEGRATORS VIA LOG APPROXIMATION

-n(i/zf

z

Rectangular

T

l-z

l-z

11 +A

A1

2 k --

2i -+z

T_ +

CLASSICAL NAME

Eular

T z i-z

z

1-

-

1 = T s

fn(l/z)

1

.1-z +

"

.

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 (C-13) 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 n-lnT

xn-i

-

0O
(-

-(n-1)T

If in Equation (C-18) Ti =

(C-19)

112

then in Equation (C-17)

sn

wT/2 wT/2

(C-20)

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(C-21)

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

1-y)cos (n

+ (2y

1(

-

-oWT)

I~T(C-22)

cos nw]j

(C-23)

and finally

tan rTwT(2y -1)

tan

(C-24)

wT/2

Therefore

-

I

tan

-12-1 ta

ta

WT/ T2(C-25)

62

or tan r :T"

,v]/(1

t1

1/21

)

1 /

(C-26)

In small angles, "i' T. - i T

-

12 1

,

(C-27)

or 1/2 + rY

For any

"L,

(C-28)

T,"

TABLE C-2 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 (C-13), Tables (C-3) and (C-4) would resi It. The only difference between the tableL is one of point of view. When small angle approximations are warranted, Tables (C-5) and (C-6) 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 (C-7), 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

E-30

V-44

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

1

T z

DIFFERENTIATION OF A COSINE WAVE, wT - 7

tn~l)-w

sin nwT

Tz (y = 1/)

0W2n i~nT-r2

-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. R-K(N) CONVOLUTION

69

Simpson's 3/8 integration rule is

X(z)

=k

3T (1 + Z)O (z) + X0 -z

(D-1)

and the ratio for a sine wave input is

(D-2)

[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 Runge-Kutta integration is X(O)

1 1+z1/3 )3.

(D-3)

+ -

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 (D-3) may be written

X(Z)

=~

[z)

(T

I.

[32

+

3 [Tz2/3]

()+XO

As with R-K(3), Equation (D-5) 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

(D-6)

R-K(4) Convolution.

Higher order convolutions, R-K(N)C, could be developed in this fashion using trapezoidol integration/convolution and mean va'ue integration/convolution.

71

REFERENCES Pennsylvania: Addison-Wesley 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.

77-12, US Army Missile Research and Development Command, Redstone Arsenal, Alabama, 5 May 1977.

9r Atens 9. Halijak C.A.. t Sanpled-Iata S and Genratmg Fun'tions. NASA CR-61349. 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 T-78-75, 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 TD-77-14. US Army Missile Research and Development Command. Redstone Arsenal, Alabama. 3 June 1977.

6. Jury. E.L., Theory and Application of the Z-Transform 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

DRSMI-LP. Mr. Voigt DROMI-T. Dr, Kobler DRDMI.TE1 DRDMI-TO. D. Cifiax DRDMI-TD1 DRDMI-TDD1 DRDMI-TDK1 DROM(-TOS ORDMI-TDR1 DRDMt-TDF. V. Grimes1 DROMI-TOD. R. Dickson DRDMI-TL DRDMI-TKI DRDMI-TB DOMI-TBD DROMI-TI (Reference Copy)I DROMI-TI (Record Copy)

I I

73

I

100 I 5

3 I

Loading...

Trapezoidal Convolution Revisited - Defense Technical Information

IN,. 0 TECHNICAL REPORT T-78-66 U.S. ARMY TRAPEZOIDAL CONVOLUTION REVISITED: R-K CONVOLUTION OR THE DIGITAL SIMULATION MISSILEOF CONTINUOUS RESEA...

2MB Sizes 1 Downloads 0 Views

Recommend Documents

Annual - Defense Technical Information Center
Sep 30, 1996 - Patients. Other nominees were: CPT Bradley F. Schwartz, Department of Surgery, Urology Service, for his r

mhmmmmmhl - Defense Technical Information Center
Nov 5, 2016 - I: 1 -'H'.S. OL. 69 291~419 23 4" 1,'. ' 0. 6c. 0-. C. ClA. EY-2345. 70. 311'. 3". OE. 4. 0I. 2 .'- 70a24d

(CCISs) and Information Technology - Defense Technical Information
The Institute for Defense Analyses (IDA) has completed a review of technical standards applicable to future information

Information in Warfare - Defense Technical Information Center
STRATEGIC. APPRAISAL. DISTRIBUTION STATEMENT A. Approved for Public Release. Distribution Unlimited. Information in Warf

a394343.tiff - Defense Technical Information Center
exploitation tools such as computer Viruses, Trojan horses, WOrms, logic bombs, and eavesdropping Sniffers that can dest

a328974.tiff - Defense Technical Information Center
ATI-136615. Visual acuity tasks in a submarine. D. Farnsworth and F. L. Dimmick. 1951. ATI-155137. Oxygen consumption an

Job Enrichment - Defense Technical Information Center
(a) a general review and evaluation of job enrichment and Its related ... between job satisfaction and career decisions,

Complex Adaptive Systems - Defense Technical Information Center
Apr 30, 2014 - Each agent has the ability to make purposeful decisions in pursuit of a desired ... Purposeful behavior i

20000405 035 5Ä - Defense Technical Information Center
Aug 30, 1977 - JPRS 69708. 30 August 1977. REFERENCE AID. ABBREVIATIONS IN THE: LATIN AMERICAN PRESS. This report consis

Family Law Guide - Defense Technical Information Center
Jun 11, 1998 - FORM OF AGREEMENT 3-19. A. Introduction and Directions for Use 3-19. B. Standard Provisions 3-19. C. Samp

720p.BluRay | La Nonne vf | Vengeance: A Love ...