Add tests for ricatti() function, thx to Sabine Gauzere

2 files changed, 107 insertions, 0 deletions

diff --git a/tests/basic_tests/slicot.dia.ref b/tests/basic_tests/slicot.dia.ref
index 03d971e..58038b5 100644
--- a/tests/basic_tests/slicot.dia.ref
+++ b/tests/basic_tests/slicot.dia.ref
@@ -1694,3 +1694,72 @@ X = linmeq(3,A,C,flag);
 if norm(A'*(X'*X)*A-X'*X+C'*C,1)>100000*%eps then bugmes();quit;end
  
  
+ 
+// ============================================================================
+ 
+// Tests for riccati() function
+ 
+// Added : 25/06/2007
+ 
+// Thanks to Sabine GAUZERE
+ 
+// ============================================================================
+ 
+ 
+n = 10;
+ 
+A = rand(n,n);
+ 
+B = rand(n,n);
+ 
+C = rand(n,n);
+ 
+C = C*C';
+ 
+R = rand(n,n);
+ 
+R = R*R'+eye();
+ 
+B = B*inv(R)*B';
+ 
+ 
+// Test de l'équation en temps continu
+ 
+X = riccati(A,B,C,'c','eigen');
+ 
+ 
+// Vérification
+ 
+C_cont = A'*X+X*A-X*B*X;
+ 
+ 
+// Test de l'équation en temps discret
+ 
+F      = A;
+ 
+B      = rand(n,n);
+ 
+G1     = B;
+ 
+G2     = R;
+ 
+G      = G1/G2*G1';
+ 
+H      = C;
+ 
+[X1,X2]= riccati(F,G,H,'d','schur');
+ 
+ 
+// Vérification
+ 
+X      = X1/X2;
+ 
+C_disc = F'*X*F-(F'*X*G1/(G2+G1'*X*G1))*(G1'*X*F)-X;
+ 
+ 
+// Comparaison des différents résultats obtenus
+ 
+if norm(C_cont+C,1) > (10000*%eps) then bugmes();quit;end
+ 
+if norm(C_disc+H,1) > (10000*%eps) then bugmes();quit;end
+ 
diff --git a/tests/basic_tests/slicot.tst b/tests/basic_tests/slicot.tst
index 0d3f529..dd3d9ed 100644
--- a/tests/basic_tests/slicot.tst
+++ b/tests/basic_tests/slicot.tst
@@ -885,3 +885,41 @@ flag=[1,2];
 X = linmeq(3,A,C,flag);
 if norm(A'*(X'*X)*A-X'*X+C'*C,1)>100000*%eps then pause,end
 
+
+// ============================================================================
+// Tests for riccati() function
+// Added : 25/06/2007
+// Thanks to Sabine GAUZERE
+// ============================================================================
+
+n = 10;
+A = rand(n,n);
+B = rand(n,n);
+C = rand(n,n);
+C = C*C';
+R = rand(n,n);
+R = R*R'+eye();
+B = B*inv(R)*B';
+
+// Test de l'équation en temps continu
+X = riccati(A,B,C,'c','eigen');
+
+// Vérification
+C_cont = A'*X+X*A-X*B*X;
+
+// Test de l'équation en temps discret
+F      = A;
+B      = rand(n,n);
+G1     = B;
+G2     = R;
+G      = G1/G2*G1';
+H      = C;
+[X1,X2]= riccati(F,G,H,'d','schur');
+
+// Vérification
+X      = X1/X2;
+C_disc = F'*X*F-(F'*X*G1/(G2+G1'*X*G1))*(G1'*X*F)-X;
+
+// Comparaison des différents résultats obtenus
+if norm(C_cont+C,1) > (10000*%eps) then pause,end
+if norm(C_disc+H,1) > (10000*%eps) then pause,end