Något om Minsta kvadratmetoden och Mathematica

Något om Minsta kvadratmetoden och Mathematica

Bertil Nilsson 2021-08-15

x1 xi xn

x y1

yi

yn

y

xi,yi

xi, fxi,͍

fi

1 2 3 4 5 6 7 x

1 2 3 4 5 y

Extrapol Interpolation Extrapol

min _{i 1} ⁿ f x _i , ͍ y _i ² min f x, ͍ y x ² x

5 10 15 20 25ts

0 1 2 3 4 5 6 ci,ctmgl

1 2 3 4 x

2 3 4 5 y

ť Förord

På följande sidor presenteras en elementär “streetwise guide” till minsta kvadratmetoden med flitig användning av Mathematica.

Framställningen är fåordig, fri från pedanteri men i någon mening fullständig. Det man väsentligen behöver veta om begrepp, terminologi, beteckningar och teori för att modellera och lösa problem i framtida kurser och yrkesliv som ingenjör, naturvetare eller lärare klarläggs och typiska exempel ges.

ť Inledning

Vi börjar med en problemställning från verkliga livet. Anders på Andersberg har köpt en ny bil med konstantfarthållare. För att undvika framtida bekymmer med lagens väktare vill nu Anders kontrollera riktigheten i utrustningen. För detta ändamål uppsöker han en lång raksträcka där han vet att avståndet mellan lyktstolparna är 100 m. Med konstantfarthållaren inställd på 50 km/h bär det iväg. Eftersom han bara har ett vanligt armbandsur till hands så blir det till att uppskatta restiderna på sekunden när. Han får följande serie från en stolpe till de fem efterföljande

tid s 0 7 14 22 29 36

Anders strategi är nu att äntligen dra nytta av den gamla mekanikformeln s vt som han tvingades exercera under skoltiden. Han börjar med att lägga in mätdata i ett diagram

data 0 7 14 22 29 36

0 100 200 300 400 500 ;

ListPlot data, PlotStyle Red, PointSize 0.03 , AxesLabel "t s ", "s m "

5 10 15 20 25 30 35 ts 100

200 300 400 500

sm

Han förstår att s vt är en rät linje genom origo, men hur ska den luta? Entusiastiskt låter han den gå genom sista mätpunkten, det vill säga han väljer modellen s vt ⁵⁰⁰₃₆t och lägger till den i diagrammet

Plot500 36

t, t, 0, 40 , PlotStyle Blue, AxesLabel "t s ", "s m " , Epilog Red, PointSize 0.03 , Point data 

10 20 30 40 ts

100 200 300 400 500

sm

Nja, riktigt bra blev det inte. Genast dyker det upp några frågeställningar. Det kanske finns någon lutning som gör att fler punkter kommer ännu närmare linjen? Å varför är första mätpunkten 0, 0 helig? Det borde ju vara lika svårt att “pricka” första stolpen som de andra? I så fall kanske han skulle prova s vt m. Vi ska hjälpa honom att utreda detta.

ť Minsta kvadratmetoden

Antag att vi har givet mätdata och vill anpassa en funktion till dessa “så bra som möjligt”. Låt denna funktion allmänt tecknas y f x, , där x, y är den oberoende respektive beroende variabeln och en vektor med parametrar som ska bestämmas så anpass- ningen blir så bra som möjligt. Naturligtvis kan vi ha flera oberoende variabler, t.ex. T f , f x, y, z, t , som modell för hur temperaturen varierar med tiden t på olika platser i rummet x, y, z . Vi kallar f för modellen och för modellparametrarna.

Exempelvis

y f x, k x m k, m y f x, c0 c1x c2x² c0, c1, c2

y f x, a ^bx a, b y f t, acosΩt bsinΩt a, b,Ω

Det typiska är att vi måste välja modellen själva, se figur nedan där tre olika modeller är anpassade till samma mätdata. Minsta kvadratmetoden (MKM) hjälper sedan till med att bestämma parametrarna . Vi får alltså olika anpassningar beroende på vilka val vi gör. Detta är viktigt att förstå! Ofta är det så att man vet vilka funktioner som ska ingå i f. Kanske har man en matematisk modell där man behöver trimma några parametrar mot försök.

I Minsta kvadratmetoden måste man välja modellen själv och man får olika anpassningar beroende på vilka val man gör!

Det är viktigt att man väljer funktioner i modellen så att den på ett kvalitativt sätt fångar det man vet om systemet som studeras, t.ex.

