|
53
|
1 function x = lyap (a, b, c) |
|
40
|
2 |
|
73
|
3 # Usage: x = lyap (a, b {,c}) |
|
53
|
4 # |
|
|
5 # If (a, b, c) are specified, then lyap returns the solution of the |
|
|
6 # Sylvester equation |
|
|
7 # |
|
|
8 # a x + x b + c = 0 |
|
|
9 # |
|
|
10 # If only (a, b) are specified, then lyap returns the solution of the |
|
|
11 # Lyapunov equation |
|
|
12 # |
|
|
13 # a' x + x a + b = 0 |
|
|
14 # |
|
|
15 # If b is not square, then lyap returns the solution of either |
|
|
16 # |
|
|
17 # a' x + x a + b' b = 0 |
|
|
18 # |
|
|
19 # or |
|
|
20 # |
|
|
21 # a x + x a' + b b' = 0 |
|
|
22 # |
|
|
23 # whichever is appropriate. |
|
73
|
24 # |
|
|
25 # Solves by using the Bartels-Stewart algorithm (1972). |
|
|
26 |
|
|
27 # Written by A. S. Hodel (scotte@eng.auburn.edu) August 1993. |
|
|
28 |
|
40
|
29 |
|
53
|
30 if (nargin != 3 && nargin != 2) |
|
57
|
31 error ("usage: lyap (a, b {,c})"); |
|
53
|
32 endif |
|
|
33 |
|
|
34 if ((n = is_square(a)) == 0) |
|
|
35 error ("lyap: a is not square"); |
|
|
36 endif |
|
|
37 |
|
|
38 if (nargin == 2) |
|
40
|
39 |
|
53
|
40 # Transform Lyapunov equation to Sylvester equation form. |
|
|
41 |
|
|
42 if ((m = is_square (b)) == 0) |
|
|
43 if ((m = rows (b)) == n) |
|
|
44 |
|
|
45 # solve a x + x a' + b b' = 0 |
|
|
46 |
|
|
47 b = b * b'; |
|
|
48 a = a'; |
|
|
49 else |
|
|
50 |
|
|
51 # Try to solve a'x + x a + b' b = 0. |
|
|
52 |
|
|
53 m = columns (b); |
|
|
54 b = b' * b; |
|
|
55 endif |
|
40
|
56 |
|
53
|
57 if (m != n) |
|
|
58 error ("lyap: a, b not conformably dimensioned"); |
|
|
59 endif |
|
40
|
60 endif |
|
53
|
61 |
|
|
62 # Set up Sylvester equation. |
|
|
63 |
|
|
64 c = b; |
|
|
65 b = a; |
|
|
66 a = b' |
|
|
67 |
|
|
68 else |
|
|
69 |
|
|
70 # Check dimensions. |
|
|
71 |
|
|
72 if ((m = is_square (b)) == 0) |
|
|
73 error ("lyap: b must be square in a sylvester equation"); |
|
40
|
74 endif |
|
53
|
75 |
|
73
|
76 [n1, m1] = size(c); |
|
53
|
77 |
|
|
78 if (n != n1 || m != m1) |
|
|
79 error("lyap: a,b,c not conformably dimensioned"); |
|
|
80 endif; |
|
40
|
81 endif |
|
|
82 |
|
53
|
83 # Call octave built-in function. |
|
40
|
84 |
|
73
|
85 x = syl (a, b, c); |
|
40
|
86 |
|
|
87 endfunction |
