3.2 Optimum Second Order Recurrence Formulas

3.2.1 Basic Assumptions

It was convenient in the treatment of first order simulation to employ a simple mnemonic subscripting scheme, using p and c to denote "past" and "current" values of the time dependent variables. However, when dealing with higher order response the basic recurrence formulas normally involve more than two time states, and it is expedient to adopt a more general notation. We will represent any time dependent variable q by an ordered sequence of discrete values

... q_i-2 , q_i-1 , q_i , q_i+1 , q_i+2 , ...

where each of these discrete values is understood to be the true value (or at least the best available estimate) of the variable q at the times

... t_i-2 , t_i-1 , t_i , t_i+1 , t_i+2 , ...

respectively. These values of time form an ordered sequence such that for all n

t_n+1 � t_n = Dt

where Dt is defined as the execution period of the simulation. By convention, the subscript "i" denotes the "current value", and variables subscripted with i-1, i-2, and so on are referred to as "past values". The subscript n is used to denote the general term of the sequence.

For purposes of a simulation algorithm, it is assumed that all values of the input sequence x_n are known up to and including the current value x_i. It is also assumed that all values of the output sequence y_n are known up to but not including the current value y_i. The specific task of a second-order simulation is to compute y_i given that the continuous variables x and y are related according to equation (3.1-1). For real-time simulations, we will assume that the time required for the computation is negligible relative to the sampling interval Dt.

Even though all past values of x and y are theoretically known, a practical algorithm will retain only a limited number of them. To determine how many past values are necessary, we note that equation (3.1-1) includes the second derivative of x, and that a minimum of three values are required to provide an estimate of the secnd derivative of a function. Therefore, in addition to x_i, the algorithm must have access to the past values x_i-1 and x_i-2. Equation (3.1-1) also includes the second derivative of y, so we anticipate that two conditions on y must be specified to determine the solution. For this purpose we will use the past values y_i-1 and y_i-2.

The solution of equation (3.1-1) requires a continuous definition of x(t), so the algorithm must interpolate between the three available discrete values. However, unlike the first-order case, the optimum interpolation here is not obvious. Since three points can only provide a single estimate of the second derivative, it might seem reasonable to assume a function x(t) wirth a single, constant, second derivative during the double interval from x_i-2 to x_i. This would imply that x(t) can be represented by a second-degree polynomial passing through the three given values of x.

However, there are two possible objections to this approach. First, it introduces an inconsistency of interpolation on successive executions of the algorithm. For example, when y_i is being computed, the input function between t_i-2 and t_i is assumed to be the second-degree polynomial that passes through x_i-2, x_i-1 and x_i. But on the next pass, when y_i+1 is being calculated, the input function between t_i-1 and t_i+1 is assumed to be the second-degree polynomial that passes through x_i-1, x_i, and x_i+1. Hence the function between t_i-1 and t_i has been represented by two different polynomials on successive passes.

The other possible objection to the use of second order interpolation concerns the response to "step" input functions, as illustrated in Figure 4.

Figure 4

Second-order interpolation over-shoots a step change in x by 1/8 of the actual change in x. This may produce a slight amplification of any "noise"that is present in the x_n sequence, particularly because in digital systems all changes are really step changes at the limit of resolution.

An alternative is to use point-to-point interpolation, which ensures a unique interpolation for each interval, and precludes overshoot. However, it also increases the tendancy to underestimate the amplitude of actual input oscillations (recall the discussions of pre-warping in Sections 2.3.3 and 2.3.5). Also, with regard to uniqueness, it must be recognized that consistency does not imply accuracy. It is likely that in most cases both second-order interpolations for a given interval are better representations of the true function than the linear interpolation. Furthermore, since the second derivative of x appears explicitly in equation (3.1-1), we must recognize that, using point-to-point interpolation, the second derivative is zero across each interval and infinite at each breakpoint.