y då x 0 eller y 0 då x . I figuren nedan ser vi också att det är mycket vanskligt att använda modellen utanför mätom- rådet. Detta kallas extrapolation till skillnad mot interpolation då man befinner sig i det område som omsluts av mätområdet.

Notera också i figuren att med ett komplicerat modellval så kan man lätt få den att gå genom alla mätpunkter, exempelvis ett polynom av grad n definieras entydigt av n 1 mätpunkter. Men detta kan vara en mycket dålig representation av den bakomlig- gande fysikaliska modellen, som kanske bättre speglas av de två “mindre nervösa” modellerna, speciellt om man tänkt sig tillåta

“lite” extrapolation. Spridningen i mätdata kan helt enkelt bero på svårigheter vid mätningen. Den kraftiga översvängningen vid ändpunkterna av mätområdet som uppstår vid approximation med polynom av höga gradtal brukar kallas Runges fenomen, (Carl Runge, 1856-1927, tysk matematiker), se den bruna kurvan.

1 2 3 4 5 6 7 x

1 2 3 4 5 y

Extrapol Interpolation Extrapol

y f x, kx m y f x, ax² bx c

y f x, c6x⁶ c5x⁵ c1x c0

Har man möjlighet att inspektera såväl mätpunkter som anpassad modell i samma graf ska detta alltid göras. Ögat är nämligen väldigt bra på att hitta en vanlig notorisk felkälla, så kallade uteliggare, vilka ofta förklaras av mätfel. Man får då olika modeller beroende på om dessa tas med eller inte. Exempelvis Anders igen med ett duktigt bokföringsfel vid en av stolparna.

10 20 30 40 ts

100 200 300 400 500 600

sm

Anders mäter ok Anders och en uteliggare

Minsta kvadratmetoden vilar på gamla kända kunskaper, nämligen minimering av en i detta fall kvadratisk funktion. Därav namnet.

Antag att vi har en uppsättning mätdata med n st mätpunkter, se fig nedan från vänster till höger. Bilda nu vid varje mätpunkt xi, yi

felet i mellan modell f xi, och det uppmätta värdet yi, det vill säga i f xi, yi. Nu är det lika illa om mätpunkten ligger över som under modellen. Välj därför istället att studera _i² och lägg samman dem alla.

x1 y1

x2 y2

xn yn

x1 xi xn

x y1

yi

yn

y

xi,yi

xi, fxi,͍

fi

S _{i 1}ⁿ _i² ₁² ₂² _n² _{i 1}ⁿ f xi, yi2

En kvadratsumma är alltid icke-negativ, det vill säga S 0, och i det ideala fallet S 0 går modellen genom alla mätpunkter.

Eftersom alla mätpunkter är applicerade blir S en funktion av parametrarna, det vill säga S . Minsta kvadratmetoden går nu ut på att minimera S med avseende på parametrarna .

min S

Man säger att f är anpassad i minsta kvadratmetodens mening. Vi drar oss nu till minnes att söka minimum till en funktion är det samma som att söka nollställe till derivatan. Eftersom vi har flera variabler i måste vi simultant söka nollställe till respektive derivata. Att derivera en funktion av flera variabler med avseende på någon speciell variabel kallas för partiell derivata och skrivs med istället för . Inget märkvärdigt, man gör som vanligt, det vill säga deriverar med avseende på den variabel som önskas och låter de andra vara konstanter vilka som helst. Vi tar ett

Exempel: Låt funktionen g x, y x³ 5y 4xy² xsin 2 y 7 vara given. Bestäm de partiella derivatorna ^g_x och ^g_y. Lösningsförslag: Vi får direkt ^g_x 3x² 4y² sin 2 y och ^g_y 5 8x y 2xcos 2 y .

Så att söka minimum av S, min S , blir då att lösa ekvationssystemet ^S

S p₁ 0

S p₂ 0

Alltså ett ekvationssystem med lika många ekvationer som vi har modellparametrar p1, p2, . Om detta blir linjärt har vi linjär minsta kvadratmetod annars olinjär. En linjär modell kan alltid skrivas som en skalärprodukt mellan ingående funktioner, samlade i en vektor, och modellparametrarna, det vill säga f x, f1, f2, . Exempelvis har vi modeller som leder till linjär MKM

y f x, kx x , k y f x, c0 c1x c2x² 1, x, x² , c0, c1, c2

y f x, kx m x, 1 , k, m y f x, ax bsin x x, sin x , a, b

Vid linjär MKM är koefficientmatrisen i ^S både kvadratisk och symmetrisk, och S har endast ett lokalt minimum, som då också är globalt, med entydig lösning som följd. Tryggt! I detta fall finns till vår hjälp i Mathematica en lättanvänd funktion som heter just

