# Coverage for pySDC/projects/DAE/problems/simpleDAE.py: 100%

## 22 statements

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

1import numpy as np

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

6class SimpleDAE(ProblemDAE):

7 r"""

8 Example implementing a smooth linear index-2 differential-algebraic equation (DAE) with known analytical solution.

9 The DAE system is given by

11 .. math::

12 \frac{d u_1 (t)}{dt} = (\alpha - \frac{1}{2 - t}) u_1 (t) + (2-t) \alpha z (t) + \frac{3 - t}{2 - t},

14 .. math::

15 \frac{d u_2 (t)}{dt} = \frac{1 - \alpha}{t - 2} u_1 (t) - u_2 (t) + (\alpha - 1) z (t) + 2 e^{t},

17 .. math::

18 0 = (t + 2) u_1 (t) + (t^{2} - 4) u_2 (t) - (t^{2} + t - 2) e^{t}.

20 The exact solution of this system is

22 .. math::

23 u_1 (t) = u_2 (t) = e^{t},

25 .. math::

26 z (t) = -\frac{e^{t}}{2 - t}.

28 This example is commonly used to test that numerical implementations are functioning correctly. See, for example,

29 page 267 of [1]_.

31 Parameters

32 ----------

33 nvars : int

34 Number of unknowns of the system of DAEs.

35 newton_tol : float

36 Tolerance for Newton solver.

38 References

39 ----------

40 .. [1] U. Ascher, L. R. Petzold. Computer method for ordinary differential equations and differential-algebraic

41 equations. Society for Industrial and Applied Mathematics (1998).

42 """

44 def __init__(self, newton_tol=1e-10):

45 """Initialization routine"""

46 super().__init__(nvars=3, newton_tol=newton_tol)

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

49 r"""

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

52 Parameters

53 ----------

54 u : dtype_u

55 Current values of the numerical solution at time t.

56 du : dtype_u

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

58 t : float

59 Current time of the numerical solution.

61 Returns

62 -------

63 f : dtype_f

64 Current value of the right-hand side of f (which includes three components).

65 """

66 # Smooth index-2 DAE pg. 267 Ascher and Petzold (also the first example in KDC Minion paper)

67 a = 10.0

68 f = self.dtype_f(self.init)

70 f.diff[:2] = (

71 -du.diff[0] + (a - 1 / (2 - t)) * u.diff[0] + (2 - t) * a * u.alg[0] + (3 - t) / (2 - t) * np.exp(t),

72 -du.diff[1] + (1 - a) / (t - 2) * u.diff[0] - u.diff[1] + (a - 1) * u.alg[0] + 2 * np.exp(t),

73 )

74 f.alg[0] = (t + 2) * u.diff[0] + (t**2 - 4) * u.diff[1] - (t**2 + t - 2) * np.exp(t)

75 self.work_counters['rhs']()

76 return f

78 def u_exact(self, t):

79 """

80 Routine for the exact solution.

82 Parameters

83 ----------

84 t : float

85 The time of the reference solution.

87 Returns

88 -------

89 me : dtype_u

90 The reference solution as mesh object containing three components.

91 """

92 me = self.dtype_u(self.init)

93 me.diff[:2] = (np.exp(t), np.exp(t))

94 me.alg[0] = -np.exp(t) / (2 - t)

95 return me

97 def du_exact(self, t):

98 """

99 Routine for the derivative of the exact solution.

101 Parameters

102 ----------

103 t : float

104 The time of the reference solution.

106 Returns

107 -------

108 me : dtype_u

109 The reference solution as mesh object containing three components.

110 """

112 me = self.dtype_u(self.init)

113 me.diff[:2] = (np.exp(t), np.exp(t))

114 me.alg[0] = (np.exp(t) * (t - 3)) / ((2 - t) ** 2)

115 return me