Based on the preceding remarks, it appears that no single interpolation scheme is best for all applications. An important consideration is the "quality" of the input signal. When the input is sampled from a smooth continuous function, as illustrated in Figure 5a, then a second-order fit will normally be more representative than point-to-point interpolation. However, when the input is ampled from a "coarse" digital representation as shown in Figure 5c, then the second derivative implied by three consecutive samples is unreliable, and point-to-point interpolation is probably preferable. In the intermediate case shown in Figure 5b the preference is less clear, but presumably the optimum representation cannot be outside the range between point-to-point interpolation and second-order interpolation. Therefore, our approach will be to derive the exact recurrence formulas for these two limiting cases.

Figure 5

It is shown in the next section that the ambiguity of interpolation results in only a single degree of freedom in the exact recurrence formula for a given differential equation, so the range of reasonable solutions is fully delimited by these two cases.

3.2.2 General Second-Order Recurrence Formulas

As discussed previously, the second-order recurrence formula computes y_i as a function of y_i-2, y_i-1, x_i-2, x_i-1, and x_i. As the general form we propose a linear combination of these five values:

where the coefficients c₁ through c₅ are all functions of Dt and of a₁, b₁, a₂, and b₂ from equation (3.1-1).

It may appear that equation (3.2.2-1) has 5 degrees of freedom, but because of the relative nature of the response, it�s clear that the five coefficients must satisfy the identity

c₁ + c₂ + c₃ + c₄ + c₅ = 1

which removes one degree of freedom. Furthermore, if we define

Dx_n = x_n � x_n-1DDx_n = Dx_n - Dx_n-1

then equation (3.2.2-1) can be written

y_i = (c₁) y_i-2 + (c₂) y_i-1 + (c₃ + c₄ + c₅) x_i-2 + (c₄ + 2c₅) Dx_i-1 + (c₅) DDx_i

This shows that among all possible recurrence formulas based on reasonable interpolation schemes for a given differential equation and Dt, there is really only one degree of freedom. This is because when x_i-2, x_i-1, and x_i are colinear versus time, the optimum interpolation is clearly linear, so we can assume that all reasonable formulas are identical for this case. Since DDxi is zero in this case, the last term of (3.2.2-2) drops out, and the remaining terms must be identical for all reasonable formulas. Therefore, c₁ and c₂ are unique determined b this requirement, as are the sums c₃ + c₄ + c₅ and c₄ + 2c₅. Thus the general recurrence formula written in the form

has only one free coefficient, L₅, the value of which depends on the type of interpolation used when the xn are not colinear. The other four coefficients are fully determined by the linear case. Therefore, comparing equations (3.2.2-2) and (3.2.2-3), we see that any reasonable set of "c-coefficients" must satisfy the equations

c₁ = L₁

c₂ = L₂

c₃ + c₄ + c₅ = L₃

c₄ + 2c₅ = L₄

Taking c₅ (which equals L₅) as the free variable, we can solve the last two of these equations for c₃ and c₄ to give

c₃ = L₃ � L₄ + c₅c₄ = L₄ - 2c₅

This shows that, within the class of reasonable coefficients sets for a given differential equation and Dt, the values of c₃, c₄, and c₅ are mutually co-linear. Hence, not only is there just a single degree of freedom among these sets, but an interpolation between any two sets varies all the coefficients in an equal proportion to their total differences.

3.2.3 Optimum Recurrence Coefficients For Second-Order Interpolation

In this section we derive the coefficients of equation (3.2.2-1) assuming that the input function x(t) is the second-degree polynomial that passes through x_i-2, x_i-1, and x_i. If we set

t_i-2 = -Dtt_i-1 = 0t_i = Dt

then a second-order curve fit of the two points gives

where

Differentiating equation (3.2.3-1) gives

Therefore, equation (3.1-1) can be written

We propose a particular solution of the form

which gives

Substituting these expressions for y and its derivatives into equation (3.2.3-2) and equating coefficients of like powers of t, we can solve for k₁, k₂, and k₃. If we define the dimensionless parameters

then we have

Now, the homogeneous form of equation (3.2.3-2) has the general solution

where

