2 where the Hamiltonian function H has the special structure Hp, q) = Tp) + V q). In this case, the Hamiltonian system )-3) is of the form p q = fq) = gp) with f = q V, g = p T. For the numerical integration of ) we consider partitioned Runge-Kutta PRK) methods. The idea of the PRK methods is to take two different Runge-Kutta RK) methods, A, b, Â,ˆb), and to treat the variables y with the first method, A, b), and the z variables with the second one, Â,ˆb). In this way, the numerical solution from t n, y n, z n ) to t n+, y n+, z n+ ) with the s-stage PRK method A, b, Â,ˆb) is given by y n+ = y n + h b t I l )fy n+, Z n+ ), z n+ = z n + h ˆb t I m )gy n+, Z n+ ), where the internal stages vectors Y n+, Z n+ are obtained from the system Y n+ = e y n + h A I l )fy n+, Z n+ ), 4) Z n+ = e z n + h Â I m)gy n+, Z n+ ). ) As usual, we have denoted by the Kronecker product and by e =,..., ) t. For implicit PRK methods, the computational effort may be dominated by the cost of solving the non-linear systems 4)-). These non-linear systems are solved with an iterative method that requires starting values Y ) n+, Z ) n+, and it is well known that these values must be as accurate as possible, because in other case, the number of iterations in each step may be too high or even worse, the convergence may fail. A common way to proceed is to take as starting values the last computed numerical solution, i.e., Y ) n+ = e y n, Z ) n+ = e z n. We will refer to them as the trivial predictors. However it is also possible to involve the last computed internal stages, Y n, Z n, and the approximation y n, and consider starting algorithms of the form Y ) n+ = b y n + B I l )Y n, ) Z ) n+ = c z n + C I m )Z n, 7) where b, c IR s, B and C are s s matrices which have to be determined. Several kinds of starting algorithms for RK methods, also called stage value predictors, have been studied by different authors [], [], [4], [7], [4], [3], [], [], [], []). The

3 starting algorithms considered for example in [], [] or [7] are of the form )-7), whereas the ones considered for example in [], [4] or [] also involve function evaluations of the internal stages. Although in this paper we only deal with Hamiltonian systems without restrictions ODEs), we are also interested in Hamiltonian systems with restrictions DAEs). For ODE systems we can use any of the stage value predictors mentioned in the previous paragraph, but this is not the case for DAEs. Predictors involving function evaluations cannot be used for this type of systems [], [8]). This is the reason why we use predictors of the form )-7). These predictors have already been studied for PRK methods in [9]. In this paper we consider the PRK methods Lobatto IIIA-IIIB, that have been proved to be suitable ones for Hamiltonian systems [], [8]) and construct good starting algorithms for them. This is done in Section. The numerical experiments done in Section 3, show the efficiency of these starting algorithms. Starting Algorithms for Lobatto IIIA-IIIB Methods We consider the PRK methods A, b, Â,ˆb), where A, b) is the Lobatto IIIA method and Â,ˆb) is the Lobatto IIB method. We will refer to this PRK method as the Lobatto IIIA-IIB method. Observe that these methods share the same nodes, that is to say c = ĉ. Furthermore, the s-stage Lobatto IIIA method satisfies Cs), and thus Ac q = c q /q, q =,..., s; whereas the s-stage Lobatto IIIB method satisfies Cs ), and therefore Âĉ q = ĉ q /q, q =,..., s Figure : The pair Lobatto IIIA-IIIB for s = Figure : The pair Lobatto IIIA-IIIB for s = 4

4 In [9] the general order conditions for predictors of the form )-7) have been obtained. Fortunately, these conditions are reduced considerably if the method PRK satisfies some simplifying assumptions, as it occurs with Lobatto IIIA-IIIB methods. For s 3, the s-stage Lobatto IIIA and IIIB methods satisfy at least C3) and C) respectively. In this case, the general order conditions in [9] for the variable y are reduced to the following ones Consistency: b + B e = e, Order : B c = e + r c, Order : B A c = eb t c + rae + rc), 8) Order 3: B Ac = eb t c + rae + rc), B AÂc = ebt Âc + raeb t c + r AÂe + rc), whereas the order conditions for the variable z are reduced to Consistency: c + C e = e, Order : C c = e + r c, Order : C Â c = eˆb t c + râe + rc), 9) Order 3: C Âc = eˆb t c + râe + rc), C Â c = eˆb t Âc + râeˆb t c + r Â e + rc). The parameter r = h n+ /h n, where h n+ denotes the current stepsize, h n+ = t n+ t n, and h n denotes the last stepsize, h n = t n t n. For s = 3, if we take into account the number of parameters available and the number of order conditions, we get that the maximum order achieved is for both variables, y and z. It is still possible to impose one of the two conditions of order 3, but not both of them. For s 4, the s-stage Lobatto IIIA and IIIB methods satisfy at least C4) and C) respectively. Both of them satisfy b t c = ˆb t c = /. In this case, in 8) and 9) the two conditions of order 3 are equivalent. Hence for s = 4, it is possible to achieve order 3 for both variables, y and z. In order to simplify the implementation, it is interesting to initialize the variables y and z with the same coefficients. This situation corresponds to the case b = c and B = C in )-7).

