Stiff3 - adaptive solver for stiff systems of differential equations

Announcements
1

stiff3 is a semi-implicit Runge Kutta solver with adaptive time-stepping originally from the book by Villadsen & Michelsen, Solutions of Differential Equation Models by Polynomial Approximation, Prentice-Hall, 1978.

The solver is now available as an fpm package and uses LAPACK for the matrix factorization and back-substitution in the semi-implicit time-stepping.

The solver is designed for autonomous problems of the form

\frac{dy}{dx} = f(y)

In case of a non-autonomous problem where x appears explicitly, i.e. y' = f(x,y), it suffices to switch to a new integration variable t and introduce the additional ODE

\frac{dx}{dt} = 1

The repository provides an example for the Van der Pol oscillator. Here are some plots of the solutions and phase cycle.

image
image1269Ă—313 31 KB

image
image1269Ă—313 29.1 KB

image

There is one remaining issue with the stiff3 solver where I’m seeking help. It has to do with a section of code that sets to zero some elements of the Jacobian matrix. I believe this section is just an outdated practice that has to do with legacy usage patterns. I’ve posted it as a question on SciComp StackExchange, but so far it is still unanswered.

7 Likes
2

Nice! Thanks for creating the library. We need more of such libraries. :slight_smile:

3

You could code the do-loops in question as:

where ( abs(df) > df_tol ) 
    df = -h * a * df
elsewhere
    df = 0.0_wp
endwhere

forall (i=1:n) 
    df(i,i) = df(i,i) + 1.0_wp
end forall

(If I remember the forall syntax properly)

Whether you actually want to do it that way, depends on your appreciation of these constructs :slight_smile:

4

Indeed, that is more elegant. However, I think the zero values should be set explicitly in the user-provided Jacobian routine, and there is no need for this chunk of code at all. It doesn’t make sense to me to hardcode a df_tol cutoff in the Jacobian matrix.

5

Is LAPACK supported in FPM now ?

6

Not exactly. But you can add it as a build dependency:

[build]
link = ["lapack", "blas"]

As long as you have a system-installed LAPACK library that the linker is able to locate with flag -llapack, then things should just work. :rocket:

7

Great stuff @ivanpribec! Nice to see the fpm ecosystem growing with more useful packages!

Regarding your final point, I’ve previously come across thresholding of small matrix elements for iterative methods - I believe it’s to avoid having very small eigenvalues which degrade convergence of the overall method or place severe limits on the region of stability.

As pointed out on the scicomp post, these small values creep in unavoidably due to floating-point arithmetic, but they are essentially zero because they are smaller than the effective numerical precision which gets smaller as matrix condition number increases. (I think I read somewhere that every order of magnitude in the condition number takes off one significant figure in numerical precision.) The value of the cut-off threshold was perhaps set based on some expected condition number.

2 Likes
8

Very nice, thank you. I suggest submitting stiff3 to GitHub - fortran-lang/fpm-registry: Centralized registry of fpm packages.

9

Just curious, how is Stiff3 compared with DVODE by Byrne and Thompson?

DVODE can solve stiff problem too.

10

I don’t have a direct comparison but I can imagine that VODE is much more robust, having originated from the group of Hindmarsh. VODE is also younger, published in 1989, while the stiff3 code I refactored was originally published in 1978.

There are several things missing in stiff3, including:

  • dense output of variables
  • support for banded or sparse Jacobians
  • more precise control over the error tolerances and time-stepping

Stiff3 was designed primarily for autonomous systems, i.e. those without an explicit “time” dependence. As a semi-implicit method it is also designed to work exclusively with analytic Jacobian matrices. Use of finite-difference approximations for the Jacobian matrix will not work.

Use VODE if you need a robust and mature ODE solver. Stiff3 is more of a teaching code, perhaps still useful for smaller problems.

1 Like