Minstakvadratproblem

I denna övning kommer Du att få ställa upp och lösa både linjära och ickelinjära minstakvadratproblem.

Greve Rumford utförde 1798 följande experiment. Ett kanonrör uppvärmdes till 130° F genom att en trubbig borr vreds runt av hästar i 30 minuter. Sedan fick röret svalna medan man då och då mätte temperaturen T i röret. Man fick följande mätserie:

Den omgivande temperaturen T_omg var 60 °F. Enligt Newtons avsvalningslag gäller att:

dT(t) / dt = -ß (T(t) - T_omg)

Om vi löser differentialekvationen får vi:

T(t) = T_omg + (T(0) - T_omg) e^{-ß t},  t >= 0

Vi vill bestämma parametern ß med hjälp av mätvärdena. Du kommer att lösa detta anpassningsproblem på tre olika sätt nedan.

Grundproblemet är dock att bestämma ß givet tider och temperaturer. Som problemet är formulerat ovan ingår inte parametern ß linjärt (problemet passar inte in i den teori vi såg på under föreläsningen). Vi kan dock formulera ett så kallat ickelinjärt minstakvadratproblem:
 
Givet tider och temperaturer (t_k, T_k), minimera kvadratsumman av:

(T₁ - T_omg - (T(0) - T_omg) e^{-ß t}₁)² + (T₂ - T_omg - (T(0) - T_omg) e^{-ß t}₂)² +...

med avseende på ß. Som tur är, finns ett färdigt Matlabprogram (lsqnonlin) som kan göra detta åt oss. Du kommer att använda detta program i andra och tredje deluppgiften. I första deluppgiften kommer Du att transformera det ickelinjära problemet till ett problem där parametern ingår linjärt. Du kan då använda teorin från föreläsningen.  

1. Parametern ß ingår inte linjärt i T(t), men genom att skriva om uttrycket kan man få ett linjärt problem (ledning flytta en term i uttrycket och applicera logaritmfunktionen). Detta linjära problem kan lösas på vanligt sätt (med \ ). (Eftersom vi bildar temperaturdifferenser i Newtons lag spelar det ingen roll att vi arbetar i °F och inte i K).
   Formulera problemet på vanligt sätt med A, x och b samt tala om vad som svarar mot dessa tre storheter i det aktuella problemet.
    

2. Använd Matlab-rutinen lsqnonlin (ge kommandot help lsqnonlin) för att lösa det ickelinjära problemet direkt (dvs. utan att skriva om problemet).
 
lsqnonlin ingår i den så kallade optimeringstoolboxen (något som inte följer med standard-Matlab).
    
   Hur man skriver Matlabfunktioner. Läs detta först om du har glömt hur Matlabfunktioner fungerar.
    
   Så här används lsqnonlin.
    
   lsqnonlin kräver ett startvärde på parametern. Vad är ett lämpligt startvärde i detta problem?
    

3. Ovanstående approximationer stämmer inte så bra (om Du plottar upp mätvärden och approximationer). Kanske har några värden blivit fel avskrivna (något kan ha tappats bort under de 200 år sedan experimentet utfördes). Kan Du komma på någon annan, fysikalisk, förklaring till diskrepansen?
   Betrakta nu även omgivningstemperaturen som okänd och lös det ickelinjära minstakvadratproblemet med avseende på både ß och T_omg. Vad är lämpliga startvärden?

Lös de tre problemen ovan, plotta sedan mätpunkter och de tre approximationerna i samma diagram. Spara parametervärdena. Kommentarer?

I verkliga experiment eller forskning får man givetvis inte justera data som i deluppgift tre. Sådana saker har dock hänt, läs (om Du vill) om Schönskandalen.

1. Bestäm a och b genom att logaritmera och sedan lösa ett linjärt minstakvadratproblem.
2. Använd lsqnonlin för att lösa det ickelinjära problemet direkt (använd bra startvärden).
3. Använd lsqnonlin för att lösa ett viktat ickelinjärt problem; använd modellen 1 = a x^b/ y och samma startvärden som i 2.

Lös de tre problemen ovan. Plotta sedan mätpunkter och de tre approximationerna i samma diagram, använd loglog-funktionen. Spara parametervärdena. Kommentarer? Tolka resultatet!

Använd följande data (första kolonnen innehåller x-värdena och andra kolonnen y-värdena).

   1.000e-03   3.155e-04
   1.584e-02   1.787e-02
   2.511e-01   3.475e-01
   3.981e+00   9.596e+00
   6.309e+01   5.388e+02
   1.000e+03   7.147e+03