[Octave-bug-tracker] [bug #66567] control package lyap function returns incorrect output

[Andrea]

URL:
  <https://savannah.gnu.org/bugs/?66567>

                 Summary: control package lyap function returns incorrect
output
                   Group: GNU Octave
               Submitter: andrea993
               Submitted: Mon 16 Dec 2024 08:22:57 PM UTC
                Category: Octave Function
                Severity: 3 - Normal
                Priority: 5 - Normal
              Item Group: Incorrect Result
                  Status: None
             Assigned to: None
         Originator Name:
        Originator Email:
             Open/Closed: Open
         Discussion Lock: Any
                 Release: stable
        Operating System: Any
           Fixed Release: None
         Planned Release: None


    _______________________________________________________

Follow-up Comments:


-------------------------------------------------------
Date: Mon 16 Dec 2024 08:22:57 PM UTC By: Andrea <andrea993>
This is a sample code involving the function lyap() from the control package.
The same following code has been executed in Octave and in MATLAB and the
Octave result is wrong because the norm(A*X + X*A' + Q) should be zero:

A = [
-0.016, 0.4, -0.705, -0.717,
0.386, -0.147, 0.933, -0.693,
0.643, -0.617, 0.634, -0.689,
0.464, -0.411, 0.364, 0.444];

Q = [
-0.016, 0.1, -0.705, -0.717,
0.3, -0.147, 0.933, -0.693,
0.643, -0.67, 0.634, -0.689,
0.464, -0.411, 0.4, 0.444];


X1 = lyap(A,Q)
X2 = lyapunov2(A,Q)

norm(A*X1 + X1*A' + Q)
norm(A*X2 + X2*A' + Q)


% Octave
X1 =

  -31.9640   23.0586    0.1309   13.4373
   23.0586  -34.0085  -31.7876  -22.8448
    0.1309  -31.7876  -57.3029  -23.6806
   13.4373  -22.8448  -23.6806  -16.2757

X2 =

    5.6744  -53.5457  -87.7349  -38.8693
   45.2712    6.0224   68.5902   10.6251
   87.8241  -56.7541   10.0535  -31.3875
   33.8899   -2.0659   40.2937    2.4127

ans = 2.2311
ans = 1.3398e-13

% Matlab

X1 =

    5.6744  -53.5457  -87.7349  -38.8693
   45.2712    6.0224   68.5902   10.6251
   87.8241  -56.7541   10.0535  -31.3875
   33.8899   -2.0659   40.2937    2.4127


X2 =

    5.6744  -53.5457  -87.7349  -38.8693
   45.2712    6.0224   68.5902   10.6251
   87.8241  -56.7541   10.0535  -31.3875
   33.8899   -2.0659   40.2937    2.4127


ans =

   1.2669e-13


ans =

   1.1840e-13



The function lyapunov2() that works in both programs is the following:

function X = lyapunov2(A,Q)
m = size(Q,1);
[U,A1] = schur(A,'complex');
Q1 = -U'*Q*U;
X = zeros(m);
X(:,m) = (A1 + conj(A1(m,m))*eye(m))\Q1(:,m);
for i = m-1:-1:1
    v = Q1(:,i) - X(:,i+1:m)*A1(i,i+1:m)';
    X(:,i) = (A1 + conj(A1(i,i))*eye(m))\v;
end
X = real(U*X*U');









    _______________________________________________________

Reply to this item at:

  <https://savannah.gnu.org/bugs/?66567>

_______________________________________________
Message sent via Savannah
https://savannah.gnu.org/

signature.asc
Description: PGP signature