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

## 40 statements

, created at 2024-09-19 09:13 +0000

1import numpy as np

3from pySDC.core.errors import ProblemError

4from pySDC.core.problem import Problem

5from pySDC.implementations.datatype_classes.mesh import mesh

8# noinspection PyUnusedLocal

9class logistics_equation(Problem):

10 r"""

11 Problem implementing a specific form of the Logistic Differential Equation

13 .. math::

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

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

17 given by

19 .. math::

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

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

38 dtype_u = mesh

39 dtype_f = mesh

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

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

55 )

57 def u_exact(self, t):

58 r"""

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

61 Parameters

62 ----------

63 t : float

64 Time of the exact solution.

66 Returns

67 -------

68 me : dtype_u

69 The exact solution.

70 """

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

76 def eval_f(self, u, t):

77 """

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

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.

87 Returns

88 -------

89 f : dtype_f

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

91 """

93 f = self.dtype_f(self.init)

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

95 return f

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

98 """

99 Simple Newton solver for the nonlinear equation.

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

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)

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

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

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

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

143 # increase iteration count

144 n += 1

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)

151 if n == self.newton_maxiter:

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

154 return u