Uppsala universitet
Institutionen för informationsteknologi
Teknisk databehandling
Lösningsförslag
Tentamen i Beräkningsvetenskap I/KF, 5.0 hp, 2015-03-17
Del A
1. Gäller att Ax = b ⇒ P Ax = P b ⇒ LU x = P b.
Beräkna P b = (6, 4, 2)T . Lös Ld = P b med framåtsubstitution vilket ger d =
(4, 6, 0)T . Lös U x = d med bakåtsubstitution vilket ger x = (−2, 3, 0)T .
2. (a) Iteration 1:
Hitta mittpunkten x1 = 0.7
Välj nytt intervall genom att hitta teckenbyte: f (0.6) > 0, f (0.7) > 0, f (0.8) <
0 ⇒ välj höger intervall, [0.7 0.8].
Felet (intervallets längd): |0.8 − 0.7| = 0.1 > toleransen
Iteration 2:
Mittpunkt: x2 = 0.75
Nytt intervall: f (0.7) > 0, f (0.75) < 0 ⇒ välj vänsterintervall [0.7 0.75].
Felet: |0.75 − 0.7| = 0.05 = toleransen (men kravet är strikt < toleransen)
Iteration nr 3:
Mittpunkt: x3 = 0.725
Nytt intervall: f (0.7) > 0, f (0.725) < 0 ⇒ välj vänsterintervall [0.7 0.725].
Felet: |0.725 − 0.7| = 0.025 < toleransen. Stoppvillkoret uppfyllt.
Lösningen blir x = x3 = 0.725 ± 0.025.
(b) Vilken rot man hittar är väldigt känsligt för vilken startgissning man väljer.
Väljer man t ex x0 = 1 så hittar N-R t ex den rot som ligger mellan 1 och 1.5
istället för den första positiva roten. Om man inte vet något om funktionen är
det sannolikt att man hamnar i fel rot.
Ett annat problem är att N-R någon gång under iterationerna hamnar i en extrempunkt eller nära en extrempunkt (där f 0 (x) är noll eller nära noll). NewtonRaphson hamnar då i ∞ eller väldigt långt bort från lösningen.
3. Då h minskar med en faktor 2 minskar felet med en faktor 23 = 8, så det nya felet
blir 0.152/8 = 0.019
1
4. (a) Metoden konvergerar mot lösningen och eftersom felet i x3 är kvadraten på felet
i x2 , som i sin tur är kvadraten på felet i x1 , så är det kvadratisk konvergens.
Rätt svar är alltså (ii) och (iv).
(b) För h < 10−3 börjar avrundningsfelen dominera och för h > 10−3 dominerar
diskretiseringsfelet.
Diskretiseringsfelet avtar med minskande h och följer en rät ”snygg” förutsägbar
linje. Avrundningsfelet beter sig mer slumpmässigt, men trenden är att det ökar
med minskande h.
5. (a)
i+j är:
2
j är:
2
i är:
3
(b) Korrekt följd är 3, 7, 2, 5, 1, 4 (3 och 7 kan växla plats). Rad 6 är onödig.
Programmet blir alltså:
x = [5 1 2 3 1];
summa = 0;
for i = 1:length(x)
summa = summa + x(i);
end
medel = summa/length(x);
2
Del B
6. När man har mätvärden måste man använda Trapetsmetoden, eftersom Simpsons
metod kräver dubbelintervall, dvs udda antal punkter. I uppgiften skulle metoden
klara av udda och jämnt antal punkter.
När man har en kontinuerlig funktion kan man använda vilken indelning som helst
och det är då ingen begränsning att använda dubbelintervall. Vi väljer då Simpsons
metod eftersom den är noggrannare och kräver färre beräkningar för att uppnå en
viss noggrannhet.
När det gäller det totala felet så består det i båda fallen av diskretiseringsfel E och
funktionsfel. Det totala felet blir då E + (b − a) · .
I fallet då man har en datamängd är E bestämd av antalet datavärden, dvs diskretiseringen är given. I fallet då man har kontinuerlig funktion så kan man diskretisera
”fritt” så att E uppfyller en given tolerans.
Funktionsfelet bestäms av längden på integrationsintervallet och , där = M då
man har en kontinuerlig funktion. Om man har en datamängd där absoluta felet i
f (xi ) ≤ 0.5 · 10−3 så blir funktionsfelet (b − a) · 0.5 · 10−3 .
7. Formulera problemet som en ickelinjär ekvation. Tryckvektorn fås ur lösningen av
ekvationsystemet, som beror av indata r, dvs p beror av r. Ekvationen som ska lösas
är
min pk (r) = c ⇒ min pk (r) − c = 0.
k=1,...,n
k=1,...,n
Det här är en ickelinjär ekvation och vi kan använda Matlabs fzero för iterera oss
fram till ett r som är ett nollställe till ekvationen.
Om vi låter g(r) beteckna vänsterledet i ovanstående ekvation och implementerar
funktionen i Matlab:
function gr = g(r)
global L U P c
b = hogerled(r);
d = L\(P*b);
p = U\d;
gr = min(p) - c
%
%
%
%
Skapa högerledet givet r
Framåtsubstitution
Bakåtsubstitution
Beräkna g(r) = min(p) - c
Här förutsätts att L, U och P beräknats i huvudprogrammet genom LU-faktorisering
av A. Eftersom Ap = b kommer att lösas en gång per iteration i fzero så är det ur
effektivitetssynpunkt lämpligt att använda LU-faktorisering.
Följande huvudprogram löser problemet enligt idén ovan. A förutsetts vara given.
global L U P c
[L, U, P] = lu(A);
% LU-faktorisera A
3
c = input(’Ge kritiskt tryck: ’); % Ge det kritiska trycket c
r = fzero(@g, r0);
% Hitta nollstället m iterativ metod
% Skriv ut lösningen
disp([’Trycket i reservoaren måste vara minst ’ , num2str(r) , ’ bar’]);
Startgissningen eller startintervallet r0 kan man t ex hitta genom att rita upp funktionen (använd t ex fplot i Matlab).
Det är ur programmeringssynpunkt bättre att inte använda globals, utan parameteröverföring. Global-deklarationerna tas då bort och funktionshuvudet ändras till
function gr = g(r, L, U, P, c). Anropet till fzero ändras till r = fzero(@(r)
g(r,L,U,P,c), r0).
Observera att problemet behöver inte alls lösas med Matlabkod, utan du kan beskriva
en algoritm med någon typ av pseudokod eller vanlig text, t ex genom att enbart
använda kommentarsatserna ovan.
4