1.5. Minsta kvadratmetoden

forts. p˚a f¨oreg˚aende f¨orel¨asning:

Som vi ser, subtraheras i den inre slingan skal¨arprodukten Xi−1

j=1

`_ijx_j = L(i, 1 : i − 1) ∗ x(1 : i − 1)

fr˚an b_i, s˚a att vi kan vektorisera programfragmentet p˚a f¨oljande s¨att:

x(1) = b(1)/L(1,1);

for i =2:n

x(i) = (b(i) - L(i,1:i-1)*x(1:i-1))/L(i,i);

end

Emedan ber¨akningen av xi kr¨aver ca 2i flops, s˚a beh¨over hela processen 2(1 + 2 + · · · n) ≈ n² flops.

Denna algoritm ¨ar radorienterad, men vi kan ocks˚a konstruera en kolumnorienterad algoritm, som bygger p˚a saxpy operationen. Vi skall visa detta utg˚aende fr˚an ett exempel med tre ekvationer:

2

41 0 0

2 3 0

4 5 6

3 5

2 4x₁

x₂ x₃

3 5 =

2 47

8 9

3 5

Efter att vi l¨ost x₁ = 7 ˚aterst˚ar ett system med tv˚a obekanta:

3 0

5 6

 x₂ x₃



=

8 9



− 7

2 4



=

 −6

−19

 ,

varav f¨oljer x₂ = −2, varp˚a systemet reduceras till6x₃ = −19+5·2, som ger x₃ = −3/2. I allm¨anhet l¨oser vi i avseende p˚a x_j vid det j:te steget, och eliminerar d¨arp˚a denna obekant fr˚an ekvationerna med ordningsnumret j + 1 t.o.m. n. F¨orst f˚ar vi x₁ = b₁/`₁₁ varefter de ¨ovriga ekvationerna transformeras till 2

66 4

`₂₂ 0 · · · 0

`<sub>32</sub> `₃₃ · · · 0 ... ... . . . ...

`<sub>n2</sub> `_n3 · · · `_nn 3 77 5

2 66 4

x₂ x₃ ...

x_n 3 77 5 =

2 66 4

b₂ − x1`<sub>21</sub> b<sub>3</sub> − x1`₃₁

...

b_n − x₁`_n1 3 77

5 = b(2 : n) − x¹L(2 : n, 1).

I allm¨anhet ber¨aknar man i det j:te steget f¨orst x_j = b_j/`_jj, och utf¨or d¨arp˚a saxpy operationen b(j + 1 : n) ← b(j + 1 : n) − x_jL(j + 1 : n, j).

Genom kombination f˚ar vi slutligen MATLAB–funktionsrutinen function x=ltri(L,b)

% Invariabler:

% L: undre triangul¨ar, icke-singul¨ar nxn matris

% b: nx1 vektor

% Utvariabel:

% x: Lx = b

n = length(b);

x = zeros(n,1);

for j=1:n-1

x(j) = b(j)/L(j,j);

b(j+1:n) = b(j+1:n) - x(j)*L(j+1:n,j);

end

x(n) = b(n)/L(n,n);

I denna kolumnorienterade version av programmet utf¨ors ocks˚a n² flops, liksom i den radorienterade versionen.

En ¨ovre triangul¨ar matris kan behandlas p˚a ett liknande s¨att. Den enda skillnaden ¨ar att de obekanta ber¨aknas i motsatt ordningsf¨oljd. F¨or att l¨osa t.ex. systemet

2

4u₁₁ u₁₂ u₁₃ 0 u₂₂ u₂₃

0 0 u₃₃

3 5

2 4x₁

x₂ x₃

3 5 =

2 4b₁

b₂ b₃

3 5

behandlar vi det nerifr˚an upp:

x₃ = b₃/u₃₃

x₂ = (b₂ − u23x₃)/u₂₂

x₁ = (b₁ − u12x₂ − u13x₃)/u₁₁. F¨or ett allm¨ant fall kan programmet skrivas

