< M A T L A B > Copyright 1984-2005 The MathWorks, Inc. Version 7.1.0.183 (R14) Service Pack 3 August 02, 2005 To get started, select MATLAB Help or Demos from the Help menu. >> % per fare delle prove ci fabbrichiamo due vettori m e s a caso >> n=5 n = 5 >> m=rand(1,n+1) m = 0.9501 0.2311 0.6068 0.4860 0.8913 0.7621 >> s=rand(1,n) s = 0.4565 0.0185 0.8214 0.4447 0.6154 >> gamma=sum(m) gamma = 3.9275 >> s/gamma ans = 0.1162 0.0047 0.2091 0.1132 0.1567 >> % siamo fortunati, s/gamma sono tutti tra 0 ed 1 e le ipotesi sono tutte soddisfatte >> % ora costruiamo a e b >> a=[m;diag(s),zeros(n,1)] a = 0.9501 0.2311 0.6068 0.4860 0.8913 0.7621 0.4565 0 0 0 0 0 0 0.0185 0 0 0 0 0 0 0.8214 0 0 0 0 0 0 0.4447 0 0 0 0 0 0 0.6154 0 >> b=a/gamma b = 0.2419 0.0589 0.1545 0.1237 0.2269 0.1940 0.1162 0 0 0 0 0 0 0.0047 0 0 0 0 0 0 0.2091 0 0 0 0 0 0 0.1132 0 0 0 0 0 0 0.1567 0 >> % ora siamo pronti a fare i calcoli sugli autovalori di b >> w=eig(b) w = 0.2700 0.0502 + 0.0524i 0.0502 - 0.0524i -0.0652 -0.0316 + 0.0570i -0.0316 - 0.0570i >> % calcoliamone i moduli >> for k=1:n+1 abs(w(k)) end ans = 0.2700 ans = 0.0726 ans = 0.0726 ans = 0.0652 ans = 0.0652 ans = 0.0652 >> % abbiamo trovato quello di modulo massimo, 0.2700 >> % disegnamo i cerchi di Gershgorin usando gershcircles.m (da A.Quarteroni,F.Saleri,P.Gervasio) >> gershcircles(b) >> % calcoliamo l'autovalore di modulo massimo con eigpower.m (da A.Quarteroni,F.Saleri,P.Gervasio) >> lambda=eigpower(b) lambda = 0.2700 >> % proviamo ad effettuare a mano alcuni passi del metodo delle potenze >> % inziamo da un vettore c di tutti 1 >> c=ones(n+1,1) c = 1 1 1 1 1 1 >> % si normalizza c, si moltiblica per b; na norma stima il massimo modulo dell'autovalore. >> c=b*(c/norm(c)); >> norm(c) ans = 0.4271 >> % passo successivo >> c=b*(c/norm(c)); >> norm(c) ans = 0.3369 >> % passo successivo >> c=b*(c/norm(c)); >> norm(c) ans = 0.2938 >> % passo successivo >> c=b*(c/norm(c)); >> norm(c) ans = 0.2763 >> % passo successivo >> c=b*(c/norm(c)); >> norm(c) ans = 0.2696 >> % passo successivo >> c=b*(c/norm(c)); >> norm(c) ans = 0.2695 >> % passo successivo >> c=b*(c/norm(c)); >> norm(c) ans = 0.2698 >> % passo successivo >> c=b*(c/norm(c)); >> norm(c) ans = 0.2699 >> % passo successivo >> c=b*(c/norm(c)); >> norm(c) ans = 0.2700 >> % vediamo a che punto siamo arrivati: quanto buono sembra c come autovettore e |c| come autovalore ? c = 0.2480 0.1067 0.0019 0.0014 0.0006 0.0004 >> b*c-norm(c)*c ans = 1.0e-06 * 0.3910 0.2725 -0.0255 -0.1212 -0.1871 -0.4489 >> % sembrano risultati ragionevoli, con 5 cifre decimali >> % a questo punto si possono scrivere delle piccole routines per fare automaticamente i vari passi