Coverage for pySDC/implementations/problem_classes/LogisticEquation.py: 50%

1import numpy as np 

2 

3from pySDC.core.errors import ProblemError 

4from pySDC.core.problem import Problem 

5from pySDC.implementations.datatype_classes.mesh import mesh 

6 

7 

8# noinspection PyUnusedLocal 

9class logistics_equation(Problem): 

10 r""" 

11 Problem implementing a specific form of the Logistic Differential Equation 

12 

13 .. math:: 

14 \frac{du}{dt} = \lambda u(t)(1-u(t)) 

15 

16 with :math:`\lambda` a given real coefficient. Its analytical solution is 

17 given by 

18 

19 .. math:: 

20 u(t) = u(0) \frac{e^{\lambda t}}{1-u(0)+u(0)e^{\lambda t}} 

21 

22 Parameters 

23 ---------- 

24 u0 : float, optional 

25 Initial condition. 

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 direct : bool, optional 

31 Indicates if the problem should be solved directly or not. If False, it will be approximated by Newton. 

32 lam : float, optional 

33 Problem parameter :math:`\lambda`. 

34 stop_at_nan : bool, optional 

35 Indicates if the Newton solver stops when nan values arise. 

36 """ 

37 

38 dtype_u = mesh 

39 dtype_f = mesh 

40 

41 def __init__(self, u0=0.5, newton_maxiter=100, newton_tol=1e-12, direct=True, lam=1, stop_at_nan=True): 

42 """Initialization routine""" 

43 nvars = 1 

44 

45 super().__init__((nvars, None, np.dtype('float64'))) 

46 self._makeAttributeAndRegister( 

47 'u0', 

48 'lam', 

49 'newton_maxiter', 

50 'newton_tol', 

51 'direct', 

52 'nvars', 

53 'stop_at_nan', 

54 localVars=locals(), 

55 readOnly=True, 

56 ) 

57 

58 def u_exact(self, t): 

59 r""" 

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

61 

62 Parameters 

63 ---------- 

64 t : float 

65 Time of the exact solution. 

66 

67 Returns 

68 ------- 

69 me : dtype_u 

70 The exact solution. 

71 """ 

72 

73 me = self.dtype_u(self.init) 

74 me[:] = self.u0 * np.exp(self.lam * t) / (1 - self.u0 + self.u0 * np.exp(self.lam * t)) 

75 return me 

76 

77 def eval_f(self, u, t): 

78 """ 

79 Routine to compute the right-hand side of the problem. 

80 

81 Parameters 

82 ---------- 

83 u : dtype_u 

84 Current values of the numerical solution. 

85 t : float 

86 Current time of the numerical solution is computed. 

87 

88 Returns 

89 ------- 

90 f : dtype_f 

91 The right-hand side of the problem (contains one component). 

92 """ 

93 

94 f = self.dtype_f(self.init) 

95 f[:] = self.lam * u * (1 - u) 

96 return f 

97 

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

99 """ 

100 Simple Newton solver for the nonlinear equation. 

101 

102 Parameters 

103 ---------- 

104 rhs : dtype_f) 

105 Right-hand side for the nonlinear system. 

106 dt : float 

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

108 u0 : dtype_u 

109 Initial guess for the iterative solver. 

110 t : float 

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

112 

113 Returns 

114 ------- 

115 u : dtype_u 

116 The solution as mesh. 

117 """ 

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

119 u = self.dtype_u(u0) 

120 

121 if self.direct: 

122 d = (1 - dt * self.lam) ** 2 + 4 * dt * self.lam * rhs 

123 u = (-(1 - dt * self.lam) + np.sqrt(d)) / (2 * dt * self.lam) 

124 return u 

125 

126 else: 

127 # start newton iteration 

128 n = 0 

129 res = 99 

130 while n < self.newton_maxiter: 

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

132 g = u - dt * self.lam * u * (1 - u) - rhs 

133 

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

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

136 if res < self.newton_tol or np.isnan(res): 

137 break 

138 

139 # assemble dg/du 

140 dg = 1 - dt * self.lam * (1 - 2 * u) 

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

142 u -= 1.0 / dg * g 

143 

144 # increase iteration count 

145 n += 1 

146 

147 if np.isnan(res) and self.stop_at_nan: 

148 raise ProblemError('Newton got nan after %i iterations, aborting...' % n) 

149 elif np.isnan(res): 

150 self.logger.warning('Newton got nan after %i iterations...' % n) 

151 

152 if n == self.newton_maxiter: 

153 raise ProblemError('Newton did not converge after %i iterations, error is %s' % (n, res)) 

154 

155 return u