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 Dt is defined as the execution period of the simulation. By convention, the subscript “i” denotes the “current value”, and variables subscripted with i1, i2, 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 secondorder simulation is to compute y_{i} given that the continuous variables x and y are related according to equation (3.11). For realtime 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.11) 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 x_{i}, the algorithm must have access to the past values x_{i1} and x_{i2}. Equation (3.11) 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_{i1} and y_{i2}. 

The solution of equation (3.11) requires a continuous definition of x(t), so the algorithm must interpolate between the three available discrete values. However, unlike the firstorder 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 x_{i2} to x_{i}. This would imply that x(t) can be represented by a seconddegree 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_{i2} and t_{i} is assumed to be the seconddegree polynomial that passes through x_{i2}, x_{i1} and x_{i}. But on the next pass, when y_{i+1} is being calculated, the input function between t_{i1} and t_{i+1} is assumed to be the seconddegree polynomial that passes through x_{i1}, x_{i}, and x_{i+1}. Hence the function between t_{i1} 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. 



Secondorder interpolation overshoots 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 pointtopoint 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 prewarping 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 secondorder 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.11), we must recognize that, using pointtopoint 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 secondorder fit will normally be more representative than pointtopoint 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 pointtopoint 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 pointtopoint interpolation and secondorder 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 SecondOrder Recurrence Formulas 

As discussed previously, the secondorder recurrence formula computes yi as a function of y_{i2}, y_{i1}, x_{i2}, x_{i1}, and x_{i}. As the general form we propose a linear combination of these five values: 

_{} 

where the coefficients c_{1} through c_{5} are all functions of Dt and of a_{1}, b_{1}, a_{2}, and b_{2} from equation (3.11). It may appear that equation (3.2.21) 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.21) can be written 

_{} 

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_{i2}, x_{i1}, 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.22) drops out, and the remaining terms must be identical for all reasonable formulas. Therefore, c_{1} and c_{2} are unique determined b this requirement, as are the sums c_{3} + c_{4} + c_{5} and c_{4} + 2c_{5}. Thus the general recurrence formula written in the form 

_{} 

has only one free coefficient, L_{5}, 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.22) and (3.2.23), we see that any reasonable set of “ccoefficients” must satisfy the equations 

_{} 

Taking c_{5} (which equals L_{5}) as the free variable, we can solve the last two of these equations for c_{3} and c_{4} to give 

_{} 

This shows that, within the class of reasonable coefficients sets for a given differential equation and Dt, the values of c_{3}, c_{4}, and c_{5} are mutually colinear. 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 SecondOrder Interpolation 

In this section we derive the coefficients of equation (3.2.21) assuming that the input function x(t) is the seconddegree polynomial that passes through x_{i2}, x_{i1}, and x_{i}. If we set 

_{} 

then a secondorder curve fit of the two points gives 

_{} 

where 
_{} 

Differentiating equation (3.2.31) gives 

_{} 

Therefore, equation (3.11) 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.32) and equating coefficients of like powers of t, we can solve for k_{1}, k_{2}, and k_{3}. If we define the dimensionless parameters 

_{} 

then we have 

_{} 

Now, the homogeneous form of equation (3.2.32) has the general solution 

_{} 

where 
_{} 

and the coefficients R_{1} and R_{2} 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_{1} ¹ r_{2}. If r_{1} = r_{2} the expressions for R_{1} and R_{2} below are singular. However, these singularities are removable when the solution is expressed in discretetime form. The solution derived from (3.2.34) can be shown to be identical whether or not r_{1} = r_{2}.) Adding equations (3.2.33) and (3.2.34) (the particular and complimentary solutions respectively) gives the general total solution of equation (3.2.32) 

_{} 

We will determine the coefficients R_{1} and R_{2} by imposing the following two conditions on y(t): y = y_{i1} at t = t_{i1}, and y = y_{i2} at t = t_{i2} = Dt. Inserting these two conditions into equation (3.2.35) gives 

_{} 

and 
_{} 

Solving these two equations simultaneously for R_{1} and R_{2} gives 

_{} 

_{} 

Substituting these expressions for R_{1} and R_{2} into equation (3.2.35), 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_{1}, k_{2}, and k_{3} into equation (3.2.36) 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_{1} through c_{5} 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_{i2}, x_{i1}, and x_{i} are colinear versus time), we can now determine the first coefficients of equation (3.2.23) using equations (3.2.24) 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.35) with the constants R_{1} and R_{2} determined by the two conditions y = y_{i1} at t = t_{i1} = 0 and dy/dt = _{} at t = t_{i1} = 0. Inserting the first of these conditions into equation (3.2.35) gives 

_{} 

To apply the second condition we first differentiate equation (3.2.35) as follows 

_{} 

Inserting the second condition, we have 

_{} 

Solving equations (3.2.39) and (3.2.311) simultaneously for R_{1} and R_{2} gives 

_{} 
(3.2.312) 

Substituting these expressions for R_{1} and R_{2} into equation (3.2.35), solving for y_{i} at t = t_{i} = Dt, and making the appropriate substitutions for k_{1}, k_{2}, and k_{3}, we arrive at the recurrence formula 

_{} 

If we define the dimensionless parameters 

_{} 

