[Maxima] Matrix Riccati code

[Donna and Dan Stanger]

This is a multi-part message in MIME format.
--------------070604040209060108090505
Content-Type: text/plain; charset=ISO-8859-1; format=flowed
Content-Transfer-Encoding: 7bit

I have written some code to solve the matrix riccati differential equation,
which is useful for control systems.
It is proof of concept, and has several bugs.
It is attached.  I would appreciate any comments.
The biggest issue is setting the solution constants,
which for control systems, the solution is zero at time T.
Dan Stanger


--------------070604040209060108090505
Content-Type: image/x-macpaint;
 name="odematriccati.mac"
Content-Transfer-Encoding: base64
Content-Disposition: inline;
 filename="odematriccati.mac"

LyogSSBoYXZlIHRyaWVkIHRvIG1ha2UgdGhpcyBjb2RlIHJ1biBvbiBib3RoIHRoZSBjb21t
ZXJjaWFsIG1hY3N5bWEsDQogICBhbmQgdGhlIEZBSUYgbWF4aW1hLiAgVGhpcyB2ZXJzaW9u
IG9mIHRoZSBjb2RlIHJlcGxhY2VzIHRoZSBpbml0aWFsDQogICBjb25kaXRpb24gdmFsdWUg
YW5kIGRlcml2YXRpdmUgdmFsdWUgYXQgdD0wIHdpdGggazAgYW5kIGsxIGJlY2F1c2UgaXQn
cw0KICAgZWFzaWVyLg0KICAgVG8gc29sdmUgdGhlIG9wdGltYWwgY29udHJvbCBwcm9ibGVt
LCB0aGUgc29sdXRpb24gaXMgMCBhdCB0PVQsDQogICBzbyB1c2Ugb2RlbWF0cmljY2F0aV9i
YyB3aXRoIHRoZSBzb2x1dGlvbiBmcm9tIG9kZW1hdHJpY2NhdGkuDQoqLyAgIA0KaWYgdmVy
c2lvbj0ndmVyc2lvbiB0aGVuICgNCiAgIGlmIGdldCgnb2RlbWF0c3lzLCAndmVyc2lvbikg
PSAnZmFsc2UgdGhlbiBsb2FkKCdvZGVtYXRzeXMpLA0KICAgaWYgZ2V0KCdkaWFnLCAndmVy
c2lvbikgPSAnZmFsc2UgdGhlbiBsb2FkKCdkaWFnKSwNCiAgIGlmIGdldCgnY2hvbGVza3lf
ZGVjb21wX3N5bWIsJ3ZlcnNpb24pID0gJ2ZhbHNlIHRoZW4gbG9hZCgnY2hvbGVza3lfZGVj
b21wX3N5bWIpLA0KICAgYWxpYXMob2RlbWF0c3lzQSwgb2RlbWF0c3lzKSkNCiAgIGVsc2Ug
KA0KICAgYWxpYXMob2RlbWF0c3lzQSwgb2RlbWF0c3lzX2ljKSwNCiAgIC8qIGJhdGNoKCJl
OlxcYXJyYXlfc3ltYm9scy5tYWMiKSwNCiAgIGFycmF5X3N5bWJvbHNfaW5pdCgpLCovDQog
ICBvZGVtYXRzeXNfaWMoYSxiLFtjXSk6PWJsb2NrKFtMOmFycmF5X3N5bWJvbHNfdW5hcnJh
eShbYSxiXSksdXNlZCxlXSwNCiAgICAgIGU6Zmlyc3QoTCksDQogICAgICB1c2VkOnNlY29u
ZChMKSwNCiAgICAgIGE6Zmlyc3QoZSksDQogICAgICBiOnNlY29uZChlKSwNCiAgICAgIGlm
IGxlbmd0aChjKSA+IDAgdGhlbiBhcnJheV9zeW1ib2xzX2F0dmFsdWUxKGMsdXNlZCx0LDAp
LA0KICAgICAgb2RlbWF0c3lzKGEsYikpKTsNCg0Kb2RlbWF0cmljY2F0aShRLEIsSix2YXJz
KTo9YmxvY2soW0gsQSxHLFMsTCxYLFhzb2wsVyxZLA0KICAgICAgaW5kZXB2YXIsbnJvd3Ms
bmNvbHNdLA0KICAgaW5kZXB2YXI6Zmlyc3QobGlzdG9mdmFycyh2YXJzKSksDQogICBIOmNo
b2xlc2t5X2RlY29tcF9zeW1iKFEpLA0KICAgRzp0cmFuc3Bvc2UoSCkuSi5ILA0KICAgQTpp
bnZlcnQoSCkuKEIuSC1kaWZmKEgsaW5kZXB2YXIpKSwNCiAgIFM6dHJhbnNwb3NlKEEpLUEs
DQogICBwcmludCgiSiIsSiwiSCIsSCwiQSIsQSwgIlMiLFMsIkciLEcpLA0KICAgbnJvd3M6
bGVuZ3RoKHZhcnMpLA0KICAgbmNvbHM6bGVuZ3RoKHZhcnNbMV0pLA0KICAgTDpibG9jayhb
YXV4ZXFuLGF1eHZhcnNdLA0KICAgICAgYXV4dmFyczpnZW5tYXRyaXgobGFtYmRhKFtpLGtd
LExbaSxrXShpbmRlcHZhcikpLG5yb3dzLG5jb2xzKSwNCiAgICAgIG9kZW1hdHN5c0EoDQog
ICAgICAgICBnZW5tYXRyaXgobGFtYmRhKFtpLGtdLCdkaWZmKExbaSxrXShpbmRlcHZhciks
aW5kZXB2YXIpKSxucm93cyxuY29scykNCiAgICAgICAgICAgIC0oYXV4dmFycy5TLzIpLA0K
ICAgICAgICAgYXV4dmFycywNCiAgICAgICAgIGlkZW50KG5yb3dzKSkpLA0KICAgTDpyaHMo
Zmlyc3QoTCkpLA0KICAgcHJpbnQoImNvbnN0YW50IiwoKEwpLihHK3RyYW5zcG9zZShBKS5B
KyhTLlMvNCkNCiAgICAgICAgICAgICsoKHRyYW5zcG9zZShkaWZmKEEsaW5kZXB2YXIpKStk
aWZmKEEsaW5kZXB2YXIpKS8yKSkudHJhbnNwb3NlKEwpKSksDQogICBYc29sOmJsb2NrKFtq
YWNlcW4samFjdmFycyxtZXFuXSwNCiAgICAgIGphY3ZhcnM6Z2VubWF0cml4KGxhbWJkYShb
aSxrXSxYW2ksa10oaW5kZXB2YXIpKSxucm93cyxuY29scyksDQogICAgICBtZXFuOmdlbm1h
dHJpeChsYW1iZGEoW2ksa10sJ2RpZmYoWFtpLGtdKGluZGVwdmFyKSxpbmRlcHZhciwyKSks
DQoJCSAgICAgbnJvd3MsbmNvbHMpDQogICAgICAgICAgICAtKChqYWN2YXJzLkwpLihHK3Ry
YW5zcG9zZShBKS5BKyhTLlMvNCkNCiAgICAgICAgICAgICsoKHRyYW5zcG9zZShkaWZmKEEs
aW5kZXB2YXIpKStkaWZmKEEsaW5kZXB2YXIpKS8yKSkudHJhbnNwb3NlKEwpKSwNCiAgIG9k
ZW1hdHN5c0EobWVxbixqYWN2YXJzLA0KICAgICAgZ2VubWF0cml4KGxhbWJkYShbaSxrXSxL
MFtpLGtdKSxucm93cyxuY29scyksDQogICAgICBnZW5tYXRyaXgobGFtYmRhKFtpLGtdLEsx
W2ksa10pLG5yb3dzLG5jb2xzKSkpLA0KICAgWHNvbDpyaHMoZmlyc3QoWHNvbCkpLA0KICAg
VzoodHJhbnNwb3NlKEwpLmludmVydChYc29sKS5kaWZmKFhzb2wsaW5kZXB2YXIpLkwpLSgo
QSt0cmFuc3Bvc2UoQSkpLzIpLA0KICAgWTppbnZlcnQodHJhbnNwb3NlKEgpKS5XLkguaW52
ZXJ0KEgpDQopOw0KLyogZGlzcGxheTJkOmZhbHNlOyAqLw0KdHJhY2Uob2RlbWF0c3lzKTsN
Clk6b2RlbWF0cmljY2F0aShtYXRyaXgoWzEsMF0sWzAsMV0pLA0KICAgbWF0cml4KFsxLDBd
LFswLDFdKSxpZGVudCgyKSxnZW5tYXRyaXgobGFtYmRhKFtpLGpdLFlbaSxqXSh0KSksMiwy
KSk7DQpzb2x2ZShzdWJzdChbazBbMSwxXT0xLGswWzEsMl09MSxrMFsyLDFdPTEsazBbMiwy
XT0xXSxhcHBseShhcHBlbmQsbWFwKGFwcGVuZCxhcmdzKFkpKSkpLFsgazFbMSwxXSxrMVsy
LDJdLGsxWzEsMl0sazFbMiwxXV0pOw0KLyoNCg0KICAgaWYgbGVuZ3RoKGluaXRjb25kcyk9
MSB0aGVuIGJsb2NrKFtpbml0cG9pbnQsWGkseGRvdF0sDQogICAgICBpbml0cG9pbnQ6aW5p
dGNvbmRzWzFdLA0KICAgICAgWGk6YXQoZGlmZihYc29sLGluZGVwdmFyKSxpbmRlcHZhcj1p
bml0cG9pbnQpLQ0KICAgICAgICAgaW52ZXJ0KHRyYW5zcG9zZShMKSkuKChBK3RyYW5zcG9z
ZShBKSkvMikuaW52ZXJ0KEwpLA0KICAgICAgeGRvdDpnZW5tYXRyaXgoDQogICAgICAgICBs
YW1iZGEoW2ksa10sYXQoJ2RpZmYoWFtpLGtdKGluZGVwdmFyKSxpbmRlcHZhciksaW5kZXB2
YXI9MCk9MCksDQogICAgICAgICBucm93cyxuY29scyksDQogICAgICBYc29sOmxyYXRzdWJz
dCggYXBwbHkoIGFwcGVuZCxhcmdzKHhkb3QpKSwgWHNvbCksDQogICAgICBYaTpscmF0c3Vi
c3QoIGFwcGx5KCBhcHBlbmQsYXJncyh4ZG90KSksIFhpKSwNCiAgICAgIFhpOmFsZ3N5cyhh
cHBseShhcHBlbmQsYXJncyhYaSkpLA0KICAgICAgICAgYXBwbHkoYXBwZW5kLGFyZ3MoDQog
ICAgICAgICAgICBnZW5tYXRyaXgobGFtYmRhKFtpLGtdLGF0KFhbaSxrXShpbmRlcHZhciks
aW5kZXB2YXI9MCkpLA0KICAgICAgICAgICAgICAgbnJvd3MsbmNvbHMpKSkpLA0KICAgICAg
WHNvbDpscmF0c3Vic3QoZmlyc3QoWGkpLFhzb2wpLA0KICAgICAgcHJpbnQoIlhzb2wgIixY
c29sKQ0KICAgKSwNCiovDQo=
--------------070604040209060108090505
Content-Type: image/x-macpaint;
 name="array_symbols.mac"
