jacobwilliams / rklib Goto Github PK

Add some code to automatically change the color and/or line style for the methods in the plots, so we don't have to manually pick colors for each one.

Could add some Runge-Kutta Nystrom methods (for 2nd order equations).

Lines 3950-3954 of file rklib_variable_steps.f90 read:

real(wp),parameter :: d7   = 7170564158.0_wp / 30263027435.0_wp    ! 0.01158056612815422_wp
real(wp),parameter :: d8   = 7206303747.0_wp / 16758195910.0_wp    ! -0.7627079959184843_wp
real(wp),parameter :: d9   = -1059739258.0_wp / 8472387467.0_wp    ! 2.046330367018225_wp
real(wp),parameter :: d10  = 16534129531.0_wp / 11550853505.0_wp   ! -4.163198889384351_wp
real(wp),parameter :: d11  = -3.0_wp / 2.0_wp                      ! 2.901200440989918_wp

The decimal value shown in the comment section of each of these lines does not agree with the rational number on the same line, sometimes even the sign is wrong. Look, for example, at d11.

Line 1561 of file rklib_variable_steps.f90 contains
call me%f(t+h, x,f1)
The correct version is
call me%f(t, x,f1)

See #26

Allow for the additional root finder options to be optional user inputs for the even finding integration.

I think it would be helpful if the method classes included a character component storing the name of the method. Among other things, we could then run tests in a more programmatic manner.

Actually, pushing the idea one step further and combing it with #22, I was wondering if it would make sense to add the properties of a method as components of the corresponding concrete class:

type, extends (rk_fixed_step_class), public :: euler_class
    !! Basic Euler method
    integer :: p =1 ! order
    integer :: ncache = 1 ! cache size
    character(:), allocatable ::  name = "Euler"
    contains
    procedure :: step => euler
end type euler_class

I believe this would also avoid the need to have $N$ "method_order" procedures, whose sole purpose is to return p.

Some of the methods have dense interpolates. Those should be added to support dense output.

Currently, step methods have an interface like so:

subroutine step_func_fixed(me,t,x,h,xf)
!! rk step function for the fixed-step methods.
import :: rk_fixed_step_class,wp
implicit none
  class(rk_fixed_step_class),intent(inout) :: me
  real(wp),intent(in)                      :: t  !! initial time
  real(wp),dimension(me%n),intent(in)      :: x  !! initial state vector
  real(wp),intent(in)                      :: h  !! time step \( |\Delta t| \)
  real(wp),dimension(me%n),intent(out)     :: xf !! final state vector
end subroutine step_func_fixed

There is no need for a separate arguments x intent(in) and xf, intent(out). They can be merged into a single argument x, intent(inout), thereby avoiding a redundant state vector.

For instance, for Euler:

module procedure euler
associate (f1 => me%funcs(:,1))
    call me%f(t,x,f1)
    x = x + h*f1
end associate
end procedure euler

The change should be more or less straightforward for all methods expressed in the classical the rk form. For SSP and LS methods, one registry will probably have to be added because the implementations were already making use of the extra x. I can help with that; just let me know what you think.

Currently the event-finding algorithm uses the brent_solver from roots-fortran. This could be an optional input, so the caller can select any of the methods in that library.

Hi,

Could you add an example to show how to select an integration method?

Best regards,

The step procedures require a number of auxiliary variables to carry out the corresponding calculations, for instance f1,...f6 in the algorithm below.

module procedure rkf45
real(wp),dimension(me%n) :: f1,f2,f3,f4,f5,f6
...

For each time step, memory to store these variables needs to be allocated and then freed up. For very large ODE sets, these vectors will be huge and the corresponding memory operations can be significantly time consuming.

For optimal efficiency, the corresponding memory should be allocated once and then reused at each step, without repeated (de)allocation. AFAIU, DifferentialEquations.js also does so as well.

In principle, the abstract rk_class could include an allocatable component real(wp), allocatable :: cache(:,:). The cache would be allocated during the instantiation of the concrete method class. Inside the step method, we could then associate:

module procedure rkf45
associate (f1=>me%cache(:,1), f2=>me%cache(:,2), ... )
...

Quite likely, I am missing some important aspects, but I hope the main idea makes some sense.

