1. No category

Oblig 2 i MAT-INF1100
Tor Hedin Brønner (torhbr)
1
Oppgave 1.
a) from numpy import ∗
def d(x ) :
””” Takes an array , x , o f l e n g t h n and p r o d u c e s
an array , dx , o f l e n g t h n−1 c o n s i s t i n g o f t he s i g n e d
i n t e r v a l l e n g t h s dx [ i ] = x [ i +1] − x [ i ] where i r a n g e s
from 0 t o n − 1 .
”””
return x [ 1 : ] − x[: −1]
d e f D( x , y ) : # u s i n g newton−q u o t i e n t s
dy = z e r o s ( l e n ( x ) )
dy [ : − 1 ] = d ( y ) / d ( x ) # f o r w a r d q u o t i e n t s a t x [ : − 1 ]
dy [ −1] = dy [ −2]
r e t u r n dy
Her får vi strengt tatt ikke en funksjon som kan ta en vilkårlig t. Vi
kan lage en slik funskjon ved å tilnærme f (t) når t ∈ (ti , ti+1 ) med yi +
yi+1 −yi
(t − ti ) (dvs. strekke rette linjer mellom punktene). En python
ti+1 −ti
implementasjon kan se slik ut:
d e f make f ( x , y ) :
dy = d ( y ) / d ( x )
def f ( t ) :
f o r i i n range ( l e n ( x ) ) :
i f t <= x [ i + 1 ] :
1
break
r e t u r n y [ i ] + dy [ i ] ∗ ( t − x [ i ] )
return f
Når vi plotter D(x, y) på gps-sporet får vi:
a
1.5
1.0
0.5
0.0
0.5
0
1000 2000 3000 4000 5000 6000
t
b) d e f I ( x , y ) : # Trapezoid method
tra = d(x )∗( y[: −1] + d(y )/2)
i = zeros ( len (x )) # i [ 0 ] = 0
f o r j i n range ( l e n ( t r a ) ) :
i [ j +1] = i [ j ] + t r a [ j ]
return i
Når vi plotter denne får vi:
2
2.0
×10
4
s
1.5
1.0
0.5
0.0 0
1000 2000 3000 4000 5000 6000
t
c) Kort skisse over koden til plottingen (se vedlegg for full kode):
from m a t p l o t l i b . pylab import ∗
f = open ( ’ running . txt ’ )
x , y = f l o a t ( [ l . s p l i t ( ’ , ’) for l in f ] ) .T # vector f l o a t
a , s = D( x , y ) , I ( x , y )
fg , ax = s u b p l o t s ( )
ax . p l o t ( x , a )
f g . s a v e f i g ( ’ a . pdf ’ )
fg , ax = s u b p l o t s ( )
ax . p l o t ( x , s )
f g . s a v e f i g ( ’ s . pdf ’ )
Som vi ser på grafen på grafene ovenfor er det to perioder, en lang og en
kort, der a-grafen har få punkter og s-grafen er slakere.
Hvis vi ser på running.txt ser vi dette:
4063 ,3.1667
4 56 3 , 0
4 56 4 , 0
4567 ,1.75
4569 ,2.1667
3
Gps-sporet stopper altså brått (ingen nedakselerasjon) og er borte i lang
tid. Her blir approksimasjonene av akselerasjon og strekning ubrukelige.
Det er snakk om en dx ≈ 500, i forhold til en median på dx = 6. Spesielt
for s-grafen får vi problemer, man får en slak linje, som blir helt feil hvis
løperen faktisk stoppet.
Alt i alt er dette et problem med data-grunnlaget. Vi må dermed fikse
det ved å analysere situasjonen manuelt. Vi ser at løperen er stille når
sporet starter opp igjen. Vi kan kanskje anta at hun stoppet gpsen for
så å stoppe selv. I så fall kan vi dele opp hele intervalet slik at vi ikke
inkluderer de intervalene der vi ikke har data. Hvis vi gjør dette med s
får vi:
20000
15000
10000
5000
00
1000 2000 3000 4000 5000 6000 7000
Det er også mulig at gpsen stoppet av seg selv. Løperen ser dermed ut til
å ha stoppet for å fikse problemet. I så fall kan vi anta at løperen bare
stoppet så lenge det tok å få i gang gpsen igjen. I dette tilfellet kunne
vi legge på noen punkter med snitthastigheten. Det ville nok vært lurt å
spørre løperen hva som skjedde, hvis man var interessert i få en nøyaktig
løpelengde.
4
2
Oppgave 2.
a) x0 − x2 = 1 ⇒
x0
=1⇒
1 + x2
Z
Z
x0
dt =
1 + x2
Z
dt ⇒
1
dx = t + C ⇒ arctan x = t + C ⇒
1 + x2
x = tan (t + C)
x(0) = tan (C) = 1 ⇒ C = π4 . Dermed er x(t) = tan (t + π4 ).
b) Her plotter vi en tilnærming til x(t) ved å bruke Eulers metode:
5
4
tan(t + π/4)
euler
3
2
10.0
0.1
0.2
0.3
0.4
0.5
Her har vi brukt denne implementasjonen av Eulers metode, der x0 = f (x)
bare avhenger x:
def euler ( f , x 0 , I , n ) :
dt = ( I [ 1 ] − I [ 0 ] ) / f l o a t ( n )
t s = arange ( I [ 0 ] , I [ 1 ] , dt )
x = z e r o s ( n+1)
x [0]= x 0
f o r i i n range ( n ) :
x [ i +1] = x [ i ] + f ( x [ i ] ) ∗ dt
return ts , x
Plottingen er gjort slik (se kode vedlegg for all kode):
5
et , ex = e u l e r ( lambda x : 1 + x ∗∗2 , 1 , ( 0 , 0 . 6 ) , 6 )
fg , a x e s = s u b p l o t s ( f i g s i z e =(5 , 2 . 5 ) )
ts = linspace (0 ,0.6)
a x e s . p l o t ( t s , tan ( p i /4 + a r r a y ( t s ) ) , l a b e l=r ’ $tan ( t+\p i / 4) $ ’ )
a x e s . p l o t ( et , ex , l a b e l=r ’ $ e u l e r $ ’ )
axes . a x i s ( ’ tight ’ )
a x e s . l e g e n d ( l o c =2)
c) Hvis vi bruker Eulers midtpunktmetode på x(t) får vi dette plottet:
5
4
3
tan(t + π/4)
euler
euler_mid
2
10.0
0.1
0.2
0.3
0.4
Her har vi brukt en lignende implementasjon som i b):
def euler mid ( f , x 0 , I , n ) :
dt = ( I [ 1 ] − I [ 0 ] ) / f l o a t ( n )
t s = [ I [ 0 ] + i ∗ dt f o r i i n range ( n +1)]
x = z e r o s ( n+1)
x [0]= x 0
f o r i i n range ( n ) :
x mid = x [ i ] + f ( x [ i ] ) ∗ dt /2
x [ i +1] = x [ i ] + f ( x mid )∗ dt
return ts , x
Plottingen er gjort på samme måte som i b).
d) Ved å plotte den alternative metoden får vi:
6
0.5
5
4
3
tan(t + π/4)
euler
euler_mid
mid
2
10.0
0.1
0.2
0.3
0.4
0.5
Eulers midtpunktmetode ser ut til å fungere litt bedre enn den alternative
metoden, men det er ganske jevnt. Eulers metode er klart dårligst. Den
alternative metoden er implementert slik:
d e f mid ( f , x 0 , I , n ) :
dt = ( I [ 1 ] − I [ 0 ] ) / f l o a t ( n )
t s = [ I [ 0 ] + i ∗ dt f o r i i n range ( n +1)]
x = z e r o s ( n+1)
x [0]= x 0
f o r i i n range ( n ) :
x [ i +1] = ( 2 − dt ∗x [ i ] − 2∗ s q r t ( 1 − dt ∗∗2 −2∗x [ i ] ∗ dt ) ) / dt
return ts , x
1 − h2 − 2xk h kan ikke være mindre enn 0. Når xk blir større må h bli
mindre. Hvis xk blir veldig stor kan det være at vi må bruke upraktiske
små h.
Løser vi denne ligningen får vi:
p
p
−2xk ± 2 x2k + 1
−2xk ± 4x2k + 4
=
h=
2
2
q
= −xk ±
0<
x2k + 1
p
x2k + 1 − xk → 0 når xk → ∞, så
7