x(n) = b(n)/U(n,n);

for i=n-1:-1:1

x(i) = (b(i) - U(i,i+1:n)*x(i+1:n))/U(i,i);

end

Liksom f¨or den undre triangul¨ara matrisen, kan vi ocks˚a i detta fall skriva en kolumnorienterad variant av programmet, som anv¨ander saxpy operationen:

function x=utri(U,b)

% Invariabler:

% U: ¨ovre triangul¨ar, icke-singul¨ar nxn matris

% b: nx1 vektor

% Utvariabel:

% x: Ux = b

n = length(b);

x = zeros(n,1);

for j=n:-1:2

x(j) = b(j)/U(j,j);

b(1:j-1) = b(1:j-1) - x(j)*U(1:j-1,j);

end

x(1) = b(1)/U(1,1);

Denna algoritm brukar man kalla f¨or bak˚atsubstitution.

En allm¨an matris kan man skriva som en produkt av en undre och en ¨ovre triangul¨ar matris: A = LU . Detta tillg˚ar med hj¨alp av Gauss elimineringsmetod, d¨ar man systematiskt eliminerar de obekanta.

F¨or att f¨orklara den skall vi f¨orst studera ett enkelt exempel med tre obekanta:

2x₁ − x₂ + 3x₃ = 13

−4x₁ + 6x₂ − 5x₃ = −28

6x₁ + 13x₂ + 16x₃ = 37

Om vi multiplicerar den f¨orsta ekvationen med −4/2 = −2, och subtraherar den fr˚an den andra ekva- tionen, s˚a kan vi eliminera x₁ fr˚an den. Likas˚a kan vi eliminera x₁ fr˚an den tredje ekvationen genom att subtrahera fr˚an den 6/2 = 3 g˚anger den f¨orsta ekvationen. Vi f˚ar d˚a ett reducerat ekvationssystem

2x₁ − x₂ + 3x₃ = 13

4x₂ + x₃ = −2

16x₂ + 7x₃ = −2

Sedan multiplicerar vi den andra ekvationen i det nya systemet med 16/4 = 4, och subtraherar den fr˚an den tredje, varigenom x₂ elimineras:

2x₁ − x₂ + 3x₃ = 13

4x₂ + x₃ = −2

3x₃ = 6

Denna elimineringsprocedur kommer allts˚a att f¨orvandla det ursprungliga kvadratiska ekvationssystemet till ett ekvivalent ¨ovre triangul¨art system, som har samma l¨osning. Enligt den metod f¨or l¨osning av triangul¨ara system som beskrevs tidigare, f˚ar vi

x₃ = 6/3 = 2

x₂ = (−2 − x₃)/4 = −1

x₁ = (13 − 3x₃ + x₂)/2 = 3

Denna beskrivning av Gauss elimineringsprocedur ¨ar enklast att uttrycka i matrisform. Med Gauss metod finner man allts˚a en undre triangul¨ar matris L och en ¨ovre triangul¨ar matris U , s˚a att A = LU . I v˚art exempel kan uppdelningen skrivas

A = 2

4 2 −1 3

−4 6 −5

6 13 16

3 5 =

2

4 1 0 0

−2 1 0

3 4 1

3 5

2

42 −1 3

0 4 1

0 0 3

3

5 ≡ LU.

Observera, att elementen under diagonalen i L ¨ar de faktorer, som anv¨ants under elimineringsprocessen.

Diagonalelementen i L ¨ar alla 1.

S¨attet att uppdela A brukar ¨aven kallas LU –faktorisering. S˚a snart man k¨anner denna uppdelning av A, kan man finna l¨osningen till det ursprungliga ekvationssystemet Ax = LU x = b genom att f¨orst l¨osa Ly = b och sedan U x = y.