we can express the coefficients f_{1} through f_{5} in equation (3.2.313) as 

_{} 

Similarly, we can solve equation (3.2.310) for (dy/dt)_{i} at t = t_{i2} = Dt with R_{1} and R_{2} given by equations (3.2.312). This leads to the recurrence formula 

_{} 

In terms of the dimensionless parameters 

_{} 

we can express the coefficients g_{1} through g_{5} in equation (3.2.314) as 

_{} 


Equations (3.2.313) and (3.2.314) can be used as a twostage recurrence scheme for generating both y_{n} and _{}. 

Another useful recurrence formula is found by solving equation (3.2.35) for y_{i2} with R_{1} and R_{2} given by the last two of equations (3.2.312) and setting t = t_{i2} = Dt. This leads to the formula 

_{} 

In terms of the dimensionless parameters 

_{} 

the coefficients h_{1} through h_{5} can be expressed as 

_{} 


3.2.4 Optimum Recurrence Coefficients Based on PointToPoint Interpolation 

In this section we derive the coefficients of equation (3.2.21) 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 t_{n}, where the slope changes discontinuously. Therefore, we will first determine the response to equation (3.11) 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.313) and (3.2.314) provide the means of solving this problem if we can infer the appropriate values of x_{i2} and x_{i} from the given information. The function x(t) is uniquely determined in the interval from t_{i1} to t_{i} by the values of x_{i1}, _{}, _{}, 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_{i2} in equations (3.2.313) and (3.2.314) is to determine the function x(t) in the interval from t_{i1} to t_{i}, which it does by assuming a constant second derivative over the entire double interval from t_{i2} to t_{i}, then it is clear that the appropriate value for x_{i2} is to be found by extrapolating the function x(t) back to t_{i2}. On this basis, it can be shown that 

_{} 

Substituting the two preceding expressions into equations (3.2.313) and (3.2.314) 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 pointtopoint interpolation of the input function. By equations (3.2.315) and (3.2.313) 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 = t_{i1} is 

_{} 

we have, by equation (3.2.41), the relation 

_{} 

Making the indicated substitutions for X_{i2}, X_{i}, and _{} in equations (3.2.42) and (3.2.43), and then eliminating _{} from the resulting equations leads to the general recurrence formula based on pointtopoint interpolation 

_{} 

where, if we define the dimensionless parameter 

_{} 

the coefficients c_{1} through c_{5} are given by 

_{} 

As required, the sum of c_{1} through c_{5} is identically equal to 1. Also, substituting these coefficients into equations (3.2.24) gives the same values for L_{1} through L_{4} and do the coefficients based on second order interpolation (see equations (3.2.38)). This is in accord with the assertion that all reasonable coefficient sets give identical results when x_{i2}, x_{i1}, 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_{1} and r_{2} 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_{1} = u + iv and r_{2} = 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 Dt = 1 sec, then by the equations in Table 2 we have 

_{} 

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

_{} 

We can now use equations (3.2.37) 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.22). setting x_{i2} = 100 and Dx_{i1} = DDx_{i} = 0, to give the formula 
_{} 

With the initial conditions y_{0} = y_{1} = 0, the subsequent values of y_{n} generated by this formula are plotted in Figure 6. 



In the second case, DDx_{i} is still zero, but Dx_{i1} is a constant equal to 1 (since the slope of the function x(t) = t is 1). Therefore, equation (3.2.22) gives the formula 

_{} 

where we have substituted t_{i2} for x_{i2}. With the initial conditions y_{0} = y_{1} = 0, the subsequent values of y_{n} 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.21). Substituting for x_{n} the corresponding value function of t_{n}, we have 

_{} 

With the initial conditions y_{0} = y_{1} = 0, the subsequent values of y_{n} generated by this formula are plotted in Figure 8. 



Figure 9 shows the response to the same input function from the initial conditions y_{0} = y_{1} = 50. 



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 seconddegree 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 nonoptimum 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.51) with the assumption of pointtopoint interpolation of the input samples. Since the pointtopoint formula is identical to the secondorderfit formula when x_{i2}, x_{i1}, 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 pointtopoint recurrence formula, we first compute the dimensionless parameters 

_{} 

using the equations in Table 2. Then by equations (3.2.44) we can compute the coefficients of the exact pointtopoint recurrence formula 

_{} 

_{} 

These coefficients differ hardly at all from those computed on the basis of secondorder interpolation, so the predicted response to the input for case (3) is virtually identical to that given by the secondorderfit method, especially since every three consecutive samples x_{n} are nearly colinear. 

To conclude this section, we will develop the “standard form” of the secondorder recurrence formula (analogous to equation (2.2.22) for firstorder formulas), in terms of which the optimum secondorder solutions can be expressed succinctly. If we define 

_{} 

then equation (3.2.22) can be written as 

_{} 

Since the coefficients sum to unity we have c_{3} + c_{4} + c_{5} = 1 – c_{1} – c_{2}, so we can make this substitution to give 

_{} 

which can be written as 

_{} 

Thus we have the “standard form” for secondorder recurrences 

_{} 

where 
_{} 

and the D coefficient is given by 

_{} 

for secondorder interpolation, and by 

_{} 

for pointtopoint 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. 
