Efter laborationen ska du
I den avslutande temauppgiften får du se hur man kan arbeta med mätdata i Matlab.
Differentialekvationen
saknar analytisk lösning.
Plotta lösningen och bestäm y.
Använd någon av Matlabs inbyggda differentialekvationslösare, förslagsvis ode45.
Om du skriver y() i Matlab (där y är vektorn av y-värden från den numeriska lösningen) får du fel värde. Varför då?
Funktionen y och vektorn y är två olika saker. Se även jämförelsen mellan Matlab och Java i laboration 1.
help end
y ≈ Lösningen divergerar kring t ≈ . y kan därför inte fastställas.
I samband med ett experiment har vi fått följande mätdata. För bekvämlighets skull anges de här i Matlab-syntax. (Kom ihåg att koden markeras automatiskt om du placerar muspilen över den, och sedan lätt kan kopieras med (Windows, Linux) ctrl-C eller (Mac) cmd--C.)
Plotta värdena som punkter (t.ex. '+'). Fundera på hur ett eventuellt underliggande samband mellan x och y kan se ut.
Skapa ett polynom av lägsta möjliga grad som passar in exakt på mätpunkterna. Plotta polynomet i samma graf som punkterna ovan.
Med sex punkter är lägsta möjliga grad fem (ett femtegradspolynom).
Tycker du att polynomet är en bra beskrivning av det undersökta sambandet? Pröva gärna att sträcka ut x-axeln till –10 ≤ x ≤ 10 och plotta polynomets värde igen.
För att sträcka ut x-axeln räcker det inte att justera grafens axlar med axis. Du måste även skapa en ny vektor av x-värden för det nya intervallet, och sedan plotta polynomets värde för dessa x-värden.
Återställ därefter axlarna:
Anpassa en rät linje till mätpunkterna enligt minsta kvadrat-metoden. Plotta linjen som en blå linje i samma graf som ovan.
Från ett annat experiment har vi fått följande mätdata:
Vi tror att vi har att göra med ett samband av typen
Uppskatta värdena på parametrarna a och b med minsta kvadrat-metoden.
Kurvanpassningen blir mycket lättare om vi skriver om sambandet som en rät linje. Det kan göras med hjälp av några variabelbyten.
Vi vill alltså skriva sambandet som wx = kx + m för lämpligt valda w, k och m.
Ett snarlikt exempel (för ett samband av typen yt = C·eat) visas i följande film.
(här har något gått fel; kontakta Patrik) båda leden i sambandet 1 ovan, det vill säga inför
där k = och m = .
Genom att göra
får du alltså ett antal nya punkter (x,w). Om vårt antagande (antagande 1 ovan) stämmer borde dessa bilda en rät linje. Pröva gärna
Kan detta tänkas vara en rät linje (med lite mätstörningar)?
Denna räta linje är inte så värst intressant i sig, men den låter oss få fram värdena på k och m. Dessa får du av Matlab när du anpassar en rät linje till punkterna (x,w). Använd polyfit.
När du nu har värden på k och m kan du räkna ut vad a och b blir.
Plotta mätvärdena x och y som punkter, tillsammans med din kurvanpassning yx som en heldragen kurva. Ser den anpassade kurvan rimlig ut?
Genom att sätta in dessa värden på a och b i det ursprungliga sambandet (1 ovan) kan du plotta den anpassade kurvan.
Definiera en funktion, exempelvis y2, för sambandet 1 ovan:
Här ska dina värden på a och b ingå. Gör sedan
a =
b =
Återstoden av laborationen ska vi ägna åt en (fiktiv) serie mätdata. Serien kommer från en provkörning med en elbil av nästa generation. En sensor mäter bilens hastighet 150 gånger per sekund. Vi ska utifrån dessa hastighetsdata skapa oss en uppfattning om bilens prestanda.
Mätdata från provkörningen finns i filen race.txt. Vi kan låta Matlab ladda ner filen åt oss:
Ladda ner race.txt enligt ovan. Därmed hamnar den rätt bland dina Matlab-filer.
(Det går lika bra att själv ladda ner filen från webbadressen ovan, bara man ser till att den hamnar på rätt plats i datorn.)
För att läsa in mätdata använder vi Matlab-funktionen csvread. Den läser tal från en textfil och passar alltså bra här. För vår fil returneras en kolonnvektor av värden.
Läs in mätvärdena och undersök resultatet:
Vi har alltså en kolonnvektor med 6000 element. Eftersom sensorn mäter 150 ggr/s omfattar mätningarna 6000 / 150 = 40 sekunder. Värdena har enheten m/s, och elbilen har alltså efter tre sekunder uppnått drygt 25 m/s. (Det är en ganska kraftig acceleration.)
Skapa en vektor med motsvarade tidsvärden, dvs för tiden 0 ≤ t ≤ 40.
(Hur många värden behövs?)
Plotta hastighetskurvan. Lägger du märke till något speciellt?
Mätserien innehåller spikar: vid några tillfällen har vi fått mätvärdet 100 (som råkar vara den högsta hastighet sensorn kan mäta) istället för den verkliga hastigheten. Detta beror troligen på glappkontakt i mätutrustningen.
Vi vill filtrera bort dessa spikar, bland annat för att vi vill uppskatta maxhastigheten. (Man ska naturligtvis vara försiktig med att filtrera mätdata, och här gör vi det endast för att kunna arbeta vidare i väntan på bättre mätningar.)
Observera att du inte ska använda någon loop i den här uppgiften.
Ersätt alla spikar med värdet 0 (noll). Plotta den filtrerade hastighetskurvan.
Vi kan se i diagrammet (från föregående uppgift) att den "riktiga" hastigheten är < 80.
Tänk på att vektorindex i Matlab kan vara vektorer i sig. Då pekas många vektorelement ut på en gång. Vad ger följande Matlab-uttryck för resultat?
find fungerade igen? Även om spikarna fortfarande syns, så påverkar de inte längre maxhastigheten. Vi kan emellertid göra något bättre: vi kan ersätta varje spik med det närmast tidigare värdet. Med andra ord: v(k) ersätts med v(k-1) om v(k) kan anses vara en spik.
Läs in mätvärdena på nytt. Ersätt alla spikar med närmast föregående värde enligt ovan.
Ovan föreslås att vi ersätter v(k) med v(k-1). Här kan k vara ett uttryck. Vilket då?
Ett lämpligt uttryck finns i en tidigare ledtråd i denna uppgift.
Vilken maxhastighet uppnåddes under mätperioden? Kontrollera att ditt svar är rimligt genom att jämföra med de plottade hastighetsvärdena.
help max
Den tillryggalagda sträckan är
Hur långt rörde sig testfordonet under de 40 sekunderna?
Notera att integralen här har alltså bygger på mätvärden, snarare än en analytisk funktion.
help trapz
Pröva
där t är den vektor med tidsvärden du skapade i uppgift D4, och v är motsvarande hastighetsvärden.
Är den beräknade sträckan rimlig? Beräkna medelhastigheten v = s / t och jämför med maxhastigheten i tidigare uppgift. Kan det stämma?
Vi ska nu undersöka en annan mätserie, uppmätt under konstant acceleration i fem sekunder. Föraren har alltså hållit fordonets acceleration så konstant som möjligt. Som tidigare är det hastigheten som mätts.
Mätresultaten finns i filen const_accel.txt. Ladda som tidigare ner filen till din dator med Matlab-raden nedan.
Läs in mätresultaten till en vektor v. Skapa också en motsvarande vektor av t-värden för tiden 0 ≤ t ≤ 5.
Plotta mätvärdena för 0 ≤ t ≤ 5. Verkar föraren ha lyckats hålla konstant acceleration?
(Vi ser att glappkontakten i sensorn tycks vara åtgärdad.)
Vi ska nu ta reda på värdet av denna konstanta acceleration. Vi vet ju att
Man skulle kunna tänka sig att derivera kurvan numeriskt (ungefär så som du gjorde i laboration 1). Varför skulle det inte fungera här?
Anpassa istället ett lämpligt polynom till hastighetskurvan. Vad är bilens (nästan) konstanta acceleration?
Om accelerationen är konstant är vt = at + v0.
Uttrycket ovan är ett (förstagrads-)polynom. En av koefficienterna i det anpassade polynomet motsvarar accelerationen.
Om du får felmeddelandet
så beror det på att du anropat polyfit med en kolonnvektor och en radvektor. Du kan se detta genom att skriva (givet att dina vektorer heter t och v)
Då ser du att vektorerna har något olika dimensioner.
Du kan lösa problemet med vektorstorlekarna genom att transponera en av vektorerna, exempelvis genom att skriva t.' istället för t i anropet till polyfit.
När du kommit hit är du klar med kursens tredje och sista Matlab-laboration. Läs igenom uppgifterna igen och se till att du förstår dina lösningar inför redovisningen.