Harj.4 teht. 2 - 6
===================
Matlabilla:


Teht. 2: (Tässä Matlab, voidaan hyvin tehdä myös Maplella, vrt. ohje)
==========

>> A=rand(50,50);
>> spektri(A)
>> title('50 x 50 -satunnaismatrisiin spektri) 

>> print
>> B=A'*A;
>> sum(sum(B-B'))

ans =

     0

>> % Kylla siis oli symmetrinen
>> % Kokeillaan kuitenkin pienemmällä ja A+A':lla, jotta mittakaava
>> % ei mene tolkuttomaksi.

>> A=rand(10,10);spektri(A+A')
>> title('10x10' symmetrinen matriisi A+A''')

>> print -dgif8 symm.gif
>> A=rand(50,50);spektri(orth(A))             
>> title('50x50 symmetrinen matriisi orth(rand(50,50))') 
>> print
>> print -dgif8 orth.gif

Tehdään vielä uudestaan subplotilla:
>> A=rand(50,50);
>> B=rand(10,10);B=B+B';
>> C=orth(A);
>> subplot(2,2,1);spektri(A);subplot(2,2,2);spektri(B);subplot(2,2,3);spektri(C)
------------------------------------------------------------------

Kaikenlaista voi kokeilla, kuvat voi myös tehdä Maplella, tutkitaan, miten
helpoimmin.

Kas näin:

A:=randmatrix(20,20);
A:=convert(op(A),float);     # jotta eigenvals laskisi numeerisesti
eig:=eigenvals(A);
c:=seq([Re(eig[i]),Im(eig[i])],i=1..20);  # [x1,y1],[x2,y2],...
plot([c],style=point);       # Piirtää listassa olevat pisteet
r:=max(op(map(abs,[eig])));  # spektraalisäde

-------------------
plot([r*cos(t),r*sin(t),t=0..2*Pi]):  # parametrikäyrien syntaksi
plot([c],x=-100..100,y=-100.100,style=point): # x-ja y-välit kannattaa antaa
                                              # ainakin Re- tai Im-akselilla
                         # sijaitsevia pistejoukkoja piirrettäessä
Kuvat saadaan samaan näin:

with(plots): # Ladataan lisägrafiikkapakkaus, jossa mm. display-funktio
p1:=plot([r*cos(t),r*sin(t),t=0..2*Pi]):  # muista kaksoispiste, eikä 
p2:=plot([c],x=-100..100,y=-100.100,style=point): # puolipistettä (tai kärsi)
display({p1,p2}):

Teht. 3
=======
>> A=[8 0 3;2 2 1;2 0 3],T=[0 1 0;1 0 0;0 0 1]

A =

     8     0     3
     2     2     1
     2     0     3


T =

     0     1     0
     1     0     0
     0     0     1

>> [V,D]=eig(A)

V =

         0    0.9045   -0.0723
    1.0000    0.3015   -0.9869
         0    0.3015    0.1446


D =

     2     0     0
     0     9     0
     0     0     2

>> [V1,D1]=eig(inv(T)*A*T)

V1 =

    1.0000    0.3015   -0.9869
         0    0.9045   -0.0723
         0    0.3015    0.1446


D1 =

     2     0     0
     0     9     0
     0     0     2

Ominaisarvot säilyvät, ominaisvektorit eivät.
---------------------------------------------

----------------------------------------------------------------------

Teht. 4
=======

A=[.5 .25 .25;.25 .5 .25;.25 .25 .5]

>> sum(A)

ans =

     1     1     1

>> u0=rand(3,1)

u0 =

    0.2478
    0.8416
    0.3070

>> u0=A*u0

u0 =

    0.4110
    0.5595
    0.4259

>> u0=A*u0

u0 =

    0.4519
    0.4890
    0.4556

>> u0=A*u0
u0 =

    0.4621
    0.4713
    0.4630

>> u0=A*u0

u0 =

    0.4646
    0.4669
    0.4648

>> u0=A*u0

u0 =

    0.4653
    0.4658
    0.4653

Uudestaan toisesta alkup:

>> u0=[-4 100 -50]'

>> u0=(A^10)*u0


u0 =

   15.3333
   15.3334
   15.3333

>> sum(A')

ans =

     1     1     1

Myös rivisummat = 1 => 1 on ominaisarvo ja [1,1,1] vastaava om. vekt.

>> [V,D]=eig(A)

V =

   -0.8054    0.1340    0.5774
    0.5187    0.6305    0.5774
    0.2867   -0.7645    0.5774


D =

    0.2500         0         0
         0    0.2500         0
         0         0    1.0000

No niinpä on

>> D^10
ans =

    0.0000         0         0
         0    0.0000         0
         0         0    1.0000