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

## 31 statements

, created at 2024-09-09 14:59 +0000

1import numpy as np

3from pySDC.core.problem import Problem

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

7# noinspection PyUnusedLocal

8class henon_heiles(Problem):

9 r"""

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

12 .. math::

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

15 .. math::

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

18 with Hamiltonian

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 """

25 dtype_u = particles

26 dtype_f = acceleration

28 def __init__(self):

29 """Initialization routine"""

30 # invoke super init, passing nparts

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

33 def eval_f(self, u, t):

34 """

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

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).

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

54 def u_exact(self, t):

55 r"""

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

58 Parameters

59 ----------

60 t : float

61 Time of the exact/initial trajectory.

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)

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

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

83 def eval_hamiltonian(self, u):

84 """

85 Routine to compute the Hamiltonian.

87 Parameters

88 ----------

89 u : dtype_u

90 Current values of The particles.

92 Returns

93 -------

94 ham : float

95 The Hamiltonian.

96 """

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