diff --git a/caotic_pendulum.f90 b/caotic_pendulum.f90 index 61d72a5..0be591f 100644 --- a/caotic_pendulum.f90 +++ b/caotic_pendulum.f90 @@ -17,6 +17,9 @@ subroutine solve(Nstep, Nperiods, x0, v0) real(dp), allocatable :: u(:), u0(:), uex(:) integer :: Ntot, Ntrans, id, id2 integer :: i + + ! w0 = 1 + ! dt = 0.01 dt = 2.0_dp*pi/(w*real(Nstep, dp)) diff --git a/functions.f90 b/functions.f90 index a2356fb..0a40eb3 100644 --- a/functions.f90 +++ b/functions.f90 @@ -6,12 +6,23 @@ module functions public :: ff public :: pendolo public :: harmonic + public :: func real(dp), public :: k1 real(dp), public :: k2 real(dp), public :: Q real(dp), public :: A - real(dp), public :: w + real(dp), public :: w + + interface + function func(t,u) result(up) + use precision + real(dp), intent(in) :: t + real(dp), intent(in) :: u(:) + real(dp), allocatable :: up(:) + end function + end interface + contains diff --git a/make.dep b/make.dep index 7e54030..d39d716 100644 --- a/make.dep +++ b/make.dep @@ -2,6 +2,9 @@ caotic_pendulum.o: precision.o caotic_pendulum.o: functions.o caotic_pendulum.o: solvers.o functions.o: precision.o +functions.o: precision.o +solvers.o: precision.o +solvers.o: functions.o solvers.o: precision.o solvers.o: precision.o test.o: precision.o diff --git a/solvers.f90 b/solvers.f90 index 7be0e7b..3c28592 100644 --- a/solvers.f90 +++ b/solvers.f90 @@ -1,22 +1,37 @@ module solvers use precision + use functions implicit none private public :: rk2 public :: rk4 + public :: dopri54 public :: dopri87 - interface - function func(t,u) result(up) - use precision - real(dp), intent(in) :: t - real(dp), intent(in) :: u(:) - real(dp), allocatable :: up(:) - end function - end interface + public :: solver + interface + subroutine solver(ff, t, dt, u0, u, err) + use precision + interface + function ff(t,u) result(up) + use precision + real(dp), intent(in) :: t + real(dp), intent(in) :: u(:) + real(dp), allocatable :: up(:) + end function + end interface + real(dp), intent(in) :: t + real(dp), intent(in) :: dt + real(dp), intent(in) :: u0(:) + real(dp), intent(inout) :: u(:) + real(dp), intent(inout) :: err + end subroutine solver + end interface + + contains subroutine rk2(f, t, dt, u0, u) diff --git a/test.f90 b/test.f90 index 994acb2..e51fd85 100644 --- a/test.f90 +++ b/test.f90 @@ -4,33 +4,16 @@ program test use solvers implicit none - interface - subroutine solver(ff, t, dt, u0, u, err) - use precision - interface - function ff(t,u) result(up) - use precision - real(dp), intent(in) :: t - real(dp), intent(in) :: u(:) - real(dp), allocatable :: up(:) - end function - end interface - real(dp), intent(in) :: t - real(dp), intent(in) :: dt - real(dp), intent(in) :: u0(:) - real(dp), intent(inout) :: u(:) - real(dp), intent(inout) :: err - end subroutine solver - end interface - real(dp) :: err, t, dt, x0, v0, AA, BB, xex real(dp), allocatable :: u(:), u0(:) - integer :: i, Nstep + integer :: i, Nstep, funit procedure(solver), pointer :: psolver character(50) :: arg + character(10) :: solname + logical :: booleano ! = .true. || .false. - if (command_argument_count() < 3) then - write(*,*) "test dt Nstep method" + if (command_argument_count() < 8) then + write(*,*) "test dt Nstep method Q A w x0 v0" stop endif @@ -40,35 +23,56 @@ end subroutine solver call get_command_argument(2,arg) read(arg,*) Nstep - call get_command_argument(3,arg) + call get_command_argument(3,solname) + + call get_command_argument(4,arg) + read(arg,*) Q + + call get_command_argument(5,arg) + read(arg,*) A + + call get_command_argument(6,arg) + read(arg,*) w + + call get_command_argument(7,arg) + read(arg,*) x0 + + call get_command_argument(8,arg) + read(arg,*) v0 + allocate(u(2)) allocate(u0(2)) - k1 = -1.0_dp - k2 = 0.0_dp - x0 = 9.0_dp - v0 = 0.0_dp + ! u0(1) = u0 + ! u0(2) = v0 + u0 = (/x0, v0/) t = 0.0_dp err = 0.0_dp - AA = x0 - !AA = (v0 - k2*x0) / (k1-k2) - !BB = (x0*k1 - v0) / (k1-k2) - select case(trim(arg)) + ! arg ='rk2 <-- trailing white spaces -->' + !if (trim(arg) == 'rk2') then + ! ... + !else if (trim(arg) = 'rk4') then + ! ... + !end if + + open(newunit=funit, file="solution.dat") + + select case(trim(solname)) case('rk2') - do i = 1, Nstep - call rk2(ff, t, dt, u0, u) + do i = 1, Nstep + write(funit,*) t, u0(1), u0(2) + call rk2(harmonic, t, dt, u0, u) t = t+dt - xex = AA*cos(-k1*t) u0 = u end do case('rk4') do i = 1, Nstep - call rk4(ff, t, dt, u0, u) + write(funit,*) t, u0(1), u0(2) + call rk4(harmonic, t, dt, u0, u) t = t+dt - xex = AA*cos(-k1*t) u0 = u end do case('dp54') @@ -86,10 +90,14 @@ end subroutine solver u0 = u end do case default - write(*,*) 'error: solver ',trim(arg),' not found' + write(*,*) 'error: solver ',trim(arg),' not found' + stop end select - write(*,*) ' time computed exact error step-error' - write(*,'(f10.4,4(es18.8))') t, u0(1), xex, abs(u0(1)-xex), err + + close(funit) + + !write(*,*) ' time computed exact error step-error' + !write(*,'(f10.4,4(es18.8))') t, u0(1), xex, abs(u0(1)-xex), err end program test