1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
|
// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
// Copyright (C) ????-2008 - INRIA
//
// This file must be used under the terms of the CeCILL.
// This source file is licensed as described in the file COPYING, which
// you should have received as part of this distribution. The terms
// are also available at
// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
function [f,d]=invr(h,flag)
//if h is a scalar, polynomial or rational fonction matrix, invr
//computes h^(-1).
//!
if argn(2)==1 then
flag="C";
end
lhs=argn(1)
select typeof(h)
case "constant" then
f=inv(h);
case "polynomial" then //POLYNOMIAL MATRIX
[m,n]=size(h);
if m<>n then error(20),end
ndeg=max(degree(h));
if ndeg==1 then //MATRIX PENCIL
E=coeff(h,1);A=-coeff(h,0);
if norm(E-eye(E),1) < 100*%eps then
// sI -A
[num,den]=coff(A,varn(h));f=num/den;
else
[Bfs,Bis,chis]=glever(E,A,varn(h));
f=Bfs/chis - Bis;
if lhs==2 then
d=lcm(f("den"));
f=f*d;f=f("num");
end
end
else // GENERAL POLYNOMIAL MATRIX
select flag
case "L"
f=eye(n,n);
for k=1:n-1,
b=h*f,
d=-sum(diag(b))/k
f=b+eye(n,n)*d,
end;
d=sum(diag(h*f))/n,
if degree(d)==0 then d=coeff(d),end,
if lhs==1 then f=f/d;end
case "C"
[f,d]=coffg(h);
if degree(d)==0 then d=coeff(d),end
if lhs==1 then f=f/d;end
else
error(msprintf(gettext("%s: Wrong value for input argument #%d: Must be in the set {%s}.\n"),..
"invr",2,"''C'',''D''"))
end;
end
case "rational" then
[m,n]=size(h(2));
if m<>n then error(20),end
select flag
case "L" // Leverrier
f=eye(n,n);
for k=1:n-1,
b=h*f,
d=0;for l=1:n,d=d+b(l,l),end,d=-d/k;
f=b+eye(n,n)*d,
end;
b=h*f;d=0;for l=1:n,d=d+b(l,l),end;d=d/n,
if lhs==1 then f=f/d;end
case "A" // lcm of all denominator entries
denh=lcm(h("den"));
Num=h*denh;Num=Num("num");
[N,d]=coffg(Num);
f=N*denh;
if lhs==1 then f=f/d;end
case "C"// default method by polynomial inverse
[Nh,Dh]=lcmdiag(h); //h=Nh*inv(Dh); Dh diagonal;
[N,d]=coffg(Nh);
f=Dh*N;
if lhs==1 then f=f/d;end
case "Cof"// cofactors method
[f,d]=coffg(h);
if lhs==1 then f=f/d;end
else
error(msprintf(gettext("%s: Wrong value for input argument #%d: Must be in the set {%s}.\n"),..
"invr",2,"''L'',''A'',''C'',''Cof''"))
end;
else
error(msprintf(gettext("%s: Wrong type for input argument #%d: A floating point number or polynomial or rational fraction array expected.\n"),"invr",1))
end;
endfunction
|
