1. No category

TNA005: Tillämpad matematik i teknik och naturvetenskap
VT 2015
UPG6—Miniprojekt 3:
Kastparabler, projektiler, och raketer
Många fysikaliska fenomen inom natur och tillämpningar kan beskrivas och modelleras genom ordinära differentialekvationer (ODE). Det är endast i de enklare fallen som en analytisk lösning till
differentialekvationerna kan bestämmas. Det blir därmed viktigt att kunna lösa ODE-ekvationer
numeriskt. I detta projekt kommer vi att undersöka metoder för numeriska/approximativa lösningar av differentialekvationer. Ekvationerna för modellproblemet vi behandlar kommer från
Newtons kraft- och rörelselagar. Projektarbetet kan delas in i tre områden:
Modellering: Vi börjar med enkla modeller, för vilka analytiska lösningar finns, och går stegvis
över till mer realistiska modeller, vars differentialekvationer måste lösas numeriskt.
Numeriska metoder: Ni kommer att implementera två metoder för ODE-lösare och undersöka
deras noggrannhet. Ni ska även undersöka en av MATLABs inbyggda ODE-lösare.
Experiment: Projektet innehåller även ett antal ingenjörsmässiga uppgifter. Det kommer att
krävas en del experimenterande för att lösa dessa frågor.
1 Tillämpningar av Newtons lagar
Vi har två realistiska(?!) tillämpningscenarier som bakgrund till projektet. I båda fallen gäller
det att undersöka rörelsen av något föremål under inverkan av olika krafter.
Lunchleverans: Företaget Lunch AB vill utveckla en ny och innovativ tjänst för lunchleveranser. Den nya tjänsten går ut på att leverera nylagad varm lunch till byggarbetare på toppen av
en skyskrapa. Lunchen vill man leverera snabbt genom att skjuta upp den med en liten raket.
Er uppgift, som nyanställda på företaget, blir att undersöka möjligheterna för detta med hjälp
av modellering och datorsimulering. Ni ska redovisa resultaten inom 3–4 veckor till företagsledningen.
Snälla grisar: Det är många snälla grisar i svenska gårdar. Det är många snälla fåglar också, men fåglarna är
något nonchalanta och tappar ofta bort sina ägg. När
grisarna hittar något ägg vill de hemskt gärna hjälpa
fåglarna och lämna tillbaka det borttappade ägget. Er
uppgift är att hjälpa grisarna lämna tillbaka äggen. Grisarna vill använda en dräkt med raketmotor, se figur 1,
för att nå höga höjder.
Figur 1: Gris med reketmotor.
Linköpings universitet, ITN, Berkant Savas
1
TNA005: Tillämpad matematik i teknik och naturvetenskap
VT 2015
Huvuduppgift
Miniprojektets huvuduppgift är att lösa rörelseekvationerna för ett föremål som drivs med en
raketmotor. Rörelsen beskrivs av differentialekvationer, och det går inte att lösa dessa analytiskt
utan ekvationerna måste lösas numeriskt. För att kunna bedöma noggrannheten i de numeriska
lösningarna skall vi studera rörelseekvationer med linjärt luftmotstånd, och jämföra analytiska
lösningar med olika numeriska lösningar.
Projektbeskrivningen (detta dokument) innehåller avsnitt med nio deluppgifter som utgör underlag för er att strukturera upp arbetet med miniprojektet. Alla dessa uppgifter skall lösas!
Rapporten och den muntliga redovisningen bygger på att ni har löst dessa uppgifter.
2 Modellering och deluppgifter
I detta avsnitt kommer vi att beskriva olika fysikaliska modeller som beskriver rörelsen av ett
föremål under inverkan av olika krafter. Vi vet att ett föremål som rör sig följer Newtons lagar.
Speciellt intressant för denna uppgift är Newtons andra lag, accelerationslagen,
d(mv)
,
dt
där kraften F, som verkar på ett föremål, är lika med förändringen i rörelsemängd. Rörelsemängden ges av massan gånger hastigheten, mv, och dess förändring blir tidsderivatan d(mv)/dt. Om
vi antar att massan m inte förändras över tiden får vi
F=
d(mv)
dv
=m
= ma,
(1)
dt
dt
där a är föremålets acceleration. Notera att storheterna F, v, och a är vektorstorheter. Vi ser
vidare att rörelseekvationen (1) är en differentialekvation då accelerationen är andraderivatan av
positionen.
F=
Vi kommer att undersöka tre olika fall av rörelsekvationen:
1. Rörelse av ett föremål där luftmotstånd försummas, och inga externa krafter verkar på
föremålet förutom gravitationen.
2. Sedan kommer vi att introducera luftmotstånd, som är proportionell mot hastigheten.
3. Slutligen kommer vi att betrakta ett variabel-mass system med kvadratiskt luftmotstånd.
Detta är ett typiskt fall vid raketuppskjutning eftersom massan av raketen ändras genom
att raketbränsle förbränns och skjuts ut från motorn med hög fart.
2.1 Rörelse utan luftmotstånd, kastparabler
I det enklaste fallet bortser vi från luftmotståndet, och låter endast tyngdkraften verka på föremålet. Vi får
F = ma = mg,
(2)
där g är tyngdaccelerationen. Förenkling av (2) ger a = g och delar vi upp den i x- och ykomponenter får vi
ax
g
0
a = g ⇐⇒
= x =
ay
gy
−g
Linköpings universitet, ITN, Berkant Savas
2
TNA005: Tillämpad matematik i teknik och naturvetenskap
VT 2015
där g ≈ 9.82 m/s2 är vertikalkomponenten av tyngdaccelerationen1 . Skriver vi ekvationerna i
termer av de sökta funktionerna x(t) och y(t) får vi
ax = x00 (t) = 0,
(3)
00
ay = y (t) = −g.
(4)
Analytiska lösningar för (3) och (4) hittas enkelt! Integrera två gånger! I x-led får vi
x00 (t) = 0
⇐⇒
x0 (t) = k1
⇐⇒
x(t) = k1 t + k2 ,
där k1 och k2 är koefficienter från integrationen som kan bestämmas med hjälp av villkor för någon
position och hastighet vid någon tidpunkt. På samma sätt får vi lösningen i y-led enligt
y 00 (t) = −g
⇐⇒
y 0 (t) = −gt + k3
⇐⇒
g
y(t) = − t2 + k3 t + k4 ,
2
där k3 och k4 är nya integrationskoefficienter. Låt oss använda begynnelsevillkoren
x(0) = 0,
y(0) = 0,
x0 (0) = v0 cos(α),
y 0 (0) = v0 sin(α).
(5)
Detta betyder att positionen för kulan vid t = 0 är origo, samt kulan har farten2 v0 m/s i en
riktning som bildar vinkeln α från x-axeln. Även detta vid t = 0. Begynnelsevillkoren (5) ger
k2 = k4 = 0 samt
k1 = v0 cos(α),
k3 = v0 sin(α).
Sätter vi in dessa i lösningarna för x(t) och y(t) fås
x(t) = v0 cos(α)t,
(6)
g
y(t) = v0 sin(α)t − t2 .
2
(7)
Vi kan nu kombinera (6) och (7) och härleda höjden y som en funktion av x. Skriv om (6) som
t=
x
v0 cos(α)
och insättning i (7) ger
y(x) = tan(α)x −
2v02
gx2
.
cos2 (α)
Vi ser att lösningarna blir andragradsfunktioner, som ju beskriver parabler, därav namnet kastparabel. Notera att initialfarten v0 och utgångsvinkeln α är givna konstanter.
2.2 Projektilbana med luftmotstånd
Det finns två förhållandevis enkla sätt att modellera luftmotståndet och därmed hitta mer realistiska lösningar. Dessa är antingen att modellera kraften som linjärt proportionell mot hastigheten
eller som kvadratiskt proportionell mot farten. I båda fallen verkar kraften från luftmotståndet
i motsatt riktning som föremålet rör sig. Vi ska undersöka fallet med linjärt luftmotstånd. Ekvationen blir
ma = mg − kv,
(8)
1 Tyngdacceleration
g varierar beroende på var någonstans man är, och i Sverige (Norrköping) är g = 9.82 en
bra approximation.
2 Vad är det för skillnad mellan fart och hastighet?
Linköpings universitet, ITN, Berkant Savas
3
TNA005: Tillämpad matematik i teknik och naturvetenskap
VT 2015
där v är föremålets hastighet och k är en konstant som modellerar övriga faktorer (t.ex. fluidens
densitet, area och ytform av föremålet) i kraften från luftmotståndet. Vi kan återigen skiva om
ekvationen komponentvis och får
k 0
x (t),
m
k
ay = y 00 (t) = − y 0 (t) − g.
m
ax = x00 (t) = −
(9)
(10)
Notera att differentialekvationerna (9) och (10) är linjära och inte kopplade till varandra. Dessa
kan lösas analytiskt för att erhålla explicita uttryck för x(t), y(t), och även y(x).
Deluppgifter 1
1. Lös ekvationerna (9) och (10) analytiskt så att du får explicita uttryck för x(t) och y(t).
Använd begynnelsevillkoren i (5). Konsultera Forsling och Neymark [2] om det behövs.
2. Implementera lösningarna x(t) och y(t) till (9) respektive (10) i en designerad MATLABfunktion. Låt t, men även modellparametrarna m, k, v0 , och α, vara inparametrar till
funktionen så att det blir enkelt att variera dessa när ni experimenterar.
3. Placera en fyrkantig måltavla med hörn i koordinaterna (40,30), (50,30), (50,40), (40,40).
Låt kulan ha vikten m = 5 kg. Ange utgångsvinkel α och begynnelsefart v0 så att kulan
träffar ovansidan av måltavlan med en negativ hastighet i y-led. Gör detta för tre olika fall
med luftmotstånd, använd k = 0.01 kg/s; k = 0.5 kg/s; och k = 3 kg/s.
4. Differentialekvationerna (9) och (10) är av andra ordningen. Skriv om dessa (på papper)
som ett system av differentialekvationer av första ordningen,
z0 = f1 (t, z).
(11)
Ange komponenterna av z och f1 . Vi kommer att använda detta för att lösa problemet
numeriskt, och jämföra de approximativa lösningarna med de analytiska lösningarna.
2.3 Till rymden! Projektil med raketmotor!!
Ja, i alla fall mot högre höjder! Låt oss montera en raketmotor på kulan! Eftersom raketmotorer
förbränner bränsle, som skjuts ut i någon riktning, kommer massan av raketen inte att vara
konstant, utan den kommer att ändras med tiden, dvs m = m(t). Vi får därmed ett variabelmass system. Genom att använda Newtons tredje lag, impulslagen, kan man visa att följande
ekvation gäller
m(t)a = F + m0 (t)u,
(12)
där m(t)a beskriver totala accelerationskraften på raketen, F anger summan av externa krafter
som verkar på raketen (gravitation och luftmotstånd), och m0 (t)u beskriver kraften raketen
utsätts för genom att partiklar (bränsle) skjuts ut från raketen. Hastighetsvektorn
ux (t)
u = u(t) =
,
uy (t)
är en funktion av t och anger riktning och fart av bränslet som skjuts ut från raketen. Vi kan
anta att u(t) är en känd funktion som piloten (du) kan använda för att styra raketen i en önskad
Linköpings universitet, ITN, Berkant Savas
4
TNA005: Tillämpad matematik i teknik och naturvetenskap
VT 2015
bana. Vi förutsätter även att m(t) och dess derivata m0 (t) är kända funktioner. Kombinerar vi
ekvation (12) med krafter för gravitation och luftmotstånd får vi
m(t)a = m(t)g − ckvkv + m0 (t)u,
(13)
där c är en ny proportionalitetskonstant. Notera att vi, i detta fall, modellerar luftmotståndet som
kvadratiskt proportionell mot farten. Division med m(t) och komponentvis uppdelning ger
m0 (t)
c p 0 2
(x ) + (y 0 )2 x0 +
ux (t),
(14)
m(t)
m(t)
m0 (t)
c p 0 2
(x ) + (y 0 )2 y 0 +
ay = y 00 (t) = −
uy (t) − g,
(15)
m(t)
m(t)
p
där vi använder kvk = (x0 )2 + (y 0 )2 . Vi ser att ekvationerna är kopplade till varandra på ett
ickelinjärt sätt. Jämför även (14) och (15) med ekvationerna i Exempel 14.7 i Jönsson [3].
ax = x00 (t) = −
Deluppgifter 2
5. Skriv om (14) och (15) som ett ickelinjärt system av första ordningens ODE
z0 = f2 (t, z).
(16)
Ange komponenterna av z och f2 . Vi kommer att använda dessa längre fram.
2.4 Numerisk lösning av differentialekvationer
I projektet ska du undersöka lösningar till differentialekvationerna som har introducerats. Du ska
speciellt jämföra exakta (analytiska) lösningar, där dessa finns, med olika numeriska lösningar.
Betrakta en ODE i formen
z0 = f (t, z),
där z(t) är en vektorvärd funktion som söks och f är en given vektorvärd funktion som kan
evalueras för givet t och z(t). Den numeriska lösningsproceduren för differentialekvationer är
följande. Bestäm ekvidistanta punkter t0 , . . . , tn som delar in aktuell tidsintervall [a b] i n stycken delintervall med längden h = (b − a)/n. Numeriska approximationer för z(t1 ), . . . , z(tn ) fås
rekursivt, där vi använder ett givet begynnelsevillkor z(t0 ) = z0 . Du skall själv implementera två
metoder: Eulers method och fjärde ordningens Runge-Kutta metod (RK-4). En beskrivning av
dessa metoder ges i Algoritm 1 och 2, se [1] för mer information.
Algoritm 1 Eulers metod
for i = 0, 1, 2, . . . , n − 1 do
zi+1 = zi + hf (ti , zi )
end for
Linköpings universitet, ITN, Berkant Savas
5
TNA005: Tillämpad matematik i teknik och naturvetenskap
VT 2015
Algoritm 2 Runge-Kutta ordning 4
for i = 0, 1, 2, . . . , n − 1 do
k1 = hf (ti , zi )
k2 = hf (ti + h/2, zi + k1 /2)
k3 = hf (ti + h/2, zi + k2 /2)
k4 = hf (ti + h, zi + k3 )
zi+1 = zi + 16 (k1 + 2k2 + 2k3 + k4 )
end for
Deluppgifter 3
6. Ta fram uttryck för f1 och f2 i ODE problemen (11) resp. (16). Du har bestämt dessa i
tidigare uppgifter. Implementera f1 (t, z) och f2 (t, z) i lämpliga MATLAB-funktioner, t.ex.
funk1.m och funk2.m. Varje funktionsfil tar en tidsvariabel t (skalär) och en kolumnvektor
z som inparametrar och returnerar evalueringen av den aktuella funktionen i t och z. Notera
att f1 , f2 , och z måste ha samma dimension. Funktionerna skall kunna anropas genom
>> dz1 = funk1 (t , z );
>> dz2 = funk2 (t , z );
Utöver t och z, behöver funktionerna ha tillgång till parametrar som massan m. Denna
gång definiera dessa parametrar internt i funktionerna. I fallet med raketmotor behöver
du komma åt en massfunktion m(t), dess derivata m0 (t), och en hastighetsfunktion u(t).
Det är lämpligt att implementera dessa funktioner i egna filer och få ut aktuella storheter
genom funktionsanrop. Mer om detta i kommande uppgift.
7. Implementera Eulers metod och RK-4 metoden i MATLAB. Låt förslagsvis funktionerna
vara odeSolverEuler.m och odeSolverRk4.m med fyra inparametrar som ges nedan:
fun
z0
tspan
n
funktionshandtag som anger funktionen för f (t, z);
vektorn z0 anger begynnelsevillkor för ODEn;
anger tstart och tslut , gränser för tidsintervallet, t.ex. tspan = [0,20];
n anger antalet diskretiseringspunkter.
Låt vidare varje funktion returnera en tidsvektor t och en 4×(n+1) matris y med approximationer av den sökta lösningen. Efter implementeringen av ODE-lösarna och initiering av
relevanta variabler skall ni kunna anropa funktionerna enligt följande för att få ut respektive
approximativ lösning.
>> [ t1 , y1 ] = odeSolverEuler ( @fun , tspan , z0 , n );
>> [ t2 , y2 ] = odeSolverRk4 ( @fun , tspan , z0 , n );
8. Du ska i denna uppgift jämföra den analytiska lösningen y(t) av (10) (du har explicit uttryck
för detta sedan tidigare) med tre stycken numeriska lösningar. Använd y(t) från den tidigare
uppgiften där parametrarna (m, k, v0 , α) är sådana att lösningen prickar måltavlan. Låt
n ∈ [20, 200]. Beräkna tre stycken approximativa lösningar för samma differentialekvation
genom att använda f1 i (11). De numeriska lösningarna skall du beräkna genom att använda
din egen Euler och RK-4 metod, samt MATLABs inbyggda ODE-lösare ode45.
9. I den sista uppgiften skall du skjuta iväg raketen och styra den så att den träffar en måltavla, som är betydligt högre upp än tidigare. Hörnen på måltavlan har nu koordinaterna:
Linköpings universitet, ITN, Berkant Savas
6
TNA005: Tillämpad matematik i teknik och naturvetenskap
VT 2015
(60,130), (70,130), (70,140), (60,140).
Vi sätter initiala massan av raketen till m(0) = 5 kg, varav 1 kg är bränsle. Under tiden
som motorn är igång förbränner den 200 g bränsle per sekund. Detta gör att motorn kan
vara igång som mest 5 sekunder. Om vi startar raketmotorn vid t = 0 och håller motorn
igång tills bränslet är slut kommer massfunktionen av raketen, och dess derivata, att se ut
enligt följande
(
(
5 − 0.2t om 0 ≤ t ≤ 5,
−0.2 om 0 ≤ t ≤ 5,
m(t) =
m0 (t) =
4
om 5 < t,
0
om 5 < t.
Det är smidigt att implementera m(t) och m0 (t) i en egen funktion som anropas i f2 (t, z) för
att få aktuell massa och derivata vid en given tidpunkt. Du behöver även specificera u(t)
i ekvation (13). Låt masspartiklarna ha en konstant fart km ut från motorn, men motorns
riktning θ(t) kan du styra fritt. Detta innebär
ux (t)
km cos(θ(t))
ku(t)k = km , samt u(t) =
=
.
uy (t)
km sin(θ(t))
Fysikaliskt innebär detta att raketmotorn går på full kraft. Utgå ifrån att θ anger vinkeln
från positiva x-axeln. För att få en kraft uppåt måste partiklarna skjutas rakt ner. Detta
ger vinkeln θ = −π/2 (eller θ = 3π/2 om du vill).
I tidigare uppgifter har du skjutit iväg kulan, genom att ange en viss begynnelsefart v0 . I
denna uppgift skall du låta raketen röra sig av egen kraft, du sätter därmed v0 = 0. Låt
c = 0.05 kg/m och km = 700 m/s.
Raketen måste lyfta vertikalt, θ(0) = −π/2, och komma upp till en viss höjd innan vinkeln
θ(t) får börja ändras. Din uppgift är att experimentellt bestämma en funktion θ(t) som tar
raketen till ovansidan av målet. Använd din egen RK-4 metod med n = 200 när du löser
ODE-ekvationen numeriskt.
3 Skriftlig och muntlig redovisning
Uppgifterna ovan skall redovisas i en skriftlig teknisk rapport om max 2 sidor. Rapporten måste
inkludera följande:
1. Uttryck för x(t) och y(t) som löser ODE-ekvationerna (9) resp. (10). Använd begynnelsevillkoren i (5).
2. Figur som visar tre stycken grafer (y som funktion x och ) där du träffar ovansidan av måltavlan. En graf för vart och ett av fallen med olika luftmotstånd (se uppgift 3). Presentera
all nödvändig information så att det går att reproducera resultaten. Diskutera resultat och
förutsättningar.
3. Jämförelse mellan den analytiska lösningen y(t) och motsvarande numeriska lösningar från
egenimplementerade ODE-lösare och MATLABs ode45. Plotta absolutbeloppet av felet
med semilogy för att tydligt kunna se skillnader mellan lösningarna. Diskutera resultat
och förutsättningar.
4. En figur som visar raketens flygbana där du återigen prickar ovansidan av den högre måltavlan. Visa även en plott av motorns riktning θ(t). Diskutera och relatera olika skeden av
θ(t) till hur raketens bana påverkas.
Linköpings universitet, ITN, Berkant Savas
7
TNA005: Tillämpad matematik i teknik och naturvetenskap
VT 2015
5. Diskussion och reflektion med utgångspunkt i ett av tillämpningsscenarierna om hur du
skulle kunna gå vidare med projektet.
OBSERVERA: Denna gång blir det ingen komplettering. Det är därför viktigt att rapporten
redan vid den första inlämningen är korrekt och fullständig. Använd en lämplig disposition likt
tidigare miniprojekt; krav på reproducerbarhet gäller fortfarande; använd källor och referera i
texten till alla källor; referera till alla figurer i texten; låt alla figurer ha en beskrivande figurtext,
osv, osv.
Den muntliga redovisningen skall innehålla resultat och diskussion av analysen och experimenten
som du har utfört. Visualisera flitigt dina resultat i olika figurer med hjälp av dator.
Referenser
[1] L. Eldén and L. Wittmeyer-Koch. Numeriska beräkningar – analys och illustrationer med
MATLAB. Studentlitteratur, 4 edition, 2001.
[2] G. Forsling and M. Neymark. Matematisk analys – en variabel. Liber, 2011.
[3] P. Jönsson. MATLAB-beräkningar inom teknik och naturvetenskap. Studentlitteratur, 2010.
Linköpings universitet, ITN, Berkant Savas
8