Transcript Problemlösning 3

Probleml¨osning 3 - Kurvanpassning
Daniel Elfverson
[email protected]
Division of Scientific Computing
Uppsala University
Sweden
1 / 18
Kommentarer p˚
a M2
Miniprojekt 3 - Kurvanpassning
Gillespies algoritm
Gillespies algoritm for stokastiska simuleringar
Algorithm 1 Gillespies algoritm
1: Initialvilkor: s¨
att tillst˚
andsvektorn x = x0, starttid, t = t0 och stopp
tid T
2: while t < T do
3:
Sampla τ , n¨ar n¨asta reaktion intr¨affar
4:
Hitta vilken reaktion som intr¨affar
5:
Uppdatera tillst˚
andsvektorn och tiden
6: end while
Sv˚
arigheter
1
Ok¨ant (varierande) antal tidssteg
2
Hitta vilken reaktion som intr¨affar
3
Plotta l¨
osningen
2 / 18
Kommentarer p˚
a M2
Miniprojekt 3 - Kurvanpassning
Gillespies algoritm
Ok¨ant (varierande) antal tidssteg
• Gissa hur m˚
anga tidssteg och allokera minne
• Gissa hur m˚
anga tidssteg och allokera minne, allokera nytt n¨ar
minnet ¨ar slut
• Best¨
am vid vilka tidpunkter du vill spara l¨
osningen och allokera
minne
3 / 18
Kommentarer p˚
a M2
Miniprojekt 3 - Kurvanpassning
Gillespies algoritm
Gissa hur m˚
anga tidssteg och allokera minne
1
2
3
4
N
t
A
R
=
=
=
=
250000
z e r o s (N, 1 ) ;
z e r o s (N, 1 ) ;
z e r o s (N, 1 ) ;
( Gissning )
4 / 18
Kommentarer p˚
a M2
Miniprojekt 3 - Kurvanpassning
Gillespies algoritm
Allokera nytt n¨ar minnet ¨ar slut
1
2
3
4
5
N = 250000;
NExpand = 2 5 0 0 0 0 ;
t P l o t = z e r o s ( 1 ,N) ;
APlot = z e r o s ( 1 ,N) ;
R P l o t = z e r o s ( 1 ,N) ;
6
7
8
9
10
11
12
13
14
15
16
17
18
19
t = 0
t F i n a l = 200;
n = 0;
while t < tFinal
...
...
n = n + 1;
if n > N
t P l o t = [ t P l o t z e r o s ( 1 , NExpand ) ] ;
...
N = N + NExpand ;
end
end
5 / 18
Kommentarer p˚
a M2
Miniprojekt 3 - Kurvanpassning
Gillespies algoritm
L¨osning sparas vid best¨amda tidpunkter
1
2
3
4
N = 100;
t P l o t = z e r o s ( 1 ,N) ;
APlot = z e r o s ( 1 ,N) ;
R P l o t = z e r o s ( 1 ,N) ;
5
6
7
8
9
10
11
12
t = 0
t F i n a l = 200;
d e l t a T = ( t F i n a l −t ) /N ;
n = 1;
while t < tFinal
...
...
13
14
15
16
17
18
19
i f n∗ d e l t a T <= t
tPlot (n) = t ;
...
n = n + 1;
end
end
6 / 18
Kommentarer p˚
a M2
Miniprojekt 3 - Kurvanpassning
Gillespies algoritm
Hitta vilken reaktion som intr¨affar
• L¨
ostes p˚
a m˚
anga olika s¨att
• Tips p˚
a hur den kan l¨
osas med ’cumsum’ och ’find’ (bra att tr¨ana
p˚
a ’find’ till M3)
1
2
3
4
5
6
7
8
...
u2 = r a n d ( ) ;
w = prop vilar (x) ;
a0 = sum (w) ;
...
F = cumsum (w) ;
r = f i n d ( u2 ∗ a0 <= F , 1 ) ;
...
7 / 18
Kommentarer p˚
a M2
Miniprojekt 3 - Kurvanpassning
Gillespies algoritm
Plotta l¨osningen
• Mystiskt str¨
ack som gick genom plotten (unkown annoying object
UAO)
60
50
40
30
20
10
0
0
1
2
3
50
100
150
200
250
t = [0.5 , 1, . . . ,
200 , 0 , 0 ]
u = [ 48 , 50 , . . . , 1800 , 0 , 0 ]
plot (t , u)
8 / 18
Kommentarer p˚
a M2
Miniprojekt 3 - Kurvanpassning
Miniprojekt 3 - Kurvanpassning
• Uppgift - Hitta ett s¨
asongs beroende av Dikv¨aveoxid, N2 0 (lust gas),
i atmosf¨aren.
• V¨
axthusgas som ¨ar 200 g˚
anger starkare ¨an koldioxid.
• Baseras p˚
a Liao T., Camp C. D., Yung. Y. L.: The seasonal cykle of
N2 0 (se kurshemsida)
9 / 18
Kommentarer p˚
a M2
Miniprojekt 3 - Kurvanpassning
Kurvanpassning - repetition
Olika s¨att att approximera med polynom
• Newton interpolation (ett polynom som g˚
ar genom m¨atpunkterna)
• Minsta kvadratanpassning (ett polynom som ¨
ar ett typ av
medelv¨arde)
• Spline (styckvisa polynom)
10 / 18
Kommentarer p˚
a M2
Miniprojekt 3 - Kurvanpassning
Minsta kvadratanpassning
• Anv¨
ands d˚
a man har en “stor” datam¨angd.
• Anpassa ett 2:a grads polynom, ansats p2 (x) = a0 + a1 x + a2 x 2 .
D¨ar x = [1, 2, 3, 4, 5] kan tex vara tiden f¨
or m¨atningarna och
y = [1, 2, 1, 2, 3] vara m¨atv¨arden, men f˚
ar s˚
a ett ¨
overbest¨amt system
a0 + a1 · 1 + a2 · 12 = 1
a0 + a1 · 2 + a2 · 2
2
a0 + a1 · 3 + a2 · 32
a0 + a1 · 4 + a2 · 42
a0 + a1 · 5 + a2 · 52