Denna eleganta elimineringsprocess fungerar dock inte alltid. Om vi t.ex. ¨andrar v¨ardet av koefficienten f¨or x₁ fr˚an 2 till 0 i det ursprungliga ekvationssystemet, s˚a f˚ar vi alldeles tydligt ett system, d¨ar eliminerings- processen inte kan starta, eftersom vi inte kan eliminera x₁. Men detta problem kan man kringg˚a, som vi senare skall se.

Vi skall nu till¨ampa Gauss elimineringsprocess p˚a ett allm¨ant system med n ekvationer. Vi ser att vi kan eliminera x₁ fr˚an den i:te ekvationen (i > 1) genom att fr˚an denna ekvation subtrahera a_i1/a₁₁ g˚anger den f¨orsta ekvationen. Om man till¨ampar detta p˚a i = 2, . . . , n s˚a kommer alla element i den f¨orsta kolumnen i matrisen A att vara 0, utom det f¨orsta. Sedan eliminerar vi x₂ fr˚an den i:te raden (i > 2) genom att multiplicera den andra ekvationen i den modifierade matrisen A med a_i2/a₂₂ och subtrahera den fr˚an den i:te raden. Om detta till¨ampas p˚a i = 3, . . . , n, s˚a kommer man att f˚a en matris, d¨ar elementen i andra kolumnen ¨ar 0, utom de tv˚a f¨orsta, och s˚a vidare. M¨onstret ¨ar nu ganska klart:

for k = 1:n-1

Ber¨akna faktorerna som beh¨ovs f¨or att eliminera x(k) fr˚an ekvationerna med index k+1 till n och lagra dem i v(k+1:n).

Uppdatera ekvationerna k+1 till n.

end

Faktorerna ber¨aknas p˚a f¨oljande s¨att:

for i=k+1:n

v(i) = A(i,k)/A(k,k) end

s˚a att allts˚a v(k+1:n) = A(k+1:n,k)/A(k,k). Multiplikationen av den k:te raden med v(i) och sub- traktionen fr˚an den i:te raden kan utf¨oras med kommandot A(i,k:n) = A(i,k:n) - v(i)*A(k,k:n).

Kolumnerna startar fr˚an k eftersom de k − 1 f¨orsta elementen b˚ade i raderna k och i ¨ar noll. Dessa id´eer implementeras i f¨oljande triangulariseringsprocedur:

for k=1:n-1

v(k+1:n) = A(k+1:n,k)/A(k,k);

for i=k+1:n

A(i,k+1:n) = A(i,k+1:n) - v(i)*A(k,k+1:n);

end end

U = triu(A)

triu ¨ar en MATLAB–funktion, som i detta fall bildar en ¨ovre triangular matris av A:s element, dvs U_ij = A_ij, om i ≤ j, men U_ij = 0, om i > j.

I i-slingan utf¨ors ett antal radorienterade saxpy operationer. Om t.ex. n = 6 och k = 3, s˚a kan de tre saxpy operationerna

A(4, 4 : 6) ← A(4, 4 : 6) − v(4) ∗ A(3, 4 : 6) A(5, 4 : 6) ← A(5, 4 : 6) − v(5) ∗ A(3, 4 : 6) A(6, 4 : 6) ← A(6, 4 : 6) − v(6) ∗ A(3, 4 : 6) kombineras till vektorekvationen

A(4 : 6, 4 : 6) = 2

4A(4, 4 : 6) A(5, 4 : 6) A(6, 4 : 6)

3 5 ←

2

4A(4, 4 : 6) A(5, 4 : 6) A(6, 4 : 6)

3 5 −

2 4v₄

v₅ v₆

3

5 A(3, 4 : 6).

D¨arf¨or kan vi ers¨atta i-slingan i programmet med en yttre produkt:

for k=1:n-1

v(k+1:n) = A(k+1:n,k)/A(k,k);

A(k+1:n,k+1:n) = A(k+1:n,k+1:n) - v(k+1:n)*A(k,k+1:n);

end

U = triu(A);

