Edit Newton Solver as well as plotter function

cbcoutinho · Dec 26, 2016 · eebf44c · eebf44c
1 parent ac9ab05
commit eebf44c
Show file tree

Hide file tree

Showing 4 changed files with 32 additions and 9 deletions.
diff --git a/src/nonlinalg.f90 b/src/nonlinalg.f90
@@ -21,17 +21,20 @@ function nonlinsolve(myfun, n, x0) result(x)
     real(wp), intent(in), dimension(n)    :: x0
     procedure(fun_interf)                 :: myfun
 
-    integer                   :: ii
+    integer                   :: ii, maxit = 100
     real(wp), dimension(n)    :: x, z, b
     real(wp), dimension(n,n)  :: jac
     real(wp), parameter       :: eps = epsilon(1e0)
 
     x = x0
+    ii = 1
 
-    do ii = 1, 100
+    do
 
       b = myfun(n,x)
-      if ( ii == 1 .or. modulo(ii, 1) == 0 ) jac = jacobian(myfun, n, x)
+      if ( ii == 1 .or. modulo(ii, 1) == 0 ) then
+        jac = jacobian(myfun, n, x)
+      end if
 
       call linsolve_quick(n, jac, 1, -b, z)
 
@@ -41,6 +44,14 @@ function nonlinsolve(myfun, n, x0) result(x)
 
       ! if ( norm2(myfun(n, x)) < eps ) exit
       if ( norm2(z) < eps ) exit
+
+      if ( ii == maxIt ) then
+        print*, " *** Error: Maximum Number of Iterations Reached *** "
+        print*, " *** Number of Iterations in multivariate Newton's Method reached *** "
+        stop
+      end if
+
+      ii = ii + 1
     end do
 
   end function nonlinsolve

diff --git a/src/test/plotter.py b/src/test/plotter.py
@@ -30,10 +30,20 @@
 y3  = df.y3.values
 
 ax1 = fig.add_subplot(211)
-ax1.plot(t, y1, marker='o', linestyle='-', label='y1')
-ax1.plot(t, y2, marker='o', linestyle='-', label='y2')
-ax1.plot(t, y3, marker='o', linestyle='-', label='y3')
-# ax1.plot(t, x3, marker='o', linestyle='None', label='x3')
+ax1.plot(t, y1, marker='o', linestyle='None', label='y1', color='b')
+ax1.plot(t, y2, marker='o', linestyle='None', label='y2', color='g')
+ax1.plot(t, y3, marker='o', linestyle='None', label='y3', color='r')
+
+t   = df_raw.t.values
+y1  = df_raw.y1.values
+y2  = df_raw.y2.values
+y3  = df_raw.y3.values
+
+ax1.plot(t, y1, marker='None', linestyle='-', color='b')
+ax1.plot(t, y2, marker='None', linestyle='-', color='g')
+ax1.plot(t, y3, marker='None', linestyle='-', color='r')
+
+
 
 ax1.set_xlabel('t')
 ax1.set_ylabel('y')

diff --git a/src/test/src/main.f90 b/src/test/src/main.f90
@@ -6,7 +6,7 @@ program main
   implicit none
 
   integer                       :: ii, ios
-  integer, parameter            :: n = 3, num_t = 2
+  integer, parameter            :: n = 3, num_t = 10
   real(wp)                      :: t0, tend
   real(wp), dimension(num_t)    :: t
   real(wp), dimension(n)        :: y0
@@ -15,15 +15,16 @@ program main
 
   ! Set time vector `t`
   t0   = 0.0_wp
-  ! tend    = 4_wp*pi
   tend = 100._wp
+  ! tend = 4_wp*pi
   call linspace(t0, tend, t)
 
   ! Set y vector `y` using y0
   ! y0 = [0._wp, 0.1_wp]
   ! y0 = [2._wp, 0._wp]
   ! y0 = [ 0._wp ]
   ! y0 = [0.1_wp, 10.0_wp, 10.0_wp]
+  ! y0      = [1._wp, 0._wp, -1._wp]
   y0      = [1._wp, 0._wp, 0._wp]
   y(1,:)  = y0
 

diff --git a/src/test/src/misc.f90 b/src/test/src/misc.f90
@@ -43,6 +43,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