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 nvars = 1 

43 

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

45 self._makeAttributeAndRegister( 

46 'u0', 

47 'lam', 

48 'newton_maxiter', 

49 'newton_tol', 

50 'direct', 

51 'nvars', 

52 'stop_at_nan', 

53 localVars=locals(), 

54 readOnly=True, 

55 ) 

56 

57 def u_exact(self, t): 

58 r""" 

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

60 

61 Parameters 

62 ---------- 

63 t : float 

64 Time of the exact solution. 

65 

66 Returns 

67 ------- 

68 me : dtype_u 

69 The exact solution. 

70 """ 

71 

72 me = self.dtype_u(self.init) 

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

74 return me 

75 

76 def eval_f(self, u, t): 

77 """ 

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

79 

80 Parameters 

81 ---------- 

82 u : dtype_u 

83 Current values of the numerical solution. 

84 t : float 

85 Current time of the numerical solution is computed. 

86 

87 Returns 

88 ------- 

89 f : dtype_f 

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

91 """ 

92 

93 f = self.dtype_f(self.init) 

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

95 return f 

96 

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

98 """ 

99 Simple Newton solver for the nonlinear equation. 

100 

101 Parameters 

102 ---------- 

103 rhs : dtype_f) 

104 Right-hand side for the nonlinear system. 

105 dt : float 

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

107 u0 : dtype_u 

108 Initial guess for the iterative solver. 

109 t : float 

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

111 

112 Returns 

113 ------- 

114 u : dtype_u 

115 The solution as mesh. 

116 """ 

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

118 u = self.dtype_u(u0) 

119 

120 if self.direct: 

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

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

123 return u 

124 

125 else: 

126 # start newton iteration 

127 n = 0 

128 res = 99 

129 while n < self.newton_maxiter: 

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

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

132 

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

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

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

136 break 

137 

138 # assemble dg/du 

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

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

141 u -= 1.0 / dg * g 

142 

143 # increase iteration count 

144 n += 1 

145 

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

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

148 elif np.isnan(res): 

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

150 

151 if n == self.newton_maxiter: 

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

153 

154 return u