Transcript Lab 2

SF1547,numd15
LABORATION 2
Ekvationer, integraler, differentialekvationer
Genomför denna laboration genom att göra de handräkningar och MATLAB-program
som begärs. Var noga med att skriva upp resultaten och svaren på frågorna. Sammanställ problem, resultat och svar på frågorna i en kort skriftlig rapport som lämnas in
vid kontrollskrivningen (KS) den 7/5 2015. Rapporten är tillåtet hjälpmedel vid KS.
Bonus, 2 poäng, ges om KS samt rapport är godkänd. Rapport inlämnad efter den 7/5
ger ej bonuspoäng. OM KS ej är godkänd eller om man uteblir, måste Lab2 redovisas
muntligt enligt särskilt schema, men då ges ingen bonus. Om KS är godkänd men ej
rapporten måste rapporten kompletteras och lämnas in före ett på rapporten angivet
datum.
1. 1. Olinjär skalär ekvation - konvergensundersökning
Man vill bestämma samtliga rötter till följande ekvation
7x − 20 cos (3x − 5) = 0
(a) Rita grafen för funktionen y(x) = 7x − 20 cos (3x − 5) med MATLAB.
Samtliga nollställen till y(x) skall synas i grafen. Pröva med olika intervalllängder. Hur många nollställen finns det och i vilket intervall ligger dom?
(b) Genomför, med papper, penna och miniräknare tre iterationer med NewtonRaphsons metod och startvärdet x0 = 1.5.
(c) Skriv ett program för Newton-Raphsons metod i MATLAB som beräknar
samtliga rötter till ekvationen med relativfel högst 0.5 · 10−9 . Låt programmet skriva ut tabeller som visar hur iteraten konvergerar mot nollställena.
Utöka tabellen så att följande storheter kan avläsas: 1) Antal iterationer
som krävs för att få önskat relativfel, 2) Konstanten K i definitionen av
kvadratisk konvergens (se NAM sid 75).
(d) Visa att ekvationen kan skrivas om till fixpunktsiterationen
xn+1 = F (xn )
där
5
7
cos (3x − 5)
F (x) = x +
8
14
Genomför med papper, penna och miniräknare, tre iterationer med startvärdet x0 = 1.5. Förklara teoretiskt vilka av ekvationens rötter just denna
fixpunktsiteration kan konvergera mot om man väljer startgissningar nära
rötterna. Ledning: Läs i NAM avsnitt 6.6.
(e) Skriv ett MATLAB-program som undersöker praktiskt mot vilka rötter just
denna fixpunktsiteration konvergerar om man väljer startgissningar nära
rötterna. Avbryt fixpunktsiterationen efter ett maximalt antal iterationer,
t ex 100, för att undvika en oändlig slinga om metoden inte konvergerar.
Stämmer konvergensen för rötterna med de som konvergerar enligt d)?
2. Numerisk integrering - noggrannhetsordning
Följande integral ska beräknas
I=
Z
1
√
x + 2 dx
−1
A. Gör först följande laborationsförberedande uppgifter a)-c):
a) Rita först en graf över integranden, dvs y =
integralens värde ur grafen. Svar: I ≈
√
x + 2. Uppskatta den sökta
b) Beräkna integralens värde analytiskt, dvs exakt. Svar med 10 korrekta decimaler: I =
I nedanstående tabell har ett antal trapetsvärden som approximerar I redan
beräknats:
h
1
0.5
0.25
0.125
0.0625
T (h)
ET = T (h) − I
Textr (h)
ET extr = Textr (h) − I
2.796336
2.797160
2.797366
c) Räkna ut och fyll i för hand de värden som ska stå i de tomma rutorna.
Antag att exakta värdet tas från b)-uppgiften ovan. Avrunda till lämpligt
antal decimaler! I kolumnerna för felen i trapetsregeln resp det extrapolerade trapetsvärdet (kolumn 3 resp 5) avtar felen ungefär proportionellt
mot en potens av steglängden h. Vilken är denna potens för trapetsvärdet
resp det extrapolerade trapetsvärdet.
B. Skriv sedan ett MATLAB-program som beräknar integralen numeriskt med trapetsregeln (se NAM 5.2). Undersök hur trunkeringsfelet, ET , beror av steglängden, h, genom att plotta felen som funktion av steglängderna
h = 1/8, 1/16, 1/32, 1/64. Använd MATLABs kommando loglog för ploten.
Visa hur vi uppskattar metodens noggrannhetsordning, p, med hjälp av denna
plot. Vad blir p? Stämmer det med uppgift c) ovan?
Utöka sedan programmet så att integralen beräknas med trapetsregeln och en
extrapolation. Plotta i samma diagram som ovan felet ET extr som funktion av
h. Vad blir noggrannhetsordningen p nu? Stämmer det med c) ovan?
3. Skärning mellan ellipser - ickelinjärt ekvationssystem
Beräkna med Newtons metod alla skärningspunkter mellan den sneda ellipsen
som definieras av 0.4x2 + y 2 − xy = 10 och den ellips med centrum i (4, 2) och
halvaxlarna a = 4 och b = 6 (parallella med koordinataxlarna) som beskrivs av
(x − 4)2 /a2 + (y − 2)2 /b2 = 1.
Bägge ellipserna ska först ritas i en graf. (Den sneda ellipsen skrivs bäst i polär
form inför uppritningen.) och de erhållna skärningspunkterna markeras med en
ring. Skriv ut tabeller som visar konvergensen mot skärningspunkterna. Hur är
det med egenskapen kvadratisk konvergens i din algoritm?
4. Differentialekvation - lösning plottad i 2D och 3D
Givet differentialekvationen
dy
d2 y
+ πyex/3 (2 sin(πx) + πycos(πx)) = y/9,
2
dx
dx
y(0) = 1,
dy
1
(0) = −
dx
3
a) Vi vill lösa differentialekvationen på intervallet 0 ≤ x ≤ L. Inför nya vady
riabler u1 = y och u2 = dx
så att differentialekvationen kan skrivas om
till ett system av två första ordningens ODEer, se NAM avsnitt 8.5. Lös
detta system med Runge-Kuttas metod, RK4 fram till x = L, där L = 2.6.
Använd steglängden h = 0.2. Jämför denna lösning med den lösning som
erhålles med h = 0.1. Hur mycket skiljer sig dessa lösningar då x-värdet är
x = L? Rita upp lösningskurvan y(x) då h = 0.1.
b) Om lösningskurvan i a) roteras kring x-axeln erhålles en 3D figur som liknar
en lur, se fram sidan av NAM. En snygg 3D lurbild ritas i MATLAB på förljande sätt: Låt x och y vara kolumnvektorer för lösningen erhållen i a) med
h = 0.1. Skapa en radvektor för rotationsvinkeln 0 ≤ ϕ ≤ 2π med lagom
steg, t ex 2π/30. Bilda matriser X, Y och Z enligt X=x*ones(size(fi));
Y=y*cos(fi); Z=y*sin(fi); Skriv mesh(X,Y,Z) som ger en nätfigur eller
välj surf(X,Y,Z) som ger en fylld 3D-figur, se även 15.14 i MATLAB 7 i
korthet.
Skriv en kort rapport som lämnas in av gruppen (2 personer) vid KS-tillfället den
7/5 2015. Den ska innehålla: 1) Problemen som löses, 2) Numeriska metoder som
används, 3) Resultat i form av text, siffror, tabeller och grafer, 4) Ev slutsatser, 5)
Referenser, 6) MATLAB-programmen som använts för att lösa uppgifterna. Som framsida till rapporten används den försättssida som finns på kursens hemsida. Förutom
MATLAB-programmen bör rapporten inte bestå av fler än ca 4-6 sidor.