“anpassa”, Fit[data,f,x]. Mata med mätdata i form av en matris med två kolonner xi; yi då vi har en oberoende variabel, eller i vid flera, samt och oberoende variabler. Den levererar sedan f på direkt användbar form .

Låt oss nu använda det typiska exemplet y f x, kx m, med modellparametrarna k, m , som modell och arbeta igenom processen. Först summan av alla kvadratfel

S _{i 1}ⁿ k xi m yi2

Med inre derivata i färskt minne får vi så ekvationssystemet, med kvadratisk och symmetrisk koefficientmatris, som bestämmer .

min S mink,mS ^S

S

k i 1n 2 k xi m yi xi 0

S

m i 1n 2 k xi m yi 0

i 1n x_i² _{i 1}ⁿ xi i 1n xi n

k m

i 1n xiyi

i 1n yi (1)

Exempel: Studera Anders problem.

Lösningsförslag: Först den modell som går genom origo, det vill säga s vt t v . Vi får S ⁶_{i 1}vti si2 S

v i 16 2 vti si ti 0 v _{i 1}⁶ t_i² _{i 1}⁶ tisi

där

i 16 t_i² 0² 7² 14² 22² 29² 36² 2866

i 16 tisi 0 0 7 100 14 200 22 300 29 400 36 500 39 700 v ^{39 700}₂₈₆₆ 13.8521 Provkör Fit[data,f,x] på Anders problem. Med s vt t v får vi direkt

s1 Fit data, t , t 13.8521 t

Vi känner igen resultatet från handräkningen ovan. Sedan över till hans fundering s vt m t, 1 v, m . Vi startar direkt med slutresultatet av minimeringen (1) ovan. Kontrollräkna gärna!

i 16 t_i² _{i 1}⁶ ti i 16 ti n

v m

i 16 tisi i 16 si

2866 108

108 6

v m

39 700 1500

v m

1 461

6350 950

13.7744 2.06074 Visst, Mathematica instämmer

s2 Fit data, t, 1 , t 13.7744 t 2.06074

En liten jämförelse visar att han har mätt väl. Modellerna är inte lika och i detta fall är skillnaden nästan av akademiskt intresse, men indikerar ändå det viktiga att olika funktionsval ger olika modeller. Vi avslutar med att visualisera data och båda modellerna i samma diagram. Detta återkommer i tid och otid, bara så du vet !

Plot s1, s2 , t, 0, 40 , PlotStyle Blue, Orange , AxesLabel "t", "si,s1,s2" , Epilog Red, PointSize 0.03 , Point data

10 20 30 40 t

100 200 300 400 500

si,s1,s2

Om vi har en linjär modell kan vi också komma fram till det önskade ekvationssystemet genom att studera det överbestämda ekvationssystemet som i fallet y f x, k x m blir

x1 1 x2 1 xn 1

k m

y1

y2

yn

Detta överbestämda ekvationssystem har i allmänhet ingen lösning eftersom vi har fler ekvationer än obekanta. Man kan visa att det kan ges en lösning i minsta kvadratmetodens mening, och börjar med att övertyga oss genom en geometrisk betraktelse, se figur.

̾1

̾2

̾3

̾k

Ν͍

͏ Ν͍

͖

͏

Kolonnerna k i spänner upp ett plan och är då en vektor i planet med som komponenter i basen k. Eftersom inte löser ligger följaktligen inte i planet. Då är residualvektorn en vektor som går från spetsen på till spetsen på i planet. Minsta kvadratmetoden motsvarar då att söka minimum av . Vi kommer ihåg att närmsta vägen från en punkt till ett plan går vinkelrät mot planet. Detta innebär att vi har minimum av då är vinkelrät mot alla basvektorerna k. Eftersom vinkelräthet innebär att skalärprodukten är noll får vi till slut

1 0

2 0

Mera formellt har vi min ² min

2 2 .

Vi kallar  Ν ͖ för det utjämnade systemet eller normalekvationerna till det överbestämda ekvationssystemet ovan.

Koefficientmatrisen är alltid både kvadratisk och symmetrisk beroende på just konstruktionen . Man kan visa att normalekva- tionerna alltid har en lösning ¹ , där ¹ kallas pseudoinversen till . Naturligtvis bestäms inte på detta kostsamma sätt med inverser utan med Gausselimination direkt på normalekvationerna.

Om Ν͍  ͖ är överbestämt kallas Ν˞Ν͍  Ν˞͖ normalekvationerna och har entydig lösning ͍ i minsta kvadratmetodens mening, det vill säga ^ÙS

