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

1import numpy as np 

2 

3from pySDC.core.errors import ParameterError 

4from pySDC.core.problem import Problem 

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

6 

7 

8# noinspection PyUnusedLocal 

9class harmonic_oscillator(Problem): 

10 r""" 

11 Example implementing the harmonic oscillator [#]_ with mass :math:`1` 

12 

13 .. math:: 

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

15 

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. 

18 

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 

32 References 

33 ---------- 

34 .. [#] https://beltoforion.de/en/harmonic_oscillator/ 

35 """ 

36 

37 dtype_u = particles 

38 dtype_f = acceleration 

39 

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

41 """Initialization routine""" 

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

43 u0 = np.asarray(u0) 

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

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

46 

47 def eval_f(self, u, t): 

48 """ 

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

50 

51 Parameters 

52 ---------- 

53 u : dtype_u 

54 Current values of the particles. 

55 t : float 

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

57 

58 Returns 

59 ------- 

60 me : dtype_f 

61 The right-hand side of the problem. 

62 """ 

63 me = self.dtype_f(self.init) 

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

65 return me 

66 

67 def u_init(self): 

68 """ 

69 Helper function to compute the initial condition for u. 

70 """ 

71 u0 = self.u0 

72 

73 u = self.dtype_u(self.init) 

74 

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

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

77 

78 return u 

79 

80 def u_exact(self, t): 

81 r""" 

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

83 

84 Parameters 

85 ---------- 

86 t : float 

87 Time of the exact trajectory. 

88 

89 Returns 

90 ------- 

91 me : dtype_u 

92 Exact position and velocity. 

93 """ 

94 me = self.dtype_u(self.init) 

95 delta = self.mu / (2) 

96 omega = np.sqrt(self.k) 

97 

98 U_0 = self.u0 

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

100 if delta > omega: 

101 """ 

102 Overdamped case 

103 """ 

104 

105 lam_1 = -delta + alpha 

106 lam_2 = -delta - alpha 

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

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

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

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

111 

112 elif delta == omega: 

113 """ 

114 Critically damped case 

115 """ 

116 

117 A = U_0[0] 

118 B = U_0[1] + delta * A 

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

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

121 

122 elif delta < omega: 

123 """ 

124 Underdamped case 

125 """ 

126 

127 lam_1 = -delta + alpha * 1j 

128 lam_2 = -delta - alpha * 1j 

129 

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

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

132 

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

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

135 

136 else: 

137 pass 

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

139 return me 

140 

141 def eval_hamiltonian(self, u): 

142 """ 

143 Routine to compute the Hamiltonian. 

144 

145 Parameters 

146 ---------- 

147 u : dtype_u 

148 Current values of the particles. 

149 

150 Returns 

151 ------- 

152 ham : float 

153 The Hamiltonian. 

154 """ 

155 

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

157 return ham