Det ˚aterst˚ar f¨or oss att ber¨akna den undre triangul¨ara matrisen L. Vi kan visa, att den uppbyggs av faktorerna v_i. Om vi p˚a nytt studerar fallet d˚a n = 6 och k = 3, s˚a finner vi att matrisen A blir multiplicerad fr˚an v¨anster med matrisen

M₃ = 2 66 66 66 64

1 0 0 0 0 0

0 1 0 0 0 0

0 0 1 0 0 0

0 0 −v₄ 1 0 0

0 0 −v₅ 0 1 0

0 0 −v₆ 0 0 1

3 77 77 77 75

som allts˚a har ettor i diagonalen, och faktorerna i k:te kolumnen under diagonalen. Efter n − 1 steg har matrisen M_n−1 · · · M2M₁A = U blivit en ¨ovre triangul¨ar matris, och s˚aledes f˚ar den ursprungliga matrisen A formen

A = (M₁⁻¹M₂⁻¹ · · · M_n−1⁻¹ )U.

Inversen av en faktormatris ¨ar mycket enkel. S˚alunda f˚ar vi t.ex.

M₃⁻¹ = 2 66 66 66 64

1 0 0 0 0 0

0 1 0 0 0 0

0 0 1 0 0 0

0 0 v₄ 1 0 0

0 0 v₅ 0 1 0

0 0 v₆ 0 0 1

3 77 77 77 75

(kontrollera t.ex. genom att ber¨akna M3M₃⁻¹). Man kan ocks˚a visa, att L = M₁⁻¹M₂⁻¹ · · · M_n−1⁻¹

¨ar en undre triangul¨ar matris sum uppfyller villkoret att L(:, k) ¨ar den k:te kolumnen av M_k⁻¹. I fallet n = 6 f˚ar vi t.ex.

L = M₁⁻¹M₂⁻¹M₃⁻¹M₄⁻¹M₅⁻¹ = 2 66 66 66 66 4

1 0 0 0 0 0

v₂⁽¹⁾ 1 0 0 0 0

v₃⁽¹⁾ v₃⁽²⁾ 1 0 0 0 v₄⁽¹⁾ v₄⁽²⁾ v⁽³⁾₄ 1 0 0 v₅⁽¹⁾ v₅⁽²⁾ v⁽³⁾₅ v₅⁽⁴⁾ 1 0 v₆⁽¹⁾ v₆⁽²⁾ v⁽³⁾₆ v₆⁽⁴⁾ v₆⁽⁵⁾ 1

3 77 77 77 77 5

d¨ar ¨ovre index anger ordningsnumret f¨or det steg, d¨ar faktorn ing˚ar. Faktorerna kan allts˚a lagras i de positioner i matrisen, som ¨ar avsedda att bli nollst¨allda. S˚alunda kommer t.ex. v₅⁽²⁾ att nollst¨alla a52

under det andra steget, och kan d¨arf¨or lagras i positionen (5, 2). Vi kan d¨arf¨or implementera Gauss elimineringsmetod p˚a f¨oljande s¨att:

function [L,U] = gauss(A)

% Invariabel:

% A: nxn matris

% Utvariabler:

% L: nxn undre triangul¨ar matris

% U: nxn ¨ovre triangul¨ar matris

%

[n,n] = size(A);

for k=1:n-1

A(k+1:n,k) = A(k+1:n,k)/A(k,k);

A(k+1:n,k+1:n) = A(k+1:n,k+1:n) - A(k+1:n,k)*A(k,k+1:n);

end

L = eye(n,n) + tril(A,-1);

U = triu(A);

Detta program anropar tril(A,-1), en MATLAB–funktion, som konstruerar en undre triangul¨ar matris vars element ¨ar Aij d˚a j ≤ i − 1, men 0 d˚a j > i − 1.

Ett linj¨art ekvationssystem Ax = b kan d¨arp˚a l¨osas p˚a f¨oljande s¨att

[L,U] = gauss(A);

y = ltri(L,b);

x = utri(U,y);

om ingen division med noll intr¨affar under Gauss’ elimineringsprocess. En matris beh¨over inte alltid ha en LU –uppdelning. Se t.ex. p˚a 

