riccati

Arguments

A, B, C

Square matrices of real numbers, of size n x n: Matricial coefficients of the equation.

H, E

Square matrices of real numbers, of size 2n x 2n: Hamiltonian matrices | pencil of the equation.

dom

Type / time domain of the Riccati equation: "c" or "continuous", or "d" or "discrete".

method

string : "eigen" for block diagonalization (default), or "schur" for Schur method.

X1, X2, X

square matrices of real numbers (X2 invertible). X is symetric when B and C are so.

residual

real number: norm₁ of the actual Equation(X) residue.

Description

Output options

[X1, X2] = riccati(..) provides X1 and X2, with X2 invertible, such that X = X1 / X2.

residual is the L1-norm of the actual equation's result. If X is an actual solution, residual should be 0. Most often its value is close to %eps*norm(X).

In the discrete case, sometimes B1 and B2 can't be retrieved from B. Then residual can't be assessed and is set to %nan.

When no solution is found, X=[] or X1=[] is returned, and a warning is yielded.

Examples

n = 3;
// [A, B, C] = (grand(n,n,"uin",-9,9), grand(n,n,"uin",-9,9), grand(n,n,"uin",-9,9))
A = [
  -62.   91.   57.
  -43.  -45.  -19.
   58.   83.  -62.
   ];
B = [
   75.  -31.  -10.
  -79.   70.   68.
  -72.  -5.    32.
   ];
C = [
  -56.   70.   58.
  -41.   54.   50.
   90.   2.   -40.
   ];

// With A, B, C
// ------------
X = riccati(A, B, C, "c")
clean(A'*X + X*A - X*B*X + C)
[x1, x2, res] = riccati(A,B,C, "c");
x = x1 / x2;
and(x==X)
res

// With the Hamiltonian
// --------------------
H = [A, -B; -C, -A'];
X = riccati(H)
clean(A'*X + X*A - X*B*X + C)
[x1, x2] = riccati(H);
x = x1 / x2;
and(x==X)

--> X = riccati(A, B, C, "c")
 X  =
  -0.1790367   0.4166084   0.2319152
  -0.4977093   0.7993495   0.3086213
   0.5595286   0.3202094  -0.103394

--> clean(A'*X + X*A - X*B*X + C)
 ans  =
   0.   0.   0.
   0.   0.   0.
   0.   0.   0.

--> [x1, x2, res] = riccati(A,B,C, "c");
--> x = x1 / x2;
--> and(x==X)
 ans  =
  T

--> res
 res  =
   3.340D-13


--> // With the Hamiltonian
--> // --------------------
--> H = [A, -B; -C, -A'];
--> X = riccati(H);
 X  =
  -0.1790367   0.4166084   0.2319152
  -0.4977093   0.7993495   0.3086213
   0.5595286   0.3202094  -0.103394

--> clean(A'*X + X*A - X*B*X + C)
 ans  =
   0.   0.   0.
   0.   0.   0.
   0.   0.   0.

--> [x1, x2] = riccati(H);
--> x = x1 / x2;
--> and(x==X)
 ans  =
  T

// Eq: A'*X*A - (A'*X*B1) / (B2 + B1'*X*B1) * (B1'*X*A) + C - X = 0
n = 4;
[A, B1, B2, C] = (rand(n,n), rand(n,n), rand(n,n), rand(n,n));
B  = B1 / B2 * B1';

X = riccati(A, B, C, 'd', 'schur')
if X <> [] then
    clean(A'*X*A - (A'*X*B1) / (B2 + B1'*X*B1) * (B1'*X*A) + C - X)
    [x1, x2, res] = riccati(A, B, C, 'd', 'schur');
    x = x1 / x2;
    and(x==X)
    res
else
    disp("Riccati: No solution found")
end

See also

• ricc — Riccati equation
• srfaur — square-root algorithm
• lindquist — Lindquist's algorithm

History