A bunch of updates

-Moved rk constants to separate module
-Changed step calculation to be function of tableau order
-Made adaptive routine call explicit routine, modularize

.gitignore

@@ -1,5 +1,6 @@
 .ipynb_checkpoints
 data.out
+plot.png
 
 main
 

Makefile

@@ -21,6 +21,7 @@ FLIBS = -lblas -llapack
 # Dependencies of main program
 objects=$(OBJ)/misc.o \
 	$(OBJ)/runge_kutta.o \
+	$(OBJ)/rk_constants.o \
 	$(OBJ)/lib_array.o \
 	$(OBJ)/lib_constants.o
 
@@ -33,7 +34,9 @@ $(OBJ)/lib_constants.o: $(FORTRANLIB_SRC)/lib_constants.f90
 # Modules
 $(OBJ)/misc.o: $(SRC)/misc.f90 $(OBJ)/lib_constants.o
 	$(FF) $(FFLAGS) -J$(OBJ) -c -o $@ $<
-$(OBJ)/runge_kutta.o: $(SRC)/runge_kutta.f90
+$(OBJ)/runge_kutta.o: $(SRC)/runge_kutta.f90 $(OBJ)/rk_constants.o
+	$(FF) $(FFLAGS) -J$(OBJ) -c -o $@ $<
+$(OBJ)/rk_constants.o: $(SRC)/rk_constants.f90
 	$(FF) $(FFLAGS) -J$(OBJ) -c -o $@ $<
 
 # Main program

plotter.py

@@ -7,19 +7,33 @@
 
 filename = 'data.out'
 
+fig = plt.figure()
+
 df = pd.read_csv(filename,
-                 names=['t', 'x', 'y', 'z'],
-                 delim_whitespace=True)
+                 delim_whitespace=True,
+                 names=['t', 'y1', 'y2'])
+                #  names=['t', 'x1', 'x2', 'x3'])
+                #  names=['t', 'x', 'y', 'z'])
 
-# df.sort_values(by='t', inplace=True)
-# df.plot('t', ['x(t)', 'y(t)', 'z(t)'], linestyle='-', marker='o')
+df.sort_values(by='t', inplace=True)
 
-x = df.x.values
-y = df.y.values
-z = df.z.values
+t   = df.t.values
+x1  = df.y1.values
+x2  = df.y2.values
+# x3  = df.x3.values
 
-fig = plt.figure()
-ax = fig.add_subplot(111, projection='3d')
-ax.plot(x, y, z, label='lorenz attractor')
+ax = fig.add_subplot(111)
+ax.plot(t, x1, marker='o', linestyle='-', label='x1')
+ax.plot(t, x2, marker='o', linestyle='-', label='x2')
+# ax.plot(t, x3, marker='o', linestyle='None', label='x3')
+
+# tt = np.linspace(t[0], t[-1], 100)
+# ax.plot(tt, np.cos(tt), linestyle='-', label='cos(t)')
+# ax.plot(tt, np.sin(tt), linestyle='-', label='sin(t)')
+# ax.plot(tt, -np.sin(tt), linestyle='-', label='-cos(t)')
+
+# ax = fig.add_subplot(111, projection='3d')
+# ax.plot(x, y, z, label='lorenz attractor')
 
+plt.savefig('plot.png')
 plt.show()

src/main.f90

@@ -6,48 +6,45 @@ program main
   use lib_constants, only: pi => pi_dp
   implicit none
 
-  integer                 :: ii, ios
+  integer                       :: ii, ios
+  integer, parameter            :: n = 2, num_t = 100
+  real(wp)                      :: t0, tend
+  real(wp), dimension(num_t)    :: t
+  real(wp), dimension(n)        :: y0
+  real(wp), dimension(num_t,n)  :: y
 
-  integer, parameter      :: n = 3
-  real(wp)                :: t0, tend, dt
-  ! real(wp), dimension(num_t)    :: t
-  real(wp), dimension(n)  :: y0
-  ! real(wp), dimension(num_t,n)  :: y
+  ! Set time vector `t`
+  t0   = 0.0_wp
+  ! tend    = 4_wp*pi
+  tend = 20._wp
+  call linspace(t0, tend, t)
 