In the variable-step method definitions, the symbols x and t are used to the denote the dependent and independent variables, respectively.

In variable-step methods the truncation error (on $x$) is denoted terr. I suppose the first letter "t" stands for "truncation", but my eyes+brain automatically link first letter "t" to the symbol t, i.e. to time. So, my reflex is to read: "time error".

Perhaps, you could rename terr as xerr, to make it more obvious that we are talking about an error on x.

See #31

Investigate the cases where this test is failing and see if there is some error.

Currently, the ones failing are:

    Method  Variable                      xf  Success
 ----------------------------------------------------
      rkr4         F   1.699041129050761E+00        F
      rkc5         F   1.555555555555555E+00        F
      rkl5         F   1.388888888888889E+00        F
    rklk5a         F   1.416666666666667E+00        F
     rko10         F   1.333333333333333E+00        F
     rkh10         F   1.333333333333333E+00        F
    rkdp65         T   1.423611111111112E+00        F
    rkv65e         T   9.999999999997726E-01        F
     rkv65         T   1.438888888888889E+00        F
   dverk65         T   1.450000000000000E+00        F
    rktp75         T   7.022258831547581E-01        F
    rktmy7         T   2.639828876496722E+00        F
    rkv76r         T   1.331980864835947E+00        F
   dverk78         T   1.000000000000023E+00        F
    rkdp85         T   1.076664947110369E+00        F
    rktp86         T   2.389109029661114E-01        F
    rkv87e         T   9.999999999999840E-01        F
    rkv87r         T   1.357836517333566E+00        F
    rkev87         T   1.351341629155199E+00        F
     rkk87         T   1.370168115647311E+00        F
    rkt98a         T   9.999999999999893E-01        F
    rkv98r         T   1.131299675056499E+00        F
     rks98         T   1.266647547720950E+00        F
    rkc108         T   1.333333333333333E+00        F
    rkb109         T   1.307444093579318E+00        F
  rks1110a         T   1.281633046469373E+00        F

You could add https://jacobwilliams.github.io/rklib/ to the website field of the About section of your github repo to make it easier for users to find the doc.

I am not 100% on this because you are using slightly different step size computation formulae. In step computation routine for each integrator (variable-step), you compute the step errors as this

xerr = sum (h*e(i)*F(i))

Then in step size computation routine (rklib_module.F90) , the max_err is computed like this

max_err = abs(xerr*h/tol)

So you are saying that step error for your integrators is O(h^2) and not O(h). Is that correct? Since you have included so many different algorithms so it may be you are using a different step size formulae than what is more typically used (esp. with DOP54 and DOP853). In case your formula is correct then you can close this issue.

Should have some more comprehensive tests of all the methods, and their performance comparisons for different problems.

I've just built the lib with fpm build --verbose and noticed some warnings about unused parameters in rklib_variable_steps.f90. Maybe it's nothing, but it could also be the consequence of some typo (e.g., one parameter is used two twice and another none).

I moved stuff into Fortran submodules because I didn't want everything all in one file. But then I didn't like all that boilerplate that had to be in the main file, so I moved those into .inc files, which isn't great either.

Is there a better way?

Also, GitHub thinks the .inc files are C++. Is there an extension it would recognize as a Fortran include file?

All fixed-step methods have the following check at the start:

if (h==zero) then
    xf = x
    return
end if

I have the impression that this check is not required and actually a bit "detrimental":

• strictly speaking, a real==real comparison is questionable.
• if such a check were necessary, wouldn't it be preferable to move it up the procedure hierarchy, to avoid repeating it in every step call and every method implementation (DRY principle)?
• all methods must anyway be consistent, in the sense that dx = xf -x must be 0 for h=0. So, there is no need to protect against a step call with h=0.

Hi Jacob, a fantastic effort!

I have a very stupid question; How does one pass non-constant variables to a function where the ODE is solved?

X' = AX + B

where A and B cannot simply be hard-coded but vary with time.

Your van der Pol oscillator has mu hardcoded as 0.2.... min ODE has varying constants..

Greetings

Jan Theron

In the example where the lib is used to solve the van der Pol oscillator, the independent variable is x. This choice of symbol is counter-intuitive because the oscillator is a function of time. Thus it would be preferable to use t, i.e. $y(t)$ or $x(t)$.