diff --git a/NumLib/ODESolver/ConvergenceCriterion.cpp b/NumLib/ODESolver/ConvergenceCriterion.cpp
index e89a854f690c813793bccbeac42b0647c18b3988..42191dc81cfc76931da4b918dc970f44f949aa20 100644
--- a/NumLib/ODESolver/ConvergenceCriterion.cpp
+++ b/NumLib/ODESolver/ConvergenceCriterion.cpp
@@ -37,4 +37,12 @@ std::unique_ptr<ConvergenceCriterion> createConvergenceCriterion(
     OGS_FATAL("There is no convergence criterion of type `%s'.", type.c_str());
 }
 
+bool checkRelativeTolerance(const double reltol, const double numerator,
+                            const double denominator)
+{
+    auto const eps = std::numeric_limits<double>::epsilon();
+    return std::abs(numerator) <
+           std::abs(reltol) * (std::abs(denominator) + eps);
+}
+
 }  // NumLib
diff --git a/NumLib/ODESolver/ConvergenceCriterion.h b/NumLib/ODESolver/ConvergenceCriterion.h
index aaf07e8c15fe6e4210c6be548c3f31cb78ab5eb7..45d473d8823067246a3efda1825f8338f6761aca 100644
--- a/NumLib/ODESolver/ConvergenceCriterion.h
+++ b/NumLib/ODESolver/ConvergenceCriterion.h
@@ -77,6 +77,13 @@ public:
 std::unique_ptr<ConvergenceCriterion> createConvergenceCriterion(
     BaseLib::ConfigTree const& config);
 
+//! Returns if |numerator/denominator| < |reltol|.
+//! This method copes with the case that denominator = 0 by always adding
+//! epsilon to the denominator.
+bool checkRelativeTolerance(double const reltol,
+                            double const numerator,
+                            double const denominator);
+
 } // namespace NumLib
 
 #endif  // NUMLIB_CONVERGENCECRITERION_H
diff --git a/NumLib/ODESolver/ConvergenceCriterionDeltaX.cpp b/NumLib/ODESolver/ConvergenceCriterionDeltaX.cpp
index 5ab4ab6edbd8a525e97f9bd33f6f68c7da97f490..2c2249c1dbda3905ddd98dc47eb0749bde18cafe 100644
--- a/NumLib/ODESolver/ConvergenceCriterionDeltaX.cpp
+++ b/NumLib/ODESolver/ConvergenceCriterionDeltaX.cpp
@@ -45,7 +45,7 @@ void ConvergenceCriterionDeltaX::checkDeltaX(const GlobalVector& minus_delta_x,
         satisfied_abs = error_dx < *_abstol;
     }
     if (_reltol) {
-        satisfied_rel = error_dx < *_reltol * norm_x;
+        satisfied_rel = checkRelativeTolerance(*_reltol, error_dx, norm_x);
     }
 
     _satisfied = _satisfied && (satisfied_abs || satisfied_rel);
diff --git a/NumLib/ODESolver/ConvergenceCriterionPerComponentDeltaX.cpp b/NumLib/ODESolver/ConvergenceCriterionPerComponentDeltaX.cpp
index fe5f73b4528eed9d0c34d26b3d33f09dfb6a233c..9615ab813d148c72e6f42fab3eaa2a23633b7ee6 100644
--- a/NumLib/ODESolver/ConvergenceCriterionPerComponentDeltaX.cpp
+++ b/NumLib/ODESolver/ConvergenceCriterionPerComponentDeltaX.cpp
@@ -58,7 +58,8 @@ void ConvergenceCriterionPerComponentDeltaX::checkDeltaX(
 
         satisfied_abs = satisfied_abs && error_dx < _abstols[global_component];
         satisfied_rel =
-            satisfied_rel && error_dx < _reltols[global_component] * norm_x;
+            satisfied_rel && checkRelativeTolerance(_reltols[global_component],
+                                                    error_dx, norm_x);
     }
 
     _satisfied = _satisfied && (satisfied_abs || satisfied_rel);
diff --git a/NumLib/ODESolver/ConvergenceCriterionPerComponentResidual.cpp b/NumLib/ODESolver/ConvergenceCriterionPerComponentResidual.cpp
index dc94fb9b6e9483e89a125e717512a503a0ad45e5..155857c71fcd28b14204a171d7e65391b12e32fa 100644
--- a/NumLib/ODESolver/ConvergenceCriterionPerComponentResidual.cpp
+++ b/NumLib/ODESolver/ConvergenceCriterionPerComponentResidual.cpp
@@ -65,9 +65,10 @@ void ConvergenceCriterionPerComponentResidual::checkResidual(
         }
 
         satisfied_abs = satisfied_abs && norm_res < _abstols[global_component];
-        satisfied_rel = satisfied_rel &&
-                        norm_res < _reltols[global_component] *
-                                       _residual_norms_0[global_component];
+        satisfied_rel =
+            satisfied_rel &&
+            checkRelativeTolerance(_reltols[global_component], norm_res,
+                                   _residual_norms_0[global_component]);
     }
 
     _satisfied = _satisfied && (satisfied_abs || satisfied_rel);
diff --git a/NumLib/ODESolver/ConvergenceCriterionResidual.cpp b/NumLib/ODESolver/ConvergenceCriterionResidual.cpp
index ff0d719c078b8aff1554489b287bfa7991f67fe9..141245a44921daee3ab4bc123d4435dc192023c7 100644
--- a/NumLib/ODESolver/ConvergenceCriterionResidual.cpp
+++ b/NumLib/ODESolver/ConvergenceCriterionResidual.cpp
@@ -48,7 +48,8 @@ void ConvergenceCriterionResidual::checkResidual(const GlobalVector& residual)
         satisfied_abs = norm_res < *_abstol;
     }
     if (_reltol && !_is_first_iteration) {
-        satisfied_rel = norm_res < *_reltol * _residual_norm_0;
+        satisfied_rel =
+            checkRelativeTolerance(*_reltol, norm_res, _residual_norm_0);
     }
 
     _satisfied = _satisfied && (satisfied_abs || satisfied_rel);