# Coverage for pySDC/implementations/problem_classes/HarmonicOscillator.py: 69%

## 54 statements

, created at 2024-09-19 09:13 +0000

1import numpy as np

3from pySDC.core.errors import ParameterError

4from pySDC.core.problem import Problem

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

8# noinspection PyUnusedLocal

9class harmonic_oscillator(Problem):

10 r"""

11 Example implementing the harmonic oscillator with mass :math:1

13 .. math::

14 \frac{d^2 x}{dt^2} = -kx - \mu \frac{d x}{dt},

16 which is a second-order problem. The unknown function :math:x denotes the position of the mass, and the

17 derivative is the velocity. :math:\mu defines the damping and :math:k is the spring constant.

19 Parameters

20 ----------

21 k : float, optional

22 Spring constant :math:k.

23 mu : float, optional

24 Damping parameter :math:\mu.

25 u0 : tuple, optional

26 Initial condition for the position, and the velocity. Should be a tuple, e.g. u0=(1, 0).

27 phase : float, optional

28 Phase of the oscillation.

29 amp : float, optional

30 Amplitude of the oscillation.

31 Source: https://beltoforion.de/en/harmonic_oscillator/

32 """

34 dtype_u = particles

35 dtype_f = acceleration

37 def __init__(self, k=1.0, mu=0.0, u0=(1, 0), phase=0.0, amp=1.0):

38 """Initialization routine"""

39 # invoke super init, passing nparts, dtype_u and dtype_f

40 u0 = np.asarray(u0)

41 super().__init__((1, None, np.dtype("float64")))

42 self._makeAttributeAndRegister('k', 'mu', 'u0', 'phase', 'amp', localVars=locals(), readOnly=True)

44 def eval_f(self, u, t):

45 """

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

48 Parameters

49 ----------

50 u : dtype_u

51 Current values of the particles.

52 t : float

53 Current time of the numerical solution is computed (not used here).

55 Returns

56 -------

57 me : dtype_f

58 The right-hand side of the problem.

59 """

60 me = self.dtype_f(self.init)

61 me[:] = -self.k * u.pos - self.mu * u.vel

62 return me

64 def u_init(self):

65 """

66 Helper function to compute the initial condition for u.

67 """

68 u0 = self.u0

70 u = self.dtype_u(self.init)

72 u.pos[0] = u0[0]

73 u.vel[0] = u0[1]

75 return u

77 def u_exact(self, t):

78 r"""

79 Routine to compute the exact trajectory at time :math:t.

81 Parameters

82 ----------

83 t : float

84 Time of the exact trajectory.

86 Returns

87 -------

88 me : dtype_u

89 Exact position and velocity.

90 """

91 me = self.dtype_u(self.init)

92 delta = self.mu / (2)

93 omega = np.sqrt(self.k)

95 U_0 = self.u0

96 alpha = np.sqrt(np.abs(delta**2 - omega**2))

97 print(self.mu)

98 if delta > omega:

99 """

100 Overdamped case

101 """

103 lam_1 = -delta + alpha

104 lam_2 = -delta - alpha

105 L = np.array([[1, 1], [lam_1, lam_2]])

106 A, B = np.linalg.solve(L, U_0)

107 me.pos[:] = A * np.exp(lam_1 * t) + B * np.exp(lam_2 * t)

108 me.vel[:] = A * lam_1 * np.exp(lam_1 * t) + B * lam_2 * np.exp(lam_2 * t)

110 elif delta == omega:

111 """

112 Critically damped case

113 """

115 A = U_0[0]

116 B = U_0[1] + delta * A

117 me.pos[:] = np.exp(-delta * t) * (A + t * B)

118 me.vel[:] = -delta * me.pos[:] + np.exp(-delta * t) * B

120 elif delta < omega:

121 """

122 Underdamped case

123 """

125 lam_1 = -delta + alpha * 1j

126 lam_2 = -delta - alpha * 1j

128 M = np.array([[1, 1], [lam_1, lam_2]], dtype=complex)

129 A, B = np.linalg.solve(M, U_0)

131 me.pos[:] = np.real(A * np.exp(lam_1 * t) + B * np.exp(lam_2 * t))

132 me.vel[:] = np.real(A * lam_1 * np.exp(lam_1 * t) + B * lam_2 * np.exp(lam_2 * t))

134 else:

135 pass

136 raise ParameterError("Exact solution is not working")

137 return me

139 def eval_hamiltonian(self, u):

140 """

141 Routine to compute the Hamiltonian.

143 Parameters

144 ----------

145 u : dtype_u

146 Current values of the particles.

148 Returns

149 -------

150 ham : float

151 The Hamiltonian.

152 """

154 ham = 0.5 * self.k * u.pos[0] ** 2 + 0.5 * u.vel[0] ** 2

155 return ham