and the coefficients R₁ and R₂ are to be determined by two conditions on y(t). (Classically this is the general solution, i.e., a "basis" for all possible solutions, only if r₁ � r₂. If r₁ = r₂ the expressions for R₁ and R₂ below are singular. However, these singularities are removeable when the solution is expressed in discrete-time form. The solution derived from (3.2.3-4) can be shown to be identical whether or not r₁ = r₂.) Adding equations (3.2.3-3) and (3.2.3-4) (the particular and complimentary solutions respectively) gives the general total solution of equation (3.2.3-2)

We will determine the coefficients R₁ and R₂ by imposing the following two conditions on y(t): y = y_i-1 at t = t_i-1, and y = y_i-2 at t = t_i-2 = -Dt. Inserting these two conditions into equation (3.2.3-5) gives

and

Solving these two equations simultaneously for R₁ and R₂ gives

Substituting these expressions for R₁ and R₂ into equation (3.2.3-5), solving for y_i at t = t_i = Dt, and simplifying, we have

where g and d are dimensionless parameters defined by

Substituting the prior expressions for k₁, k₂, and k₃ into equation (3.2.3-6) and combining terms by x_n and y_n we arrive at the final recurrence formula

For a given discrete time increment Dt, the recurrence coefficients are given by

As required, the sum of c₁ through c₅ is identically equal to 1. Also, since this solution qualifies as a "reasonable" recurrence formula (in the sense that it gives x(t) as a linear function when x_i-2, x_i-1, and x_i are co-linear versus time), we can now determine the first coefficients of equation (3.2.2-3) using equations (3.2.2-4) as follows

(3.2.3-8)

Before proceding to the next section we will derive three more recurrence formulas that will prove useful. They are based on equation (3.2.3-5) with the constants R₁ and R₂ determined by the two conditions y = y_i-1 at t = t_i-1 = 0 and dy/dt = at t = t_i-1 = 0. Inserting the first of these conditions into equation (3.2.3-5) gives

To apply the second condition we first differentiate equation (3.2.3-5) as follows

Inserting the second condition, we have

Solving equations (3.2.3-9) and (3.2.3-11) simultaneously for R₁ and R₂ gives

Substituting these expressions for R₁ and R₂ into equation (3.2.3-5), solving for y_i at t = t_i = Dt, and making the appropriate substitutions for k₁, k₂, and k₃, we arrive at the recurrence formula

If we define the dimensionless parameters

we can express the coefficients f₁ through f₅ in equation (3.2.3-13) as

Similarly, we can solve equation (3.2.3-10) for (dy/dt)_i at t = t_i-2 = Dt with R₁ and R₂ given by equations (3.2.3-12). This leads to the recurrence formula

In terms of the dimensionless parameters

we can express the coefficients g₁ through g₅ in equation (3.2.3-14) as

Equations (3.2.3-13) and (3.2.3-14) can be used as a two-stage recurrence scheme for generating both y_n and .

Another useful recurrence formula is found by solving equation (3.2.3-5) for y_i-2 with R₁ and R₂ given by the last two of equations (3.2.3-12) and setting t = t_i-2 = -Dt. This leads to the formula

(3.2.3-15)

In terms of the dimensionless parameters

the coefficients h₁ through h₅ can be expressed as

3.2.4 Optimum Recurrence Coefficients Based on Point-To-Point Interpolation

In this section we derive the coefficients of equation (3.2.2-1) with the assumption that the input function x(t) is linear between successive values of x_n. As mentioned previously, this implies that the second derivative is zero during each interval and infinite at each instant tn, where the slope changes discontinuously. Therefore, we will first determine the response to equation (3.1-1) to a discontinuous change in the derivative of x(t). Our approach will be based on the solution of the following problem:

Given the values of x_i-1, , y_i-1, , and Dt, find y_i and assuming the function x(t) has a constant second derivative in the interval from t_i-1 to t_i.

Equations (3.2.3-13) and (3.2.3-14) provide the means of solving this problem if we can infer the appropriate values of x_i-2 and x_i from the given information. The function x(t) is uniquely determined in the interval from t_i-1 to t_i by the values of x_i-1, , , and Dt, and the requirement of a constant second derivative. Therefore, the value of x_i is also uniquely determined, and is given by