Content-Transfer-Encoding: base64
Content-Disposition: inline;
 filename="array_symbols.mac"

YXJyYXlfc3ltYm9sc19pbml0KCk6PWFycmF5X3N5bWJvbHNfdW51c2VkOm1ha2VsaXN0KGV2
YWwoP2dlbnN5bSgpKSxpLDAsMTAwKSQNCg0KYXJyYXlfc3ltYm9sc19nZXRfZ2Vuc3ltKCk6
PWJsb2NrKFtnXSwNCiAgIGlmIGxlbmd0aChhcnJheV9zeW1ib2xzX3VudXNlZCkgPiAwIHRo
ZW4NCiAgICAgIGc6cG9wKGFycmF5X3N5bWJvbHNfdW51c2VkKQ0KICAgZWxzZSBlcnJvcigi
YXJyYXlfc3ltYm9sc191bmFycmF5OiBydW4gb3V0IG9mIGdlbnN5bXMsIG1vZGlmeSBhcnJh
eV9zeW1ib2xzX2luaXQiKSwNCiAgIGcpJA0KDQphcnJheV9zeW1ib2xzX3VuYXJyYXkoZXhw
KTo9YmxvY2soW3N2OltdLHM6W10sdXNlZDpbXV0sDQogICBzY2FubWFwKGxhbWJkYShbeF0s
DQogICAgICBpZiBub3QgbWFwYXRvbSh4KSBhbmQgc3VidmFycChvcCh4KSkgdGhlbiBibG9j
ayhbbzpvcCh4KV0sDQogICAgICAgICBpZiBub3QgbWVtYmVyKG8sc3YpIHRoZW4gcHVzaChv
LHN2KSksDQogICAgICB4KSwNCiAgICAgIGV4cCksDQogICBmb3IgYSBpbiBzdiBkbyBibG9j
ayhbZ10sDQogICAgICAgZzphcnJheV9zeW1ib2xzX2dldF9nZW5zeW0oKSwNCiAgICAgICBw
dXNoKGc9YSx1c2VkKSwNCiAgICAgICBwdXNoKGE9ZyxzKSksDQogICBbc3Vic3QocyxleHAp
LHVzZWRdDQopJA0KDQphcnJheV9zeW1ib2xzX3JlYXJyYXkoZXhwLHVzZWQpOj1ibG9jaygN
CiAgIHN1YnN0KHVzZWQsZXhwKSwNCiAgIGZvciBlIGluIHVzZWQgZG8gcHVzaChyaHMoZSks
YXJyYXlfc3ltYm9sc191bnVzZWQpDQopJA0KDQphcnJheV9zeW1ib2xzX2F0dmFsdWUoZXhw
LHVzZWQsdmFyLHBvaW50KTo9YmxvY2soDQogICBmb3IgdSBpbiB1c2VkIGRvIGJsb2NrKFts
Omxocyh1KSxyOnJocyh1KV0sDQogICAgICBhdHZhbHVlKGZ1bm1ha2UobCxbdmFyXSksW3Zh
cj1wb2ludF0sDQogICAgICAgIGFwcGx5KGxhbWJkYShbdSx2XSxleHBbdSx2XSksYXJncyhy
KSkpKQ0KKSQNCg0KYXJyYXlfc3ltYm9sc19hdHZhbHVlMShleHAxLHVzZWQsdmFyLHBvaW50
KTo9YmxvY2soW2o6MF0sDQogIGZvciBleHAgaW4gZXhwMSBkbyBibG9jaygNCiAgIGZvciB1
IGluIHVzZWQgZG8gYmxvY2soW2w6bGhzKHUpLHI6cmhzKHUpXSwNCiAgICAgIGF0dmFsdWUo
ZGlmZihmdW5tYWtlKGwsW3Zhcl0pLHZhcixqKSxbdmFyPXBvaW50XSwNCiAgICAgICAgYXBw
bHkobGFtYmRhKFt1LHZdLGV4cFt1LHZdKSxhcmdzKHIpKSkpLA0KICAgajpqKzEpDQopJA0K
YXJyYXlfc3ltYm9sc19pbml0KCk7DQo=
--------------070604040209060108090505--