first commit

Makefile

@@ -0,0 +1,18 @@
+
+FC = gfortran
+FLAGS = -O3  # -fbounds-check
+TARGET = test 
+
+SOURCES = $(wildcard *.f90)
+OBJS = $(SOURCES:.f90=.o)
+
+%.o: %.f90 
+	$(FC) $(FLAGS) -c  $<
+
+all: $(OBJS)
+	$(FC) -o $(TARGET) $(OBJS) $(LIBS)
+
+clean:
+	rm *.o *.mod $(TARGET)
+
+include make.dep

README.md

@@ -0,0 +1 @@
+# ode

caotic_pendulum.f90

@@ -0,0 +1,69 @@
+module caotic_pendulum
+  use precision
+  use functions
+  use solvers
+  implicit none
+
+  real(dp), parameter :: pi = 3.14159265358979323846264338327_dp 
+  real(dp), parameter :: eps = 1.d-10
+
+contains
+
+  subroutine solve(Nstep, Nperiods, x0, v0)
+    integer, intent(in) :: Nstep, Nperiods    
+    real(dp), intent(in) :: x0, v0
+
+    real(dp) :: t, tfin, ttrans, dt, err
+    real(dp), allocatable :: u(:), u0(:), uex(:)
+    integer :: Ntot, Ntrans, id, id2
+    integer :: i
+
+    dt = 2.0_dp*pi/(w*real(Nstep, dp))
+
+    Ntot = Nstep*Nperiods  
+    tfin = Ntot*dt
+    Ntrans = int(Nperiods*0.0)
+    ttrans = Nstep*Ntrans*dt
+
+    write(*,*) 'dt=',dt
+
+    allocate(u0(2))
+    allocate(uex(2))
+    u0 = (/x0, v0/)
+    u = u0
+    t = 0.0_dp
+
+    do while (t<ttrans)
+       call dopri54(pendolo, t, dt, u0, u, err)
+       !call rk4(pendolo, t, dt, u0, u)
+       t = t + dt
+       u0 = u
+    end do
+
+
+    open(newunit=id, file='sol.dat')
+    !open(newunit=id2, file='poincare.dat', position='APPEND')
+
+    i = 0
+    do while (t<tfin)
+       write(id,'(F20.8,4(ES20.8))') t, u, u(2)*u(2)/2.0_dp + 1.0_dp - cos(u(1)), err
+       write(*,'(F20.8,3(ES20.8))') t, u, err
+
+       !if (mod(i,Nstep) == 0 ) then
+       !   write(id2,*) Q, u(2)
+       !end if
+
+       call dopri54(pendolo, t, dt, u0, u, err)
+       !call rk4(pendolo, t, dt, u0, u)
+       t = t + dt
+       i = i + 1
+       u0 = u
+    end do
+
+    close(id)
+    !close(id2)
+
+  end subroutine solve
+
+end module caotic_pendulum
+

functions.f90

@@ -0,0 +1,81 @@
+module functions
+  use precision
+  implicit none
+  private
+
+  public :: ff
+  public :: pendolo
+  public :: harmonic
+
+  real(dp), public :: k1 
+  real(dp), public :: k2
+  real(dp), public :: Q
+  real(dp), public :: A
+  real(dp), public :: w  
+
+  contains
+
+  subroutine fsub(t,u,up)
+    real(dp), intent(in) :: t    
+    real(dp), intent(in) :: u(:)     
+    real(dp), intent(out) :: up(:)       
+
+    up(1) = u(2)
+    up(2) = k1*u(1) + k2*u(2)
+
+  end subroutine fsub
+
+  function f1(t,u) result(up)
+    real(dp), intent(in) :: t    
+    real(dp), intent(in) :: u(:)
+    real(dp), allocatable :: up(:)
+
+    allocate(up(size(u)))
+
+    up(1) = k1*u(1) 
+
+  end function f1
+
+
+  function ff(t,u) result(up)
+    real(dp), intent(in) :: t    
+    real(dp), intent(in) :: u(:)
+    real(dp), allocatable :: up(:)
+
+    allocate(up(size(u)))
+
+    up(1) = u(2)
+    up(2) = k1*u(1) + k2*u(2)
+
+  end function ff
+
+  function pendolo(t,u) result(up)
+    real(dp), intent(in) :: t    
+    real(dp), intent(in) :: u(:)
+    real(dp), allocatable :: up(:)
+
+    allocate(up(size(u)))
+
+    up(1) = u(2)
+    up(2) = -sin(u(1)) - u(2)/Q + A*cos(w*t)
+
+  end function pendolo
+
+    function harmonic(t,u) result(up)
+    real(dp), intent(in) :: t    
+    real(dp), intent(in) :: u(:)
+    real(dp), allocatable :: up(:)
+
+    allocate(up(size(u)))
+
+    up(1) = u(2)
+    up(2) = -u(1) - u(2)/Q + A*cos(w*t)
+
+  end function harmonic
+
+
+
+
+end module functions
+
+

job

@@ -0,0 +1,26 @@
+#!/bin/bash
+
+
+echo ./test 0.01 100 rk2
+./test 0.01 100 rk2
+echo ./test 0.01 100 rk4
+./test 0.01 100 rk4
+echo ./test 0.01 100 dp54
+./test 0.01 100 dp54
+echo ./test 0.01 100 dp87
+./test 0.01 100 dp87
+
+
+echo ./test 0.001 1000 rk2
+./test 0.001 1000 rk2
+echo ./test 0.001 1000 rk4
+./test 0.001 1000 rk4
+echo ./test 0.001 1000 dp54
+./test 0.001 1000 dp54
+echo ./test 0.001 1000 dp87
+./test 0.001 1000 dp87
+
+echo ./test 0.001 100000 dp54
+./test 0.001 100000 dp54
+echo ./test 0.001 100000 dp87
+./test 0.001 100000 dp87

make.dep

@@ -0,0 +1,9 @@
+caotic_pendulum.o: precision.o
+caotic_pendulum.o: functions.o
+caotic_pendulum.o: solvers.o
+functions.o: precision.o
+solvers.o: precision.o
+solvers.o: precision.o
+test.o: precision.o
+test.o: functions.o
+test.o: solvers.o

makedep

@@ -0,0 +1,5 @@
+#!/bin/bash
+
+grep use *.f90 | awk '{print gensub(".f90:",".o:",1,$1), $3".o"}' 
+
+

precision.f90

@@ -0,0 +1,6 @@
+module precision
+  ! Double precision
+  integer, parameter, public :: dp = selected_real_kind(14) 
+  ! Quadrupole precision
+  !integer, parameter, public :: dp = selected_real_kind(30) 
+end module precision