Ù͍  σ Ó Ν˞Ν͍  Ν˞͖.

Vi provar på Anders modell s vt m. Vi ser att koefficientmatrisen är både kvadratisk och symmetrisk som sig bör. Kontroll- räkna gärna matrisgodiset och ekvationslösningen som följer, vilken åter igen levererar parametrarna i MKM:s mening

, 1 & data All, 1 ; . .data All, 2 0 1

7 1 14 1 22 1 29 1 36 1

0 7 14 22 29 36

1 1 1 1 1 1

2866 108

108 6 39 700, 1500

NSolve . . v, m .data All, 2 m 2.06074, v 13.7744

Så här via normalekvationerna kan man alltid gå då man har linjär minsta kvadratmetod. Vid olinjär modell är man placerad i osäker terräng och har man väl hittat ett minimum kan det alltid ifrågasättas om det är ett globalt. I Mathematica finns till vår direkta hjälp FindFit,(N)Minimize och(N)Maximize som utger sig för att kunna leverera globalt optimum. De två sistnämnda kan även matas med linjära och olinjära bivillkor. Fördelen med FindFit är att man inte behöver formulera kvadratsumman explicit utan kan ange mätdata och den olinjära modellen som indata. Utöver dessa finns naturligtvis de gamla kända FindMinimum och FindMaximum som då i någon mening får anses vara något omsprungna av de tre tidigare nämnda.

Hittills har vi behandlat diskret minsta kvadratmetod eftersom vi anpassat en modell till ett ändligt antal mätpunkter. Om vi i någon mening låter antalet mätpunkter växa obegränsat får vi kontinuerlig minsta kvadratmetod.

min S min _{i 1}ⁿ f xi, yi2

n min f x, y x ² x

Som sinnebild kan man se detta som att man minimerar den målade arean som innestängs mellan (modell)funktionen f x, och funktionen y x i intervallet , fast inte riktigt eftersom vi kvadrerar integranden, så bättre, minimerar den rotationsvolym som uppstår då y x roterar runt f x, , sånär som ett Π kommer vi ihåg från rotationsvolymer. Detta kan komma till användning då man i ett intervall vill approximera en besvärlig funktion med någon enklare. Även här har vi de två fallen med linjär och olinjär MKM.

Exempel: Som illustration till kontinuerlig minsta kvadratmetod väljer vi att i intervallet 0, 1 approximera y ^x 0.1sin 10x med den betydligt enklare linjära modellen f x, ax² bx c x², x, 1 a, b, c , där x², x, 1 är funktionsvalen och a, b, c de sökta modellparametrarna.

Lösningsförslag: Med ett snabbt hack och ointressanta utdataceller undertryckta får vi så en inblick i skådespelet

a, b, c ; f x², x , 1. ; S 

0 1

f ^x 0.1 Sin 10 x ² x;

Plot ^x 0.1 Sin 10 x , f . Solve D S, 0 & First , x, 0, 1 , PlotStyle Red, Blue , AxesLabel "x", "y x ,f x, "

0.2 0.4 0.6 0.8 1.0 x 1.5

2.0 2.5 yx,fx,͍

Som läsaren säkert upptäckt nu så är minsta kvadratmetoden inget som räknas för hand. Eftersom den är mycket vanlig i praktiken och då kanske matas med flera tusen mätpunkter och modeller med många parametrar och komplicerade funktioner överlämnas gärna råräknandet till datorerna. Det viktiga är att man vet vad man gör innan man lämnar över till Mathematica som har effektiva lösare för alla förekommande minsta kvadratmodeller. Som avslutning ska man rita modell och mätdata i samma diagram och göra en visuell betraktelse av hur noggrant modellen ansluter till mätdata. Detta är inte helt lätt då antalet oberoende variabler i modellen är större än två, då får man studera olika modeller och deras kvadratfel. Med andra ord så ligger det lite “godtycke” med i metoden och vid komplicerade problem bör någon erfaren person konsulteras så man inte sprider illäror med en modell som är helt uppåt väggarna! Så

Se upp för uteliggare. Var noga med modellval, “man får den anpassning man förtjänar”.

Undvik interpolation i mätområdets utkanter och befatta dig inte med extrapolation!

Resten av skriften är en bukett exempel, de flesta utan bakgrund och med få mätpunkter för att hålla beräkningsvolymen på en beskedlig nivå. Utanför skolans värld krävs definitivt information om vilket system som det är mätt på och fler mätpunkter för att kunna göra ett välunderbyggt val av modell! Vi avslutar med några sådana exempel.