0 1

1 1



=

 1 0

`₂₁ 1

 u₁₁ u₁₂ 0 u₂₂



F¨or att identitet skall g¨alla i position (1, 1) b¨or u₁₁ = 0. Men d˚a kan identitet inte g¨alla i position (2, 1), eftersom `₂₁u₁₁ borde vara 1. Som vi redan n¨amnt, kan a11 = 0 ocks˚a g¨alla, t.ex. i

0 1

1 1

 x₁ x₂



=

1 2



som har l¨osningen x = [1 1]⁰.

Men vi kan ocks˚a f˚a problem, om faktorerna ¨ar stora. Ett enkelt s¨att att l¨osa problemen ¨ar att kasta om rader. Detta tillg˚ar s˚a att man f¨orst s¨oker upp det st¨orsta elementet i kolumnvektorn A(k : n, k). Antag, att detta ¨ar Aqk (pivotelementet), och allts˚a befinner sig i q:te raden. D¨arp˚a kastar man om raderna k och q, innan man dividerar med A_kk. N¨ar man l¨oser ekvationssystemet med pivotering, b¨or man ocks˚a minnas att permutera raderna i b. Denna metod stabiliserar algoritmen, och minimerar ¨aven avrundningsfelet.

1.5. Minsta kvadratmetoden

Den f¨orsta praktiska till¨ampningen av minsta kvadratmetoden som blivit k¨and ¨ar Gauss’ ber¨akning av sm˚a- planeten Ceres’ bana. Denna sm˚aplanet, som var den f¨orsta som uppt¨acktes, s˚ags av Giuseppe Piazzi¹ f¨orsta g˚angen ny˚arsnatten 1801. Han f¨oljde den under 40 dagars tid, men tappade sedan bort den, n¨ar den kom alltf¨or n¨ara solen. Ingen s¨aker efemerid hade hunnit ber¨aknas, men med hj¨alp av minsta kvadratmetoden kunde Gauss (som d˚a var 24 ˚ar gammal) ber¨akna planetens banelement p˚a basen av endast tre observa- tioner, s˚a att den kunde ˚aterfinnas vid ˚arets slut. Metoden beskrevs i detalj i Gauss’ bok: ”Theoria Motus Corporum Coelestium in sectionibus conicis solem ambientium” (en teori f¨or r¨orelsen hos de himlakroppar, som r¨or sig i k¨agelsnitt runt solen), som utkom ˚ar 1809. Han s¨ager d¨ar sj¨alv, att han anv¨ant metoden s˚a tidigt som 1795. Samma metod hade redan tidigare publicerats av Legendre ˚ar 1806 i en avhandling om best¨amningen av kometbanor, d¨ar han bifogat ett till¨agg ”sur la m´ethode des moindres carr´es”. Det ¨ar allts˚a egentligen Legendre som givit metoden dess namn.

Antag, att vi ¨onskar anpassa en line¨ar modell till givna (fysikaliska) m¨atv¨arden. Vanligen ¨ar antalet m¨at- ningar st¨orre ¨an antalet obekanta. Det g¨aller d˚a att l¨osa ett s k ¨overdeterminerat system, som i matrisform kan uttryckas Ax = b.

1italiensk astronom (1746-1826), som grundade ett observatorium i Palermo och publicerade tv˚a stj¨arnkataloger.

Problem: A ¨ar en given n × m matris, d¨ar n (= antalet observationer) ≥ m (= antalet obekanta), och b ¨ar en kolumnvektor med n element. Best¨am vektorn x (som har m element) s˚a, att Ax ¨ar den ”b¨asta”

approximationen till b.

Emedan ifr˚agavarande ekvationssystem ¨ar ¨overdeterminerat, kan det inte l¨osas exakt, och man kan d¨arf¨or inte entydigt definiera en ”b¨asta” l¨osning. En metod, som ¨ar statistiskt motiverad, och dessutom leder till relativt enkla r¨akningar, ¨ar den minsta kvadratmetoden.

