Linjära ekvationssystem och minsta-kvadratmetoden i Matlab
Matematik 2 MVE340 Sjöingenjörsprogrammet
14 maj 2019
I bl.a. kursen i Reglerteknik kommer ni att använda beräkningsprogram- met Matlab. Vi visar här hur man använder detta program för att lösa lin- jära ekvationssytem samt anpassa en rät linje till givna mätdata med minsta- kvadratmetoden.
Linjärt ekvationssystem.
Anta att vi vill beräkna skärningspunkten mellan de tre planen x + 2y + 3z = 4, x + 3y − 5z = −2, 2x − 5z = 1. Detta innebär att lösa ekvationssystemet
x +2y +3z = 4 x +3y −5z = −2
2x −5z = 1.
I Matlab löser vi detta system på följande sätt. Vi bildar en så kallad matris av talen framför variablerna i vänsterleden och en kolumn av högerleden
>> A=[1 2 3;1 3 -5;2 0 -5] % Koefficientmatrisen
>> b=[4 ;-2; 1] % Högerledet
>> x=A\b % Notera det omvända divisionstecknet
Variabeln x innehåller nu lösningen x, y, z i turordning. Får man en varning kan man använda kommandot rref([A b]) för att undersöka vidare.
Minsta kvadratmetoden
Vi skall nu beskriva hur man kan anpassa en rät linje till mätdata i minsta- kvadratmetodens mening. Vi har en serie med n stycken mätpunkter (t1, y1), (t2, y2), . . . , (tn, yn) numrerade från 1 till n, som verkar ligga ungefär utefter en rät lin- je om man ritar dem i ett koordinatsystem. För att hitta den “bästa” linjen y = k · t + m enligt minsta-kvadratmetoden, skall man minimera summan av
kvadraterna på avvikelserna mellan linjens y-värden och mätvärdena yi. Mate- matiskt innebär det att minimera uttrycket
(y1− k · t1− m)2+ (y2− k · t2− m)2+ · · · + (yn− k · tn− m)2
Om vi sätter derivatorna av detta uttryck med avseende på k resp. m till noll fås efter lite omskrivning följande ekvationssystem
k ·X
t2 +m ·X
t =X
ty k ·X
t +m · n =X
y
Här löser vi ut k och m på vanligt sätt. I Matlab löser man detta på föjande smidiga sätt. Man skriver upp det ekvationssystem som innebär att den räta linjen går genom alla mätpunkterna (detta system är förstås olösbart eftersom det normalt inte finns någon sådan linje) och “löser” systemet som ovan med
“\”. Vi visar med ett exempel.
Exempel minsta-kvadratmetoden.
Anpassa en rät linje till till mätpunkterna t 1 2 3 4 5 y 2.1 3.3 4.0 4.9 5.7 . Lösning.
Vi löser först uppgiften för hand.
Xt = 1 + 2 + 3 + 4 + 5 = 15, X
t2= 12+ 22+ 32+ 42+ 52= 55, Xty = 1 · 2.1 + 2 · 3.3 + 3 · 4.0 + 4 · 4.7 + 5 · 5.9 = 69.0,
Xy = 2.1 + 3.3 + 4.0 + 4.7 + 5.9 = 20.0 vilket ger oss ekvationssystemet
(55k +15m = 69 15k +5m = 20
som har lösningen k = 0.9, m = 1.3 (t.ex. kan man lösa ut m ur den andra ekvationen och sätta in i den första).
För att lösa uppgiften i Matlab så ställer vi upp ekvationssystemet:
1k +m = 2.1 2k +m = 3.3 3k +m = 4.0 4k +m = 4.7 5k +m = 5.9
>> t=[1; 2; 3; 4; 5]
>> y=[2.1; 3.3; 4.0; 4.7; 5.9] % Högerledet
>> A=[1 1; 2 1; 3 1; 4 1; 5 1]; % Koefficientmatrisen, alternativt A=[t ones(size(t)]
>> losn=A\y % Beräkna k och m
>> k=losn(1);m=losn(2);
>> % Plotta mätpunkter ’*’ och räta linjen i samma bild:
>> plot(t,y,’*’)
>> hold on
>> x=linspace(t(1),t(end))’;
>> plot(x,k*x+m) ans =
0.9000 % k 1.3000 % m
t y
0 1 2 3 4 5 6
0 1 2 3 4 5 6
Övningsuppgifter
1. Bestäm mha av Matlab det andragradspolynom (dvs y = ax2+ bx + c) som passerar genom punkterna (−1, −4), (1, 3), (3, −7). (Svar: y =
−2.125x2+ 3.5x + 1.625.)
2. Bestäm för hand och mha Matlab den räta linje som i minsta-kvadratmetodens mening är bäst anpassad till mätdata (Svar: y = −1.63t + 6.73)
t 1 2 3 4 5
y 5.1 3.4 1.9. 0.3 -1.5
3. Strålningsintensiteten I hos en blandning av två radioaktiva ämnen avtar med tiden enligt
I(t) = A0e−λt+ B0e−µt
Man känner halveringstiderna för ämnena och därav vet man att λ = 0.29 och µ = 0.17. Genom mätningar har föjande data erhållits
t 1 2 3 4 5 6
I 8.01 6.18 4.71 3.68 2.86 2.20
Bestäm halterna A0och B0av ämnena i blandningen med minstakvadrat- metoden.
4. För bestämning av fjäderkonstanten k för en fjäder gjordes ett försök, där fjädern belastades med en kraft F och dess längd l mättes. Vid försöket fick man följande mätdata
F 1 2 3 4 5
l 7.97 10.2 14.2 16.0 21.2
Hookes lag säger att fjäderns längd beror linjärt på den kraft som belastar fjädern enligt
l = e +F k
där e är fjäderns längd då den är obelastad. Bestäm k med minstakvadrat- metoden.
Lösning
Uppgift 3. Vi får ett överbestämt linjärt ekvationssystem A0e−λtk+ B0e−µtk= Ik, k = 1, 2, . . . , 6 Lösning med Matlab:
>> t=[1 2 3 4 5 6]’;
>> y=[8.01 6.18 4.71 3.68 2.86 2.20]’;
>> lambda=0.29; mu=0.17;
>> A=[exp(-lambda*t) exp(-mu*t)];
>> x=A\y;
>> A0=x(1), B0=x(2) A0 =
8.4199 B0 = 2.0301
>> ym=@(t)A0*exp(-lambda*t)+B0*exp(-mu*t);
>> tm=linspace(t(1),t(end),100);
Uppgift 4. Vi får ett överbestämt linjärt ekvationssystem e +Fi
k = li, i = 1, . . . , 5 Lösning med Matlab:
>> F=[1 2 3 4 5]’; l=[7.97 10.2 14.2 16.0 21.2]’;
>> A=[ones(size(F)) F];
>> x=A\l;
>> e=x(1) e = 4.2360
>> k=1/x(2) k =
0.3100
