# Stiff3 - adaptive solver for stiff systems of differential equations **URL:** \ **Category:** Announcements\ **Created:** [August 12, 2021, 1:27pm UTC](https://fortran-lang.discourse.group/t/stiff3-adaptive-solver-for-stiff-systems-of-differential-equations/1687 "2021-08-12T13:27:32Z")\ **Posts on this page:** 10\ **Page:** 1
**Author:** ![ivanpribec](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/ivanpribec/32/3290_2.png) [@ivanpribec](https://fortran-lang.discourse.group/u/ivanpribec)\ **Post date:** [August 12, 2021, 1:27pm UTC](https://fortran-lang.discourse.group/t/stiff3-adaptive-solver-for-stiff-systems-of-differential-equations/1687/1 "2021-08-12T13:27:32Z")
[`stiff3`](https://github.com/ivan-pi/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. > **[GitHub - ivan-pi/stiff3: Adaptive solver for stiff systems of ODEs using...](https://github.com/ivan-pi/stiff3)** > > Adaptive solver for stiff systems of ODEs using three-step semi-implicit Runge-Kutta method - GitHub - ivan-pi/stiff3: Adaptive solver for stiff systems of ODEs using three-step semi-implicit Runge... 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](https://en.wikipedia.org/wiki/Van_der_Pol_oscillator). Here are some plots of the solutions and phase cycle. ![image](https://global.discourse-cdn.com/free1/uploads/fortran_lang/original/1X/090d22ec1a7baf0ce3c16b48b78574a68fe33cc3.png) ![image](https://global.discourse-cdn.com/free1/uploads/fortran_lang/original/1X/e1df40208939096ac71b6de84f3db2358225d45f.png) ![image](https://global.discourse-cdn.com/free1/uploads/fortran_lang/original/1X/c9780f9664d66bbc976917ac5b87b3c2e186288a.png) 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](https://github.com/ivan-pi/stiff3/blob/eeaeefce09777029e2a36e751dd9f0716646805a/src/stiff3_solver.f90#L264) 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](https://scicomp.stackexchange.com/questions/37412/jacobian-matrix-cutoff-in-ode-solver), but so far it is still unanswered. ---
**Author:** ![certik](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/certik/32/4_2.png) [@certik](https://fortran-lang.discourse.group/u/certik)\ **Post date:** [August 12, 2021, 2:06pm UTC](https://fortran-lang.discourse.group/t/stiff3-adaptive-solver-for-stiff-systems-of-differential-equations/1687/2 "2021-08-12T14:06:14Z")
Nice! Thanks for creating the library. We need more of such libraries. 🙂 ---
**Author:** ![Arjen](https://avatars.discourse-cdn.com/v4/letter/a/b9bd4f/32.png) [@Arjen](https://fortran-lang.discourse.group/u/Arjen)\ **Post date:** [August 12, 2021, 2:14pm UTC](https://fortran-lang.discourse.group/t/stiff3-adaptive-solver-for-stiff-systems-of-differential-equations/1687/3 "2021-08-12T14:14:57Z")
You could code the do-loops in question as: ```auto 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 🙂 ---
**Author:** ![ivanpribec](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/ivanpribec/32/3290_2.png) [@ivanpribec](https://fortran-lang.discourse.group/u/ivanpribec)\ **Post date:** [August 12, 2021, 2:37pm UTC](https://fortran-lang.discourse.group/t/stiff3-adaptive-solver-for-stiff-systems-of-differential-equations/1687/4 "2021-08-12T14:37:55Z")
> [@Arjen](#): > > ```auto > where ( abs(df) > df_tol ) > df = -h * a * df > elsewhere > df = 0.0_wp > endwhere > > ``` 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. ---
**Author:** ![Ashok](https://avatars.discourse-cdn.com/v4/letter/a/ed655f/32.png) [@Ashok](https://fortran-lang.discourse.group/u/Ashok)\ **Post date:** [August 12, 2021, 3:08pm UTC](https://fortran-lang.discourse.group/t/stiff3-adaptive-solver-for-stiff-systems-of-differential-equations/1687/5 "2021-08-12T15:08:40Z")
Is LAPACK supported in FPM now ? ---
**Author:** ![ivanpribec](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/ivanpribec/32/3290_2.png) [@ivanpribec](https://fortran-lang.discourse.group/u/ivanpribec)\ **Post date:** [August 12, 2021, 3:16pm UTC](https://fortran-lang.discourse.group/t/stiff3-adaptive-solver-for-stiff-systems-of-differential-equations/1687/6 "2021-08-12T15:16:41Z")
Not exactly. But you can add it as a build dependency: ```auto [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. 🚀 ---
**Author:** ![lkedward](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/lkedward/32/72_2.png) [@lkedward](https://fortran-lang.discourse.group/u/lkedward)\ **Post date:** [August 12, 2021, 4:05pm UTC](https://fortran-lang.discourse.group/t/stiff3-adaptive-solver-for-stiff-systems-of-differential-equations/1687/7 "2021-08-12T16:05:04Z")
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. ---
**Author:** ![milancurcic](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/milancurcic/32/2_2.png) [@milancurcic](https://fortran-lang.discourse.group/u/milancurcic)\ **Post date:** [August 12, 2021, 5:00pm UTC](https://fortran-lang.discourse.group/t/stiff3-adaptive-solver-for-stiff-systems-of-differential-equations/1687/8 "2021-08-12T17:00:21Z")
Very nice, thank you. I suggest submitting stiff3 to [GitHub - fortran-lang/fpm-registry: Centralized registry of fpm packages](https://github.com/fortran-lang/fpm-registry). ---
**Author:** ![CRquantum](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/crquantum/32/730_2.png) [@CRquantum](https://fortran-lang.discourse.group/u/CRquantum)\ **Post date:** [January 8, 2022, 8:34am UTC](https://fortran-lang.discourse.group/t/stiff3-adaptive-solver-for-stiff-systems-of-differential-equations/1687/9 "2022-01-08T08:34:11Z")
Just curious, how is Stiff3 compared with DVODE by Byrne and Thompson? > [@Has anyone used DVODE and is it good?](https://fortran-lang.discourse.group/t/has-anyone-used-dvode-and-is-it-good/2281/11): > > Thank you very much for the message @nicholaswogan , true, strangely, the original link does not work any more. Fortunately, I have downloaded their whole website (the latest version) of dvode before. Now I pushed it to a new gitlab repo, you could download it from the link below, In the folder VODE\_F90 Support Page , double click index.html, then you should see the whole original website, and you can download all the corresponding files locally. Hope that helps slight_smile DVODE can solve stiff problem too. ---
**Author:** ![ivanpribec](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/ivanpribec/32/3290_2.png) [@ivanpribec](https://fortran-lang.discourse.group/u/ivanpribec)\ **Post date:** [January 9, 2022, 6:57pm UTC](https://fortran-lang.discourse.group/t/stiff3-adaptive-solver-for-stiff-systems-of-differential-equations/1687/10 "2022-01-09T18:57:23Z")
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.