ť Blandade exempel för resten

Exempel: Använd MKM för att anpassa en rät linje p1x kx m, ett andragradspolynom p2 x ax² bx c och ett tredjegrads- polynom p3 x c3x³ c2x² c1x c0 till mätpunkterna

x 1 2 3.5 5.8 y 1.1 1.8 4 5.5

Lösningsförslag: Typisk råräkning. Börja med räta linjen och möblera det överbestämda ekvationssystemet x1 1

x2 1 xn 1

k m

y1

y2

yn

Bestäm sedan k och m enligt MKM receptet och normalekvationerna k

m . Vi tar lite matrisgodis på vägen

data 1 2 3.5 5.8

1.1 1.8 4 5.5 ; , 1 & data All, 1 ;

. .data All, 2 1 1

2 1 3.5 1 5.8 1

1 1 2 1 3.5 1 5.8 1

1 2 3.5 5.8

1 1 1 1

50.89 12.3

12.3 4. 50.6, 12.4

Slutligen de efterlängtade.

Solve . . k, m .data All, 2 k 0.954276, m 0.165602

Å så en liten bild

Plot k x m . , x, 0, 7 , PlotStyle Blue, AxesLabel "x", "yi,f x, " , Epilog Red, PointSize 0.03 , Point data

1 2 3 4 5 6 7 x

1 2 3 4 5 6 7 yi,fx,͍

Det går även bra att direkt använda funktionen Fit. Här väljer man vilka funktioner man vill använda, i retur kommer en linjärkom- bination av dessa i MKM:s mening.

p1 Fit data, 1, x , x 0.954276 x 0.165602

Med samma resultat som ovan. Vi provar nu med ett andragradspolynom, p2x ax² bx c. Detta går naturligtvis att lösa genom att möblera det överbestämda ekvationssystemet

x₁² x1 1 x₂² x2 1

x_n² xn 1 a b c

y1

y₂

yn

Först , sedan lite matrisgodis och slutligen parametrarna ur normalekvationerna .

 ², , 1 & data All, 1 ; . .data All, 2

1 1 1

4 2 1

12.25 3.5 1 33.64 5.8 1

1 4 12.25 33.64 1 2 3.5 5.8

1 1 1 1

1298.71 246.987 50.89 246.987 50.89 12.3

50.89 12.3 4.

242.32, 50.6, 12.4

Solve . . a, b, c .data All, 2 a 0.0735314, b 1.46352, c 0.464835

Man ska kunna genomföra kalkylen den hårda vägen via normalekvationerna ovan, men att för hand multiplicera ihop och lösa ett ekvationssystem med tre obekanta är direkt olämpligt i datorernas tidevarv! Det viktiga är att förstå hur det fungerar, sedan överläm- nar man gärna råräknandet! I praktiken hade man valt den lätta vägen med Fit[data,f,x]och istället lagt kraft på att prova några olika modeller

p2 Fitdata, 1, x, x², x

0.0735314 x² 1.46352 x 0.464835

Avslutningsvis de två modellerna tillsammans med mätdata. Jag har varnat Plot p1, p2 , x, 0, 7 , PlotStyle Blue, Orange ,

AxesLabel "x", "yi,p1,p2" , Epilog Red, PointSize 0.03 , Point data

1 2 3 4 5 6 7 x

2 4 6 yi,p1,p2

Avslutningsvis modellen med tredjegradspolynom direkt med Fit[data,f,x].

p3 Fitdata, 1, x, x², x³, x

0.108543 x³ 1.0122 x² 1.57679 x 1.77314

Det är nu bara att välja den modell som önskas beware

Plot p₁, p₂, p₃ , x, 0, 7 , PlotStyle Blue, Orange, Green ,

AxesLabel "x", "yi,p1,p2,p3" , Epilog Red, PointSize 0.03 , Point data

1 2 3 4 5 6 7 x

2 4 6 yi,p1,p2,p3

Exempel: Anpassa med minsta kvadratmetoden y x a bx² till mätvärdena

x 5 7.5 12 15 25

y 13.1 18.1 70.2 109 301

Lösningsförslag: De fem mätvärdena möblerar och i det överbestämda ekvationssystemet

1 x₁² 1 x₂² a

b y1

y2 a

b

för de efterlängtade parametrarna a och b, som levereras av normalekvationerna a

b .

data 5 7.5 12 15 25

13.1 18.1 70.2 109 301 ; 1, ² & data All, 1 ; . .data All, 2 1 25

1 56.25 1 144 1 225 1 625

