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
|
// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
// Copyright (C) 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 [flag]=%choose(x)
// Utility function for use with schur
// [U,dim]=schur(A,choose) returns orth. basis
// for eigenspace associated to selected polynomials
//Needs two global variables :
// %sel = list of selected polynomials (user defined)
// eps = threshold for polynomials selection (eps= 0.0001 as default value)
// see below
// Example:
// A=...
// chis=poly(A,'s'); //Characteristic polynomial
// w=factors(chis); //Factors of chis in a list
// %sel=list(w(2),w(3)); // two selected polynomials
// eps=0.01; //Threshold (see almosteq below)
// [U,dim]=schur(A,%choose); //Ordered Schur form
// U1=U(:,1:dim);chi1=poly(U1'*A*U1,'s') //Check
// w1=factors(chi1) // w1 = %sel ?
//
// Copyright INRIA
eps=0.0001; //modify eps here !
//
flag=0;s=poly(0,"s");
select x(1)
case 1 then
// ASSUME x(3) NOT ZERO (for gev pb. x(3)=0 => eval @ infty)
vp=x(2)/x(3);pol=s-vp; //disp(pol);
for p=%sel; if almosteq(pol,p,eps) then flag=1;end;end
case 2 then
pol=s^2-x(2)*s+x(3); //disp(pol);
for p=%sel; if almosteq(pol,p,eps) then flag=1;end;end
end
endfunction
function trfa=almosteq(pol,p,eps)
// returns %T if pol ~ p %F if not
if degree(pol)<>degree(p) then trfa=%F;return;end
if norm((coeff(p)-coeff(pol)),1)<=eps then trfa=%T;return;end
trfa=%F;
endfunction
|
