1. No category

TAIU07 Föreläsning 5
• Linjära ekvationssystem.
• Minsta kvadrat problem. Tillämpning: Cirkelpassning.
• Geometriska objekt. Translationer. Rotationer.
• Funktioner som inargument. Tillämpning: Derivata.
15 februari 2016 Sida 1 / 32
Linjära Ekvationssystem
Division mellan matriser skall tolkas som att A/B betyder AB−1
och A\B betyder A−1B.
Exempel Lös ekvationssystemet




1
2 1 3
 4 −2 5  x =  −2  .
3
−2 2 −3
med hjälp av MATLAB.
15 februari 2016 Sida 2 / 32
I Matlab skapar vi först matrisen och högerledet med
>> A = [ 2 1 3 ; 4 -2 5 ; -2 2 -3 ];
>> b = [ 1 -2 3]’;
% OBS kolumnvektor!
Vi löser sedan ekvationssystemet med x = A−1b som beräknas med
\-operatorn. Alltså
>> x = A\b;
ans =
1.8333
x’
1.3333
-1.3333
Kontrollera genom att beräkna residualen
>> r = b - A*x;
ans =
1.0e-15 *
0.8882
r’
0
-0.4441
15 februari 2016 Sida 3 / 32
Funktionen inv beräknar inversen till en kvadratisk matris.
Exempel Beräkna lösningen till Ax = b med hjälp av inversen
>> A
>> b
>> x
ans
= [ 2 1 3 ; 4 -2 5 ; -2 2 -3 ];
= [ 1 -2 3]’;
= inv(A)*b; x ’
=
1.8333
1.3333
-1.3333
Ger samma lösning som tidigare. Beräkna nu A−1A med
>> C = inv(A)*A
C =
1.0000
0
0
1.0000
-0.0000
0
0
-0.0000
1.0000
15 februari 2016 Sida 4 / 32
Exempel En linje i R3 kan skrivas på formen

  
1
2
P(t) = ~p1 + ~p2t =  −1  +  2  t.
3
2
Ett plan i R3 spänns upp av de två vektorerna




1
0
~v1 =  −1  och ~v2 =  1  .
0
−2
Var skär linjen planet? Formulera problemet som ett linjärt ekvationssystem.
15 februari 2016 Sida 5 / 32
Lösning Alla punkter i planet kan skrivas som en linjär kombination av ~v1 och
~v2. Vi söker alltså α1, α2 och t sådana att
  





1
2
0
1
 −1  α1 +  1  α2 =  −1  +  2  t.
2
3
−2
0
I Matlab:
>> v1=[1 -1 0]’;v2=[0 1 -2]’;p1=[2 -1 3]’;p2=[1 2 2]’
>> x=[ v1 v2 -p2]\p1; x’
ans =
1.3750
-0.8750
-0.6250
Skärningspunkten ges av
>> P = v1*x(1)+v2*x(2); P’
ans =
1.3750
-2.2500
1.7500
15 februari 2016 Sida 6 / 32
Minsta Kvadratproblemet
Exempel Vi vet att en fysikalisk storhet ges av
y = c1 x1 + c2 x2 ,
där c1 och c2 är konstanter vi vill bestämma och x1 och x2 är parametrar vi
enkelt kan variera. Hur skall vi bestämma c1 och c2?
Lösning Genomför experiment och mät y för olika värden på (x1, x2).
x1
3.0
2.7
1.8
4.2
x2
1.60
1.35
0.70
1.20
y
12.6
11.2
6.3
13.2
=⇒

3.0
 2.7

 1.8
4.2



1.60 12.6
 11.2 
1.35 
c
1



=
 6.3 
0.70  c2
1.20
13.2
Nöjer vi oss med två mätningar får vi ett linjärt ekvationssystem. Fler
mätningar bör gå att utnyttja för att få ett bättre resultat.
15 februari 2016 Sida 7 / 32
Euklidisk Längd
Den Euklidiska normen av en vektor x är
1/2

