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

## 22 statements

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

1import numpy as np

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

4from pySDC.implementations.datatype_classes.mesh import mesh

7class ProblematicF(ProblemDAE):

8 r"""

9 Standard example of a very simple fully implicit index-2 differential algebraic equation (DAE) that is not

10 numerically solvable for certain choices of the parameter :math:\eta. The DAE system is given by

12 .. math::

13 \frac{d y(t)}{dt} + \eta t \frac{d z(t)}{dt} + (1 + \eta) z (t) = g (t).

15 .. math::

16 y (t) + \eta t z (t) = f(t),

18 See, for example, page 264 of [1]_.

20 Parameters

21 ----------

22 nvars : int

23 Number of unknowns of the system of DAEs.

24 newton_tol : float

25 Tolerance for Newton solver.

27 Attributes

28 ----------

29 eta : float

30 Specific parameter of the problem.

32 References

33 ----------

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

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

36 """

38 dtype_u = mesh

39 dtype_f = mesh

41 def __init__(self, newton_tol, eta=1):

42 """Initialization routine"""

43 super().__init__(nvars=2, newton_tol=newton_tol)

44 self._makeAttributeAndRegister('eta', localVars=locals())

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

47 r"""

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

50 Parameters

51 ----------

52 u : dtype_u

53 Current values of the numerical solution at time t.

54 du : dtype_u

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

56 t : float

57 Current time of the numerical solution.

59 Returns

60 -------

61 f : dtype_f

62 Current value of the right-hand side of f (which includes two components).

63 """

64 f = self.dtype_f(self.init)

65 f[:] = (

66 u[0] + self.eta * t * u[1] - np.sin(t),

67 du[0] + self.eta * t * du[1] + (1 + self.eta) * u[1] - np.cos(t),

68 )

69 self.work_counters['rhs']()

70 return f

72 def u_exact(self, t):

73 """

74 Routine for the exact solution.

76 Parameters

77 ----------

78 t : float

79 The time of the reference solution.

81 Returns

82 -------

83 me : dtype_u

84 The reference solution as mesh object containing two components.

85 """

86 me = self.dtype_u(self.init)

87 me[:] = (np.sin(t), 0)

88 return me

90 def du_exact(self, t):

91 """

92 Routine for the derivative of the exact solution.

94 Parameters

95 ----------

96 t : float

97 The time of the reference solution.

99 Returns

100 -------

101 me : dtype_u

102 The reference solution as mesh object containing two components.

103 """

105 me = self.dtype_u(self.init)

106 me[:] = (np.cos(t), 0)

107 return me