+  ! Set y vector `y` using y0
+  ! y0 = [0._wp, 1._wp]
+  y0 = [1.5_wp, 3._wp]
+  ! y0 = [0.1_wp, 0.0_wp, 0.0_wp]
+  ! y0      = [0._wp, 1._wp, 0._wp]
+  y(1,:)  = y0
+
+  call rk_wrapper(n, num_t, t, y, mysub)
 
 
   open(unit=21, file='data.out', iostat=ios, status="replace", action="write")
   if ( ios /= 0 ) stop "Error opening file 21"
 
-  t0 = 0.0_wp
-  dt = 0.1_wp
-  ! tend = 10_wp*pi
-  tend = 100._wp
-  y0 = [0.1_wp, 0.0_wp, 0.0_wp]
-
-  ! call linspace(t0, tend, t)
-  ! dt = t(2) - t(1)
-  ! do ii = 1, num_t-1
-  !   y0 = y(ii, :)
-  !   call rk_explicit(n, t(ii), dt, y0, mysub)
-  !   y(ii+1, :) = y0
-  ! end do
-  ! do ii = 1, num_t
-  !   write(21,*) t(ii), y(ii, :)
-  ! end do
-
-  write(21,*) t0, y0
-
-  call rk_wrapper(n, t0, tend, dt, y0, mysub)
-
-
-
-
+  open(unit=22, file='raw.out', iostat=ios, status="replace", action="write")
+  if ( ios /= 0 ) stop "Error opening file 22"
 
+  do ii = 1, num_t
+      write(21,*) t(ii), y(ii, :)
+  end do
 
   close(unit=21, iostat=ios)
   if ( ios /= 0 ) stop "Error closing file unit 21"
 
+  close(unit=22, iostat=ios, status="delete")
+  if ( ios /= 0 ) stop "Error closing file unit 22"
+
 
 
 end program main

src/misc.f90

@@ -14,7 +14,8 @@ subroutine mysub(n, t, y, dy)
     ! call simple_trig(n, t, y, dy)
     ! call simple_ode(n, t, y, dy)
     ! call vanderpol(n, t, y, dy)
-    call lorenz(n, t, y, dy)
+    ! call lorenz(n, t, y, dy)
+    call brusselator(n, t, y, dy)
 
 
     return
@@ -40,6 +41,7 @@ subroutine simple_ode(n, t, y, dy)
 
     dy(1) = y(2)
     dy(2) = -y(1)
+    dy(3) = -y(2)
 
     return
   end subroutine simple_ode
@@ -50,12 +52,12 @@ subroutine vanderpol(n, t, y, dy)
     real(wp), intent(in), dimension(n)  :: y
     real(wp), intent(out), dimension(n) :: dy
 
-    real(wp), parameter                 :: mu = 50.53_wp
+    real(wp), parameter                 :: mu = 30.53_wp
     real(wp), parameter                 :: A = 3.2_wp
     real(wp), parameter                 :: omega = 2._wp*pi/10_wp
 
     dy(1) = y(2)
-    dy(2) = mu * (1-y(1)*y(1)) * y(2) - y(1) + A*dsin(t*omega)
+    dy(2) = mu * (1._wp-y(1)*y(1)) * y(2) - y(1) + A*dsin(t*omega)
 
     return
   end subroutine vanderpol
@@ -77,4 +79,16 @@ subroutine lorenz(n, t, y, dy)
     return
   end subroutine lorenz
 
+  subroutine brusselator(n, t, y, dy)
+    integer, intent(in)                 :: n
+    real(wp), intent(in)                :: t
+    real(wp), intent(in), dimension(n)  :: y
+    real(wp), intent(out), dimension(n) :: dy
+
+    dy(1) = 1._wp + y(1)*y(1)*y(2) - 4._wp*y(1)
+    dy(2) = 3._wp*y(1) - y(1)*y(1)*y(2)
+
+    return
+  end subroutine brusselator
+
 end module misc