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

 

 

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

 

 

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

 

 

where Δt 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 xn are known up to and including the current value xi. It is also assumed that all values of the output sequence yn are known up to but not including the current value yi. The specific task of a second-order simulation is to compute yi 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 Δt.

 

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 second derivative of a function. Therefore, in addition to xi, the algorithm must have access to the past values xi-1 and xi-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 yi-1 and yi-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) with a single, constant, second derivative during the double interval from xi-2 to xi. 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 yi is being computed, the input function between ti-2 and ti is assumed to be the second-degree polynomial that passes through xi-2, xi-1 and xi. But on the next pass, when yi+1 is being calculated, the input function between ti-1 and ti+1 is assumed to be the second-degree polynomial that passes through xi-1, xi, and xi+1. Hence the function between ti-1 and ti 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.

 

 

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 xn 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 tendency 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’s clear 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 sampled 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.

 

 

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 yi as a function of yi-2, yi-1, xi-2, xi-1, and xi. As the general form we propose a linear combination of these five values:

 

 

where the coefficients c1 through c5 are all functions of Δt and of a1, b1, a2, and b2 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

 

 

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

 

 

then equation (3.2.2-1) can be written

 

 

This shows that among all possible recurrence formulas based on reasonable interpolation schemes for a given differential equation and Δt, there is really only one degree of freedom. This is because when xi-2, xi-1, and xi are co-linear versus time, the optimum interpolation is clearly linear, so we can assume that all reasonable formulas are identical for this case. Since ΔΔxi 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, c1 and c2 are uniquely determined by this requirement, as are the sums c3 + c4 + c5 and c4 + 2c5. Thus the general recurrence formula written in the form

 

 

has only one free coefficient, L5, the value of which depends on the type of interpolation used when the xn are not co-linear. 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

 

 

Taking c5 (which equals L5) as the free variable, we can solve the last two of these equations for c3 and c4 to give

 

 

This shows that, within the class of reasonable coefficients sets for a given differential equation and Δt, the values of c3, c4, and c5 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 xi-2, xi-1, and xi. If we set

 

 

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 k1, k2, and k3. 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 R1 and R2 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 r1 ≠ r2. If r1 = r2 the expressions for R1 and R2 below are singular. However, these singularities are removable 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 r1 = r2.) 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 R1 and R2 by imposing the following two conditions on y(t): y = yi-1  at  t = ti-1, and y = yi-2  at  t = ti-2  =  –Δt.  Inserting these two conditions into equation (3.2.3-5) gives

 

 

and

 

Solving these two equations simultaneously for R1 and R2 gives

 

 

 

Substituting these expressions for R1 and R2 into equation (3.2.3-5), solving for yi  at  t = ti = Δt, and simplifying, we have

 

 

where g and d are dimensionless parameters defined by

 

 

Substituting the prior expressions for k1, k2, and k3 into equation (3.2.3-6) and combining terms by xn and yn we arrive at the final recurrence formula

 

 

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

 

 

As required, the sum of c1 through c5 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 xi-2, xi-1, and xi 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

 

 

Before proceeding 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 R1 and R2 determined by the two conditions  y = yi-1  at  t = ti-1 = 0  and  dy/dt  =    at  t = ti-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 R1 and R2 gives

 

 

(3.2.3-12)

 

Substituting these expressions for R1 and R2 into equation (3.2.3-5), solving for yi at t = ti = Δt, and making the appropriate substitutions for k1, k2, and k3, we arrive at the recurrence formula

 

 

If we define the dimensionless parameters

 

 

we can express the coefficients f1 through f5 in equation (3.2.3-13) as

 

 

Similarly, we can solve equation (3.2.3-10) for (dy/dt)i  at  t = ti-2 = Δt  with R1 and R2 given by equations (3.2.3-12). This leads to the recurrence formula

 

 

In terms of the dimensionless parameters

 

 

we can express the coefficients g1 through g5 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 yn and .

 

Another useful recurrence formula is found by solving equation (3.2.3-5) for yi-2 with R1 and R2 given by the last two of equations (3.2.3-12) and setting t = ti-2 = –Δt. This leads to the formula

 

 

In terms of the dimensionless parameters

 

 

the coefficients h1 through h5 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 xn. 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:

 

 

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 xi-2 and xi from the given information. The function x(t) is uniquely determined in the interval from ti-1 to ti by the values of xi-1, , , and Δt, and the requirement of a constant second derivative. Therefore, the value of xi is also uniquely determined, and is given by

 

 

Further, if we recognize that the only purpose of xi-2 in equations (3.2.3-13) and (3.2.3-14) is to determine the function x(t) in the interval from ti-1 to ti, which it does by assuming a constant second derivative over the entire double interval from ti-2 to ti, then it is clear that the appropriate value for xi-2 is to be found by extrapolating the function x(t) back to ti-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 Δt. In the limit as Δt 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

 

 

 

where

 

The synthesized values of x can be expressed as

 

 

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

 

 

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

 

 

Making the indicated substitutions for Xi-2, Xi, 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

 

 

where, if we define the dimensionless parameter

 

 

the coefficients c1 through c5 are given by

 

 

As required, the sum of c1 through c5 is identically equal to 1. Also, substituting these coefficients into equations (3.2.2-4) gives the same values for L1 through L4 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 xi-2, xi-1, and xi are co-linear versus time.

 

 

3.2.5  Discussion

 

In the preceding sections no distinction was made between the cases when the roots r1 and r2 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 r1 = u + iv and r2 = u – iv.

 

Table 2

 

 

 

 

 

 

 

 

 

 

 

 

 

                       

 

 

 

 

 

The use of the recurrence formulas developed in the preceding 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

 

 

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

 

 

Since the function is under-damped, we use the trigonometric form of the equation for d to compute

 

 

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:

 

 

 

We can verify that these coefficients sum to 1, as required. To demonstrate 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 xi-2 = 100 and  Δxi-1 = ΔΔxi = 0, to give the formula

 

 

With the initial conditions y0 = y1 = 0, the subsequent values of yn generated by this formula are plotted in Figure 6.

 

 

In the second case,  ΔΔxi  is still zero, but  Δxi-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

 

 

where we have substituted ti-2 for xi-2. With the initial conditions y0 = y1 = 0, the subsequent values of yn generated by this formula are plotted in Figure 7.

 

 

In the third case, neither the 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 xn the corresponding value function of tn, we have

 

 

With the initial conditions y0 = y1 = 0, the subsequent values of yn generated by this formula are plotted in Figure 8.

 

 

Figure 9 shows the response to the same input function from the initial conditions y0 = y1 = 50.

 

 

In these examples the computed values of yn 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 Δt = 5 sec, the coefficients would have been different but the computed values of yn would still have fallen exactly on the theoretical function y(t). The only consideration limiting the size of Δt in such cases is that the output of yn 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 Δt be small enough so that our interpolation of the xn samples is sufficiently representative of the true continuous input function f(t). This restriction on Δt 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 xi-2, xi-1, and xi are co-linear 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

 

 

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

 

 

 

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 xn are nearly co-linear.

 

To conclude this section, we will develop the “standard form” of the second-order recurrence formula (analogous 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

 

 

then equation (3.2.2-2) can be written as

 

 

Since the coefficients sum to unity we have c3 + c4 + c5 = 1 – c1 – c2, so we can make this substitution to give

 

 

which can be written as

 

 

Thus we have the “standard form” for second-order recurrences

 

 

where

 

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 by the D coefficient of the standard form.

 

Return to Table of Contents