> restart: Zagadnienie własne macierzy
> with(LinearAlgebra): # CharacteristicMatrix, CharacteristicPolynomial, Eigenvalues, Eigenvectors,JordanBlockMatrix, JordanForm...
Idea i terminologia zagadnienia własnego macierzy
> n:=5;
:=
n
5
> A:=BandMatrix([-1,2,-1],1,n);
2 -1
0
0
0
-1
2 -1
0
0
:=
A
0 -1
2 -1
0
0 0 -1 2 -1
0
0
0 -1
2
> IsDefinite(A); # wszystkie wartości własne sa dodatnie true
A I
> C:=CharacteristicMatrix(A,alpha); # macierz charakterystyczna
2
1
0
0
0
1
2
1
0
0
:=
C
0
1
2
1
0
0
0
1
2
1
0
0
0
1
2
> A-alpha*IdentityMatrix(n); # to samo inaczej
2
-1
0
0
0
-1
2
-1
0
0
0
-1
2
-1
0
0
0
-1
2
-1
0
0
0
-1
2
> w:=CharacteristicPolynomial(A,alpha); # wielomian charakterystyczny
:= 5
w
10 436 356 235 6
> Determinant(C); # to samo inaczej
510 436 356 235 6
> alpha:=Eigenvalues(A); # wartości własne macierzy A
1
2
:= 3
2
3
2
3
> alpha:='alpha':
> solve(w=0,alpha); # sekwencja wartości własnych macierzy A
,
1 ,
2 ,
3
2
3 ,
2
3
> alpha,U:=Eigenvectors(A); # alpha - wartości własne, U - macierz utworzona z wektorów własnych (kolumnowych)
1 -1
1
1
-1
1
2
3 -1 3
3
1
0
, :=
U
2
3 ,
0
2
2
0 -1
3
1 3
3
-1
0
2 1
1
1
1
1
A u u
> alpha[1],U[1..n,1]; # wartość własna i odpowiadający jej wektor własny
-1
-1
,
1
0
1
1
> L:=simplify(A.U[1..n,1]);
-1
-1
L :=
0
1
1
> P:=simplify(alpha[1]*U[1..n,1]);
-1
-1
P :=
0
1
1
> for i to n do is(L[i]=P[i]) end do; true
true
true
true
true
> evalf(alpha);
1.
3.732050808
0.267949192
3.
2.
> sort(%);
0.267949192
1.
2.
3.
3.732050808
>
Ortogonalność wektorów wlasnych
> for i to n do
for j from i+1 to n do
print(simplify(Transpose(U[1..n,i]).U[1..n,j])); end do;
end do;
0
0
0
0
0
0
0
0
0
0
> for i to n do evalf(sqrt(Transpose(U[1..n,i]).U[1..n,i])); end do; # długosc wektorów własnych
2.
3.464101616
2.
1.732050808
Własności wartości i wektorów własnych
> A;
2 -1
0
0
0
-1
2 -1
0
0
0 -1
2 -1
0
0 0 -1 2 -1
0
0
0 -1
2
> alpha,U;
1 -1
1
1
-1
1
2
3 -1 3
3
1
0
2
3 ,
0
2
2
0 -1
3
1 3
3
-1
0
2 1
1
1
1
1
> add(alpha[i],i=1..5), Trace(A);
,
10 10
> expand(mul(alpha[i],i=1..5)), Determinant(A);
,
6 6
> simplify(U^(-1).A.U);
1
0
0
0
0
3 3 5
0
0
0
0
3 1
0
0
2
3
0
0
0
0
0
3 0
0
0
0
0
2
Normalizacja wektorów
> v:=<3,5,-4>;
3
:=
v
5
-4
> l:=sqrt(add(v[i]^2,i=1..3)); # długość wektora l := 5
2
> l:=sqrt(v.v);
l := 5
2
> l:=VectorNorm(v,Euclidean); l := 5
2
> l:=VectorNorm(v,2);
l := 5
2
1
n
p
p
x
x
, p 1,2,...
i
p
i 1
> p:=1;
:=
p
1
> nv1:=add(abs(v[i])^p,i=1..3)^(1/p); # norma pierwsza
:=
nv1
12
> nv1:=VectorNorm(v,1); # norma pierwsza
:=
nv1
12
> vu1:=v/nv1; # wektor unormowany
1
4
5
:=
vu1
12
-1
3
> sqrt(add(vu1[i]^2,i=1..3)); # dlugosc wektora unormowanego 5 2
12
> sqrt(vu1.vu1); # to samo inaczej 5 2
12
> p:=2;
:=
p
2
> nv2:=add(abs(v[i])^p,i=1..3)^(1/p); # norma druga nv2 := 50
> nv2:=VectorNorm(v,2); # norma druga
:=
nv2
5 2
> vu2:=v/nv2; # wektor unormowany
3 2
10
2
vu2 :=
2
2 2
5
> sqrt(vu2.vu2); # dlugosc wektora unormowanego 1
> # p-> infinity
x
max( x , i 1.. n)
i
> nmax:=max(seq(abs(v[i]),i=1..Dimension(v))); # norma maksimum nmax := 5
> nmax:=VectorNorm(v); # norma maksimum nmax := 5
> vum:=v/nmax; # wektor unormowany
3
5
vum := 1
-4
5
> sqrt(vum.vum); # dlugosc wektora unormowanego 2
> nv1;evalf(nv2);nmax; # porównanie norm 12
7.071067810
5
Jak szybko norma zmierza do normy maksimum?
> for p to 85 do
nv||p:=evalf(add(abs(v[i])^p,i=1..3)^(1/p)); end do:
>