Further, if we recognize that the only purpose of x_i-2 in equations (3.2.3-13) and (3.2.3-14) is to determine the function x(t) in the interval from t_i-1 to t_i, which it does by assuming a constant second derivative over the entire double interval from t_i-2 to t_i, then it is clear that the appropriate value for x_i-2 is to be found by extrapolating the function x(t) back to t_i-2. On this basis, it can be shown that

Substituting the two preceding expressions into equations (3.2.3-13) and (3.2.3-14) gives the solution to the problem:

These equations allow us to compute the effect on y and of a change in from to occurring uniformly during an interval Dt. In the limit as Dt goes to zero these equations reduce to

With this result we can proceed to determine the general recurrence formula based on point-to-point interpolation of the input function. By equations (3.2.3-15) and (3.2.3-13) we can write

(3.2.4-2)

and

(3.2.4-3)

where

The synthesized values of x can be expressed as

X_i-2 = 2x_i-1 � x_iX_i = 2x_i-1 - x_i-2

Also, since the change in the derivative of x(t) at t = t_i-1 is

we have, by equation (3.2.4-1), the relation

Making the indicated substitutions for X_i-2, X_i, and in equations (3.2.4-2) and (3.2.4-3), and then eliminating from the resulting equations leads to the general recurrence formula based on point-to-point interpolation

y_i = c₁ y_i-2 + c₂ y_i-1 + c₃ x_i-2 + c₄ x_i-1 + c₅ x_i

where, if we define the dimensionless parameter

the coefficients c₁ through c₅ are given by

As required, the sum of c₁ through c₅ is identically equal to 1. Also, substituting these coefficients into equations (3.2.2-4) gives the same values for L₁ through L₄ and do the coefficients based on second order interpolation (see equations (3.2.3-8)). This is in accord with the assertion that all reasonable coefficient sets give identical results when x_i-2, x_i-1, and x_i are colinear versus time.

3.2.5 Discussion

In the preceding sections no distinction was made between the cases when the roots r₁ and r₂ of the characteristic equation are real and when they are complex. In either case the dimensionless parameters in terms of which the recurrence coefficients have been expressed are always purely real, as is easily verified using the basic exponential and trigonometric identities for complex variables. Table 2 presents a summary of these dimensionless parameters, including the equivalent trigonometric definitions where appropriate. In this table the parameter u and v are defined as

Notice that r₁ = u + iv and r₂ = u - iv.

Table 2

The use of the recurrence formulas developed in thepreceding sections can best be illustrated by an example. Consider a function with a dynamic characteristic described by the differential equation

In this example we have

a₁ = 9 sec²b₁ = 1 seca₂ = 8 sec²b₂ = 3 sec

If we arbitrarily select an execution period of Dt = 1 sec, then by the equations in Table 2 we have

a = 3b = -2g = 0.8948393168

Since the function is underdamped, we use the trigonometric form of the equation for d to compute

d = 1.790648545

We can now use equations (3.2.3-7) to compute the coefficients of the general recurrence formula based on second order interpolation of the input function:

y_i = c₁ y_i-2 + c₂ y_i-1 + c₃ x_i-2 + c₄ x_i-1 + c₅ x_i

We can verify that these coefficients sum to 1, as required. To demmonstrate the operation of this recurrence formula we will compute the response to three different input functions:

In the first case it�s convenient to use the form of the recurrence given by equation (3.2.2-2). setting x_i-2 = 100 and Dx_i-1 = DDx_i = 0, to give the formula

y_i = (-0.894839) y_i-1 + (1.790648) y_i-1 + (0.10419) 100

With the initial conditions y₀ = y₁ = 0, the subsequent values of y_n generated by this formula are plotted in Figure 6.

Figure 6

In the second case, DDx_i is still zero, but Dx_i-1 is a constant equal to 1 (since the slope of the function x(t) = t is 1). Therefore, equation (3.2.2-2) gives the formula

y_i = (-0.894839) y_i-2 + (1.790648) y_i-1 + (0.10419) t_i-2 + (0.417733)

where we have substituted t_i-2 for x_i-2. With the initial conditions y₀ = y₁ = 0, the subsequent values of y_n generated by this formula are plotted in Figure 7.

Figure 7