L˚at oss till att b¨orja med anta, att vi vill anpassa n datapunkter (x_i, y_i) till en linj¨ar modellfunktion y = ax + b. Vi kan g¨ora detta genom att minimera summan av kvadraterna p˚a observationspunkternas avst˚and fr˚an den r¨ata linjen (godhetsfunktionen): S(a, b) =

Xn i=1

(y_i − ax_i − b)².

I minimum b¨or de partiella derivatorna av godhetsfunktionen S(a, b) i avseende p˚a parametrarna a och b f¨orsvinna:

0 = ∂S(a, b)

∂a = −2

Xn i=1

x_i(y_i − ax_i − b)

0 = ∂S(a, b)

∂b = −2

Xn i=1

(y_i − ax_i − b)

Dessa villkor leder till ett ekvationssystem med tv˚a obekanta (normalekvationerna), som kan skrivas:

[xx]a + [x]b = [xy]

[x]a + nb = [y],

d¨ar klamrarna ¨ar beteckningar, som inf¨ordes av Gauss: [x] = P

i x_i, [xx] = P

i x²_i, [xy] = P

i x_iy_i etc.

I det allm¨anna fallet modellerar vi de n observerade punkterna med en line¨ar funktion y(x) = Pm

k=1a_kX_k(x), d¨ar X_k(x) ¨ar godtyckliga funktioner av x (t.ex. potenser av x, vid polynomapproxima- tion), och a_k ¨ar parametrar, som skall best¨ammas. Vi f˚ar d˚a ett system av n ekvationer med m obekanta (tillst˚andsekvationerna):

Xm k=1

a_kX_k(x_i) = y_i, i = 1, . . . , n

Koefficientmatrisens element ¨ar allts˚a i detta fall A_ik = X_k(x_i). I det allm¨anna fallet kan vi anv¨anda godhetsfunktionen

S(a) = Xn

i=1

w_i

"

y_i − Xm

k=1

a_kX_k(x_i)

#2

,

d¨ar w_i betecknar observationernas vikter (om de inte alla ¨ar uppm¨atta med samma noggrannhet), och a = [a₁, a₂, . . . , a_m] ¨ar parametervektorn.

Eftersom ^∂S

∂ak = 0 g¨aller i minimet, s˚a f˚ar vi normalekvationerna ur villkoren Xn

i=1

w_iX_k(x_i) 2

4y_i − Xm

j=1

a_jX_j(x_i) 3

5 = 0, k = 1, . . . , m.

Detta ekvationssystem kan skrivas i formen (efter omkastning av summeringsordningen) Xm

j=1

a_j Xn

i=1

w_iX_j(x_i)X_k(x_i) =

Xn i=1

w_iy_iX_k(x_i).

Med hj¨alp av koefficientmatrisen A, som definierades ovan, kan detta ekvationssystem ocks˚a skrivas i den enklare formen

A⁰W Aa = A⁰W y,

dvs normalekvationerna f˚as direkt av tillst˚andsekvationerna Aa = y genom att multiplicera fr˚an v¨anster med den transponerade koefficientmatrisen A⁰ och en diagonal viktsmatris W (W_ii = w_i, W_ij = 0, i 6= j).

Som vi redan sett, ¨ar det m¨ojligt att l¨osa ett s˚adant ¨overdeterminerat system i MATLAB med operatorn

’\’. Vi skall studera n˚agra exempel p˚a anv¨andningen av denna operator, och b¨orja med ett ekvationssystem, som studerades av Gauss²: Antag, att vi ¨onskar l¨osa f¨oljande ¨overdeterminerade system:

8>

><

>>

:

p −q +2r −3 = 0

3p +2q −5r −5 = 0

4p +q +4r −21 = 0

−p +3q +3r −14 = 0

