File: ieigen.en

package info (click to toggle)
euler 1.61.0-12
  • links: PTS, VCS
  • area: main
  • in suites: bookworm, bullseye, forky, sid, trixie
  • size: 5,164 kB
  • sloc: ansic: 24,761; sh: 8,314; makefile: 141; cpp: 47; php: 1
file content (41 lines) | stat: -rw-r--r-- 1,227 bytes parent folder | download | duplicates (8)
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
% Let us demonstrate the use of interval methods for the
% inclusion of a Eigenvalue.
% 
% We first generate a random postive definite matrix.
>A=random(20,20); A=A'.A;
% Then we compute its Eigenvalues in the usual way.
>l=sort(re(eigenvalues(A)));
% We take the first Eigenvalue.
>longformat; l=l[1]
      0.003318934192275 
>{l,x}=xeigenvalue(A,l); l,
      0.003318934192276 
% We normalize it in a special way.
>x=x/x[1];
% We now set up the function (l,x2,...,xn) -> Ax-lx with
% x1=1.
>function f(x,A)
$n=cols(x); return ((A-x[1]*id(n)).(1|x[2:n])')';
$endfunction
% This is the start of our algorithm.
>lx=l_x[2:cols(A)];
% We could use the Broyden method to determine the zero
% of this function.
>yb=broyden("f",lx';A); yb[1],
      0.003318934192276 
% We set up the derivative of f.
>function f1(x,A)
$n=cols(x); B=A-x[1]*id(n); B[:,1]=-(1|x[2:n])'; return B;
$endfunction
% We expand lx, so that it probably contains a solution
% of f(x)=0.
>ilx=expand(~lx~,1000000000);
% The interval Newton method now proves a guaranteed
% inclusion of the Eigenvalue.
>{y,v}=inewton2("f","f1",ilx';A); y[1]
    ~0.003318934192235,0.003318934192317~ 
% v=1 means, that the inclusion is verified.
>v
                      1 
>
>