Problem with lyap(Q,R)

General Usage
3

You may also want to try some alternative solvers. First, let me - just for convenience - rename your matrices so that they correspond to the common conventions used by solvers. Your A will become C, your B will become A after transposition:

C = [1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0]
A = [0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; -1252.0 1250.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -1.0 0.0 0.0 0.0; 1250.0 -1252.0 2.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -1.0 0.0 0.0; 0.0 2.0 -1252.0 1250.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -1.0 0.0; 0.0 0.0 1250.0 -1252.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -1.0; -1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1252.0 -1250.0 0.0 0.0; 0.0 -1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -1250.0 1252.0 -2.0 0.0; 0.0 0.0 -1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -2.0 1252.0 -1250.0; 0.0 0.0 0.0 -1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -1250.0 1252.0; 0.0 0.0 0.0 0.0 -1.0 0.0 0.0 0.0 -1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 -1.0 0.0 0.0 0.0 -1.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 -1.0 0.0 0.0 0.0 -1.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -1.0 0.0 0.0 0.0 -1.0 0.0 0.0 0.0 0.0]'

Now, I am aware of two alternatives for solving Lyapunov equation:

julia> using MatrixEquations

julia> X = lyapc(A,C)
16×16 Matrix{Float64}:
  15279.3          -15293.4           12471.6          …  -223.391        -119.321         115.15
 -15293.4           15299.2          -12479.8              227.77          114.793        -119.321
  12471.6          -12479.8           15299.2              114.793         227.77         -223.391
 -12455.4           12471.6          -15293.4             -119.321        -223.391         227.961
     -2.78533e-12       2.2311e-11       -2.34454e-8         1.21925e-11    -9.91488e-10     9.90487e-10
      1.06155e-11      -2.97575e-11       2.34695e-8   …     0.5             9.91008e-10    -9.90498e-10
      ⋮                                                ⋱                                     ⋮
      9.78416e-10      -9.78896e-10      -3.18567e-11        9.85301e-12     2.73115e-14     1.40998e-14
    227.961          -223.391          -119.321              7.10242        -6.05262         3.9942
   -223.391           227.77            114.793             -4.9365          4.00557        -6.05262
   -119.321           114.793           227.77               4.00557        -4.9365          7.10242
    115.15           -119.321          -223.391        …    -6.05262         7.10242        -4.92681

julia> A*X+X*A'+C
16×16 Matrix{Float64}:
  3.3527e-8    -3.49999e-8  -1.3459e-9   -9.35604e-10  …   2.41243e-9   -2.72769e-9    4.49009e-9   -2.63853e-9
 -3.49999e-8    3.64698e-8   1.17352e-9   1.10988e-9      -2.49065e-9    2.80954e-9   -4.42813e-9    2.57208e-9
 -1.3459e-9     1.17352e-9  -4.86263e-8   4.56729e-8       3.24606e-9   -3.6765e-9     4.3874e-9    -2.1422e-9
 -9.35604e-10   1.10988e-9   4.56729e-8  -4.27212e-8      -3.22284e-9    3.65187e-9   -4.42378e-9    2.18074e-9
 -2.49747e-9    2.48474e-9  -1.66074e-9   1.65346e-9      -7.63123e-11   8.03766e-11   3.14344e-11  -3.66214e-11
  5.33328e-9   -5.326e-9     4.77667e-9  -4.77121e-9   …   7.03722e-11  -7.77618e-11   3.45324e-12   5.57066e-12
  ⋮                                                    ⋱                                             ⋮
  1.88513e-9   -1.88631e-9   1.96059e-9  -1.95914e-9       2.40963e-11  -2.31539e-11  -1.49818e-11   1.37126e-11
  2.41243e-9   -2.49065e-9   3.24606e-9  -3.22284e-9      -7.75162e-11   8.55078e-11   1.25947e-10  -1.4731e-10
 -2.72769e-9    2.80954e-9  -3.6765e-9    3.65187e-9       8.55078e-11  -9.43823e-11  -1.29694e-10   1.51797e-10
  4.49009e-9   -4.42813e-9   4.3874e-9   -4.42378e-9       1.25947e-10  -1.29694e-10  -1.10189e-10   1.43357e-10
 -2.63853e-9    2.57208e-9  -2.1422e-9    2.18074e-9   …  -1.4731e-10    1.51797e-10   1.43357e-10  -1.77527e-10

and

julia> using ControlMatrixEquations

julia> X = lyapc(A,C)
16×16 Matrix{Float64}:  
15279.3          -15293.3           12471.6          …  -223.387        -119.323         115.147 -15293.3           15299.2          -12479.8              227.773         114.79         -119.323  12471.6          -12479.8           15299.2              114.79          227.773        -223.387 -12455.4           12471.6          -15293.3             -119.323        -223.387         227.964     -2.61548e-12       2.21302e-11      -2.34452e-8         1.21753e-11    -9.91495e-10     9.90495e-10      1.04541e-11      -2.9604e-11        2.34692e-8   …     0.5             9.91015e-10    -9.90506e-10      ⋮                                                ⋱                                     ⋮      9.78423e-10      -9.78905e-10      -3.18727e-11        9.853e-12       2.7605e-14      1.39602e-14
    227.964          -223.387          -119.323              7.10262        -6.0528          3.99402
   -223.387           227.773           114.79              -4.93629         4.00538        -6.0528
   -119.323           114.79            227.773              4.00538        -4.93629         7.10262
    115.147          -119.323          -223.387        …    -6.0528          7.10262        -4.92661

julia> A*X+X*A'+C16×16 Matrix{Float64}:  3.36389e-8   -3.51083e-8  -1.50165e-9  -8.36706e-10  …   2.40959e-9   -2.72439e-9    4.48801e-9   -2.64033e-9 -3.50854e-8    3.65519e-8   1.28263e-9   1.05764e-9      -2.48701e-9    2.80511e-9   -4.42472e-9    2.57256e-9 -1.57833e-9    1.30631e-9  -4.85399e-8   4.55455e-8       3.24075e-9   -3.67374e-9    4.38195e-9   -2.14322e-9 -7.61101e-10   1.03502e-9   4.56016e-8  -4.2609e-8       -3.21617e-9    3.64777e-9   -4.42071e-9    2.18392e-9 -2.49929e-9    2.49383e-9  -1.66438e-9   1.6571e-9       -7.63976e-11   8.0405e-11    3.1477e-11   -3.67208e-11  5.3351e-9    -5.32782e-9   4.7803e-9   -4.77485e-9   …   7.05427e-11  -7.78471e-11   3.31113e-12   5.65592e-12  ⋮                                                    ⋱                                             ⋮  1.88534e-9   -1.88663e-9   1.96076e-9  -1.95922e-9       2.41025e-11  -2.31735e-11  -1.49392e-11   1.36744e-11
  2.41596e-9   -2.49111e-9   3.24749e-9  -3.21498e-9      -7.74358e-11   8.5367e-11    1.25663e-10  -1.47177e-10
 -2.72917e-9    2.80772e-9  -3.67852e-9   3.64461e-9       8.53111e-11  -9.41252e-11  -1.29703e-10   1.51957e-10
  4.49071e-9   -4.42183e-9   4.38627e-9  -4.41719e-9       1.25676e-10  -1.2959e-10   -1.10448e-10   1.43313e-10
 -2.63916e-9    2.5658e-9   -2.14402e-9   2.17699e-9   …  -1.47204e-10   1.5186e-10    1.4371e-10   -1.77583e-10
2 Likes