W. E. Milne

Bio: W. E. Milne is an academic researcher from Oregon State University. The author has contributed to research in topics: Exponential integrator & Differential algebraic equation. The author has an hindex of 10, co-authored 22 publications receiving 262 citations. Previous affiliations of W. E. Milne include University of Oregon & Bowdoin College.

Numerical Integration of Ordinary Differential Equations

Abstract: (1926). Numerical Integration of Ordinary Differential Equations. The American Mathematical Monthly: Vol. 33, No. 9, pp. 455-460.

Stability of a Numerical Solution of Differential Equations

Abstract: In 1926 Milne [1] published a numerical method for the solution of ordinary differential equations. This method turns out to be unstable, as shown by Muhin [2], Hildebrand [3], Liniger [4], and others. Instability was not too serious in the day of desk calculators but is fatal in the modern era of high speed computers. The basic cause of the instability in this particular method is the use of Simpson's rule to perform the final integration. Simpson's rule integrates over two intervals, and under certain conditions can produce an error which alternates in sign from step to step and which increases in magnitude exponentially. It is the purpose of this paper to show that the occasional application of Newton's “three eighths” quadrature formula over three intervals can effectively damp out the unwanted oscillation without harm to the desired solution.Let the given differential equation be dy/dx = ƒ(x, y), and let the step length for the independent variable x be denoted by h. The quantity s = h ∂ƒ/∂y plays a basic role in the analysis, for it may be shown that when Simpson's rule is used an error E0 at x = x0 is propagated through subsequent steps according to the second order difference equation (1 - sn+1/3) En+1 - (4sn/3) En - (1 + sn-1/3) En-1 = 0. (See e.g., Hildebrand [3], p. 206, Milne [5], p. 68.)While in general s is a variable, the special case where s is constant not only permits a simple analysis but also serves to explain the behavior in other cases. Cf. Hildebrand [3], p. 202. Accordingly we treat s as a constant and assume that our differential equation is dy/dx = Gy, the general solution of which, after n steps, is yn = Aens, where A is an arbitrary constant, s = hG, and x = nh.When Simpson's rule is used, the corresponding difference equation is (1 - s/3) yn+1 - (4s/3) yn - (1 + s/3) yn-1 = 0, (1) with the general solution yn = Ar1n + Br2n, (2) in which r1 and r2 are the roots of the quadratic equation (1 - s/3) r2 - (4s/3) r - (1 + s/3) = 0. From this equation we obtain the derivative dr/ds = (r2 + 4r + 1)2(12r2 + 12r + 12)-1, which is never negative. Hence the roots r1 = [2s/3 + (1 + s2/3)1/2] (1 - s/3)-1, r2 = [2s/3 - (1 + s2/3)1/2] (1 - s/3)-1, are monotone increasing functions for all real values of s, except for a discontinuity in r1 at s = 3.Moreover, the roots r1 and r2 are analytic within a circle of radius 31/2 with center at the origin in the complex s plane. Through terms of degree five in s the power series for r1 and r2 are respectively r1 = 1 + s + s2/2! + s3/3! + s4/4! + s5/72 + ··· = es + s5/180 + ··· , (3) and r2 = -1 + s/3 - s2/18 - s3/54 + 5s4/648 + 5s5/1944 + ··· (4) Obviously r1 is the desired root and r2 is the unwanted root that produces the oscillation.Quite apart from questions of stability the process of numerical integration with Simpson's rule requires that the quantity s/3 must be numerically less than unity and in practical computation should be considerably less than unity. Cf. Milne [5], p. 67. We shall therefore assume that |s| \lt 1. Table I shows to six decimal places the value of r1 and r2 for s ranging from -1 to +1 at steps of 0.1. It is evident that in this range r2 is numerically less than one if s is positive, greater than one if s is negative. Thus the oscillating term increases exponentially with n if G is negative, decreases if G is positive.Now suppose that after k steps of the process we recompute yk from the values already found, using Newton's “three eighths rule”, to obtain yk* = yk-3 + (3s/8)(yk + 3yk-1 + 3yk-2 + yk-3). Then we replace the originally computed yk by the arithmetic mean yk = (yk + yk*)/2. (5) From (2) and (5) we find that yk = Ark-31K(r1) + Brk-32K(r2), (6) in which K(r) is defined by the equation K(r) = [r3 + 1 + (3s/8) (r + 1)3]/2. (7) This function K(r) is the key to the problem.For by means of the series for r1 and r2 it can be shown that K(r1) = r13 + s5/96 + ··· (8) while K(r2) = s/2 - s2/4 + ··· . Hence equation (3) becomes yk = Ar1k + terms of degree 5 and higher + Brk-32(s/2) + terms of degree 2 and higher. (9) Comparing yk with yk we note that the desired solution is substantially unchanged, and agrees with the true solution eks through terms of degree 4 in s, while the unwanted solution has been decreased roughly by a factor of magnitude s/2.Table I shows to six decimal places the values of K(r1) and K(r2) in the range from s = -1 to s = +1. It is seen that in this interval the absolute value of K(r2) is always less than unity.Consider now the propagation of a single error starting at n = 0 and modified after every group of k steps by means of formula (5). Since En is a solution of equation (1), in the mth group of k steps the error En can be expressed by formula (2) in the form En = amr1j + bmr2j, (10) provided n = mk + j and j 1, and Simpson's rule, if uncorrected, produces instability.Table III shows the difference E = ens - yn between the true solution ens and the computed solution yn after n steps of the computation. Six values of k are used in table III, k = ∞ (that is, no stabilization), k = 169, 39, 19, 5, and 3. Four values of s are used, namely s = -0.10, -0.07, -0.04, and -0.01. The number n of steps in the computations varies from 300 for larger values of -s to 2000 for the smallest. Not all computations were carried to the full number of steps shown at the left, hence some columns are partially blank.The number of decimal places is indicated for each division of the table. For example at s = -0.10, k = 169, and n = 300, the entry -31 means -0.000031, while for s = -0.04, k = 169, n = 500, the entry -6 means -0.00000006.From table II we obtain the integral parts of q corresponding to the given values of s and find that according to theory the solution should be stable for s = -0.10 if k

