Bug #10574 fixed - Runge-Kutta-Method failed for vector [x, 1] with x > 9. 29/6329/4
Allan CORNET [Fri, 17 Feb 2012 11:54:53 +0000 (12:54 +0100)]
Change-Id: I1e34feb810742c1b78c2d501dbf4e944aca70b6f

12 files changed:
scilab/CHANGES_5.4.X
scilab/modules/differential_equations/sci_gateway/fortran/sci_f_ode.f
scilab/modules/differential_equations/src/c/DllmainDifferential_equations.c
scilab/modules/differential_equations/src/c/rk4.c
scilab/modules/differential_equations/src/c/rk4.h
scilab/modules/differential_equations/src/fortran/lsrgk.f
scilab/modules/differential_equations/src/fortran/odeint.f
scilab/modules/differential_equations/src/fortran/rkqc.f
scilab/modules/differential_equations/tests/nonreg_tests/bug_10574.dia.ref [new file with mode: 0644]
scilab/modules/differential_equations/tests/nonreg_tests/bug_10574.tst [new file with mode: 0644]
scilab/modules/differential_equations/tests/unit_tests/ode.dia.ref
scilab/modules/differential_equations/tests/unit_tests/ode.tst

index 8ab1cb1..ffaef38 100644 (file)
@@ -372,6 +372,8 @@ Bug fixes
 * Bug #10565 fixed - demo simulation --> bicycle simulation --> unstable
                      trajectory failed.
 
+* Bug #10574 fixed - Runge-Kutta-Method failed for vector [x, 1] with x > 9.
+
 * Bug #10576 fixed - getdate (on Windows) did not manage dates after January 19, 2038
 
 * Bug #10577 fixed - getdate did not check input argument.
index c5ca391..7ffe3f2 100644 (file)
@@ -579,7 +579,7 @@ c     lsdisc
          liw=1
       elseif(meth.eq.5) then
 c     lsrgk
-         lrw=3*ny
+         lrw=9*ny
          liw=1
       elseif(meth.eq.6) then
 c     rkf45
index edc485b..50124c3 100644 (file)
@@ -142,10 +142,6 @@ DIFFERENTIAL_EQUATIONS_IMPEXP struct {
     int jupbnd;
 } C2F(dqa001);
 
