# Coverage for pySDC/projects/DAE/problems/pendulum2D.py: 88%

## 25 statements

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

1from pySDC.projects.DAE.misc.problemDAE import ProblemDAE

4class Pendulum2D(ProblemDAE):

5 r"""

6 Example implementing the well known 2D pendulum as a first order differential-algebraic equation (DAE) of index 3.

7 The DAE system is given by the equations

9 .. math::

10 \frac{dp}{dt} = u,

12 .. math::

13 \frac{dq}{dt} = v,

15 .. math::

16 m\frac{du}{dt} = -p \lambda,

18 .. math::

19 m\frac{dv}{dt} = -q \lambda - g,

21 .. math::

22 0 = p^2 + q^2 - l^2

24 for :math:l=1 and :math:m=1. The pendulum is used in most introductory literature on DAEs, for example on page 8

25 of [1]_.

27 Parameters

28 ----------

29 nvars : int

30 Number of unknowns of the system of DAEs.

31 newton_tol : float

32 Tolerance for Newton solver.

34 Attributes

35 ----------

36 t_end: float

37 The end time at which the reference solution is determined.

39 References

40 ----------

41 .. [1] E. Hairer, C. Lubich, M. Roche. The numerical solution of differential-algebraic systems by Runge-Kutta methods.

42 Lect. Notes Math. (1989).

43 """

45 def __init__(self, newton_tol):

46 """Initialization routine"""

47 super().__init__(nvars=5, newton_tol=newton_tol)

49 # data file must be generated and stored under misc/data and self.t_end = t[-1]

51 # t = data[:, 0]

52 # solution = data[:, 1:]

53 # self.u_ref = interp1d(t, solution, kind='cubic', axis=0, fill_value='extrapolate')

54 self.t_end = 0.0

56 def eval_f(self, u, du, t):

57 r"""

58 Routine to evaluate the implicit representation of the problem, i.e., :math:F(u, u', t).

60 Parameters

61 ----------

62 u : dtype_u

63 Current values of the numerical solution at time t.

64 du : dtype_u

65 Current values of the derivative of the numerical solution at time t.

66 t : float

67 Current time of the numerical solution.

69 Returns

70 -------

71 f : dtype_f

72 Current value of the right-hand side of f (which includes five components).

73 """

74 g = 9.8

75 # The last element of u is a Lagrange multiplier. Not sure if this needs to be time dependent, but must model the

76 # weight somehow

77 f = self.dtype_f(self.init)

78 f.diff[:4] = (

79 du.diff[0] - u.diff[2],

80 du.diff[1] - u.diff[3],

81 du.diff[2] + u.alg[0] * u.diff[0],

82 du.diff[3] + u.alg[0] * u.diff[1] + g,

83 )

84 f.alg[0] = u.diff[0] ** 2 + u.diff[1] ** 2 - 1

85 self.work_counters['rhs']()

86 return f

88 def u_exact(self, t):

89 """

90 Approximation of the exact solution generated by spline interpolation of an extremely accurate numerical reference solution.

92 Parameters

93 ----------

94 t : float

95 The time of the reference solution.

97 Returns

98 -------

99 me : dtype_u

100 The reference solution as mesh object. It contains fixed initial conditions at initial time.

101 """

102 me = self.dtype_u(self.init)

103 if t == 0:

104 me.diff[:4] = (-1, 0, 0, 0)

105 me.alg[0] = 0

106 elif t < self.t_end:

107 u_ref = self.u_ref(t)

108 me.diff[:4] = u_ref[:4]

109 me.alg[0] = u_ref[5]

110 else:

111 self.logger.warning("Requested time exceeds domain of the reference solution. Returning zero.")

112 me.diff[:4] = (0, 0, 0, 0)

113 me.alg[0] = 0

114 return me