1 1 1 1 1

25 56.25 144 225 625

5. 1075.25

1075.25 465 775. 511.4, 224 104.

Solve . . a, b .data All, 2 a 2.36283, b 0.486598

Eller direkt med Fit.

f Fitdata, 1, x², x

0.486598 x² 2.36283

En bild piggar alltid upp

Plot f, x, 0, 25 , PlotStyle Blue, AxesLabel "x", "yi,f x, " , Epilog Red, PointSize 0.03 , Point data

5 10 15 20 25 x

50 100 150 200 250 300 yi,fx,͍

Exempel: Anpassa med minsta kvadratmetoden y x ax bx² till mätvärdena

x 2 3 5 8

y 1.9 0.1 11 39

Lösningsförslag: De fyra mätvärdena möblerar och i det överbestämda ekvationssystemet x1 x₁²

x2 x₂² a b

y1

y2 a

b

för de sökta parametrarna a och b. Bestäm nu a och b enligt MKM receptet och normalekvationerna a

b .

data 2 3 5 8

1.9 0.1 11 39 ;  , ² & data All, 1 ; . .data All, 2 2 4

3 9 5 25 8 64

2 3 5 8 4 9 25 64

102 672

672 4818 362.9, 2762.5

Solve . . a, b .data All, 2 a 2.70872, b 0.951174

Eller direkt med Fit.

f Fitdata, x, x², x

2.70872 x 0.951174 x²

Javisst vill vi se en bild

Plot f, x, 2, 8 , PlotStyle Blue, AxesLabel "x", "yi,f x, " , Epilog Red, PointSize 0.03 , Point data

3 4 5 6 7 8 x

40 30 20 10

yi,fx,͍

Exempel: Kväve och syre har atommassorna N 14 och O 16. I tabellen nedan finns experimentellt uppmätta molekylmassor för sex olika kväveoxider. Använd minsta kvadratmetoden för att bestämma kvävets och syrets atommassor.

NO N₂O NO₂ N₂O₃ N₂O₄ N₂O₅ 30.01 44.01 46.03 75.98 92.02 108.01

Lösningsförslag: Lite olika blandningar av “äpplen och päron” på vågen. De sex vägningarna möblerar och i det överbestämda ekvationssystemet mN

mO för de sökta atommassorna. Bestäm nu mN och mO med normalekvationerna mN

mO .

1 2 1 2 2 2

1 1 2 3 4 5 ; 30.01, 44.01, 46.03, 75.98, 92.02, 108.01 ; . . 1 1

2 1 1 2 2 3 2 4 2 5

1 2 1 2 2 2 1 1 2 3 4 5

18 29

29 56 716.08, 1302.15

Solve . . m_N, m_O . m_N 14.0008, m_O 16.0023

Exempel: Anpassa med minsta kvadratmetoden y x _{ax b}^x till mätvärdena x 1 2.5 4 6 8 y 1 1.5 2 2.5 3

Lösningsförslag: Här har vi en olinjär modell, men kan efter omskrivning lineariseras och möblera ett överbestämt ekvationssystem.

y ^x

ax b x 0 1 y

ax b x

1

y a ^b_x

1 _x¹

1

1 _x¹

2

a b

1 y₁ 1 y₂

a b

Bestäm nu a och b enligt MKM. Först matrisgodis sedan lösning av normalekvationerna

data 1 2.5 4 6 8

1 1.5 2 2.5 3 ; 1, 1

 & data All, 1 ;  . . ¹

data All,2  1 1 1 1 1

1 0.4 ¹₄ ¹₆ ¹₈

5. 1.94167

1.94167 1.2659 2.9, 1.5

Solve . . a, b .

1 data All, 2  a 0.2964, b 0.730302

Här kan det vara lämpligt att introducera den olinjära hjälpredan FindFit som sväljer modellen med hull och hår FindFitdata, x

a x b

, a, b , x

a 0.212872, b 1.0605 Ok, dé får väl bli en bild då

Plot x a x b

. , x, 1, 9 , PlotStyle Blue, AxesLabel "x", "y_i,f x, " , PlotRange 0, 4 , Epilog Red, PointSize 0.03 , Point data 

2 4 6 8 x

1 2 3 4 yi,fx,͍

Exempel: Sambandet mellan trycket p och volymen v i en gas följer lagen pvⁿ c, där n och c är konstanter. Bestäm dessa med minsta kvadratmetoden och mätserien

v 4.60 7.20 10.1 15.3 20.4 30.0 p 14.2 7.59 4.74 2.66 1.78 1.04

