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

1import numpy as np 

2 

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

4 

5 

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 

10 

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}, 

13 

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}, 

16 

17 .. math:: 

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

19 

20 The exact solution of this system is 

21 

22 .. math:: 

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

24 

25 .. math:: 

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

27 

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

29 page 267 of [1]_. 

30 

31 Parameters 

32 ---------- 

33 nvars : int 

34 Number of unknowns of the system of DAEs. 

35 newton_tol : float 

36 Tolerance for Newton solver. 

37 

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

43 

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

45 """Initialization routine""" 

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

47 

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

51 

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. 

60 

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) 

69 

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 

77 

78 def u_exact(self, t): 

79 """ 

80 Routine for the exact solution. 

81 

82 Parameters 

83 ---------- 

84 t : float 

85 The time of the reference solution. 

86 

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 

96 

97 def du_exact(self, t): 

98 """ 

99 Routine for the derivative of the exact solution. 

100 

101 Parameters 

102 ---------- 

103 t : float 

104 The time of the reference solution. 

105 

106 Returns 

107 ------- 

108 me : dtype_u 

109 The reference solution as mesh object containing three components. 

110 """ 

111 

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