Math. 128: Taylor Series for Differential Equation Solvers’ Local Error Prof. W. Kahan page 1/2 August 25, 2004 5:53 am The Taylor Series of an infinitely differentiable vector-valued function y(t) of a scalar t isy(t) = y(0) + t·y ' (0) + t 2 ·y " (0)/2 + t 3 ·y "' (0)/6 + t 4 ·y "" (0)/24 + … .These derivatives can be computed for a solution of the Initial Value Problem “ y(0) = y 0 is given, and y ' (t) = ƒ(y(t)) for all t ≥ 0 ”from the derivatives of the given vector-valued function ƒ(y) . In fact, from the Chain Rule , y ' = ƒ , y " = ƒ ' ·ƒ , y '" = ƒ " ·ƒ·ƒ + ƒ ' 2 ·ƒ , y "" = ƒ '" ·ƒ·ƒ·ƒ + 3ƒ " ·ƒ ' ·ƒ·ƒ + ƒ ' ·ƒ " ·ƒ·ƒ + ƒ ' 3 ·ƒ , … .Then each derivative of y(t) can be evaluated at t = 0 by evaluating each derivative of ƒ(y) at y = y(0) = y 0 . Note that the higher derivatives of ƒ are symmetric multilinear operators ; for instance, ƒ " (y) is a symmetric bilinear operator: ƒ " ·u·v = ƒ " ·v·u is a vector-valued linear function of each vector u and v separately. Because linear operators do not necessarily commute, ƒ " ·ƒ ' ·ƒ·ƒ ≠ ƒ ' ·ƒ " ·ƒ·ƒ in general, though they are equal if y’s vector space is one-dimensional. If y’s vector space is N-dimensional, then y and ƒ can be represented by column vectors each with N components; ƒ ' by a matrix with N 2 components; ƒ " by an array with N 3 components of which at most (N+1)N 2 / 2 can be distinct; ƒ "' by an array of N 4 components … . Higher derivatives’ arrays become huge when N is large.Normally the Taylor series would be used to obtain y(t+h) from y(t) for any sufficiently small stepsize h :y(t+h) = y(t) + h·y ' (t) + h 2 ·y " (t)/2 + h 3 ·y "' (t)/6 + h 4 ·y "" (t)/24 + … ,in which the derivatives dy(t+h)/dh etc. are computed at h = 0 from the same formulas y ' = ƒ , y " = ƒ ' ·ƒ , y '" = ƒ " ·ƒ·ƒ + ƒ ' 2 ·ƒ , y "" = ƒ '" ·ƒ·ƒ·ƒ + 3ƒ " ·ƒ ' ·ƒ·ƒ + ƒ ' ·ƒ " ·ƒ·ƒ + ƒ ' 3 ·ƒ , … as before except that now ƒ(y) and its derivatives are computed at y = y(t) .A similar process generates a formal series for any one-step numerical method’s formula that advances an approximate solution y = y(t) through one step h to Y = Y(t+h) , but now we differentiate with respect to h instead of t to getY(t+h) = y + h·Y ' + h 2 ·Y " /2 + h 3 ·Y "' /6 + h 4 ·Y "" /24 + …in which y = y(t) and the derivatives Y ' etc. are derivatives of Y(t+h) with respect to h evaluated at h = 0 . These derivatives depend upon the numerical method’s formula. For example, take the (implicit) Trapezoidal Rule Y = y + h·(ƒ(y) + ƒ(Y))/2 . Now, at t+h ,Y ' = (ƒ(y) + ƒ(Y))/2 + h·ƒ ' (Y)·Y ' /2 ,Y " = ƒ ' (Y)·Y ' + h·( ƒ " (Y)·Y ' ·Y ' + ƒ ' (Y)·Y " )/2 ,Y "' = 3ƒ " (Y)·Y ' ·Y ' /2 + ƒ ' (Y)·Y "/2 + h·( ƒ"'(Y)·Y'·Y'·Y' + … )/2 , etc.Here every instance of Y or its derivatives is evaluated at t+h . For the Taylor series we set h = 0 in the foregoing formulas and substitute for derivatives of Y in succession to getY' = ƒ , Y" = ƒ'·ƒ , Y"' = ( 3ƒ"·ƒ + ƒ'2 )·ƒ/2 , …in which now the derivatives of Y are evaluated at t , and ƒ and its derivatives are evaluated at y(t) . Hence, the computed approximationY(t+h) = y(t) + h·ƒ + h2·ƒ'·ƒ/2 + h3·( 3ƒ"·ƒ + ƒ'2 )·ƒ/12 + …can be compared with the local solutiony(t+h) = y(t) + h·ƒ + h2·ƒ'·ƒ/2 + h3·( ƒ"·ƒ + ƒ'2 )·ƒ/6 + … to reveal the 2nd-order Trapezoidal Rule’s local truncation errory(t+h) – Y(t+h) = h3·( ƒ'2 – ƒ"·ƒ)·ƒ/12 + … .Math. 128: Taylor Series for Differential Equation Solvers’ Local ErrorProf. W. Kahan page 2/2 August 25, 2004 5:53 amAnother example is the (implicit) Midpoint Rule Y = y + h·ƒ((y+Y)/2) . At t+h , Y' = ƒ((y+Y)/2) + h·ƒ'((y+Y)/2)·Y'/2 ,Y" = ƒ'((y+Y)/2)·Y' + h·( ƒ"((y+Y)/2)·Y'·Y'/4 + ƒ'((y+Y)/2)·Y"/2 ) ,Y"' = 3ƒ"((y+Y)/2)·Y'·Y'/4 + 3ƒ'((y+Y)/2)·Y"/2 + h·( ƒ"' … ) , etc.Setting h = 0 and substituting for derivatives of Y in succession yieldsY' = ƒ , Y" = ƒ'·ƒ , Y"' = ( 3ƒ"·ƒ + 6ƒ'2 )·ƒ/4 , …in which now the derivatives of Y are evaluated at t , and ƒ and its derivatives are evaluated at y(t) . Hence, by comparison with the Taylor series for y(t+h) we find the 2nd-order Midpoint Rule’s local truncation error to bey(t+h) – Y(t+h) = h3·( ƒ"·ƒ – 2ƒ'2 )·ƒ/24 + … .The foregoing manipulations are tedious enough that only a computerized algebra system should perform them. However, programming Maple or Mathematica or Derive to perform them has been more difficult than it should be. At least some of the difficulty arises, I think, because these languages disallow declaration of a variable’s linguistic Type. Besides the derivatives’ non-commutative partially associative multiplication, their multi-linear symmetry has to be taken into account in order to achieve correct simplifications of expressions like the fifth derivativeyv = ƒ""·ƒ·ƒ·ƒ·ƒ + 6ƒ"'·ƒ'·ƒ·ƒ·ƒ + 4ƒ"·ƒ"·ƒ·ƒ·ƒ + 4ƒ"·ƒ'2·ƒ·ƒ ++ 3ƒ"·ƒ'·ƒ·ƒ'·ƒ + ƒ'·ƒ"'·ƒ·ƒ·ƒ + 3ƒ'·ƒ"·ƒ'·ƒ·ƒ + ƒ'2·ƒ"·ƒ·ƒ + ƒ'4·ƒ .( I hope I’ve gotten it right.). . . . . .The foregoing formal (because their convergence is undetermined) Taylor series expansions do not reveal an important property possessed by the computed solution Y of the Initial Value Problem when it is obtained from a Ref lexive formula, which is a formula in which Y(t+h) and y(t) are merely swapped when the sign of h is reversed. The Midpoint and Trapezoidal Rules’ formulas are reflexive. The composition of T/h steps of a reflexive formula to approximate the true solution y(T) at a fixed T , but using any sufficiently small stepsize h so long as T/h is an integer, can be proved to produce a computed approximation Y(T) that depends upon h and differs from y(T) by an errory(T) – Y(T) = c2h2 + c4h4 + c6h6 + …whose formal expansion in powers of h contains only even