Lösningsförslag: Här har vi en olinjär modell, men kan lineariseras efter logaritmering av båda sidor pvⁿ c ln p nln v ln c ln c nln v ln p

1 ln v1

1 ln v2 ln c n

ln p1

ln p2 ln c

n

Bestäm nu ln c och n enligt MKM receptet och normalekvationerna ln c

n .

data 4.60 7.20 10.1 15.3 20.4 30.0 14.2 7.59 4.74 2.66 1.78 1.04 ; 1, Log & data All, 1 ;

. .Log data All, 2 1 1.52606

1 1.97408 1 2.31254 1 2.72785 1 3.01553 1 3.4012

1 1 1 1 1 1

1.52606 1.97408 2.31254 2.72785 3.01553 3.4012

6. 14.9573

14.9573 39.6764 7.83027, 16.1894

lncÅn Solve . . lnc, n .Log data All, 2 First lnc 4.77908, n 1.39359

Så modellen

p

lnc

vⁿ

. lncÅn 118.995

v^1.39359

Eller direkt med den olinjära hjälpredan FindFitdata, c

vⁿ

, c, n , v

c 119.337, n 1.39505

Detta eviga ritande

Plot p, v, 0, 25 , PlotStyle Blue, AxesLabel "v", "pi,p v " , Epilog Red, PointSize 0.03 , Point data

5 10 15 20 25 v

5 10 15 20 25 30 35 pi,pv

Exempel:Efter att en sjö försetts med ett antal fiskyngel uppskattades, bl.a. från åtgången av mat, hur antalet fiskar f utvecklades över tid t mätt i år. Man fick följande tabell

t 1 2 3 4 5 6 7 8 9 10

f 100 210 360 550 720 840 910 940 950 960 . I många biologiska sammanhang visar det sig att den så kallade logistiska modellen fungerar mycket bra. Denna har formen f t ^f

1 ^f_f

0 1 ^kt,

där de tre parametrarna f , f0och k styr dess utseende. Anpassa denna

0 2 4 6 8 10 tår

200 400 600 800 1000

fiskar

modell till data i tabellen och uppskatta hur många fiskyngel som planterades ut från början f0samt hur många fiskar det kommer att finnas vid jämvikt f , det vill säga efter lång tid.

Lösningsförslag: Vi har en typisk ickelinjär modell, med data och föreslagen modell.

data 1 2 3 4 5 6 7 8 9 10

100 210 360 550 720 840 910 940 950 960 ; modell f 1 ^f

f0

1 ^{k t}

;

Vi tar direkt hjälp av Mathematicas funktion för olinjär MKM. Vi kan därefter direkt avläsa de önskade parametrarna och rita en liten figur ända från år noll då utplanteringen skedde.

FindFit data, modell, f , f₀, k , t f 967.146, f₀ 49.6372, k 0.799629

Plot modell . , t, 0, 10 , AxesLabel "t år ", "fiskar" , PlotStyle Blue, PlotRange Automatic, 0, 1000 , Epilog Red, PointSize 0.03 , Point data

0 2 4 6 8 10 tår

200 400 600 800 1000

fiskar

Exempel: Inte sällan har man nödvändigtvis en modell som genererar ett icke linjärt system. Ta som exempel tomtegröt som ju avsvalnar enligt Newtons avsvalningslag ^T_t k T0 T . Lösningen får den välkända formen T t a b ^kt, där parametrarna a, b och k ingår olinjärt. Antag att tomten nedtecknat följande mätserie.

t min 1 3 5 10 15 20 T C 75 55 40 28 25 23

Lösningsförslag: På grund av olinjäriteten måste vi angripa problemet med en direkt minimering av felsumman S _{i 1}ⁿ _i², eller direkt med den olinjära hjälpredan FindFit. Vi provar båda! Men först mätdata och modell.

data 1 3 5 10 15 20

75 55 40 28 25 23 ; T t : a b ^{k t} Felet vid varje mätpunkt

data All, 2 T data All, 1

 a b ^k 75, a b ^{3 k} 55, a b ^{5 k} 40, a b ^{10 k} 28, a b ^{15 k} 25, a b ^{20 k} 23

Felsumman S _{i 1}ⁿ _i² får vi nu enklast om vi betraktar som en vektor. En enkel skalärprodukt gör jobbet.

S .

 a b ^{20 k} 23²  a b ^{15 k} 25²  a b ^{10 k} 28²  a b ^{5 k} 40²  a b ^{3 k} 55²  a b ^k 75²

Sedan minimering.

