Filename: Ellipsoi Ellipsoidal Error-Bounds for Trajectory Calculations November 14, 1999 11:12 pmProf. W. Kahan Page 1/5 Abbreviated Lecture Notes on Ellipsoidal Error Bounds for Trajectory Calculations Prof. W. KahanMathematics Dept., andElect. Eng. & Computer Science Dept.University of California at BerkeleyFor Presentation Oct. 12, 1993 Abstract: A practical way is outlined to bound error accrued during numerical calculations of trajectories. The algorithm accommodates uncertainties in the governing differential equation as well as error due to the numerical process. The accrued error is constrained to lie within an ellipsoid that can be proved to grow, as time τ increases to + ∞ , bigger than the worst possible accrual by a factor no worse than 1 + ß √τ for some constant ß , rather than exponentially bigger, until nonlinearity in the differential equation forces a singularity to manifest itself. Introduction: A trajectory is the graph of y ( τ ) , the solution of an Autonomous Initial Value Problem ( dy/d τ = ) y ' ( τ ) = f ( y ( τ )) for τ ≥ 0 , y ( 0 ) = y 0 . (AIVP)Here y 0 is a given vector and f a given vector-valued function of its vector argument. Perturbing y 0 to y 0 + δ y 0 and f to f + δ f perturbs the trajectory to the solution y + δ y of a perturbed AIVP(y + δ y) ' ( τ ) = (f + δ f) ( (y + δ y) ( τ ) ) for τ ≥ 0 , (y + δ y) ( 0 ) = y 0 + δ y 0 .Given only some kind of bound for δ y 0 and δ f , how can we infer a bound upon δ y ?Actually, in practice we know y + δ y but not y ; this will be overlooked now to simplify the exposition. To the same end, we shall not discuss how the truncation and rounding errors incurred by the numerical process that solves the AIVP can be incorporated into δ f along with errors due to idealizations that model a complex physical situation by a simplified expression f .Were all perturbations δ … infinitesimal, the accrued error δ y would satisfy the Adjoint or Variational Initial Value Problem associated with the given AIVP : δ y ' ( τ ) = J ( τ ) · δ y ( τ ) + δ f , δ y ( τ ) = δ y 0 . (VIVP)Here Jacobian matrix J ( τ ) := f ' ( y ) = ∂ f ( y ) / ∂ y at y = y ( τ ) ; this J ( τ ) is the matrix of first partial derivatives of f evaluated at y ( τ ) , the presumed-to-be-known trajectory. Perturbation δ f = δ f ( y ( τ )) ; however, whereas J ( τ ) can be computed from f ' and y ( τ ) , all we know about δ f is an upper bound. Likewise for δy0 . Presumably computable, perhaps functions of y(τ) , these bounds upon δf and δy0 are the data from which we wish to infer a bound upon δy(τ) . Computing this bound for accrued error along with y(τ) is the problem addressed in these notes.In practice the perturbations δ… are not infinitesimal. The effect in practice of their finiteness is to increase the bound upon δf by an amount roughly proportional to the product of the square of the computed bound upon δy(τ) and a bound upon the second derivative of f . This increase turns the linear VIVP into something nonlinear like a Riccati equation, whose solution mayFilename: Ellipsoi Ellipsoidal Error-Bounds for Trajectory Calculations November 14, 1999 11:12 pmProf. W. Kahan Page 2/5become infinite at a finite time τ even in some cases when δy(τ) is known to stay bounded for all finite τ . Despite its importance, this nonlinear contribution is not discussed in these notes. Instead, to simplify the exposition of a subject that is already complicated enough, δy is assumed to stay so small that its second-order contributions may safely be neglected.We assume J(τ) and δf to be continuous for all τ ; in practice δf may be only piecewise continuous. This is a technicality we could dispatch by converting differential equations into equivalent Volterra integral equations. For simplicity’s sake we won’t do that either.Despite all our simplifications, the computation of a bound upon δy remains so challenging that every previously published scheme I know about is prone to producing bounds too big by a factor that grows like an exponential function of τ when J(τ) behaves in a way that the scheme dislikes. The algorithm described here is far less pessimistic; its bound cannot grow too big by a factor bigger than 1 + ß√τ for some constant ß that depends upon moduli of continuity of J(τ) and of the bound for δf . The proof of this claim is too long to fit here. Space and time barely suffice for an outline of my algorithm.My algorithm adjoins to the given AIVP another differential equation whose solution, intended to be computed simultaneously with y(τ) , is a symmetric matrix A(τ) that describes an ellipsoid centered at y(τ) and surely big enough to enclose δy(τ) , but not too much bigger. The adjoined differential equation’s dimension is about half the square of y ’s dimension, so A(τ) may well cost enormously more to compute than y(τ) . Since previously published schemes typically cost twice as much as mine, I need not apologize for it.The Reachable Set:Let us now write z in place of δy and u in place of δf in the VIVP above; it becomesz'(τ) = J(τ)·z(τ) + u(τ) for τ ≥ 0 , z(0) = z0 . (VIVP)The matrix J(τ) is assumed computable; but for u(τ) and z0 only bounds are available. We construe these bounds as constraints that restrict u(τ) and z0 to certain small regions; sayu(τ) ∈ ÛÛ(τ) and z0 ∈ ÅÅ for given centrally symmetric convex bodies ÛÛ(τ) and ÅÅ characterized by parameters to be discussed later. For instance these bodies could be spheres characterized by their radii, which are then upper bounds for the lengths of u(τ) and z0 . Parallelepipeds have been used too, characterized by the matrices that map a unit hypercube onto them. But we use ellipsoids for reasons to be discussed