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

## 38 statements

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

1import numpy as np

3from pySDC.core.problem import Problem

4from pySDC.implementations.datatype_classes.mesh import mesh

7# noinspection PyUnusedLocal

8class auzinger(Problem):

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

13 .. math::

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

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

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

21 .. math::

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

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.

31 References

32 ----------

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

34 """

36 dtype_u = mesh

37 dtype_f = mesh

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

40 """Initialization routine"""

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

46 def u_exact(self, t):

47 r"""

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

50 Parameters

51 ----------

52 t : float

53 Time of the exact solution.

55 Returns

56 -------

57 me : dtype_u

58 The exact solution.

59 """

61 me = self.dtype_u(self.init)

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

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

64 return me

66 def eval_f(self, u, t):

67 """

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

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

77 Returns

78 -------

79 f : dtype_f

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

81 """

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

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

91 """

92 Simple Newton solver for the nonlinear system.

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

105 Returns

106 -------

107 u : dtype_u

108 The solution as mesh.

109 """

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]

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 )

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

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

130 if res < self.newton_tol:

131 break

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 )

141 idg = np.linalg.inv(dg)

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

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

146 # set new values and increase iteration count

147 x1 = u[0]

148 x2 = u[1]

149 n += 1

151 return u

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