Harjoitus 3, Matlabia

Järkevät tehtävänratkaisut ovat

>>>>>>>>>>>>>
>>>>>>>>>>>>>

tällaisissa laatikoissa,

>>>>>>>>>>>>
>>>>>>>>>>>>

tässä tiedostossa on myös muuta suttua.





>>>>>>>>>>
>>>>>>>>>>

h3t1

close
c=0;				%nämä voidaan arpoa
b=[1 2 3]';

A=[3 -2 c;-1 4 2;1 -1 5];	%tehdään iterointiin tarvittavat C ja e
D=diag(diag(A));
B=A-D;
C=-inv(D)*B;
e=inv(D)*b;

x(:,1)=[1 1 1]';		%arvotaan alkuarvaus

for k=1:10			%iteroidaan
x(:,k+1)=C*x(:,k)+e;
end

z=x;				%piirretään iterointipolku
r=z(1,:);t=z(2,:);y=z(3,:);plot3(r,t,y)
hold on

y=A\b;				%oikea ratkaisu

z=y;				%piirretään oikea vastaus iterointipolulle
r=z(1,:);t=z(2,:);y=z(3,:);plot3(r,t,y)


Ylimääräistä pohdintaa:

Kynällä ja paperilla nähdään, että iteraatio suppenee, kun |c|<3/2. Ratkaisu on tällöin olemassa kaikille vektoreille b, mutta miten tämä näkyy A:n determinantista? 
Piirretään C=det(A) muuttujan c funktiona välillä -2:2


c=-2:0.1:2;
for i=1:length(c)
C(i)=det([3 -2 c(i);-1 4 2; 1 -1 5]);
end
plot(c,C)

Kuvasta nähdään, että det(A)>>0 näillä c:n arvoilla eli ratkaisu on olemassa

Tallensin kuvan nimellä "h3t1onkoratkaisua.jpg" komennolla

print -djpeg 'h3t1onkoratkaisua.jpg'

>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>


>>>>>>>>>>>>>
>>>>>>>>>>>>

h3t2

Alla alfa=a.

a=0:0.01:1;
Tnormi1=max(abs(1-7*a)+5*abs(a),2*abs(a)+abs(1-4*a));
Tnormioo=max(abs(1-7*a)+2*abs(a),5*abs(a)+abs(1-4*a));
plot(a,Tnormi1,'r')

hold on
plot(a,Tnormioo,'b')
yks=ones(1,length(a));
plot(a,1)

Kuvasta nähdään, että sopivalla a:n arvolla ||T||_1<1, mutta että ||T||_oo on vähintään 1 kaikilla a:n arvoilla. Millä a:n arvoilla normit saavuttavat miniminsä?
Kun kuva on piirrettynä, voidaan käyttää komentoa ginput:

Tnormi1min=ginput

Tnormi1min =

    0.1401    0.7262



%Siis Tnormi1 minimoituu arvolla a=0.1401 eli a=1/7.
%kuvasta nähdään, että Tnormioo on aina suurempi kuin 1, mutta että Tnormioo = 1, kun a on 0. Jos a =/= 0, niin Tnormioo>1.

>>>>>>>>>>>>>>
>>>>>>>>>>>>>>

Muunnelma:

close
c=0;				
b=[1 2 3]';
A=[3 -2 c;-1 4 2;1 -1 5];	
D=diag(diag(A));
B=A-D;
C=-inv(D)*B;
e=inv(D)*b;
x(:,1)=[1 1 1]';		
for k=1:100			
x(:,k+1)=C*x(:,k)+e+randn(3,1)/k; %lisätään satunnaisuutta
end
z=x;				
r=z(1,:);t=z(2,:);y=z(3,:);plot3(r,t,y)


Satunnaisliikettä:

close
x(:,1)=[1 1 1]';		
for k=1:100			
x(:,k+1)=randn(3,1)/k; %lisätään satunnaisuutta
end
z=x;				
r=z(1,:);t=z(2,:);y=z(3,:);plot3(r,t,y)


Arvotaan c,b ja x(1):

close
c=randn(1);
b=randn(3,1);
A=[3 -2 c;-1 4 2;1 -1 5];
D=diag(diag(A));
B=A-D;
C=-inv(D)*B;
e=inv(D)*b;
x(:,1)=randn(3,1);
for k=1:10
x(:,k+1)=C*x(:,k)+e;
end
z=x;
r=z(1,:);t=z(2,:);y=z(3,:);plot3(r,t,y)
hold on
y=A\b;
z=y;
r=z(1,:);t=z(2,:);y=z(3,:);plot3(r,t,y)


Tehdään video aiheesta, jolloin voidaan katsoa arvottujen tilanteiden kuvia.
Video f voidaan tämän jälkeen katsoa komennolla movie(f,1,2). Tällöin
katsotaan video 1 kerran nopeudella 2 kuvaa sekunnissa


for i=1:20		%tehdään 20 kuvan video
close
c=randn(1);
b=randn(3,1);
A=[3 -2 c;-1 4 2;1 -1 5];
D=diag(diag(A));
B=A-D;
C=-inv(D)*B;
e=inv(D)*b;
x(:,1)=randn(3,1);
for k=1:10
x(:,k+1)=C*x(:,k)+e;
end
z=x;
r=z(1,:);t=z(2,:);y=z(3,:);plot3(r,t,y)
hold on
y=A\b;
z=y;
r=z(1,:);t=z(2,:);y=z(3,:);plot3(r,t,y)
f(i)=getframe;
end

movie(f,1,2)		%katsotaan video 1 kerran nopeudella 2 kuvaa sekunnissa


Arvotaan A,b ja x(1) ja tehdään näistä tilanteista video:

