Coverage for pySDC/implementations/problem_classes/Auzinger_implicit.py: 100%

1import numpy as np 

2 

3from pySDC.core.Problem import ptype 

4from pySDC.implementations.datatype_classes.mesh import mesh 

5 

6 

7# noinspection PyUnusedLocal 

8class auzinger(ptype): 

9 r""" 

10 This class implements the Auzinger equation [1]_ as initial value problem. The system of two ordinary differential equations 

11 (ODEs) is given by 

12 

13 .. math:: 

14 \frac{d y_1 (t)}{dt} = -y_2 (t) + y_1 (t) (1 - y^2_1 (t) - y^2_2 (t)), 

15 

16 .. math:: 

17 \frac{d y_2 (t)}{dt} = y_1 (t) + 3 y_2 (t) (1 - y^2_1 (t) - y^2_2 (t)) 

18 

19 with initial condition :math:`(y_1(t), y_2(t))^T = (1, 0)^T` for :math:`t \in [0, 10]`. The exact solution of this problem is 

20 

21 .. math:: 

22 (y_1(t), y_2(t))^T = (\cos(t), \sin(t))^T. 

23 

24 Attributes 

25 ---------- 

26 newton_maxiter : int, optional 

27 Maximum number of iterations for Newton's method. 

28 newton_tol : float, optional 

29 Tolerance for Newton's method to terminate. 

30 

31 References 

32 ---------- 

33 .. [1] doi.org/10.2140/camcos.2015.10.1 

34 """ 

35 

36 dtype_u = mesh 

37 dtype_f = mesh 

38 

39 def __init__(self, newton_maxiter=1e-12, newton_tol=100): 

40 """Initialization routine""" 

41 

42 # invoke super init, passing dtype_u and dtype_f, plus setting number of elements to 2 

43 super().__init__((2, None, np.dtype('float64'))) 

44 self._makeAttributeAndRegister('newton_maxiter', 'newton_tol', localVars=locals(), readOnly=True) 

45 

46 def u_exact(self, t): 

47 r""" 

48 Routine to compute the exact solution at time :math:`t`. 

49 

50 Parameters 

51 ---------- 

52 t : float 

53 Time of the exact solution. 

54 

55 Returns 

56 ------- 

57 me : dtype_u 

58 The exact solution. 

59 """ 

60 

61 me = self.dtype_u(self.init) 

62 me[0] = np.cos(t) 

63 me[1] = np.sin(t) 

64 return me 

65 

66 def eval_f(self, u, t): 

67 """ 

68 Routine to compute the right-hand side of the problem for both components simultaneously. 

69 

70 Parameters 

71 ---------- 

72 u : dtype_u 

73 Current values of the numerical solution. 

74 t : float 

75 Current time of the numerical solution is computed (not used here). 

76 

77 Returns 

78 ------- 

79 f : dtype_f 

80 The right-hand side of the problem (contains two components). 

81 """ 

82 

83 x1 = u[0] 

84 x2 = u[1] 

85 f = self.dtype_f(self.init) 

86 f[0] = -x2 + x1 * (1 - x1**2 - x2**2) 

87 f[1] = x1 + 3 * x2 * (1 - x1**2 - x2**2) 

88 return f 

89 

90 def solve_system(self, rhs, dt, u0, t): 

91 """ 

92 Simple Newton solver for the nonlinear system. 

93 

94 Parameters 

95 ---------- 

96 rhs : dtype_f 

97 Right-hand side for the nonlinear system. 

98 dt : float 

99 Abbrev. for the node-to-node stepsize (or any other factor required). 

100 u0 : dtype_u 

101 Initial guess for the iterative solver. 

102 t : float 

103 Current time (e.g. for time-dependent BCs.) 

104 

105 Returns 

106 ------- 

107 u : dtype_u 

108 The solution as mesh. 

109 """ 

110 

111 # create new mesh object from u0 and set initial values for iteration 

112 u = self.dtype_u(u0) 

113 x1 = u[0] 

114 x2 = u[1] 

115 

116 # start newton iteration 

117 n = 0 

118 while n < self.newton_maxiter: 

119 # form the function g with g(u) = 0 

120 g = np.array( 

121 [ 

122 x1 - dt * (-x2 + x1 * (1 - x1**2 - x2**2)) - rhs[0], 

123 x2 - dt * (x1 + 3 * x2 * (1 - x1**2 - x2**2)) - rhs[1], 

124 ] 

125 ) 

126 

127 # if g is close to 0, then we are done 

128 res = np.linalg.norm(g, np.inf) 

129 

130 if res < self.newton_tol: 

131 break 

132 

133 # assemble dg and invert the matrix (yeah, I know) 

134 dg = np.array( 

135 [ 

136 [1 - dt * (1 - 3 * x1**2 - x2**2), -dt * (-1 - 2 * x1 * x2)], 

137 [-dt * (1 - 6 * x1 * x2), 1 - dt * (3 - 3 * x1**2 - 9 * x2**2)], 

138 ] 

139 ) 

140 

141 idg = np.linalg.inv(dg) 

142 

143 # newton update: u1 = u0 - g/dg 

144 u -= np.dot(idg, g) 

145 

146 # set new values and increase iteration count 

147 x1 = u[0] 

148 x2 = u[1] 

149 n += 1 

150 

151 return u 

152 

153 # def eval_jacobian(self, u): 

154 # 

155 # x1 = u[0] 

156 # x2 = u[1] 

157 # 

158 # dfdu = np.array([[1-3*x1**2-x2**2, -1-x1], [1+6*x2*x1, 3+3*x1**2-9*x2**2]]) 

159 # 

160 # return dfdu 

161 # 

162 # 

163 # def solve_system_jacobian(self, dfdu, rhs, factor, u0, t): 

164 # 

165 # me = mesh(2) 

166 # me = LA.spsolve(sp.eye(2) - factor * dfdu, rhs) 

167 # return me