7 The problems have been integrated for different values of T OL and for different stepsizes. In Tables - we show the average number of iterations per step needed when the trivial predictor is used, and the ones when the optimum predictor given by ) is used. For each value of TOL and h, the values on the left correspond to the trivial predictor, whereas the values on the right correspond to the optimum predictor given by ). 3. Problem The problem is given by y = 4z + t) + t y) = 3) z = y t z+t) z) = The exact solution is yt) = sint + t, zt) = cost t. We have integrated this problem for t [, ]. h TOL.E-3.E-.E-7.E-./../.9 3./..E-3./../../..E-3./../../..E-3.898/../../.9 Table : Iterations per step. Problem It can be seen that except in one case, the number of iterations is lower with the predictor constructed in this paper. The exception is for h =.E-3 and TOL =.E-7 where the number of iterations coincides. We remark that in this case, the numerical solution obtained with the optimum predictor is more accurate. 3. The restricted three-body problem We have considered the restricted three-body problem from [3] x = v x x) = x y = v y y) = y z = v z z) = z v x = v y + x [ ] x+µ µ x µ + µ r 3 v r 3 x ) = v x v y = v x + y [ ] µ + µ r 3 r y 3 vy ) = v y v z = [ ] µ + µ r 3 r z 3 vz ) = v z 4)

2.6 Iterative Solutions of Linear Systems 143 2.6 Iterative Solutions of Linear Systems Consistent linear systems in real life are solved in one of two ways: by direct calculation (using a matrix factorization,

Lecture 9 Numerical methods for American options Lecture Notes by Andrzej Palczewski Computational Finance p. 1 American options The holder of an American option has the right to exercise it at any moment

Joint models for classification and comparison of mortality in different countries. Viani D. Biatat 1 and Iain D. Currie 1 1 Department of Actuarial Mathematics and Statistics, and the Maxwell Institute

Stability of the LMS Adaptive Filter by Means of a State Equation Vítor H. Nascimento and Ali H. Sayed Electrical Engineering Department University of California Los Angeles, CA 90095 Abstract This work

Roots of Polynomials (Com S 477/577 Notes) Yan-Bin Jia Sep 24, 2015 A direct corollary of the fundamental theorem of algebra is that p(x) can be factorized over the complex domain into a product a n (x

578 CHAPTER 1 NUMERICAL METHODS 1. ITERATIVE METHODS FOR SOLVING LINEAR SYSTEMS As a numerical technique, Gaussian elimination is rather unusual because it is direct. That is, a solution is obtained after

ON GALOIS REALIZATIONS OF THE 2-COVERABLE SYMMETRIC AND ALTERNATING GROUPS DANIEL RABAYEV AND JACK SONN Abstract. Let f(x) be a monic polynomial in Z[x] with no rational roots but with roots in Q p for

September 29, 201 9-1 9. Particular Solutions of Non-homogeneous second order equations Undetermined Coefficients We have seen that in order to find the general solution to the second order differential

Lecture Notes on Polynomials Arne Jensen Department of Mathematical Sciences Aalborg University c 008 Introduction These lecture notes give a very short introduction to polynomials with real and complex

1 Heriot-Watt University M.Sc. in Actuarial Science Life Insurance Mathematics I Tutorial 5 1. Consider the illness-death model in Figure 1. A life age takes out a policy with a term of n years that pays

NOTES ON LINEAR TRANSFORMATIONS Definition 1. Let V and W be vector spaces. A function T : V W is a linear transformation from V to W if the following two properties hold. i T v + v = T v + T v for all

Corrections to the First Printing Chapter 2 (i) Page 48, Paragraph 1: cells/µ l should be cells/µl without the space. (ii) Page 48, Paragraph 2: Uninfected cells T i should not have the asterisk. Chapter

Overview of Violations of the Basic Assumptions in the Classical Normal Linear Regression Model 1 September 004 A. Introduction and assumptions The classical normal linear regression model can be written

24. The Branch and Bound Method It has serious practical consequences if it is known that a combinatorial problem is NP-complete. Then one can conclude according to the present state of science that no

A Learning Based Method for Super-Resolution of Low Resolution Images Emre Ugur June 1, 2004 emre.ugur@ceng.metu.edu.tr Abstract The main objective of this project is the study of a learning based method

Continued Fractions and the Euclidean Algorithm Lecture notes prepared for MATH 326, Spring 997 Department of Mathematics and Statistics University at Albany William F Hammond Table of Contents Introduction

International Conference on Mathematical and Statistical Modeling in Honor of Enrique Castillo. June 28-30, 2006 Generalized Inverse Computation Based on an Orthogonal Decomposition Methodology. Patricia

THE SCHEDULING OF MAINTENANCE SERVICE Shoshana Anily Celia A. Glass Refael Hassin Abstract We study a discrete problem of scheduling activities of several types under the constraint that at most a single

Chapter 3 Chebyshev Expansions The best is the cheapest. Benjamin Franklin 3.1 Introduction In Chapter, approximations were considered consisting of expansions around a specific value of the variable (finite

Exam Introduction Mathematical Finance and Insurance Date: January 8, 2013. Duration: 3 hours. This is a closed-book exam. The exam does not use scrap cards. Simple calculators are allowed. The questions

ONDERZOEKSRAPPORT NR 8904 OPTIMAl PREMIUM CONTROl IN A NON-liFE INSURANCE BUSINESS BY M. VANDEBROEK & J. DHAENE D/1989/2376/5 1 IN A OPTIMAl PREMIUM CONTROl NON-liFE INSURANCE BUSINESS By Martina Vandebroek

Appendix H H.Calculating Normal Vectors This appendix describes how to calculate normal vectors for surfaces. You need to define normals to use the OpenGL lighting facility, which is described in Chapter