1
=2
1

=1⇒ 
1
1
=2
1
=3
1
2
3
4
5

 
12  
1
2
22 
a
0

 
a1  = 1
32 
 

2
42  a2
2
3
5
11 / 18
Kommentarer p˚
a M2
Miniprojekt 3 - Kurvanpassning
• L¨
oses med hj¨alp av normalekvationen AT Aa = b.
• Normalekvationen hittar v s˚
a att kb − Av kL2 minimeras.
L¨
oses i matlab med hj¨alp av
1
p = pol yfit (x , y ,2)
12 / 18
Kommentarer p˚
a M2
Miniprojekt 3 - Kurvanpassning
Avtrendifiering
314
312
310
308
306
304
302
300
298
1976
1978
1980
1982
1984
1986
1988
1990
1992
1994
1996
• Uppgift - Hitta ett s¨
asongs beroende av Dikv¨aveozid, N2 0 (lust gas),
i atmosf¨aren.
• Hitta den ˚
arliga trenden (genom minsta kvadratanpasning) och dra
bort den fr˚
an datan.
13 / 18
Kommentarer p˚
a M2
Miniprojekt 3 - Kurvanpassning
314
312
310
308
306
304
302
300
298
1976
1978
1980
1982
1984
1986
1988
1990
1992
1994
1996
• Data (bl˚
a) med ˚
arlig trend (gr¨
ont).
14 / 18
Kommentarer p˚
a M2
Miniprojekt 3 - Kurvanpassning
Programmerings tips
• load: F¨
or txt filer s˚
a m˚
aste all onumerisk data bort, d˚
a kan
load(’namn.txt’) anv¨
andas.
15 / 18
Kommentarer p˚
a M2
Miniprojekt 3 - Kurvanpassning
A bort o¨onskade data
1
2
x = [1
2 3
4
y = [ 1 0 0 120 0 140
5
6];
0 135];
3
4
i = f i n d ( y ˜= 0 ) ;
5
6
7
xnz = x ( i ) ;
ynz = y ( i ) ;
• Vad ¨
ar xnz och ynz?
16 / 18
Kommentarer p˚
a M2
Miniprojekt 3 - Kurvanpassning
Index- och logiska vektorer
Indexvektorer
Standard vektorn som nu ¨ar vana att jobba med, kan se ut som
x = [1 2 0 1 2 0 9 1 0];
och x(4) = 1.
Logiskavektorer
Element som ¨ar 1 eller 0. Operationen
y = (x ∼= 0) ger en logisk vektor om ¨ar
y = [1 1 0 1 1 0 1 1 0];
Bra f¨
or h¨
ogprestanda ber¨akningar
17 / 18
Kommentarer p˚
a M2
Miniprojekt 3 - Kurvanpassning
Deadline: 7 juni (samma datum som tenta, s˚
a det rekommenderas att
l¨amna in tidigare).
18 / 18