Coverage for pySDC/implementations/problem_classes/HenonHeiles.py: 100%

1import numpy as np 

2 

3from pySDC.core.problem import Problem 

4from pySDC.implementations.datatype_classes.particles import particles, acceleration 

5 

6 

7# noinspection PyUnusedLocal 

8class henon_heiles(Problem): 

9 r""" 

10 This class implements the second-order Hénon-Heiles system 

11 

12 .. math:: 

13 \frac{d^2 x}{dt^2} = - x - 2 x y, 

14 

15 .. math:: 

16 \frac{d^2 y}{dt} = - y - x^2 + y^2 

17 

18 with Hamiltonian 

19 

20 .. math:: 

21 H = 0.5 \left[\left(\frac{d x}{d t}\right)^2 + \left(\frac{d y}{d t}\right)^2\right] + 0.5 \left(x^2 + y^2\right) 

22 + x^2 y - \frac{y^3}{3}. 

23 """ 

24 

25 dtype_u = particles 

26 dtype_f = acceleration 

27 

28 def __init__(self): 

29 """Initialization routine""" 

30 # invoke super init, passing nparts 

31 super().__init__((2, None, np.dtype('float64'))) 

32 

33 def eval_f(self, u, t): 

34 """ 

35 Routine to compute the right-hand side of the problem. 

36 

37 Parameters 

38 ---------- 

39 u : dtype_u 

40 Current values of the particles. 

41 t : float 

42 Current time of the particles is computed (not used here). 

43 

44 Returns 

45 ------- 

46 me : dtype_f 

47 The right-hand side of the problem. 

48 """ 

49 me = self.dtype_f(self.init) 

50 me[0] = -u.pos[0] - 2 * u.pos[0] * u.pos[1] 

51 me[1] = -u.pos[1] - u.pos[0] ** 2 + u.pos[1] ** 2 

52 return me 

53 

54 def u_exact(self, t): 

55 r""" 

56 Routine to compute the exact/initial trajectory at time :math:`t`. 

57 

58 Parameters 

59 ---------- 

60 t : float 

61 Time of the exact/initial trajectory. 

62 

63 Returns 

64 ------- 

65 me : dtype_u 

66 Exact/initial position and velocity. 

67 """ 

68 assert t == 0.0, 'error, u_exact only works for the initial time t0=0' 

69 me = self.dtype_u(self.init) 

70 

71 q1 = 0.0 

72 q2 = 0.2 

73 p2 = 0.2 

74 U0 = 0.5 * (q1 * q1 + q2 * q2) + q1 * q1 * q2 - q2 * q2 * q2 / 3.0 

75 H0 = 0.125 

76 

77 me.pos[0] = q1 

78 me.pos[1] = q2 

79 me.vel[0] = np.sqrt(2.0 * (H0 - U0) - p2 * p2) 

80 me.vel[1] = p2 

81 return me 

82 

83 def eval_hamiltonian(self, u): 

84 """ 

85 Routine to compute the Hamiltonian. 

86 

87 Parameters 

88 ---------- 

89 u : dtype_u 

90 Current values of The particles. 

91 

92 Returns 

93 ------- 

94 ham : float 

95 The Hamiltonian. 

96 """ 

97 

98 ham = 0.5 * (u.vel[0] ** 2 + u.vel[1] ** 2) 

99 ham += 0.5 * (u.pos[0] ** 2 + u.pos[1] ** 2) 

100 ham += u.pos[0] ** 2 * u.pos[1] - u.pos[1] ** 3 / 3.0 

101 return ham