Smin, NMinimize S, a 0, b 0, k 0 , a, b, k 3.4631, a 23.0881, b 68.0623, k 0.264998

Alltså modellen T t .

23.0881 68.0623 ^{0.264998 t}

Visst blir man glad av upprepningar

Plot T t . , t, 0, 30 , PlotStyle Blue, AxesLabel "t min ", "Ti,T t C " , PlotRange 0, 100 , Epilog Red, PointSize 0.03 , Point data

0 5 10 15 20 25 30 tmin 20

40 60 80 100 Ti,Tt^éC

Eller varför inte den i särklass enklaste varianten med direkt anpassning av en olinjär modell.

FindFit data, T t , a, b, k , t a 23.0881, b 68.0623, k 0.264998

Exempel: För att studera en människas rörelser under gång har man utvecklat en datormodell. För att generera indata till denna utgår man från fotografier som tagits vid åtta tillfällen under två steg, en period. Tiden för en period är normerad till 1 för att enklare kunna simulera olika gånghastigheter. I nedanstående tabell är knävinkeln i grader angivet för de olika tidpunkterna. För att undvika minst sagt ryckig gång på datorskärmen vill man givetvis anpassa en funktion som levererar en vinkel vid godtycklig tidpunkt t.

data 0 0.125 0.25 0.375 0.5 0.625 0.75 0.875 1

180 160 165 175 170 130 110 150 180 ;

Lösningsförslag: Här är det viktigt att använda en periodisk funktion så att vi får en kontinuerlig gånghastighet (derivata) över periodgränsen. Använd därför en harmonisk funktion som automatiskt tar hand om periodiciteten, en så kallad Fourierserie, här

t _{k 0}³ akcos 2Πkt bksin 2Πkt , med modellparametrarna a0, a1, b1, a2, b2, a3, b3 eftersom b0 0. Så modellen.

f 1, Cos 2 Π t , Sin 2 Π t & 1, 2, 3 Flatten 1, cos 2Πt , sin 2Πt , cos 4Πt , sin 4Πt , cos 6Πt , sin 6Πt

Fit data, _f, t

23.4727 sin 2Πt 8.75 sin 4Πt 4.02728 sin 6Πt 3.55055 cos 2Πt 18.9167 cos 4Πt 1.78278 cos 6Πt 155.083

Varför inte rita lite

Plot , t, 0, 1 , PlotStyle Blue, AxesLabel "t", " i, t, " , PlotRange 100, 200 , Epilog Red, PointSize 0.03 , Point data

0.0 0.2 0.4 0.6 0.8 1.0 t

120 140 160 180 200 wi,wt,͍

Eftersom det i båda ändpunkterna är uppmätt 180 bör man kan kanske kontrollera att ingen översvängning skett. Knävinklar över 180 är klart ohälsosamma

FindMaximum , t, 1 179.523, t 0.99025

Exempel: Vid medicinsk behandling behöver man ibland mäta kapaciteten (volym och blodflöde) i ett inre organ. En känd mängd av ett spårämne injiceras i ingående artär strax intill organet. I utgående ven mäter man sedan hur spårämnets koncentration c [mg/l]

i blodet varierar med tiden. I tabellen nedan ses resultatet av en sådan undersökning sedan 5 mg av ett spårämne injicerats.

Uppgiften är nu att ta fram en modell där såväl blodflödet i [l/s] som organets volym i [l] kan bestämmas.

tid s c mg l tid s c mg l tid s c mg l tid s c mg l tid s c mg l

0 0 5 3.0 10 5.4 15 3.0 20 0.6

1 0 6 3.8 11 5.6 16 2.1 21 0.4

2 0 7 4.4 12 5.7 17 1.5 22 0.4

3 0.1 8 4.9 13 5.8 18 1.1 23 0.2

4 1.8 9 5.2 14 4.2 19 0.8 24 0.1

Lösningsförslag: Mätdata enligt uppgift. Tid i kolonn ett och koncentration i kolonn två. En bild förtydligar tabellen.

data SortPartitionFlatten

0 0 5 3.0 10 5.4 15 3.0 20 0.6 1 0 6 3.8 11 5.6 16 2.1 21 0.4 2 0 7 4.4 12 5.7 17 1.5 22 0.4 3 0.1 8 4.9 13 5.8 18 1.1 23 0.2 4 1.8 9 5.2 14 4.2 19 0.8 24 0.1

, 2;

ListPlot data, Joined True, AxesLabel "t s ", "ci mg l " , PlotStyle Red, Epilog Red, PointSize 0.03 , Point data

5 10 15 20 ts

1 2 3 4 5 6 cimgl