Om endast de tre f¨orsta ekvationerna vore givna, skulle vi f˚a p = 18/7, q = 23/7 och r = 13/7, och problemet vore d¨armed l¨ost. Dessa v¨arden satisfierar dock inte den fj¨arde ekvationen, vars v¨anstra membrum blir −8/7 ist¨allet f¨or 0 vid substitutionen. Normalekvationerna j¨amte sina l¨osningar blir i detta

fall 8

<

:

27p +6q = 88

6p +15q +r = 70

q +54r = 107

8<

:

p = 49154/19899 ≈ 2.47017

q = 2617/737 ≈ 3.55088

r = 12707/6633 ≈ 1.91572 Med MATLAB kan problemet l¨osas p˚a f¨oljande s¨att:

2Gauss: Theoria Motus ..., Art. 184

>> A=[1 -1 2; 3 2 -5; 4 1 4; -1 3 3];

>> b =[3; 5; 21; 14];

>> A\b ans =

2.4702 3.5509 1.9157

Ett annat exempel behandlar approximation av kvadratrotsfunktionen f (x) = √

x inom ett litet intervall, t.ex. [0.25, 1] med en r¨at linje y(x) = α + βx. Funktionen som minimeras ¨ar i detta fall

S_n(α, β) =

Xn i=1

[(α + βx_i) − √

x_i]².

Funktionen S kan ocks˚a uttryckas i matrisform: S = r⁰r, d¨ar r ¨ar kolumnvektorn

r = 2 66 4

1 x₁ 1 x₂ ... ...

1 x_n 3 77 5

α β



− 2 66 4

y₁ y₂ ...

y_n 3 77 5 ,

och y_i = √ x_i.

Nedan visas ett MATLAB–program, som g¨or en minsta kvadratanpassning f¨or n = 2 och n = 100, och visar kvadratrotsfunktionen, samt den linj¨ara approximationen i en graf.

% Program: lsfit

% Visar minstakvadratanpassning av f(x)=sqrt(x) f¨or 0.25<=x<=1 close

x1 = linspace(0.25,1);

y1 = sqrt(x1);

for n = [2 100]

x = linspace(0.25,1,n)’;

A = [ones(n,1) x];

b = sqrt(x);

xls = A\b;

alfa = xls(1);

beta = xls(2);

figure

plot(x1,y1,’y-’,x1,alfa+beta*x1,’g--’)

title(sprintf(’n = %2.0f, alfa = %10.6f, beta = %10.6f’,n,alfa,beta)) end

F¨or n = 2 f˚ar vi en r¨at linje mellan punkterna x = 0.25 och x = 1, men n = 100 ger en b¨attre approximation f¨or kvadratroten, som synes av nedanst˚aende graf.

Om n ¨ar stort, kan man approximera minsta kvadratsumman med en integral:

0.75 n

Xn i=1

[(α + βx_i) − √

x_i]² ≈ Z 1

0.25

[(α + βx) − √

x]²dx ≡ S_∞(α, β),

d¨ar xi = 0.25 + 0.75(i − 1)/(n − 1) d˚a i = 1 : n.

D˚a n → ∞ kommer funktionen S_n(α, β) att konvergera mot S_∞(α, β). Av ekvationerna

1 2

∂S_∞

∂α = α Z 1

0.25

dx + β Z 1

0.25

xdx − Z 1

0.25

√xdx = 0

1 2

∂S_∞

∂β = α Z 1

0.25

xdx + β Z 1

0.25

x²dx − Z 1

0.25

x^3/2dx = 0 f¨oljer d˚a det linj¨ara systemet

 3/4 15/32

15/32 21/64

 α_∗ β_∗



=

 7/12 31/80

 .

L¨osningen ¨ar 

α_∗ β_∗



=

0.370370370 0.651851851

 .

Som vi n¨amnt, kan funktionerna Xk(x) i det allm¨anna fallet vara godtyckliga funktioner av x, t.ex.

potenser. I detta fall ¨ar modellfunktionen ett polynom, och l¨osningen kan utf¨oras med MATLAB–funktionen polyfit.