n
X
xi2 .
kxk2 = 
i=1
I Matlab finns funktionen norm som beräknar normen kxk2.
I Matlab
>> x = [ 1 3 -2 ]
>> E = norm( x )
E =
3.7417
15 februari 2016 Sida 8 / 32
Definition Låt A ∈ Rm×n m > n. Ekvationssystemet Ax = b är då
överbestämt.
Fler ekvationer än obekanta gör att det oftast inte finns en lösning. Istället
minimerar vi residualen och löser
minn kAx − bk.
x∈R
Detta kallas Minsta kvadratproblemet. I statistiken kallas det Linjär
Regression.
Sats Minsta kvadrat lösningen x satisfierar normalekvationerna
AT Ax = AT b.
Om kolumnerna i A är linjärt oberoende så är AT A icke-singulär.
15 februari 2016 Sida 9 / 32
I MATLAB tolkas x=A\b som ett minsta kvadrat problem om A
är rektangulär. Det x som minimerar kAx − bk2 beräknas.
Exempel Lös Minsta kvadratproblemet med A och b som tidigare
>> A = [3.0 1.60 ; 2.7 1.35 ; 1.8 0.70 ; 4.2 1.20];
>> b = [12.6 11.2 6.3 13.2]’;
>> x_MK = A\b
% Eller (A’*A)\( A’*b )
x_MK =
1.8734
4.4208
Vi får minsta kvadrat lösningen x = ( 1.873 , 4.420 )T .
Hade valt x = ( 2 , 4 )T och b med ett slumpmässigt fel.
15 februari 2016 Sida 10 / 32
Modell Passning
Exempel Vi will anpassa ett polynom p(t) = c0 + c1t + c2t2 till en uppsättning
mätningar (ti, yi), 1 ≤ i ≤ m.
1
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
−0.2
0
0.2
0.4
0.6
0.8
1
1.2
Formulera detta som ett överbestämt ekvationssystem och lös med minsta
kvadratmetoden.
15 februari 2016 Sida 11 / 32
Lösning Antag att data (ti, yi) lagrats som två vektorer t och x. Då fås
>> A=[ t.^0 t.^1 t.^2]; b=y; c= A\b;x
>> tt=0:0.1:2;yy=c(1)+c(2)*tt+c(3)*tt.^2;
>> plot(t,y,’r+’,tt,yy);
1.2
1
0.8
0.6
0.4
0.2
0
−0.2
−0.2
0
0.2
0.4
0.6
0.8
1
1.2
1.4
Vi får en vektor
>> c’
ans =
0.0400
-0.1073
0.7112
15 februari 2016 Sida 12 / 32
Exempel Vi har ett antal punkter (xk , yk ) som borde ligga på en cirkel.
6
5.5
5
4.5
4
3.5
3
2.5
2
1.5
−2
−1
0
1
2
3
4
Hur kan vi använda Minstakvadrat metoden för att uppskatta centrum och radie
för cirkeln?
15 februari 2016 Sida 13 / 32
Lösning Varje punkt (x, y)T på cirkeln uppfyller ekvationen
F(u) = a(x2 + y2) + b1x + b2y + 1 = 0,
där u = (a, b1, b2)T är okända parametrar.
Centrum och radie för cirkeln ges av
b1 b2 T
T
z = (x0, y0) = −( , ) ,
2a 2a
2 + b2
b
1
2
1
2
− .
och, r =
2
a
4a
Vi behöver två funktioner
1. Givet x och y beräkna u:
function [u]=CircleFit( x , y )
2. Givet u beräkna z och r samt plotta cirkeln:
function [z,r]=DisplayCircle( u )
15 februari 2016 Sida 14 / 32
På filen CircleFit.m skriver vi
% Hitta parameter vektor u som beskriver en cirkel
% vektorerna x och y.
%
function [u]=CircleFit( x , y )
A=[x.^2+y.^2 x y];b=-ones(length(x),1);
u=A\b;
end
Vi använder funktioner på våra data vektorer x och y och får
>> u=CircleFit( x , y)
u =
0.1729
-0.4113
-1.1799
15 februari 2016 Sida 15 / 32
På filen DisplayCircle.m skriver vi
% Beräkna centrum och radie samt rita upp
% cirkeln givet parameter vektorn u.
%
function [z,r]=DisplayCircle( u )
z=-u(2:3)/2/u(1)
r=sqrt( (u(2)^2+u(3)^2)/4/u(1)^2-1/u(1));
theta=2*pi*(0:99)’/100;
x=z(1)+r*cos(theta);y=z(2)+r*sin(theta);
x=[x;x(1)];y=[y;y(1)];
plot(x,y);
end
15 februari 2016 Sida 16 / 32
Vi använder nu funktionen för att rita cirkeln
>> hold on, [z,r]=DisplayCircle( u ); hold off
6
5.5
5
4.5
4
3.5
3
2.5
2
1.5
1
−2
−1
0
1
2
3
4
Vi fick z = (1.1898, 3.4131)T och r = 2.6981. Hade valt z0 = (1.2, 3.4)T och
r = 2.7.
15 februari 2016 Sida 17 / 32
Geometriska Objekt och Avbildningar
Exempel En månghörning består av hörpunkter och kantlinjer. Varje punkt
Pk = (xk , yk )T kan ses som en 2 × 1 vektor.
Vi vill rita upp en femhörning med givna hörnpunkter
1
3
5
6
3
,
,
,
och
.
3
6
5
2
1
Vi vill sedan förflytta alla punkterna med en konstant vektor T = (1 , 2)T och
rita upp månghörningen igen. Hur skall vi göra?
15 februari 2016 Sida 18 / 32
Spara punkterna i en 2 × n matris P. Varje punkt representeras av en kolumn
P(:,k). I Matlab
P = [
1 3 5 6 3 ;
3 6 5 2 1 ];
plot( P(1,:),P(2,:),’r+’);
hold on
ind = [1 2 3 4 5 1 ];
plot( P(1,ind),P(2,ind),’b-’);
För att translatera punkterna skapar vi en ny 2 × n matris att addera till P.
T=ones(2,5); T(2,:)=2; P2 = P+T;
plot( P2(1,ind),P2(2,ind),’k--’)
hold off
15 februari 2016 Sida 19 / 32
Sätt koordinataxlarna så att vi ser hela bilden och skriv ut
axis([ 0 8 0 9])
print -depsc F2-5horning.eps
Den färdiga bilden blir
9
8
7
6
5
4
3
2
1
0
0
1
2
3
4
5
6
7
8
Kommentar Många olika operationer inom datorgrafik består av att enkla
operationer utförs på punktmängder. I Matlab kan man göra rätt mycket om
man är lite kreativ.
15 februari 2016 Sida 20 / 32
Linjära Avbildningar
Definition En n × m matris A kan ses som en avbildning
A : x ∈ Rm 7→ y ∈ Rn.
där
y = Ax.
Exempel Enhetsmatrisen I representerar en avbildning I : x 7→ x. Den bildas i
Matlab med kommandot
>> n = 10;
>> I = eye( n );
Andra viktiga avbildningar är rotationer, speglingar, och projektioner.
15 februari 2016 Sida 21 / 32
Rotationer
Lemma En punkt P = (x1 , x2)T roteras kring origo (0 , 0)T , med
vinkel α, med hjälp av
cos(α) sin(α)
′
P = G · P,
G=
.
− sin(α) cos(α)
Exempel Vi vill illustrera rotationer i Matlab genom att rotera en punkt
P = (4 , 2)T kring ett origo som satts till O = (2 , 3)T . Vridningsvinkeln skall
vara α = 2π/3.
Rita upp relevanta vektorer både före och efter vridningen. Samt rita upp en
cirkel med centrum i punkten O för att tydligare visa att vi fått en vridning.
15 februari 2016 Sida 22 / 32
Lösning Skapa först vektorer för att representera Origo O = (2 , 3)T och
punkten P. Beräkna även skillnadsvektorn R.
Origo=[ 2 ; 3 ]; P=[ 4 ; 2 ]; R = P - Origo;
Plotta nu vektorerna i blått
clf,hold on
plot( [0 P(1)],[0,P(2)] ,’b’ )
plot( [Origo(1) P(1)],[Origo(2),P(2)] ,’b’ )
15 februari 2016 Sida 23 / 32
Skapa en mängd punkter på cirkeln med origo i O = (2 , 3)T .
Radie=sqrt( sum( R.^2 ) );
alpha=0:0.01:2*pi;
x=Origo(1)+Radie*cos( alpha );
y=Origo(2)+Radie*sin( alpha );
plot(x,y,’r--’)
15 februari 2016 Sida 24 / 32
Skapa sedan vridningspatrisen G och applicera den på skillnadsvektorn R.
alpha=2*pi/3;
G = [cos(alpha),-sin(alpha);sin(alpha),cos(alpha
R2=G*R; P2=Origo+R2;
Rita nu upp de nya punkterna i svart.
plot( [0 P2(1)],[0,P2(2)] ,’k’ )
plot( [Origo(1) P2(1)],[Origo(2),P2(2)] ,’k’ )
axis equal
hold off
Kommentar Kommandot axis equal gör koordinataxlarna lika långa i
fönstret. Annars ser inte cirkeln rund ut.
15 februari 2016 Sida 25 / 32
5
4.5
4
3.5
3
2.5
2
1.5
1
0.5
0
-1
0
1
2
3
4
5
Den färdiga grafen innehåller de objekt som behövs för att illustrera
vridningen. Punkten blir
>> P2’
ans =
1.8660
5.2321
15 februari 2016 Sida 26 / 32
Funktioner som inargument
Kommandot feval kan användas för att anropa en funktion med
givna inparametrar.
Exempel Har vi skrivit en fil funk.m som innehåller
function [z]=funk( x , y )
z=x^2+3*y;
end
så kan vi anropa funktionen med inargument x = 2 och y = 2.5 genom att
skriva
>> feval( @funk , 2 , 2.5 )
ans =
11.5000
15 februari 2016 Sida 27 / 32
Funktioner som inargument
Exempel Definitionen av derivata är
f (x + h) − f (x)
′
f (x) = lim
.
h
h→0
Välj ett fixt h > 0 för att approximera derivatan av f (x).
Skriv en funktion Derivata.m som kan anropas med
>> Df = Derivata( @funk , x , h );
Använd funktionen för att plotta derivatan av sin(x) på intervallet [0, 2π].
15 februari 2016 Sida 28 / 32
Lösning På filen Derivata.m skriver vi
function [Df]=Derivata( f , x , h )
f1=feval( f , x + h );
f0=feval( f , x );
Df=(f1-f0)/h;
end
Vi kan sedan plotta derivatan av sin(x) genom att skriva
>> x = 0:0.01:2*pi;
>> plot( x , Derivata( @sin , x , 10^-3 ) );
15 februari 2016 Sida 29 / 32
Enkla funktioner
Enkla funktioner kan skapas direkt i Matlab. Man skriver
>> f = @( <variabler> )
<uttryck>;
Man får då ett handtag till en funktion.
Exempel Skriv
>> f = @(x) exp(2*x).*cos(x.^2)+4*x
f =
@(x)exp(2*x).*cos(x.^2)+4*x
Funktionen kan beräknas med feval eller direkt f(4).
15 februari 2016 Sida 30 / 32
Exempel Skapa funktionerna
p
f (x) = 1 − x2e−2x,
och
g(x) = (x12 − 1 ,
p
1 + x2)T .
Tänk på att g(x) måste ha en vektor som inagrument och returnera en vektor i
R2.
Försök få f (x) att acceptera vektor argument så att det går att skriva
>> x = 0:0.01:1;
>> plot( x , f(x) )
15 februari 2016 Sida 31 / 32
Funktionen integral beräknar integraler. Den anropas
>> I = integral( fun , a , b )
där f är ett funktionshandtag.
Exempel Beräkna samma integral som tidigare genom
>> f = @(x) exp(2*x).*cos(x.^2)+4*x;
>> I = integral( f , 0 , 1)
I =
4.6736
Det finns även integral2 för dubbelintegraler.
15 februari 2016 Sida 32 / 32