A snippet showing a combination of fpm libraries put together to replicate a Lorentz system
program main
use ogpf, only: gpf
use rklib_module
use iso_fortran_env, only: output_unit, wp => real64
implicit none
integer,parameter :: n = 3 !! dimension of the system
real(wp),parameter :: tol = 1.0e-12_wp !! integration tolerance
real(wp),parameter :: t0 = 0.0_wp !! initial t value
real(wp),parameter :: dt = 0.01_wp !! initial step size
integer, parameter :: num_steps = 10000 !! number of steps
real(wp),parameter :: x0(n) = [0.0_wp,0.1_wp,1.05_wp] !! initial x value
integer :: incr !! increment number
real(wp) :: time(0:num_steps)
real(wp) :: x(n,0:num_steps)
type(rkf45_class) :: prop !! rk integrator
type(gpf) :: gp !! gnuplot interface
! ---------------------------------------------
! Initialize buffers
time = [(t0+incr*dt,incr=0,num_steps)]
x(:,0) = x0; x(:,1:) = 0._wp
! ---------------------------------------------
! Initialize method
call prop%initialize(n=n,f=fvpol,rtol=[tol],atol=[tol])
do incr = 1, num_steps
call prop%integrate(time(incr-1),x(:,incr-1),dt,time(incr),x(:,incr))
end do
! static 2D multiplot
call gp%multiplot(2,2)
call gp%plot(x(1,:), x(2,:), 'with lines lt 5 lc rgb "blue"')
call gp%plot(x(1,:), x(3,:), 'with lines lt 5 lc rgb "green"')
call gp%plot(x(2,:), x(3,:), 'with lines lt 5 lc rgb "red"')
call gp%run_script()
contains
subroutine fvpol(me,t,x,f)
!! Right-hand side of Lorentz system
class(rk_class),intent(inout) :: me
real(wp),intent(in) :: t
real(wp),intent(in) :: x(:)
real(wp),intent(out) :: f(:)
real(wp),parameter :: s=10._wp, r=28._wp, b=2.667_wp
f(1) = s*(x(2) - x(1))
f(2) = r*x(1) - x(2) - x(1)*x(3)
f(3) = x(1)*x(2) - b*x(3)
end subroutine fvpol
end program main
Which produces:
To get here:
- create a fpm project
> fpm new lorentz
- Add dependencies in fpm.toml
[dependencies]
rklib = { git="https://github.com/jacobwilliams/rklib.git" }
ogpf = { git = "https://github.com/kookma/ogpf.git" }
(it also requires having gnuplot installed)