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

1import numpy as np 

2 

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

4from pySDC.implementations.datatype_classes.mesh import mesh 

5 

6 

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 

11 

12 .. math:: 

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

14 

15 .. math:: 

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

17 

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

19 

20 Parameters 

21 ---------- 

22 nvars : int 

23 Number of unknowns of the system of DAEs. 

24 newton_tol : float 

25 Tolerance for Newton solver. 

26 

27 Attributes 

28 ---------- 

29 eta : float 

30 Specific parameter of the problem. 

31 

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 """ 

37 

38 dtype_u = mesh 

39 dtype_f = mesh 

40 

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()) 

45 

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)`. 

49 

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. 

58 

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 

71 

72 def u_exact(self, t): 

73 """ 

74 Routine for the exact solution. 

75 

76 Parameters 

77 ---------- 

78 t : float 

79 The time of the reference solution. 

80 

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 

89 

90 def du_exact(self, t): 

91 """ 

92 Routine for the derivative of the exact solution. 

93 

94 Parameters 

95 ---------- 

96 t : float 

97 The time of the reference solution. 

98 

99 Returns 

100 ------- 

101 me : dtype_u 

102 The reference solution as mesh object containing two components. 

103 """ 

104 

105 me = self.dtype_u(self.init) 

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

107 return me