New numerical approximation of fractional derivative with non-local and non-singular kernel: Application to chaotic models

Abstract: Recently a new concept of fractional differentiation with non-local and non-singular kernel was introduced in order to extend the limitations of the conventional Riemann-Liouville and Caputo fractional derivatives. A new numerical scheme has been developed, in this paper, for the newly established fractional differentiation. We present in general the error analysis. The new numerical scheme was applied to solve linear and non-linear fractional differential equations. We do not need a predictor-corrector to have an efficient algorithm, in this method. The comparison of approximate and exact solutions leaves no doubt believing that, the new numerical scheme is very efficient and converges toward exact solution very rapidly.

A Polyalgorithm for the Numerical Solution of Ordinary Differential Equations

Abstract: Two variable-order, varlable-step size methods for the numerical solution of the initial value problem for ordinary differential equations are presented. These methods share a common philosophy and have been combined in a single program. The two integrators are for stiff and nonstiff ordinary differential equations, respectively. The former integrator is based on backward differentiation formulas of orders one through five, each of which is stiffly stable, while the latter is based on formulas of Adams-Moulton type of orders one through twelve. Both use a Nordsieck history array; generating polynomials to compute coefficients; and Milne-like estimates of the principal part of the local truncation error to control error, step size, and order. Some numerical results are given.

Performance degradation of OFDM systems due to Doppler spreading

Abstract: The focus of this paper is on the performance of orthogonal frequency division multiplexing (OFDM) signals in mobile radio applications, such as 802.11a and digital video broadcasting (DVB) systems, e.g., DVB-CS2. The paper considers the evaluation of the error probability of an OFDM system transmitting over channels characterized by frequency selectivity and Rayleigh fading. The time variations of the channel during one OFDM symbol interval destroy the orthogonality of the different subcarriers and generate power leakage among the subcarriers, known as inter-carrier interference (ICI). For conventional modulation methods such as phase-shift keying (PSK) and quadrature-amplitude modulation (QAM), the bivariate probability density function (pdf) of the ICI is shown to be a weighted Gaussian mixture. The large computational complexity involved in using the weighted Gaussian mixture pdf to evaluate the error probability serves as the motivation for developing a two-dimensional Gram-Charlier representation for the bivariate pdf of the ICI. It is demonstrated that its truncated version of order 4 or 6 provides a very good approximation in the evaluation of the error probability for PSK and QAM in the presence of ICI. Based on Jakes' model for the Doppler effects, and an exponential multipath intensity profile, numerical results for the error probability are illustrated for several mobile speeds