In the third case, neither ther first nor the second derivative of x(t) is zero, so we use the general form of the recurrence formula (3.2.2-1). Substituting for x_n the corresponding value function of t_n, we have

y_i = (-0.894839) y_i-2 + (1.790648) y_i-1 + (0.688398) [(1/4)t_i-2² + (1/2)t_i-2]

+ (-1.586147) [(1/4)t_i-1² + (1/2)t_i-1] + (1.001940) [(1/4)t_i² + (1/2)t_i]

With the initial conditions y₀ = y₁ = 0, the subsequent values of y_n generated by this formula are plotted in Figure 8.

Figure 8

Figure 9 shows the response to the same input function from the initial conditions y₀ = y₁ = 50.

Figure 9

In these examples the computed values of y_n fall exactly on the theoretical function y(t), because all the input functions are exactly representable by a second-degree polynomial. Any choice of execution period will yield the same accuracy. For example, if we had chosen Dt = 5 sec, the coefficients would have been different but the computed values of y_n would still have fallen exactly on the theoretical function y(t). The only consideration limiting the size of Dt in such cases is that the output of y_n must occur frequently enough to adequately represent the function y(t) for external purposes.

In practice, the continuous input function x(t) is usually specified by a sequence of discrete values from some external source, rather than computed from a given equation as in the preceding examples. It is the lack of a precise functional definition of x(t) that requires us to assume an interpolation function for the xn sequence. This, in turn, requires that Dt be small enough so that our interpolation of the x_n samples is sufficiently representative of the true continuous input function f(t). This restriction on Dt should not be confused with a restriction arising from the use of a non-optimum solution method, as discussed in Sections 2.3.5 and 3.3.5.)

We can also determine the optimum recurrence formulas for the dynamical system specified by equation (3.2.5-1) with the assumption of point-to-point interpolation of the input samples. Since the point-to-point formula is identical to the second-order-fit formula when x_i-2, x_i-1, and x_i are colinear versus time, they give identical results for the constant and ramp inputs i.e., cases (1) and (2) above. It is only when applied to the third, more general, case that the two assumptions give different results. To determine the exact point-to-point recurrence formula, we first compute the dimensionless parameters

e = 0.1111111111

m = 0.9290200587

s = -0.9469364981

l = 0.8437120471

using the equations in Table 2. Then by equations (3.2.4-4) we can compute the coefficients of the exact point-to-point recurrence formula

y_i = c₁ y_i-2 + c₂ y_i-1 + c₃ x_i-2 + c₄ x_i-1 + c₅ x_i

These coefficients differ hardly at all from those computed on the basis of second-order interpolation, so the predicted response to the input for case (3) is virtually identical to that given by the second-order-fit method, especially since every three consecutive samples x_n are nearly co-linear.

To conclude this section, we will develop the "standard form" of the second-order recurrence formula (analagous to equation (2.2.2-2) for first-order formulas), in terms of which the optimum second-order solutions can be expressed succinctly. If we define

Dy_n = y_n - y_n-1

then equation (3.2.2-2) can be written as

y_i = (c₁ + c₂) y_i-2 + c₂ Dy_i-1 + (c₃ + c₄ + c₅) x_i-2 + (c₄ + 2c₅) Dx_i-1 + c₅ DDx_i

Since the coefficients sum to unity we have c₃ + c₄ + c₅ = 1 � c₁ � c₂, so we can make this substitution to give

y_i = (c₁ + c₂) y_i-2 + c₂ Dy_i-1 + (1 � c₁ � c₂) x_i-2 + (c₄ + 2c₅) Dx_i-1 + c₅ DDx_i

which can be written as

y_i = y_i-2 + c₂ Dy_i-1 + (1 � c₁ � c₂) (x_i-2 � y_i-2) + (c₄ + 2c₅) Dx_i-1 + c₅ DDx_i

Thus we have the "standard form" for second-order recurrences

where

(3.2.5-3)

and the D coefficient is given by

for second-order interpolation, and by

for point-to-point interpolation. It�s worth recalling that the set of all reasonable second order recurrence schemes has just a single degree of freedom, which can be linearly parameterized as the D coefficient of the standard form.