Non-constant variables

[jantheron69]

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

[marin-sapunar]

I had a similar issue while trying to interface with the library and the documentation is not helpful in this regard. The solution is to extend the solver class with a class which will also hold the coefficients:

  !> Extended class holding coefficients
  type, extends(rktp86_class) :: my_rk_class
     real(wp) :: c1, c2
  end type my_rk_class

  contains

  !> RHS of van der Pol system with configurable coefficients
  subroutine fvpol(me,t,x,f)
    class(rk_class),intent(inout)     :: me
    real(wp),intent(in)               :: t
    real(wp),dimension(:),intent(in)  :: x
    real(wp),dimension(:),intent(out) :: f

    select type(self => me)
    class is (my_rk_class)
        write(*,*) 'Inside fvpol', self%c1, self%c2
        f(1) = x(2)
        f(2) = self%c1 * (self%c2 - x(1)**2) * x(2) - x(1)
    class default
        error stop "fvpol: wrong class type passed in"
    end select
  end subroutine fvpol

The prop%initialize method creates a new instance and destroys any data already stored in it so the code in the example needs to look like:

  ...
  type(my_rk_class) :: prop
  ...

  call prop%initialize(n=n,f=fvpol,rtol=[tol],atol=[tol])
  ! Set coefficients for this instance only after initializing it
  prop%c1 = 0.2_wp
  prop%c2 = 1.0_wp
  call prop%integrate(t0,x0,dt,tf,xf)
  call prop%status(message=message)
  ...

This seems like it would be a quite common problem so perhaps it would be useful to update the documentation to address it?