Johan Helsing, 19 november 2015
FMN140 HT15: Beräkningsprogrammering
Numerisk Analys, Matematikcentrum
Inlämningsuppgift 2
Sista dag för inlämning: onsdag den 25 november.
Syfte: träning på att skriva egna underfunktioner, på att hitta nollställen till
ickelinjära funktioner, och på att utföra snabb matris-vektormultiplikation.
I denna inlämningsuppgift är det bra att förstå programmen bisection.m,
newton.m, tabex.m, myvecsave.m, myown.m, och simpmat.m på kursens webbsida. Observera att nollställen till en funktion f (x) samtidigt är lösningar
till ekvationen f (x) = 0, och att dessa lösningar ibland kallas rötter.
Problem 1. En egen underfunktion. Skriv ett kort Matlabprogram
som består av ett huvudprogram och en egen underfunktion. Underfunktionen skall ha två inparameterar och en utparameter men får i övrigt heta vad
som helst och göra vad som helst. Underfunktionen skall anropas från huvudprogrammet i en for-loop som går tio varv. Svaren från de olika anropen
skall sparas i en vektor som heter resultvec.
Redovisning: lämna in programkod.
Problem 2. Nollställen till en ickelinjär funktion. Skriv ett Matlabprogram som beräknar och tabellerar de fem första nollställena till funktionen
h(x) = tan(x) − x , x > 0 .
Du har redan hittat nollställena grafiskt i Inlämningsuppgift 1. Nu skall du
beräkna dem så noggrant som möjligt med hjälp av iterativa metoder. Problemet att hitta nollställen till h(x) dyker upp i Eulers tredje knäckningsfall.
Lösningarna är relaterade till de lägsta knäckningsmoderna.
Detta skall göras: Beräkna och tabellera de fem första nollställena med
tre olika metoder och med alla siffror du tror är korrekta. Kontrollera även
om dina beräknade nollställen är riktiga nollställen genom att utvärdera
h(x) i dem och se om det blir (nästan) noll. Börja med intervallhalveringsmetoden. Tag sedan Newtons metod. Tabellera här även hur många iterationer
som behövs för att möta ditt avbrottskriterium. Beräkna och tabellera till
sist nollställena med Matlabs egen lösare fzero, som använder en kombination av metoder. Skriv help fzero eller doc fzero för mer information.
Tycker du att fzero var en besvikelse? Då blir slutsatsen att det ibland är
bättre att programmera själv än att förlita sig på “black-box lösare”.
sida 1 av 3
Programdetaljer: Kalla ditt program nonlin.m och använd flera egna
underfunktioner. Att skriva h(x) och dess derivata som underfunktioner är
naturligt. Men även intervallhalveringsmetoden och Newtons metod skall
skrivas som underfunktioner. (Jo, det blir faktiskt enklast så.) Du kan kalla
den underfunktion som utför intervallhalveringsmetoden på intervallet [a, b]
för mybisect och inleda med raden
function myzero=mybisect(a,b)
Här är a och b inparametrar och myzero är utparameter (svaret som kommer ut). Du kan kalla den underfunktion som utför Newtons metod med
startgissning xk för mynewton och inleda med raden
function [myzero,iter]=mynewton(xk)
Här är xk inparameter och myzero och iter är utparametrar (svaret som
kommer ut samt antalet iterationer som behövdes). Andra parameter- och
funktionsnamn går naturligvis också bra.
Redovisning: Lämna in programkod och tabeller. När du gör tabellerna
kan du låta Matlab skriva ut alla uträknade siffror. Sedan kan du gå in i
tabellen för hand och stryka de siffror du inte tror på.
Problem 3. Snabb matris-vektormultiplikation för gles matris. Programmet nedan gör följande: Det skapar först en tridiagonal matris A av dimension n × n och en slumpmässig kolonnvektor x av längd n. I Algoritm 1
multipliceras sedan A med x och resultatet kallas b1. Exekveringstiden, det
vill säga tiden som det tar att utföra själva matris-vektormultiplikationen,
mäts upp av tidtagaruret “tic-toc” och kallas t1.
function matvec
format compact
n=input(’Give system size: ’);
cen=rand(n,1);
sub=rand(n-1,1);
sup=rand(n-1,1);
A=diag(sub,-1)+diag(cen,0)+diag(sup,1);
x=rand(n,1);
disp(’Setup done’)
% *** Algoritm 1 ***
tic
b1=A*x;
t1=toc
Skriv av denna kod och provkör den för några olika värden på n.
I exemplet ovan används Matlabs egen matris-vektor multiplikation
“*”. Nu lägger vi till följande rader
sida 2 av 3
% *** Algoritm 2 ***
tic
b2=zeros(n,1);
for i=1:n
for j=1:n
b2(i)=b2(i)+A(i,j)*x(j);
end
end
t2=toc
check2=norm(b2-b1)
Här har vi ersatt Matlabs matris-vektormultiplikation “*”med en dubbel
for-loop. Den sista raden kollar att resultatet blir riktigt. Prova även det
utökade programmet. Jämför exekveringstiderna t1 och t2 för för några olika
värden på n. Tidtagaruret tic-toc har ofta bara en upplösning på ett par
millisekunder. Om exekveringen går fortare än så anges tiden till noll. För
att få mätbara tider bör du därför välja n > 1500. (Väljer du n alltför stort
tar dock slumptalsgenereringen tid, och i värsta fall tar även minnet slut).
Vilken algoritm går snabbast?
Ett problem med den dubbla for-loopen i Algoritm 2 är att de flesta
elementen i matrisen A är noll. Många multiplikationer görs i onödan. Skriv
nu en egen “Algoritm 3”, där for-loopen är modifierad så att den endast
går över de element i A som är skilda från noll. Kalla resultatet b3 och
beräkningstiden t3. Kontrollera att b3=b1. Observera att den första och
den sista komponenten av b3 kan behöva brytas ut ur for-loopen. Jämför på
nytt exekveringstiderna för några olika (stora) värden på n. Vilken algoritm
går snabbast? Försök att förklara resultaten. Några hjälpande fakta:
• Den programkod som döljer sig bakom “*”-operatorn i Algoritm 1 är
samma kod som i vår Algoritm 2.
• for-loopar man skriver själv kan gå långsammare än de for-loopar
som döljer sig inuti Matlabs fördefinierade operatorer.
• Beräkningsarbetet i Algoritm 2 växer proportionellt mot n2 , medan
beräkningsarbetet i Algoritm 3 växer proportionellt mot n.
Redovisning: Lämna in programkod för beräkning av b1, b2, och b3
med angiven tidsåtgång för några olika n och en diskussion av resultatet.
Överkurs: Skriv en algoritm som beräknar Ax ännu snabbare och som
klarar större matriser. Tips: Förbättra Algoritm 1 genom att skapa matrisen
A med kommandot spdiags i stället för med diag. Då förstår Matlab att
matrisen A är tridiagonal. När du sedan skriver b=A*x så väljer Matlab en
algoritm som liknar vår Algoritm 3. Gör help spdiags för mer info.
sida 3 av 3