-DIFFERENTIAL_EQUATIONS_IMPEXP struct {
-    int kmax, kount;
-    double dxsav, xp[200], yp[2000]    /* was [10][200] */;
-} C2F(path);
 
 DIFFERENTIAL_EQUATIONS_IMPEXP struct {
     int iero;
index 71c3b2f..6abc13b 100644 (file)
 /*--------------------------------------------------------------------------*/
 /* Runge-Kutta (RK4) method */
 /* http://media4.obspm.fr/public/DEA/cours/chapitre3/souschapitre2/section4/page5/section3_2_4_5.html */
+/*     The original version has been modified to replace statically
+ *    allocated arrays yt, dym and dyt by rwork arguments parts
+ *    array + blas use. Serge Steer INRIA- feb 2012*/
+
 /*--------------------------------------------------------------------------*/
-int C2F(rk4)(double *y, double *dydx, int *n,double *x, double *h, double *yout,void (*derivs)())
+int C2F(rk4)(double *y, double *dydx, int *n,double *x, double *h, double *yout,void (*derivs)(), double *rwork)
 {
     double d = 0.0;
     int i;
-    double h6 = 0.0, hh = 0.0, xh = 0.0, yt[10], dym[10], dyt[10];
-
-    /* Parameter adjustments (fortran)*/
-    --yout;
-    --dydx;
-    --y;
+    double h6 = 0.0, hh = 0.0, xh = 0.0;
 
+    double *yt = rwork;
+    double *dym = rwork+*n;
+    double *dyt = rwork+2*(*n);
     C2F(ierode).iero = 0;
     hh = *h * 0.5;
     h6 = *h / 6.0;
     xh = *x + hh;
-
-    for (i = 1; i <= *n; ++i) yt[i - 1] = y[i] + hh * dydx[i];
-
+    for (i = 0; i < *n; ++i) yt[i] = y[i] + hh * dydx[i];
     (*derivs)(n, &xh, yt, dyt);
 
     if (C2F(ierode).iero > 0) return 0;
 
-    for (i = 1; i <= *n; ++i) yt[i - 1] = y[i] + hh * dyt[i - 1];
-
+    for (i = 0; i < *n; ++i) yt[i] = y[i] + hh * dyt[i];
     (*derivs)(n, &xh, yt, dym);
 
     if (C2F(ierode).iero > 0) return 0;
 
-    for (i = 1; i <= *n; ++i) 
+    for (i = 0; i < *n; ++i) 
     {
-               yt[i - 1] = y[i] + *h * dym[i - 1];
-               dym[i - 1] = dyt[i - 1] + dym[i - 1];
+               yt[i] = y[i] + *h * dym[i];
+               dym[i] = dyt[i] + dym[i];
     }
     d = *x + *h;
-
     (*derivs)(n, &d, yt, dyt);
 
     if (C2F(ierode).iero > 0) return 0;
 
-    for (i = 1; i <= *n; ++i) yout[i] = y[i] + h6 * (dydx[i] + dyt[i - 1] + dym[i - 1] * 2.0);
-
+    for (i = 0; i < *n; ++i) yout[i] = y[i] + h6 * (dydx[i] + dyt[i] + dym[i] * 2.0);
     return 0;
 }
 /*--------------------------------------------------------------------------*/
index ce1046f..f020fa3 100644 (file)
 * @param h
 * @param yout
 * @param derivs
+* @param rwork
 * @return 0
 */
-DIFFERENTIAL_EQUATIONS_IMPEXP int C2F(rk4)(double *y, double *dydx, int *n,double *x, double *h, double *yout,void (*derivs)());
+DIFFERENTIAL_EQUATIONS_IMPEXP int C2F(rk4)(double *y, double *dydx, int *n,double *x, double *h, double *yout,void (*derivs)(), double *rwork);
 /*--------------------------------------------------------------------------*/
 #endif /* __RK4_H__ */
 /*--------------------------------------------------------------------------*/
index f195bc1..67c0563 100644 (file)
@@ -1,7 +1,10 @@
 c     ====================================
-c     runge kutta d'ordre 4 adaptatif 
-c     l'interface lsrgk a ete fait en s'inspirant de lsode 
-c     voir lsode.f pour comprendre le sens des variables 
+c     ode Gateway for Adaptative fourth order Runge Kutta 
+C
+c     The original version has been modified to pass the  rwork
+c     argument to odeint
+c     array + blas use. Serge Steer INRIA- feb 2012
+
 c     ====================================
       subroutine lsrgk (f, neq, y, t, tout, itol, rtol, atol, itask,
      1            istate, iopt, rwork, lrw, iwork, liw, jac, mf)
@@ -14,7 +17,8 @@ c     ====================================
       integer iero
       common/ierode/iero
       iero=0
-      call odeint(y,neq,t,tout,atol(1),1.0d-4,0.0d0,nok,nbad,f,rkqc)
+      call odeint(y, neq, t,tout,atol(1),1.0d-4,0.0d0,nok,nbad,f,rkqc,
+     $     rwork)
       t=tout
       if (iero.gt.0) istate=-1
       return
index b746f73..3bf2d07 100644 (file)
@@ -1,17 +1,42 @@
 c     ====================================
-      subroutine odeint(ystart,nvar,x1,x2,eps,h1,hmin,nok,nbad,derivs,rk
-     *qc)
-C      implicit undefined (a-z)
+c .   Runge-Kutta driver with adaptive stepsize control.  Integrate
+c .   the starting values ystart(1:nvar) from x1 to x2 with accuracy 
+c .   eps, storing intermediate results in the common block /path/.
+c .   h1 should be set as a guessed first stepsize, hmin as the 
+c .   minimum allowed stepsize (can be zero).
+c .   On output, nok and nbad are the number of good and bad (but
+c .   retried and fixed) steps taken, and ystart is replaced by 
+c .   values at the end of the integration interval.
+c .   derivs is the user-supplied subroutine for calculating the
+c .   right-hand side derivatives, while rkqs is the name of 
+c .   the stepper routine to be used.  
+c .   /path/ contains its own information about how often an
+c .   intermediate value is to be stored.
+c .
+c .   The original version has been modified to replace statically
+c .   allocated arrays y, yscal and dxdy by parts of rwork argument
+c     array + blas use. Serge Steer INRIA- feb 2012
+
+c     ====================================
+      subroutine odeint(ystart,nvar,x1,x2,eps,h1,hmin,nok,nbad,derivs,
+     $     rkqc,rwork)
       external derivs,rkqc
-      integer maxstp,nmax,kmax,kount,nvar,i,nok,nbad,nstp
-      double precision two,zero,tiny,dxsav,xp(200),yp(10,200),x,h
-      parameter (maxstp=10000,nmax=10,two=2.0,zero=0.0,tiny=1.e-30)
-      double precision x1,x2,eps,h1,hmin,xsav,hdid,hnext
-      double precision ystart(nvar),yscal(nmax),y(nmax),dydx(nmax)
+      integer maxstp,kount,nvar,i,nok,nbad,nstp
+      double precision two,zero,tiny,x,h
+      parameter (maxstp=10000,two=2.0d0,zero=0.0d0,tiny=1.d-30)
+      double precision x1,x2,eps,h1,hmin,hdid,hnext
+      double precision ystart(nvar)
+      double precision rwork(*)
       character*80 messag
-      common /path/ kmax,kount,dxsav,xp,yp
+      integer ly,lyscal,ldydx,lwork
       integer iero
       common/ierode/iero
+c     
+      ly=1
+      lyscal=ly+nvar
+      ldydx=lyscal+nvar
+      lwork=ldydx+nvar
+
       iero=0
       if ( abs(x2-x1).le.tiny) return
       x=x1
@@ -19,60 +44,39 @@ C      implicit undefined (a-z)
       nok=0
       nbad=0
       kount=0
-      do 11 i=1,nvar
-        y(i)=ystart(i)
-11    continue
-      xsav=x-dxsav*two
+      call dcopy(nvar,ystart,1,rwork(ly),1)
+
       do 16 nstp=1,maxstp
-        call derivs(nvar,x,y,dydx)
-        if (iero.gt.0) return 
-        do 12 i=1,nvar
-          yscal(i)=abs(y(i))+abs(h*dydx(i))+tiny
-12      continue
-        if(kmax.gt.0)then
-          if(abs(x-xsav).gt.abs(dxsav)) then
-            if(kount.lt.kmax-1)then
-              kount=kount+1
-              xp(kount)=x
-              do 13 i=1,nvar
-                yp(i,kount)=y(i)
-13            continue
-              xsav=x
-            endif
-          endif
-        endif
-        if((x+h-x2)*(x+h-x1).gt.zero) h=x2-x
-        call rkqc(y,dydx,nvar,x,h,eps,yscal,hdid,hnext,derivs)
-        if(iero.gt.0) return
-        if(hdid.eq.h)then
-          nok=nok+1
-        else
-          nbad=nbad+1
-        endif
-        if((x-x2)*(x2-x1).ge.zero)then
-          do 14 i=1,nvar
-            ystart(i)=y(i)
-14        continue
-          if(kmax.ne.0)then
-            kount=kount+1
-            xp(kount)=x
-            do 15 i=1,nvar
-              yp(i,kount)=y(i)
-15          continue
-          endif
-          return
-        endif
-        if(abs(hnext).lt.hmin) then
-           write(messag, 17) hnext,hmin
-           hnext=hmin
-        endif
-        h=hnext
-16    continue
+         call derivs(nvar,x,rwork(ly),rwork(ldydx))
+         if (iero.gt.0) return 
+         do 12 i=0,nvar-1
+            rwork(lyscal+i)=abs(rwork(ly+i))+abs(h*rwork(ldydx+i))+tiny
+ 12      continue
+
+         if((x+h-x2)*(x+h-x1).gt.zero) h=x2-x
+         call rkqc(rwork(ly),rwork(ldydx),nvar,x,h,eps,rwork(lyscal),
+     $        hdid,hnext,derivs,rwork(lwork))
+         if(iero.gt.0) return
+         if(hdid.eq.h)then
+            nok=nok+1
+         else
+            nbad=nbad+1
+         endif
+         if((x-x2)*(x2-x1).ge.zero)then
+            call dcopy(nvar,rwork(ly),1,ystart,1)
+            return
+         endif
+         if(abs(hnext).lt.hmin) then
+            write(messag, 17) hnext,hmin
+            hnext=hmin
+         endif
+         h=hnext
+ 16   continue
       iero=-1
-c      print *, 'Trop d''iterations a faire pour la precision demandee.'
+c     print *, 'Trop d''iterations a faire pour la precision demandee.'
       return
  17   format('stepsize ',e10.3,' smaller than minimum ',e10.3)
       end
 c     ====================================
       
-      
\ No newline at end of file
+      
index c34fc3d..8fc719a 100644 (file)
@@ -1,57 +1,90 @@
 c     ====================================
-      subroutine rkqc(y,dydx,n,x,htry,eps,yscal,hdid,hnext,derivs)
+      subroutine rkqc(y,dydx,n,x,htry,eps,yscal,hdid,hnext,derivs,rwork)
 C      implicit undefined (a-z)
-      integer nmax,n,i
+c ................................................................
+c .   Subroutine rkqc is a fifth-order Runge-Kutta step with 
+c .   monitoring of local truncation error to ensure accuracy
+c .   and adjust stepsize.   
+c .
+c .   Inputs:
+c .   y(1:n) = the dependent variable vector
+c .   dydx(1:n) = deriv of y wrt x at starting value of indep var x
+c .   htry = attempted step size
+c .   eps = required accuracy
+c .   yscal(1:n) = scaling vector for the error
+c .   
+c .   Outputs:
+c .   y = new value of y
+c .   x = new value of x
+c .   hdid = actual stepsize accomplished
+c .   hnext = estimated next stepsize   
+c .
+c .   derivs is the user-supplied subroutine that computes the 
+c .    right-hand side derivatives
+c ................................................................
+c
+c     The original version has been modified to replace statically
+c     allocated arrays dysav, ytemp and ysav by rwork arguments parts
+c     array + blas use. Serge Steer INRIA- feb 2012
+
+      integer n,i
       double precision fcor,one,safety,errcon
-      parameter (nmax=10,fcor=.0666666667,
-     $     one=1.0,safety=0.9,errcon=6.e-4)
+      parameter (fcor=.0666666667,one=1.0,safety=0.9,errcon=6.e-4)
       double precision x,htry,eps,hdid,hnext,pgrow,pshrnk,xsav,hh
-      double precision errmax,h,dysav(nmax)
-      double precision y(n),dydx(n),yscal(n),ytemp(nmax),ysav(nmax)
-     
+      double precision errmax,h
+      double precision y(n),dydx(n),yscal(n)
+      double precision rwork(*)
+
       external derivs
       integer iero
       common/ierode/iero
+
+      lysav=1
+      ldysav=lysav+n
+      lytemp=ldysav+n
+      lwork=lytemp+n
+
       iero=0
-      pgrow=-0.20
-      pshrnk=-0.25
+      pgrow=-0.20d0
+      pshrnk=-0.25d0
       xsav=x
-      do 11 i=1,n
-        ysav(i)=y(i)
-        dysav(i)=dydx(i)
-11    continue
+
+      call dcopy(n,y,1,rwork(lysav),1)
+      call dcopy(n,dydx,1,rwork(ldysav),1)
       h=htry
 1     hh=0.5*h
-      call rk4(ysav,dysav,n,xsav,hh,ytemp,derivs)
+      call rk4(rwork(lysav),rwork(ldysav),n,xsav,hh,rwork(lytemp),
+     $     derivs,rwork(lwork))
       x=xsav+hh
-      call derivs(n,x,ytemp,dydx)
+      call derivs(n,x,rwork(lytemp),dydx)
       if (iero.gt.0) return 
-      call rk4(ytemp,dydx,n,x,hh,y,derivs)
+      call rk4(rwork(lytemp),dydx,n,x,hh,y,derivs,rwork(lwork))
       x=xsav+h
       if(x.eq.xsav) then 
          iero=1
          return
       endif
-      call rk4(ysav,dysav,n,xsav,h,ytemp,derivs)
-      errmax=0.
-      do 12 i=1,n
-        ytemp(i)=y(i)-ytemp(i)
-        errmax=max(errmax,abs(ytemp(i)/(yscal(i)*eps)))
+      call rk4(rwork(lysav),rwork(ldysav),n,xsav,h,rwork(lytemp),
+     $     derivs,rwork(lwork))
+      errmax=0.0d0
+       do 12 i=0,n-1
+        rwork(lytemp+i)=y(i+1)-rwork(lytemp+i)
+        errmax=max(errmax,abs(rwork(lytemp+i)/(yscal(i+1)*eps)))
 12    continue
+       
       if(errmax.gt.one) then
         h=safety*h*(errmax**pshrnk)
         goto 1
       else
         hdid=h
-        if(errmax.gt.errcon)then
+        if(errmax.gt.errcon) then
           hnext=safety*h*(errmax**pgrow)
         else
-          hnext=4.*h
+          hnext=4.0d0*h
         endif
       endif
-      do 13 i=1,n
-        y(i)=y(i)+ytemp(i)*fcor
-13    continue
+      call daxpy(n,fcor,rwork(lytemp),1,y,1)
+
       return
       end
 c     ====================================
diff --git a/scilab/modules/differential_equations/tests/nonreg_tests/bug_10574.dia.ref b/scilab/modules/differential_equations/tests/nonreg_tests/bug_10574.dia.ref
new file mode 100644 (file)
index 0000000..83a04e8
--- /dev/null
@@ -0,0 +1,62 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2012 - DIGITEO - Allan CORNET
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+//
+// <-- Non-regression test for bug 10574 -->
+//
+// <-- Bugzilla URL -->
+//http://bugzilla.scilab.org/show_bug.cgi?id=10574
+//
+// <-- Short Description -->
+// Runge-Kutta-Method failed for vector [x, 1] with x > 9
+dt = 0.25;
+a = 1;
+b = a/dt;
+t = 0:dt:a;
+t0 = t(1);
+St = size(t);
+Dim = 18;
+x = ones(Dim, 1);
+X = eye(Dim,Dim);
+function xdot = Test(t, x)
+  xdot = X * x;
+endfunction
+x = ode('rk', x, t0, t, Test);
+assert_checkequal(size(x), [18 5]);
+REF_X = [    1.    1.2840254    1.6487213    2.117    2.7182818; ..
+    1.    1.2840254    1.6487213    2.117    2.7182818; ..
+    1.    1.2840254    1.6487213    2.117    2.7182818; ..
+    1.    1.2840254    1.6487213    2.117    2.7182818; ..
+    1.    1.2840254    1.6487213    2.117    2.7182818; ..
+    1.    1.2840254    1.6487213    2.117    2.7182818; ..
+    1.    1.2840254    1.6487213    2.117    2.7182818; ..
+    1.    1.2840254    1.6487213    2.117    2.7182818; ..
+    1.    1.2840254    1.6487213    2.117    2.7182818; ..
+    1.    1.2840254    1.6487213    2.117    2.7182818; ..
+    1.    1.2840254    1.6487213    2.117    2.7182818; ..
+    1.    1.2840254    1.6487213    2.117    2.7182818; ..
+    1.    1.2840254    1.6487213    2.117    2.7182818; ..
+    1.    1.2840254    1.6487213    2.117    2.7182818; ..
+    1.    1.2840254    1.6487213    2.117    2.7182818; ..
+    1.    1.2840254    1.6487213    2.117    2.7182818; ..
+    1.    1.2840254    1.6487213    2.117    2.7182818; ..
+    1.    1.2840254    1.6487213    2.117    2.7182818];
+assert_checkalmostequal(x, REF_X, 0, 1e-7);
+Dim = 9;
+x = ones(Dim, 1);
+X = eye(Dim, Dim);
+x = ode('rk', x, t0, t, Test);
+assert_checkequal(size(x), [9 5]);
+REF_X = [    1.    1.2840254    1.6487213    2.117    2.7182818; ..
+    1.    1.2840254    1.6487213    2.117    2.7182818; ..
+    1.    1.2840254    1.6487213    2.117    2.7182818; ..
+    1.    1.2840254    1.6487213    2.117    2.7182818; ..
+    1.    1.2840254    1.6487213    2.117    2.7182818; ..
+    1.    1.2840254    1.6487213    2.117    2.7182818; ..
+    1.    1.2840254    1.6487213    2.117    2.7182818; ..
+    1.    1.2840254    1.6487213    2.117    2.7182818; ..
+    1.    1.2840254    1.6487213    2.117    2.7182818];
+assert_checkalmostequal(x, REF_X, 0, 1e-7);
diff --git a/scilab/modules/differential_equations/tests/nonreg_tests/bug_10574.tst b/scilab/modules/differential_equations/tests/nonreg_tests/bug_10574.tst
new file mode 100644 (file)
index 0000000..7628fff
--- /dev/null
@@ -0,0 +1,72 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2012 - DIGITEO - Allan CORNET
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+//
+// <-- Non-regression test for bug 10574 -->
+//
+// <-- Bugzilla URL -->
+//http://bugzilla.scilab.org/show_bug.cgi?id=10574
+//
+// <-- Short Description -->
+// Runge-Kutta-Method failed for vector [x, 1] with x > 9
+
+dt = 0.25;
+a = 1;
+b = a/dt;
+t = 0:dt:a;
+t0 = t(1);
+St = size(t);
+Dim = 18;
+x = ones(Dim, 1);
+X = eye(Dim,Dim);
+
+function xdot = Test(t, x)
+  xdot = X * x;
+endfunction
+
+x = ode('rk', x, t0, t, Test);
+assert_checkequal(size(x), [18 5]);
+
+REF_X = [    1.    1.2840254    1.6487213    2.117    2.7182818; ..
+    1.    1.2840254    1.6487213    2.117    2.7182818; ..
+    1.    1.2840254    1.6487213    2.117    2.7182818; ..
+    1.    1.2840254    1.6487213    2.117    2.7182818; ..
+    1.    1.2840254    1.6487213    2.117    2.7182818; ..
+    1.    1.2840254    1.6487213    2.117    2.7182818; ..
+    1.    1.2840254    1.6487213    2.117    2.7182818; ..
+    1.    1.2840254    1.6487213    2.117    2.7182818; ..
+    1.    1.2840254    1.6487213    2.117    2.7182818; ..
+    1.    1.2840254    1.6487213    2.117    2.7182818; ..
+    1.    1.2840254    1.6487213    2.117    2.7182818; ..
+    1.    1.2840254    1.6487213    2.117    2.7182818; ..
+    1.    1.2840254    1.6487213    2.117    2.7182818; ..
+    1.    1.2840254    1.6487213    2.117    2.7182818; ..
+    1.    1.2840254    1.6487213    2.117    2.7182818; ..
+    1.    1.2840254    1.6487213    2.117    2.7182818; ..
+    1.    1.2840254    1.6487213    2.117    2.7182818; ..
+    1.    1.2840254    1.6487213    2.117    2.7182818];
+
+assert_checkalmostequal(x, REF_X, 0, 1e-7);
+
+
+Dim = 9;
+x = ones(Dim, 1);
+X = eye(Dim, Dim);
+x = ode('rk', x, t0, t, Test);
+assert_checkequal(size(x), [9 5]);
+
+
+REF_X = [    1.    1.2840254    1.6487213    2.117    2.7182818; ..
+    1.    1.2840254    1.6487213    2.117    2.7182818; ..
+    1.    1.2840254    1.6487213    2.117    2.7182818; ..
+    1.    1.2840254    1.6487213    2.117    2.7182818; ..
+    1.    1.2840254    1.6487213    2.117    2.7182818; ..
+    1.    1.2840254    1.6487213    2.117    2.7182818; ..
+    1.    1.2840254    1.6487213    2.117    2.7182818; ..
+    1.    1.2840254    1.6487213    2.117    2.7182818; ..
+    1.    1.2840254    1.6487213    2.117    2.7182818];
+
+assert_checkalmostequal(x, REF_X, 0, 1e-7);
index 7825af9..e416ffe 100644 (file)
@@ -334,6 +334,6 @@ rkRes = [0.    0.09983341664683527    0.19866933079512300    0.29552020666153711
 rkRes(17:32) = [0.99957360305287668    0.99166481046564392    0.97384763089330029    0.94630008770455387 0.90929742684492987    0.86320936667026738    0.80849640384312127    0.74570521220233155 0.67546318057873300    0.59847214413334937    0.51550137185246248    0.42737988026620155 0.33498815018939587    0.23924932924832210    0.14112000809477554    0.04158066246848785];
 rkfRes = [0.    0.09983341667063099    0.19866933087782307    0.2955202068043328    0.38941834246046680     0.47942553864900261    0.56464247314940919    0.64421768645637856    0.71735608928892303      0.78332690686480611    0.84147098056308622    0.89120735401875306    0.93203907784361117      0.96355817497519225    0.98544971704249074    0.99749497101988072    0.99957358472991464];
 rkfRes(18:32) = [    0.99166478935902302    0.97384760697150774    0.94630006094857066    0.90929739724116376     0.86320933420871493    0.80849636852157181    0.74570517403636893    0.67546313961625104      0.59847210047141619    0.51550132565379336    0.42737983177229999    0.33498809972769594      0.23924927723128114    0.14111995500988161    0.04158060885936282];
-assert_checkalmostequal(rkRes, rk, %eps, [], "matrix");
+assert_checkalmostequal(rkRes, rk, %eps * 20, [], "matrix");
 assert_checkalmostequal(rkfRes, rkf, %eps, [], "matrix");
 assert_checkalmostequal(rkf, fixx, %eps, [], "matrix");
index 04f4828..e62702f 100644 (file)
@@ -380,7 +380,7 @@ rkRes(17:32) = [0.99957360305287668    0.99166481046564392    0.9738476308933002
 rkfRes = [0.    0.09983341667063099    0.19866933087782307    0.2955202068043328    0.38941834246046680     0.47942553864900261    0.56464247314940919    0.64421768645637856    0.71735608928892303      0.78332690686480611    0.84147098056308622    0.89120735401875306    0.93203907784361117      0.96355817497519225    0.98544971704249074    0.99749497101988072    0.99957358472991464];
 rkfRes(18:32) = [    0.99166478935902302    0.97384760697150774    0.94630006094857066    0.90929739724116376     0.86320933420871493    0.80849636852157181    0.74570517403636893    0.67546313961625104      0.59847210047141619    0.51550132565379336    0.42737983177229999    0.33498809972769594      0.23924927723128114    0.14111995500988161    0.04158060885936282];
 
-assert_checkalmostequal(rkRes, rk, %eps, [], "matrix");
+assert_checkalmostequal(rkRes, rk, %eps * 20, [], "matrix");
 assert_checkalmostequal(rkfRes, rkf, %eps, [], "matrix");
 assert_checkalmostequal(rkf, fixx, %eps, [], "matrix");