* Bug #13121 fixed - Differential_equations: ode crashed Scilab 74/13274/3
Paul Bignier [Fri, 29 Nov 2013 11:22:31 +0000 (12:22 +0100)]
ode "rk" option crashed Scilab when the user derivative function failed.

Change-Id: I89a5a1f8aa70447a0457560c9e7c35c236acb093

scilab/CHANGES_5.5.X
scilab/modules/differential_equations/src/fortran/rkqc.f
scilab/modules/differential_equations/tests/nonreg_tests/bug_13121.dia.ref [new file with mode: 0644]
scilab/modules/differential_equations/tests/nonreg_tests/bug_13121.tst [new file with mode: 0644]

index e68207c..5c292de 100644 (file)
@@ -277,6 +277,8 @@ Scilab Bug Fixes
 
 * Bug #13116 fixed - qpsolve now respects upper-bounds constraints.
 
+* Bug #13121 fixed - ode "rk" option crashed Scilab when the user derivative function failed.
+
 
 Xcos Bug Fixes
 ==============
index 8fc719a..e44bce7 100644 (file)
@@ -2,9 +2,9 @@ c     ====================================
       subroutine rkqc(y,dydx,n,x,htry,eps,yscal,hdid,hnext,derivs,rwork)
 C      implicit undefined (a-z)
 c ................................................................
-c .   Subroutine rkqc is a fifth-order Runge-Kutta step with 
+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 .   and adjust stepsize.
 c .
 c .   Inputs:
 c .   y(1:n) = the dependent variable vector
@@ -12,14 +12,14 @@ 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 .
 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 .   hnext = estimated next stepsize
 c .
-c .   derivs is the user-supplied subroutine that computes the 
+c .   derivs is the user-supplied subroutine that computes the
 c .    right-hand side derivatives
 c ................................................................
 c
@@ -55,23 +55,25 @@ c     array + blas use. Serge Steer INRIA- feb 2012
 1     hh=0.5*h
       call rk4(rwork(lysav),rwork(ldysav),n,xsav,hh,rwork(lytemp),
      $     derivs,rwork(lwork))
+      if (iero.gt.0) return
       x=xsav+hh
       call derivs(n,x,rwork(lytemp),dydx)
-      if (iero.gt.0) return 
+      if (iero.gt.0) return
       call rk4(rwork(lytemp),dydx,n,x,hh,y,derivs,rwork(lwork))
       x=xsav+h
-      if(x.eq.xsav) then 
+      if(x.eq.xsav) then
          iero=1
          return
       endif
       call rk4(rwork(lysav),rwork(ldysav),n,xsav,h,rwork(lytemp),
      $     derivs,rwork(lwork))
+      if (iero.gt.0) return
       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
diff --git a/scilab/modules/differential_equations/tests/nonreg_tests/bug_13121.dia.ref b/scilab/modules/differential_equations/tests/nonreg_tests/bug_13121.dia.ref
new file mode 100644 (file)
index 0000000..3c47a56
--- /dev/null
@@ -0,0 +1,69 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2013 - Scilab Enterprises - Paul Bignier
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+//
+// <-- CLI SHELL MODE -->
+//
+// <-- Non-regression test for bug 13121 -->
+//
+// <-- Bugzilla URL -->
+//http://bugzilla.scilab.org/show_bug.cgi?id=13121
+//
+// <-- Short Description -->
+// ode "rk" option crashed Scilab when the user derivative function failed.
+t0=0;
+x0=0; y0=0;
+// снаряд
+c = 0.2;
+s = 0.07;
+m0 = 15;
+tT=-1; mT=m0;
+v0 = 50; teta0 = 0.6;
+// ракета
+c = 0.2;
+s = 0.25;
+m0 = 500;
+tT = 1; mT = 250;
+v0 = 50; teta0 = 0.6;
+T0=10000;
+r = 1.29; g = 9.8;
+function y=m(t)
+    if t <= tT then
+        y = m0*(1-t/tT) + mT*t/tT;
+    else
+        y = mT;
+    end
+endfunction
+function y=der_m(t)
+    if t <= tT then
+        y = 1 - (m0 - mT)/tT;
+    else
+        y = 0;
+    end
+endfunction
+function y=T(t)
+    if t <= tT then
+        y = T0(1-t/tT);
+    else
+        y = 0;
+    end
+endfunction
+function ydot=right(t,pl)
+    x = pl;
+    x(1) = max(x(1),0);
+    x(2) = max(x(2),0);
+    //if x(2) <= 0 & t > 0 then
+    //    disp(t,"tf");
+    //end;
+      ydot(1) = x(3) * cos(x(4));
+      ydot(2) = x(3) * sin(x(4));
+      ydot(3) = 1/m(t)*(T(t)-0.5*c*r*s*x(3)*x(3)) - der_m(t)/m(t)*x(3) - g*sin(x(4));
+      ydot(4) = -g/x(3)*cos(x(4));
+//    end
+endfunction
+t=0:0.01:40;
+refMsg = msprintf(_("Invalid index."));
+assert_checkerror("res=ode(""rk"",[x0;y0;v0;teta0;],t0,t,right)", refMsg);
diff --git a/scilab/modules/differential_equations/tests/nonreg_tests/bug_13121.tst b/scilab/modules/differential_equations/tests/nonreg_tests/bug_13121.tst
new file mode 100644 (file)
index 0000000..aa75cd8
--- /dev/null
@@ -0,0 +1,78 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2013 - Scilab Enterprises - Paul Bignier
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+//
+// <-- CLI SHELL MODE -->
+//
+// <-- Non-regression test for bug 13121 -->
+//
+// <-- Bugzilla URL -->
+//http://bugzilla.scilab.org/show_bug.cgi?id=13121
+//
+// <-- Short Description -->
+// ode "rk" option crashed Scilab when the user derivative function failed.
+
+t0=0;
+x0=0; y0=0;
+
+// снаряд
+c = 0.2;
+s = 0.07;
+m0 = 15;
+tT=-1; mT=m0;
+v0 = 50; teta0 = 0.6;
+// ракета
+c = 0.2;
+s = 0.25;
+m0 = 500;
+tT = 1; mT = 250;
+v0 = 50; teta0 = 0.6;
+T0=10000;
+
+r = 1.29; g = 9.8;
+
+function y=m(t)
+    if t <= tT then
+        y = m0*(1-t/tT) + mT*t/tT;
+    else
+        y = mT;
+    end
+endfunction
+
+function y=der_m(t)
+    if t <= tT then
+        y = 1 - (m0 - mT)/tT;
+    else
+        y = 0;
+    end
+endfunction
+
+function y=T(t)
+    if t <= tT then
+        y = T0(1-t/tT);
+    else
+        y = 0;
+    end
+endfunction
+
+function ydot=right(t,pl)
+    x = pl;
+    x(1) = max(x(1),0);
+    x(2) = max(x(2),0);
+    //if x(2) <= 0 & t > 0 then
+    //    disp(t,"tf");
+    //end;
+    ydot(1) = x(3) * cos(x(4));
+    ydot(2) = x(3) * sin(x(4));
+    ydot(3) = 1/m(t)*(T(t)-0.5*c*r*s*x(3)*x(3)) - der_m(t)/m(t)*x(3) - g*sin(x(4));
+    ydot(4) = -g/x(3)*cos(x(4));
+    //    end
+endfunction
+
+t=0:0.01:40;
+
+refMsg = msprintf(_("Invalid index."));
+assert_checkerror("res=ode(""rk"",[x0;y0;v0;teta0;],t0,t,right)", refMsg);