I am studying an implementation of a semi-implicit Runge Kutta method of 3. order (SIRK3) from the book by Villadsen & Michelsen (1978), Solution of differential equation models by polynomial approximation, Prentice-Hall.

The code expects an analytically evaluated Jacobian matrix $J = \partial F / \partial y$. In the calculation of the new function values $y^{n+1}$ at time $t + h$, the matrix $I -h a J$ is assembled by overwriting the storage used by the Jacobian matrix (array df):

    do i = 1, n
      do j = 1, n
        df(i,j) = -h * a * df(i,j)
        if (abs(df(i,j)) < DF_TOL) df(i,j) = 0.0_dp
      end do
      df(i,i) = df(i,i) + 1.0_dp
    end do

In the original code the value of DF_TOL = 1.0d-12. Afterwards an LU factorization of the matrix df is performed and used to calculate the new function values by back-substitution.

What is the purpose of truncating the Jacobian matrix values, and why choose this particular cutoff value? Is this type of truncation needed at all?

  • $\begingroup$ I think it is to prevent false non-zeros (which may occur due to floating point arithmetic), but it should be considered in a relative sense. If all the jacobian entries are small (<1e-12) then it does not make sense to compare to 1e-12; it should be replaced by 1e-12*max(abs(df)). $\endgroup$ – Abdullah Ali Sivas May 13 at 19:01
  • $\begingroup$ The value 1e-12 still feels a bit arbitrary. Shouldn't it change with the target precision (i.e. single vs double), or be left as an optional user input? If it is known all entries are small, wouldn't good modelling practice require the original problem should be non-dimensionalized/scaled. $\endgroup$ – IPribec May 15 at 13:45
  • $\begingroup$ You are right, but remember that the book was published in 1978. The double precision standard was not readily available back then (officially adopted in 1985, even though there were some vendors offering it earlier and not fully adopted until a few years later), so the arbitrary (I agree with you on that point) tolerance 1e-12 is for single precision. Also, the book is on the solution of ODEs, so probably they did not care much about the intricacies of the floating point arithmetic. $\endgroup$ – Abdullah Ali Sivas May 15 at 19:56
  • $\begingroup$ Good point. It would be nice to learn if any modern ODE libraries use this type of cutoff or is this an outdated practice that I can simply remove in the refactored code. $\endgroup$ – IPribec May 15 at 20:16

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.