1. No category

SF1518,SF1519,numpbd15
Laboration 3 i Numeriska Metoder och Programmering, ht 2015
DIFFERENTIALEKVATIONER
OBS! Rapportanvisningar på sista sidan
Deluppgifterna 1a, 1b, 2a, 2b, 2c redovisas i form av en rapport som gruppen lämnar in utskriven på papper, i brevlådan vid studentexpeditionen på Matematik senast
torsdag 22/10 2015.
Anvisningar för rapportens utformning, finns här på sista sidan.
Rapport som är godkänd och inlämnad i tid ger 1.5 bonuspoäng till tentamen.
1. Pendelrörelse
Rörelsen hos en pendel beskrivs av Newtons kraftekvation:
mL
d2 ϕ
+ mg sin (ϕ) = 0
dt2
I denna andra ordningens icke–linjära differentialekvation är m massan, L pendelns
längd och g tyngdkraftaccelerationen.
1a. Svängningstiden – numerisk integration
Här ska svängningstiden för pendeln beräknas.
Svängningstiden T beror av utslagsvinkeln ϕ0 enligt formeln:
s
T =4
L
I(ϕ0 )
g
där L är pendelns längd, g är tyngdkraftaccelerationen och
Z π/2
I(ϕ0 ) =
0
p
dϕ
,
1 − k 2 (sin ϕ)2
k = sin
ϕ0
2
Låt L = 0.25 och g = 9.81 m/s2 . Skriv ett program som beräknar svängningstiden
T för ϕ0 -värden motsvarande 5, 10, . . . , 90 grader. (Observera att integralen ovan är
uttryckt i radianer och att argumentet till funktionen sin i MATLAB uttrycks i
radianer.) För att beräkna integralen I(ϕ0 ) används MATLABs funktion integral.
Plotta resultatet, dvs T som funktion av ϕ0 i en graf (plotbild 1).
En ofta använd approximation av T är svängningstiden för små svängningar:
s
T̃ = 2π
L
g
Denna approximation är bra för små utslagsvinklar, men relativfelet ökar med ökande
utslagsvinkel. Plotta T̃ -värdet i samma graf som ovan, dvs plotbild 1.
Gör även en graf som visar det procentuella relativfelet i T̃ , dvs 100 ∗ (T − T̃ )/T som
funktion av utslagsvinkeln då 5 ≤ ϕ0 ≤ 90 grader, plotbild 2.
1b. Simulering av pendelrörelsen – begynnelsevärdesproblem
Om pendeln startar från stillastående med en given utslagsvinkel ϕ0 , ges begynnelsevillkoren av
dϕ
ϕ(0) = ϕ0 ,
(0) = 0
dt
Skriv om systemet på vektorform, dvs
du
= f (u)
dt
Skriv den MATLAB–funktion som svarar mot högerledet.
Simulera pendelns rörelse genom att anropa MATLAB-funktionen ode45. Tänk
vid funktionsdefinitionen på att ode45 internt kan använda er funktion på vektorer.
Välj själv ett lämpligt tidsintervall och några startvinklar ϕ0 och skriv ut en graf över
vinkeln ϕ som funktion av tiden t. Låt parametrarna i problemet ha samma värden
som i deluppgift 1: L = 0.25 och g = 9.81. Avläs i grafen svängningstiden T och jämför
med det värde du får i deluppgift 1.
2. Stationär värmeledning – randvärdesproblem
Randvärdesproblemet
d
− dx
k(x) dT
dx = Q(x),
T (0) = T0 ,
0<x<L
(1)
T (L) = TL
(2)
beskriver temperaturfördelningen T (x) i en cylindrisk stav av längden L. Här är k(x)
stavens värmeledningsförmåga och Q(x) den värmemängd som per tidsenhet och volymsenhet genereras i staven, t.ex. genom radioaktivitet. Vänster och höger ändpunkt
på staven hålls vid konstant temperatur T0 respektive TL .
Antag att L = 3.4 [m], T0 = 373 [K], TL = 295 [K], samt att k(x) [J/(K · m · s)]
ges av det styckvisa polynomet
k(x) = 5,
k(x) = c0 + c1 (x − 2.4) + c3 (x −
0 < x < 2.3
2.4)3
+ c5 (x −
2.4)5 ,
k(x) = 1,
c0 = 3,
c1 = −3.75 · 101 ,
2.3 ≤ x ≤ 2.5
2.5 < x < L
c3 = 2.5 · 103 ,
c5 = −7.5 · 104 ,
Visa att k(x) och k 0 (x) är kontinuerliga!
Q(x) [J/(s · m3 )] ges av funktionen
2
Q(x) = 300e−3(x−1) ,
0 < x < L.
Differentialekvation (1) med randvillkoren (2) kan lösas numeriskt med hjälp av
finita differensmetoden. Om vi diskretiserar intervallet [0, L] enligt xi = ih, i =
0, 1, 2, . . . , N − 1, N , där h · N = L, sätter y(xi ) ≈ yi och approximerar andraderivatan med centraldifferens erhålles ett linjärt ekvationssystem
AT = b
(3)
där A är en matris, T är en vektor som innehåller temperaturvärden, och b är en
vektor som beror av Q(xi )-värdena samt temperaturen i stavens ändpunkter.
a) Skriv ner matrisen A och vektorerna T och b för N = 4 med papper och
penna uttryckt i k, k 0 , och
Q. Vilken struktur har
matrisen A?
df (x)
dg(x)
d
Ledning: dx (f (x)g(x)) = dx g(x) + f (x) dx
b) Skriv en MATLAB-funktion som löser randvärdesproblemet och räknar ut ett
närmevärde för temperaturen T i punkten x = 1. Låt funktionen ta antalet delinterval,
N , som indata, och returnera närmevärdet till T (1) som utdata. Testa funktionen med
t.ex. N = 17, N = 34, och N = 68. Tar programmet lång tid att köra? Här är två
tips för snabbare program:
• Matrisen A kommer endast att ha ett fåtal nollskilda element. Kommandot
sparse är ett kommando för att skapa glesa matriser eller för att tala om för
MATLAB att matrisen kan lagras glest. MATLAB-satsen A=sparse(A) ändrar
till gles lagring i MATLAB.
• Hur mycket skall felet minska då steglängden halveras? Gör det det i ditt program?
c) Skriv ett MATLAB-program som använder funktionen i b) till att räkna ut
temperaturen T i punkten x = 1 med 4 säkra decimaler. (Glöm inte att visa vad du
gjort som gör att du vågar lita på dina 4 decimaler.) Plotta också temperaturen som
funktion av x på hela intervallet 0 ≤ x ≤ L.
ANVISNINGAR för rapporten om deluppgifterna 1a, 1b, 2a, 2b, 2c
Deluppgifterna i Lab 3 redovisas i form av en rapport, utskriven på papper. En
rapport per grupp lämnas in! Rapporten kan författas på svenska eller engelska.
Målsättningen med en teknisk rapport är att ge läsaren en välskriven, väldisponerad
och sammanfattande beskrivning av t.ex. en teknisk utredning, undersökning eller experiment. Tekniska rapporter följer vanligen ett visst mönster vad gäller dispositionen.
Man ska tänka på vilken målgrupp som rapporten vänder sig till. I detta fall kan man
tänka sig en civilingenjör med till hälften glömda kunskaper i numeriska metoder,
d.v.s. du om några år.
Rapporten ska vara en skrift, där avsnitten 4-5 nedan upptar 4-7 sidor.
Rapporten ska innehålla:
1. En speciell försättssida som Matematikinstitutionen vill ha den. Länk finns på
kursens webbplats.
2. Ett omslag som innehåller rapportens titel, författarnas namn och personnummer, program och inskrivningsår på KTH, kursens namn och kurskod och gärna
en figur.
3. En sida med kort sammanfattning på svenska (ca 5-8 rader) och samma sammanfattning på engelska.
4. En kort beskrivning av de bakomliggande problemen (totalt 2-3 sidor):
För deluppgift 1 och 2, den matematiska pendeln, pendelns ekvation, några
rader om vilka fysikaliska förutsättningar som gäller för att ekvationen ska gälla,
exempel på beräkningar som kan utföras för en pendel
För deluppgift 3 motsvarande för värmeledningsproblemet.
5. Det viktigaste: Beskrivning av de beräkningar som gjorts i uppgifterna.
Namnge metoder som använts, hur algoritmerna utformats, resultat presenterade med numeriska värden och grafer. Grafer ska vara försedda med rubrik,
variabelnamn vid axlarna, eventuell förklarande text i grafen. Grafer kan ingå
insprängda i texten eller med referens till sidor i appendix.
Kommentarer om resultatens rimlighet och noggrannhet. (totalt 2-4 sidor)
6. Eventuella referenser
7. Appendix: MATLAB-program och eventuella grafer.
OBS! INGA MATLAB-program i själva rapporten. MATLAB-program SKA
bifogas som appendix i slutet av rapporten. Använd gärna MATLAB-funktionen
publish.