for i=1:20		%tehdään 20 kuvan video
close

b=randn(3,1);
A=randn(3);		%arvotaan A suoraan tekemättä muuttujaa c
D=diag(diag(A));
B=A-D;
C=-inv(D)*B;
e=inv(D)*b;
x(:,1)=randn(3,1);
for k=1:10
x(:,k+1)=C*x(:,k)+e;
end
z=x;
r=z(1,:);t=z(2,:);y=z(3,:);plot3(r,t,y)
hold on
y=A\b;
z=y;
r=z(1,:);t=z(2,:);y=z(3,:);plot3(r,t,y)
g(i)=getframe;
end


Tehdään pupuhyppely. Supetaan pisteeseen, jonka jälkeen muutetaan systeemiä ja supetaan toiseen pisteeseen jne.

n=0;				%alussa on 0 pistettä muistissa
x(:,1)=randn(3,1);		%arvotaan ensimmäinen alkuarvo

for i=1:20			%pupuhypellään 20 kertaa

b=randn(3,1);			%arvotaan, minne suuntaan mennään tekemällä b ja A
A=randn(3);
D=diag(diag(A));		%tehdään iterointiin tarvittavat matriisi C ja vektori e
B=A-D;
C=-inv(D)*B;
e=inv(D)*b;

yksnormiC=max(sum(C));		%lasketaan matriisin C 1-normi

if abs(yksnormiC)<1,		%jos iterointi niin iteroidaan, muuten ei

for k=1:10
x(:,k+1+n)=C*x(:,k+n)+e;	%aluksi iterointi alkaa arvosta x(:,2) ja käytetään arvoa x(:,1) ok
				%yhden iteroinnin jälkeen x(:,12) ja käytetään arvoa x(:,11) ok
end				%lopetetaan iterointi
n=n+10;				%laitetaan pisteet talteen
end				%lopetetaan if
end				%lopetetaan




n=0;				%alussa on 0 pistettä muistissa
x(:,1)=randn(3,1);		%arvotaan ensimmäinen alkuarvo

for i=1:20			%pupuhypellään 20 kertaa

C=2*eye(3);
while max(sum(C))>=1
b=randn(3,1);			%arvotaan, minne suuntaan mennään tekemällä b ja A
A=randn(3);
D=diag(diag(A));		%tehdään iterointiin tarvittavat matriisi C ja vektori e
B=A-D;
C=-inv(D)*B;
e=inv(D)*b;
end

for k=1:10
x(:,k+1+n)=C*x(:,k+n)+e;	%aluksi iterointi alkaa arvosta x(:,2) ja käytetään arvoa x(:,1) ok
				%yhden iteroinnin jälkeen x(:,12) ja käytetään arvoa x(:,11) ok
end				%lopetetaan iterointi
n=n+10;				%laitetaan pisteet talteen
end				%lopetetaan

>>>>>>>>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>>>>>>>>>

h3t6

Tehdään, mitä tehtävänannossa pyydetään ja piirretään aiheesta kuva. Kuvassa iteraatin z-koordinaatti vastaa iteraatin
järjestyslukua, esim. x(2)=(x^2_1,x^2_2,2)

a=3*(rand(1)-0.5);
x(:,1)=[a,a];
for k=1:10
b=x(:,k);
N=[2*b(2) 2;2 2*b(1)]/(4*(b(1)*b(2)-1));
x(:,k+1)=N*x(:,k);
end

z=x;r=z(1,:);t=z(2,:);y=1:length(z);plot3(r,t,y)

Nähdään, että iteraatit ovat pisteiden (0,0) ja x(1) kautta kulkevalla suoralla ja että ne lähestyvät
ratkaisua (0,0) niin, että virhe pienenee monotonisesti, oskilloiden ratkaisun kummallakin puolella.

Miten virhe pienenee?

Olkoon

y=[0 0]';		%asetetaan jompi kumpi, jotta saadaan tarkka arvo
y=x(:,end);

q=max(abs(x-y*ones(1,length(x))));	%vähennetään x:n jokaisesta sarakkeesta oikea ratkaisuvektori, jolloin saadaan virhevektoreista koostuva matriisi
					%käytetään max-normia, se saadaan yksinkertaisesti

%piirretään kuva virheestä ja virheen logaritmista iteraattien järjestyslukujen funktiona
%nähdään, että virhe pienenee "lineaarisesti"

plot(q,1:length(q))			
%print -djpeg h3t6_1.jpg

Q=log(q);
figure
plot(Q,1:length(Q))	
%print -djpeg h3t6_2.jpg

>>>>>>>>>>>>>>
>>>>>>>>>>>>>>


x(:,1)=3*(rand(2,1)-0.5);
for k=1:10
b=x(:,k);
N=[2*b(2) 2;2 2*b(1)]/(4*(b(1)*b(2)-1));
x(:,k+1)=N*x(:,k);
end


Satunnaiskävely

N=1000;

x(:,1)=[0 0];
for k=1:N
x(:,k+1)=x(:,k)+randn(2,1);
end
a=min(x(1,:));b=max(x(1,:));c=min(x(2,:));d=max(x(2,:));
figure
hold on
axis([a b c d])
axis tight
for i=1:N+1
plot(x(1,1:i),x(2,1:i))
f(i)=getframe;
end


Jotain ihan ylimääräistä

A=[3 4;-2 -1]/5.1;
det(A)
eig(A)
x(:,1)=[3 3]';
for k=1:10
x(:,k+1)=A*x(:,k);
end
a=x(1,:);
b=x(2,:);
plot(a,b)
for k=1:10
plot(a(1:k),b(1:k))
axis([-1 5 -2 3])
f(k)=getframe;
end
save mallivideo.mat