1. No category

Mer om texter i MATLAB och om iterativ lösning av linjära ekvationssystem
Texter (strängar) i MATLAB skrivs omgivna av '' och behandlas
som vektorer, med samma operationer:
text = 'iss';
disp(['M' text text 'ipp' text(1)]) ger utskriften Mississippi
Ett diagonaltungt linjärt ekvationssystem Ax=b kan skrivas om
x = c + M * x, genom division av raderna och högerledet med
diagonalelementen.
M, iterationsmatrisen, består av de icke-diagonala elementen och
har max-normen || M || < 1 (definitionen av diagonaltung).
Iterationmetoden x(n+1) = c + M x(n), med x(0) = c som start,
konvergerar med feluppskattningen ||x(n+1) - x(n)|| ≈ || M ||n • ||c||
Yngve Sundblad
Föreläsning 3 sid.1
SF 1518/19 ht 2015
ACPU - 29 January 2004
7 sept.
Yngve Sundblad
Iterationsmetoden enligt Jacobi och Gauss-Seidel,
I MATLAB
x(n+1) = c + M x(n), med x(0) = c som start, n ekvationer
Jacobi: Komponenterna i x
(n+1)
beräknas direkt ur komponenterna i x
(n)
(n+1)
Gauss-Seidel: Komponenterna i x
beräknas ur de komponenter i (n+1)
x
som redan beräknats, i övrigt ur komponenter i x(n)
MATLAB, utgående från matrisen M, vektorn c, felgränsen err :
% Jacobi
% Gauss-Seidel
diffx = norm(c,Inf)
diffx = norm(c,Inf);
x=c; iter=0;
x=c; iter=0; n=rank(M); while diffx >= err
while diffx >= err
xold = x; iter = iter+1;
xold = x; iter = iter+1;
x = c + M*x;
for i=1:n
x(i) = c(i) + M(i,:)*x;
end
diffx=norm(x-xold,Inf)
diffx=norm(x-xold,Inf) end
end
disp(’Lösning:’); disp(x);
disp(’Lösning:’); disp(x);
disp(’Varv:’); disp(iter);
disp(’Varv:’); disp(iter);
OBS! norm(c,Inf) ger maxnorm. norm(c,inf) går lika bra.
Yngve Sundblad
Föreläsning 3 sid.2
SF 1518/19 ht 2015
ACPU - 29 January 2004
7 sept.
Yngve Sundblad
Jacobi och Gauss-Seidel på exempel 3.6
A=[2 10 0 -1;0 -1 1 5;5 1 0 0;-1 0 10 0], b=[30 25 10 70]'
Efter byte av rader till diagonaltung och division med
diagonalelementen och överflyttning av icke-diagonala:
M=[0 0.2 0 0; 0.2 0 0 -0.1; -0.1 0 0 0; 0 -0.2 0.2 0]
c=[2 3 7 5]'
Resulterande normer med err=0.00005:
Iteration
Jacobi
Gauss-Seidel
1
0,8
0,74 2
0,06
0,074
3
0,02
0,015
4
0,0036
0,00088
5
0,00088
0,000059
6
0,00023
0,0000039
7
0,000061
8
0,000015
Gauss-Seidel vinner klart i längden (men t.o.m. sämre i början),
extra decimal efter 6, Jacobi kräver 8 iterationer
Yngve Sundblad
Föreläsning 3 sid.3
SF 1518/19 ht 2015
ACPU - 29 January 2004
7 sept.
Yngve Sundblad
Minstakvadratmetoden
Ett vanligt problem i vetenskap och teknik är att man har ett
antal mätdata och vill anpassa en funktion / kurva till dessa.
Om man har lika många parametrar som data kan man
interpolera, vi återkommer till detta, men oftast har man fler
data än parametrar och vill anpassa med funktionen / kurvan.
Man ska alltså hantera ett överbestämt ekvationssystem.
Med minstakvadratmetoden anpassar man så att summan av
kvadraterna av avvikelserna mellan funktionsvärdena och
mätpunkterna minimeras.
Om mätvärdena är (y1 för x1), (y2 för x2), … (yn för xn)
ska alltså funktionen y = f(x) anpassas så att
∑k rk2 = ∑k (yk - f(xk))2 minimeras
Yngve Sundblad
Föreläsning 3 sid.4
SF 1518/19 ht 2015
ACPU - 29 January 2004
7 sept.
Yngve Sundblad
Linjär minstakvadratanpassning 1
Ofta vill man anpassa till en linjär kombination av
basfunktioner som en rät linje, ett högregradspolynom eller ett
uttryck som kan omformas till detta (t.ex. genom logaritmering).
Låt oss titta på fallet med två basfunktioner, f1(x), f2(x) :
y= F(x) = c1f1(x)+ c2f2(x) ska anpassas till (y1,x1),(y2,x2),… (yn,xn)
så att ∑k rk2 = ∑k (c1f1(xk) + c2f2(xk) - yk)2 minimeras
Då gäller om F = [ f1 f2 ] c att F ≈ y och att vi vill finna c så att || y – F ||2 (√kvadratsumman) minimeras
Geometriskt:
f1
F
f2
För att längden av y–F ska vara så liten
som möjligt måste den vara ortogonal
(vinkelrät) mot F, dvs både mot f1och f2,
som är kolumnerna i systemmatrisen:
T
y
Yngve Sundblad
T
AT(y – F)=0, AT(y – Ac)=0, A Ac=A y
Normalekvationerna
Föreläsning 3 sid.5
SF 1518/19 ht 2015
ACPU - 29 January 2004
7 sept.
Yngve Sundblad
Linjär minstakvadratanpassning 2 – rät linje
Analogt kan man, med mer skrivarbete, visa att det allmänt gäller
för linjär minstakvadratanpassning med godtyckligt många
basfunktioner, Ac = y, att ”bästa lösningen (koefficenterna)” ges
av normalekvationerna ATAc = ATy
I MATLAB finns inbyggt att för ett överbestämt linjärt
ekvationssystem Ac = y ger c=A\y minstakvadratlösningen.
Exempel 1(Pohl): Upphettad stång från T=20 grader, uppmätt
längdökning (y):
T=[20 30 40 50 60 70 80]'
y=[0.0 1.1 1.5 2.2 3.3 3.8 4.7]'
Anpassa till linjen y = c1 + c2T:
A=[1 1 1 1 1 1 1; 20 30 40 50 60 70 80]'
N=A’*A=[7 350; 350 20300]; z=A’*y=[16.6 1043]'
c=N\z ger [-1.4321 0.0761]; c=A\y ger samma.
Residualvektorn: r=y-Ac: 0.01*[-9 25 -11 -17 17 -9 4]
Yngve Sundblad
Föreläsning 3 sid.6
SF 1518/19 ht 2015
ACPU - 29 January 2004
7 sept.
Yngve Sundblad
Linjär minstakvadratanpassning 3 sönderfall
Exempel 2 (NAM (Eriksson)): Radioaktivt sönderfall:
t=[0:500:3000]'
x=[0.100 0.0892 0.0776 0.0705 0.0603 0.0542 0.0471]'
Anpassa till produktfunktionen x(t) = x0 exp(-kt)
Logaritmera: y = ln x = ln x0 – kt = c1 + c2t
MATLAB (med vektorerna t och x ovan): y = log(x); A=[ones(size(t)) t];
c=A\y; x0 = exp(c(1)) ger 0,1006
k = -c(2) ger 0,00025 (halveringstid ln2/k = 2772 sek, rimligt) xanp=x0*exp(-k*t) ger [0.1006 0.0888 0.0783 0.0691 0.0610 0.0538 0.0474]'
residual = x-xanp ger 1e-4*[-6 4 -7 14 7 4 3]'
Yngve Sundblad
Föreläsning 3 sid.7
SF 1518/19 ht 2015
ACPU - 29 January 2004
7 sept.
Yngve Sundblad
Plottning av anpassad kurva och avvikelser
subplot(2,1,1); plot(t,x,’*’,t,xanp);
subplot(2,1,2); plot(t,residual);
Yngve Sundblad
Föreläsning 3 sid.8
SF 1518/19 ht 2015
ACPU - 29 January 2004
7 sept.
Yngve Sundblad
Linjär minstakvadratanpassning 4 – planetbanor (Kepler)
Exempel 4.15 (EXS): Planetbanor (solavstånd R i jordradier, omloppstid T i år):
Merkurius (0,39 0,24) Venus (0,72 0,62) Mars (1,52 1,88)
Jupiter(5,20 11.9)
Saturnus (9,50 29,5) Uranus (19,0 84,0)
Anpassa T0 och q i T = T0 * Rq,
log(T) = log(T0) + q* log(R)
MATLAB: R=[0.39 0.72 1.52 5.2 9.5 19]'
T=[0.24 0.62 1.88 11.9 29.5 84]'
y = log(T); A=[ones(6,1) R];
c=A\y; T0 = exp(c(1)); q = c(2); ger T0=1,0008, q=1,5040 r=T-T0*R.^q ger 1e-3*[-3 -9 1.5 -45 -67 141]' r=T-R.^1.5 ger 1e-3*[-4 9 0.6 42 219 1181]' Utan Uranus (upptäckt 1781), med Jorden (1, 1): T0=1,0006, q=1,5036
Kepler (1619): Tre lagar för planetrörelse i ellipser, den tredje: T2 = R3 eller T = R1,5 (baserad på beräkningar från Tycho Brahes observationer)
Yngve Sundblad
Föreläsning 3 sid.9
SF 1518/19 ht 2015
ACPU - 29 January 2004
7 sept.
Yngve Sundblad
Polynomanpassning med minstakvadrat
Vid anpassning med polynom av grad n (n≥2) till data kan man
råka ut för stort konditionstal, dvs små fel i matrisen ger stora
fel i resulterande parametrarna.
Detta kan kraftigt förbättras genom centrering, dvs välja 1, x-xmedel , (x-xmedel) 2, (x-xmedel) 3, … som basfunktioner istället
för 1, x, x2, x3, …
Exempel 3 (NAM): x=[15:1:19]'; y=[15 8 5 3 3]';
Anpassa till y = F(x) = c1 + c2x + c3x2
A=[1 15 225; 1 16 256; 1 17 289; 1 18 324; 1 19 361]
A’A = [5 85 1455; 85 1455 25075; 1455 25075 434979]
A’y = [34 549 8929]; A\y ger [363.6 -39.3286 1.0714]
cond(A’*A) ger 3e9
Anpassa till y = F(x) = c1 + c2(x-17) + c3(x-17)2
B=[1 -2 4; 1 -1 1; 1 0 0; 1 1 1; 1 2 4]; B’B=[5 0 10; 0 10 0; 1 0 34]; B'y=[34 -29 83]';
B\y ger [4.6571 -2.9000 1.0714] cond(B'*B) ger 20
Yngve Sundblad
Föreläsning 3 sid.10
SF 1518/19 ht 2015
ACPU - 29 January 2004
7 sept.
Yngve Sundblad
Polynomanpassning i MATLAB
Funktionsanropet c=polyfit(x,y,n) ger koefficiemterna i
ett n:e-gradspolynom som anpassar sig till y-värdena för x-värdena med lämplig metod.
yv=polyval(c,xv) ger värdet av polynomet för xv.
I vårt exempel 3 med ursprungliga andragradspolynomet
c=polyfit(x,y,2)ger[1.0714 -39.3286 363.6000]'
OBS! Omvänd ordning, högsta koefficienten först, också i c för
polyval
polyval(c,17) ger 4.6571, stämmer med 4.6571-2.9000(x-17)+1.0714(x-17)^2
Yngve Sundblad
Föreläsning 3 sid.11
SF 1518/19 ht 2015
ACPU - 29 January 2004
7 sept.
Yngve Sundblad
Överkurs
I NAM (Eriksson), avsnitt 2.5-2.7 beskrivs för
minstakvadratmetoden
2.5 Residualanalys
2.6 Experimentell störningsräkning
2.7 Praktisk statistik
Det ingå inte i kursen och kommer inte på tentan, men är
läsvärt för speciellt intresserade.
Yngve Sundblad
Föreläsning 3 sid.12
SF 1518/19 ht 2015
ACPU - 29 January 2004
7 sept.
Yngve Sundblad