L¨ osningar till Tentamen i Ber¨ akningsvetenskap II, 5.0 hp, 2010-03-15
Del A
1. (a) Till¨ampa Eulers fram˚atdifferensmetod p˚a ekvationen y0(t) − cos y(t)(y(t) + 1) = 0, y(0) = 0.
Ta ett steg med stegl¨angd h = 0.1. Vi b¨orjar med att st¨alla upp ekvationen p˚a formen
y0(t) = f (t, y) = cos y(t)(y(t) + 1), y(0) = 0.
Eulers fram˚atdifferensmetod ges d˚a av
yi+1 = yi+ hf (ti, yi) = yi+ h cos yi(yi+ 1), y0 = 0.
Ett steg:
y1 = y0+ h cos y0(y0+ 1) = 0 + 0.1 cos(0)(0 + 1) = 0.1.
(b) Explicita metoder kr¨aver mindre arbete i varje tidssteg ¨an implicita metoder. Om stegl¨angden f¨or att uppfylla ett visst noggrannhetskrav
¨
ar mindre ¨an stegl¨angden f¨or att uppfylla stabilitetskravet kan det d¨arf¨or l¨ona sig med en explicit metod. Detta vore typiskt fallet f¨or icke-styva problem, d¨ar stabilitetskravet inte ¨ar lika strikt som f¨or styva problem.
2. Stabilitetsanalys f¨or trapetsmetoden. Stabilitet analyseras med hj¨alp av testekvationen, y0(t) = λy(t). Trapetsmetoden :
yi+1= yi+ h
2(f (ti, yi) + f (ti+1, yi+1)). (1) Ins¨attning av testekvationen i (1) ger
yi+1 = yi+h
2(λyi+ λyi+1) ⇒ (1 − λh
2 )yi+1 = (1 +λh 2 )yi. Stabilitetsvillkoret blir d˚a
|1 + λh2 |
|1 −λh2 | ≤ 1,
⇒ (1 + h
2<(λ))2+ (h
2=(λ))2 ≤ (1 − h
2<(λ))2+ (h
2=(λ))2.
Eftersom vi betraktar l¨osningar till testekvationen som ¨ar avtagande, dvs vi har <(λ) < 0, f˚ar vi
⇒ (1 + h
2<(λ))2+ (h
2=(λ))2 = 1 + h<(λ) + h2
4 <(λ)2+h2 4 =(λ)2
≤ 1 − h<(λ) + h2
4<(λ)2+ h2
4 =(λ)2 = (1 − h
2<(λ))2+ (h
2=(λ))2.
Stabilitetskravet ¨ar uppfyllt f¨or alla h och alla λ med <(λ) < 0, dvs trapetsme- toden ¨ar ovillkorligt stabil.
3. Formulera en Monte Carlo-metod i Matlab, given en funktion sim(x) som utf¨or en stokastisk simulering.
N = input(’Antal simuleringar: ’); % L˚at anv¨andaren ange antal simuleringar x = input(’Ange indata: ’); % L˚at anv¨andaren ange indata
Y = zeros(1,N); % Skapa vektor f¨or att spara
% l¨osningen for i=1:N
y = sim(x); % utf¨or N st stokastiska simuleringar Y(i) = y; % spara resultatet till varje simulering end
medel = mean(Y); % ber¨akna medelv¨ardet
4. err(i) uppskattade felet n¨ar N(i) ber¨akningspunkter anv¨andes.
(a) Formulera om ansatsen f el(N ) = cNp, sa att c och p kan best¨ammas med minstakvadratanpassning med r¨at linje.
Logaritmering av ansatsen ger
log(f el(N )) = log(cNp) = log(c) + log(Np) = log(c) + plog(N ).
c och p kan d˚a best¨ammas genom att till¨ampa minstakvadratanpass- ning med r¨at linje p˚a
f el = ˜˜ c + p ˜N ,
d¨ar ˜f el = log(f el), ˜c = log(c) och ˜N = log(N ).
(b) Matlab-kommandot som genomf¨or minstakvadratanpassningen enligt ovan ¨ar polyfit och anropas:
px = polyfit(log(N),log(err),1); % px = [p log(c)]
I det allm¨anna fallet ges minstakvadratanpassning med r¨at linje till givna data x och y av
px = polyfit(x,y,1);
5.
(a) Initialv¨ardesproblem/Begynnelsev¨ardesproblem.
(b) Att metoden ¨ar stabil.
(c) u ska vara slumpad ur sannolikhetsf¨ordelningen som ¨ar uniformt f¨ordelad i intervallet [0, 1].
(d) Normalekvationerna.
Del B
6. Kemiskt experiment upprepat fem g˚anger, volymen V och antal mol n uppm¨att.
V n
1 0.067 0.0027 2 0.098 0.0038 3 0.051 0.0021 4 0.056 0.0024 5 0.062 0.0024
Temperaturen var T = 294 K, trycket p = 103.8 kPa. Anv¨and data och allm¨anna gaslagen, pV = nRT , till att best¨amma gaskonstanten, R. Alla data ska utnyttjas. Detta utesluter interpolation eftersom det finns tv˚a v¨arden p˚a V f¨or n = 0.0024, och dessa inte kan uppfyllas samtidigt. Med minstakvadratanpassning med r¨at linje kan v¨ardet p˚a R best¨ammas ur
pV = nRT, (2)
med T = 294, p = 103.8, dvs ett linj¨art samband med avseende p˚a R
d¨ar konstanten ska vara noll. Problemuppst¨allningen med det ¨overbest¨amda systemet (2) blir
AR = b, d¨ar b = pV =
6.9546 10.1724
5.2938 5.8128 6.4356
, A = nT =
0.7861 1.1247 0.6289 0.7041 0.7041
.
R best¨ams ur normalekvationerna
ATAR = ATb, som h¨ar blir
ATA = 3.2625, ATb = 28.8288.
D˚a normalekvationerna l¨oses f˚as ett approximativt v¨arde p˚a gaskonstan- ten R = 8.84 (8.8365).
7. Ett steg i Gillespies algoritm slumpar fram numret p˚a n¨asta reaktion.
(a) Visa med handr¨akning hur det steget g˚ar till, f¨or exemplet u = 0.371, w = [3 8 2 14 5].
Den kumulativa f¨ordelningsfunktionen ges av F (w) = [
1
X
i=1
wi
2
X
i=1
wi
3
X
i=1
wi
4
X
i=1
wi
5
X
i=1
wi] = [3 11 13 27 32].
Den normaliserade kumulativa f¨ordelningsfunktionen ˜F = P5F
i=1wi blir F (w) = [0.0938 0.3438 0.4062 0.8438 1.0000].˜
N¨asta reaktion r ges av
F (w = w(r − 1)) < u ≤ ˜˜ F (w = w(r)), och i detta fall
0.3438 = ˜F (w = w(2)) < 0.371 ≤ ˜F (w = w(3)) = 0.4062,
⇒ r = 3.
(b) Skriv Matlabfunktion r = next reaction(w,u) som tar fram n¨asta reaktion med Gillespies algoritm.
function r = next_reaction(w,u) C = cumsum(w);
r = find(C >= u*C(end),1);
8. Best¨am v¨ardena f¨or Heuns metod utg˚aende fr˚an v¨ardena f¨or Eulers metod i tabellen nedan och skriv in dessa i tabellen. h ¨ar stegl¨angd, te total ber¨akningstid, n antal tidssteg och ts ber¨akningstid per steg.
Utr¨aknade v¨arden ¨ar insatta i tabellen med fet stil.
h te n ts
Eulers metod: 0.0001604 0.1075398 62364 0.0000017 Heuns Metod: 0.0421138 0.0008058 237 0.0000034
Sluttiden T best¨ams som antalet tidssteg g˚anger stegl¨angd f¨or Eulers metod till
T = h ∗ n = 0.0001604 ∗ 62364 ≈ 10.0032,
dvs simuleringer k¨ors till T = 10. D˚a stegl¨angden f¨or Heuns metod ¨ar k¨and kan antalet tidssteg best¨ammas enligt
nHeun = T
hHeun = 10
0.0421138 ≈ 237.
Vidare kan man anta att Heun kr¨aver dubbelt s˚a mycket arbete som Euler i varje steg, eftersom tv˚a funktionsanrop g¨ors i varje steg j¨amf¨ort med Eulers enda funktionsanrop. Detta leder till att
ts,Heun= 2 ∗ ts,Euler = 2 ∗ 0.0000017 = 0.0000034.
Slutligen kommer den totala ber¨akningstiden uppg˚a till antal tidssteg multi- plicerat med tids˚atgangen f¨or varje steg
te,Heun= nHeun∗ ts,Heun= 237 ∗ 0.0000034 ≈ 0.0008058.
