issue #18

Enhancement: ode45reg accepts reverse time tspan argument now
SystemAnalysisDpt-CMC-MSU · Jan 16, 2016 · a0f5c29 · a0f5c29
1 parent 65a7f44
commit a0f5c29
Show file tree

Hide file tree

Showing 2 changed files with 74 additions and 6 deletions.
diff --git a/products/+gras/+ode/+test/+mlunit/SuiteOde45Reg.m b/products/+gras/+ode/+test/+mlunit/SuiteOde45Reg.m
@@ -41,10 +41,10 @@ function compare(yMat,yyMat,yRegMat,yyRegMat)
                     [~,yyMat,yyRegMat] = interpObj.evaluate(tspanVec);
                 end
 
-                 tBeginVec = linspace(0,1,nPoints);
-                 tBeginVec(2) = 1e-6;
-                 interpObj = helpForComparision(tBeginVec,tBeginVec);
-                 compare(yMat,yyMat,yRegMat,yyRegMat);
+                tBeginVec = linspace(0,1,nPoints);
+                tBeginVec(2) = 1e-6;
+                interpObj = helpForComparision(tBeginVec,tBeginVec);
+                compare(yMat,yyMat,yRegMat,yyRegMat);
 
                 tVaryVec = tBeginVec.^2;
                 tVaryVec(2) = tBeginVec(2);
@@ -86,5 +86,59 @@ function compare(yMat,yyMat,yRegMat,yyRegMat)
 
             end
         end
+
+        function self = testInverseTspan(self)
+            %TODO: common compare and fRegDummy function for both test
+
+            function [isStrictViolVec,yRegMat] = fRegDummy(~,yMat)
+            	isStrictViolVec=false(1,size(yMat,2));
+                yRegMat=yMat;
+            end
+
+            function compare(yMat,yyMat,yRegMat,yyRegMat, t1, t2, COMPARE_TOL)
+                if ~exist('COMPARE_TOL', 'var')
+                    COMPARE_TOL = 1e-8;
+                end
+
+                [isEqual,~,~,~,~] = modgen.common.absrelcompare(...
+                    yMat,yyMat,COMPARE_TOL,COMPARE_TOL,@norm);
+                mlunitext.assert_equals(true,isEqual,...
+                    'matrix yMat and yyMat are not equal')
+                [isEqual,~,~,~,~] = modgen.common.absrelcompare(...
+                    yRegMat,yyRegMat,COMPARE_TOL,COMPARE_TOL,@norm);
+                mlunitext.assert_equals(true,isEqual,...
+                    'matrix yRegMat and yyRegMat are not equal');
+                [isEqual,~,~,~,~] = modgen.common.absrelcompare(...
+                    t1,t2,COMPARE_TOL,COMPARE_TOL,@norm);
+                mlunitext.assert_equals(true,isEqual,...
+                    'time vectors t1 and t2 are not equal');
+            end
+
+            ABS_TOL=1e-8;
+            odePropList={'NormControl','on','RelTol',ABS_TOL,...
+                    'AbsTol',ABS_TOL};
+
+            function check(fOdeDeriv1, tspan1, y0Vec)
+                [t1, y1, yr1] = feval(self.odeSolver, fOdeDeriv1,...
+                                     @fRegDummy, tspan1, y0Vec,...
+                                     odePropList{:});
+
+                t0 = tspan1(end);
+                fOdeDeriv2 = @(t, y) -fOdeDeriv1(t0 - t, y);
+                tspan2 = t0 - tspan1;
+                [t2, y2, yr2] = feval(self.odeSolver, fOdeDeriv2,...
+                                     @fRegDummy, tspan2, y0Vec,...
+                                     odePropList{:});
+                t2 = t0 - t2;
+
+                compare(y1, y2, yr1, yr2, t1, t2);
+            end
+
+            check(@(t,y) cos(y), [0 1], 0);
+            check(@(t,y) sin(y), linspace(0, 1, 9), 0);
+
+            check(@(t, y) t*cos(y), [0, 1], zeros(1, 10));
+            check(@(t, y) t*sin(y), linspace(0, 1, 9), ones(1, 10));
+        end
     end
 end
diff --git a/products/+gras/+ode/ode45reg.m b/products/+gras/+ode/ode45reg.m
@@ -41,8 +41,9 @@
 %   interpObj: gras.ode.VecOde45RegInterp[1,1] - all the data nessecary 
 %       for calculation on an arbitrary time grid is stored in  this object
 %
-% $Author: Peter Gagarinov  <pgagarinov@gmail.com> $	$Date: 2011$
-% $Author: Vadim Danilov <vadimdanilov93@gmail.com> $   $Data: 2013$
+% $Author: Peter Gagarinov  <pgagarinov@gmail.com> $	      $Date: 2011$
+% $Author: Vadim Danilov <vadimdanilov93@gmail.com> $         $Data: 2013$
+% $Author: Arthur Lukyanenko <arthurlukyanenko@gmail.com> $   $Data: 2016$
 % $Copyright: Moscow State University,
 %            Faculty of Computational Mathematics and Computer Science,
 %            System Analysis Department 2011 $
@@ -77,6 +78,13 @@
 checkgen(fOdeDeriv,'isfunction(x)');
 checkgen(fOdeReg,'isfunction(x)');
 %
+%% Handle inverse time direction
+isInvTspan = (tspan(1) > tspan(end));
+if isInvTspan
+    fOdeDeriv = @(t, y) -fOdeDeriv(-t, y);
+    fOdeReg = @(t, y) fOdeReg(-t, y);
+    tspan = -tspan;
+end
 %% Handle solver arguments
 [neq, tspan, ntspan, next, t0, tfinal, y0Vec, f0, ...
     options, threshold, rtol, normcontrol, normy, hmax, htry, htspan,...
@@ -391,10 +399,16 @@
 end
 shrinkResults();
 if nargout == 4
+    if isInvTspan
+        throwerror('notImplemented', 'Reverse time is not supported');
+    end
     SInterp.tVec = [SInterp.tBegin SInterp.tVec];
     SInterp.yCVec = [{SInterp.yBeginVec} SInterp.yCVec];
     interpObj = gras.ode.VecOde45RegInterp(SInterp);
 end;
+if isInvTspan
+	tout = -tout;
+end
 prDispObj.finish();
     function shrinkResults()
         tout = tout(1:nout).';