En slumpvandringsmetod för värmeledningsekvationen med rörlig rand

Institutionen för naturvetenskap och teknik

En slumpvandringsmetod för

värmeledningsekvationen med

rörlig rand

Andreas Lockby

Örebro universitet

Institutionen för naturvetenskap och teknik

Självständigt arbete för kandidatexamen i matematik, 15 hp

En slumpvandringsmetod för

värmeledningsekvationen med rörlig rand

Andreas Lockby 8 juni 2016

Handledare: Magnus Ögren Examinator: Mårten Gulliksson

Sammanfattning

I detta arbete har målet varit att numeriskt beräkna en lösning till värmeled-ningsekvationen med rörlig rand med en slumpvandringsmetod. Analytiska lösningar samt lösningar med andra numeriska metoder nns det gott om, det stora problemet ligger i den rörliga randen. Första halvan behandlar en del teori kring värmeledningsekvationen, i synnerhet lösningar med Di-richletvillkor. I den andra halvan behandlas slumpvandringsmetoden för den homogena värmeledningsekvationen. Metoden appliceras först på enkla ex-empel, för att sedan lösa Stefanproblemet. I slutet av arbetet görs en kort felanalys för den numeriska lösningen av Stefanproblemet.

Detta arbete har, som sagt, behandlat en lösning via en slumpprocess till ett exempel på värmeledningsekvationen med rörlig rand. Numeriska meto-der genom slumpprocesser för att lösa partiella dierentialekvationer med just rörlig rand är något som det idag forskas på, varpå detta arbete ligger i linje med vad som är aktuellt för ämnet. Exempelvis så behandlar [1] en del kring vad som är aktuellt för s.k. free boundary problems och i [2] diskuteras det aktuella metoder för att utvidga slumpprocesser till ickelinjära problem.

Innehåll

1 Värmeledningsekvationen 5

1.1 Formulering med Dirichletvillkor . . . 5

1.2 Det homogena fallet . . . 6

1.2.1 Homogena randvillkor . . . 6

1.2.2 Inhomogena randvillkor . . . 8

1.3 Det inhomogena fallet . . . 10

1.3.1 Denitioner och satser . . . 10

1.3.2 Lösningar med Greenfunktion, inhomogena fallet . . . 13

2 Slumpvandring 16 2.1 Metod . . . 16

2.1.1 Begreppet slumpvandring . . . 16

2.1.2 Koppling till värmeledningsekvationen . . . 17

2.2 Numerisk lösning, homogena randvillkor . . . 19

2.3 Resultat . . . 21

3 Utvidgning till rörlig rand 23 3.1 Formulering . . . 23

3.2 Slumpvandring . . . 24

3.2.1 Metod för rörlig rand . . . 24

3.2.2 Stefanproblemet, något att jämföra med . . . 24

3.2.3 Numerisk beräkning av Stefanproblemet . . . 25

3.3 Felanalys för Stefanproblemet . . . 27

3.4 Resultat, testproblem . . . 29

3.5 Resultat, Stefanproblemet . . . 30

Kapitel 1

Värmeledningsekvationen

1.1 Formulering med Dirichletvillkor

Detta arbete syftar till att undersöka hur man numeriskt kan lösa värme-ledningsekvationen då området förändras med tiden. Så vad är då värmeled-ningsekvationen? Vi kommer här endast att betrakta fallet för en rumsdi-mension, så vi låter T (x, t) beteckna temperaturen i punkten x vid tiden t. I det homogena fallet, då ingen värme tillförs till systemet, gäller följande samband

∂T ∂t − D ·

∂2T

∂x2 = 0, 0 < x < L, t > 0 (1.1)

där D är en konstant som bestäms utifrån det gällande materialets fysikaliska egenskaper, ofta kallad diusionskonstant.

Till detta hör vanligtvis ett begynnelsevillkor

T (x, 0) = f (x) (1.2)

och randvillkor

T (0, t) = h(t)

T (L, t) = g(t). (1.3)

Dessa villkor utgör s.k. Dirichletvillkor. När värme tillförs till systemet (det inhomogena fallet), kan ekvation (1.1) skrivas som

∂T ∂t − D ·

∂2T

∂x2 = v(x, t), 0 < x < L, t > 0 (1.4)

där v(x, t) betecknar den tillförda värmen i punkten x vid tiden t, mätt i [K/s](Kelvin per sekund).

1.2 Det homogena fallet

1.2.1 Homogena randvillkor

För att kunna avgöra huruvida den numeriskt beräknade lösningen är bra, vill vi ha en analytisk lösning att jämföra med. Vi betraktar därför ett enkelt fall av den homogena värmeledningsekvationen (1.1) med begynnelsevillkor (1.2) och randvillkor (1.3) enligt följande

f (x) = 1, h(t) = 0, g(t) = 0 (1.5) vilket leder till att vi vill lösa följande partiella dierentialekvation med tillhörande begynnelse-och randvillkor

∂T ∂t = D · ∂2T ∂x2, 0 < x < L, t > 0 T (x, 0) = 1 T (0, t) = T (L, t) = 0. (1.6)

För att lösa detta, betraktar vi en andra ordningens dierentialoperator A enligt

A = D ∂

2

∂x2 (1.7)

där D är diusionskonstanten. Vi kan då skriva vår värmeledningsekvation från (1.6) som

∂T (x, t)

∂t = AT (x, t). (1.8)

Vi kan sedan betrakta detta som ett egenvärdesproblem, ∂T (x, t)

∂t = AT (x, t) = λT (x, t), (1.9) som har lösningen

T (x, t) =

∞

X

n=1

cnξn(x)eλnt, (1.10)

där cn är konstanter, ξn(x) och λn är egenfunktioner respektive egenvärden

till operatorn A. Med lite beräkningar kan man visa att dessa egenvärden och egenfunktioner är λn= − n2π2D L2 , ξn(x) = sin nπx L  .

Det som återstår är att bestämma konstanterna cn. Om vi nu utnyttjar begynnelsevillkoret T (x, 0) = 1 får vi att T (x, 0) = 1 = ∞ X n=1 cn· sin nπx L  sinmπx L  = ∞ X n=1 cn· sin nπx L  sinmπx L  Z L 0 sinmπx L  dx = Z L 0 ∞ X n=1 cn· sin nπx L  sinmπx L  dx (1.11)

Integralen i högerledet kommer att vara 0 när n 6= m, p.g.a. att de två sinus-funktionerna är ortogonala på intervallet [0, L]. Det kan enkelt kontrolleras genom att veriera att

Z L 0 sinnπx L  sinmπx L  dx = 0, n 6= m. (1.12) Då har vi att ekvation (1.11) blir

Z L 0 sinmπx L  dx = Z L 0 cm· sin2 mπx L  dx. (1.13)

Vidare, genom att använda likheten

sin2mπx L  = 1 − cos 2mπx L
2 (1.14) kan vi beräkna Z L 0 sin2mπx L  dx = cm· L 2, (1.15)

så att vi, från (1.13) får att

cm = 2 L Z L 0 sinmπx L  dx = ( 0 om m är jämn. 4 mπ om m är udda.

Då får vi från (1.10) att (vi byter tillbaka från m till n av notationsskäl) vår lösning ges av T (x, t) = 4 π X n 1 nsin nπx L  e −n2π2D L2 t, n = 1, 3, 5, . . . (1.16)

I praktiken kommer vi att approximera den exakta lösningen med ett be-gränsat antal termer i summan, när vi sedan jämför med lösningen vi får med slumpvandringsmetoden.

1.2.2 Inhomogena randvillkor

Med inhomogena randvillkor menas, från (1.3), att

h(t) 6= 0 och/eller g(t) 6= 0. (1.17) Om vi sätter h(t) = A, g(t) = B där A, B är konstanter, får vi följande,

∂T ∂t = D · ∂2T ∂x2, 0 < x < L, t > 0 T (x, 0) = f (x) T (0, t) = A, T (L, t) = B. (1.18)

För att lösa detta med metod hämtad ur [10], inför vi en tidsoberoende (stationär) funktion Tstat(x) som uppfyller

T_stat00 (x) = 0 Tstat(0) = A

Tstat(L) = B,

(1.19)

villket har den enkla lösningen

Tstat(x) = A + B − A L x. (1.20) Om vi nu bildar u(x, t) = T (x, t) − Tstat(x) (1.21) får vi att ∂u ∂t = ∂T ∂t, ∂2u ∂x2 = ∂2T ∂x2

u(0, t) = T (0, t) − Tstat(0) = 0 = T (L, t) − Tstat(L) = u(L, t)

u(x, 0) = T (x, 0) − Tstat(x) = f (x) − Tstat(x) = ˜f (x).

(1.22)

För u(x, t) gäller nu alltså ∂u ∂t = D · ∂2u ∂x2, 0 < x < L, t > 0 u(x, 0) = ˜f (x) u(0, t) = 0 u(L, t) = 0. (1.23)

Hur man löser denna typ av system har vi behandlat i avsnitt 1.2.1 med ˜

f (x) = 1och vi kan nu bilda lösningen till T (x, t) som

I fallet då randvillkoren inte är konstanta, d.v.s. A = A(t), B = B(t) kan man på liknande sätt bilda u(x, t) som i (1.2.2), men i detta fallet kommer det att bli en inhomogen ekvation i (1.23). Hur man löser en sådan behandlar vi i nästa avsnitt. Med metoden beskriven i detta avsnitt har man något att jämföra numeriska lösningar av problem med inhomogena Dirichletvillkor med. Detta ligger dock utanför vad denna rapport täcker.

1.3 Det inhomogena fallet

1.3.1 Denitioner och satser

Ett kraftfullt verktyg för att lösa partiella dierentialekvationer är integral-transformer, exempelvis Laplacetransformen och Fouriertransformen. Samt-liga denitioner och satser i detta stycke är tagna ur [4], om inte annat är angivet.

Denition 1.3.1. [9] Laplacetransformen för en funktion f(t), denierad för alla t ≥ 0 ges som

L(f (t)) = ˜f (s) = Z ∞

0

e−stf (t)dt, s ∈ C.

Denition 1.3.2. [9] Den inversa Laplacetransformen för en funktion ˜f (s), ges som L−1( ˜f (s)) = f (t) = 1 2πiT →∞lim Z γ+iT γ−iT estf (s)ds˜

där γ är realdelen av s och väljs så att den är större än den reella delen av alla singulariteter till F (s). Detta försäkrar oss om att vi integrerar inom konvergensområdet.

Sats 1.3.1. [9] Låt ˜f (s) vara Laplacetransformen av en funktion f(t) och antag att e−st_{f (t) → 0 då t → ∞. Då gäller att}

L(f0(t)) = sL(f (t)) − f (0) = s ˜f (s) − f (0).

Denition 1.3.3. Fouriertransformen för en funktion f(x) ges som F (f (x)) = ˆf (k) = √1

2π Z ∞

−∞

e−ikxf (x)dx.

Denition 1.3.4. Den inversa Fouriertransformen för en funktion ˆf (k)ges som F−1( ˆf (k)) = f (x) = √1 2π Z ∞ −∞ eikxf (k)dk.ˆ

Sats 1.3.2. Om f är en kontinuerlig, styckvis deriverbar funktion, f, f0 _∈

L1(R), och lim|x|→∞f (x) = 0, då gäller att

Ff0 _{= ikF {f } .}

Bevis. Partiell integration ger 1 √ 2π Z ∞ −∞ f0(x)e−ikxdx = √1 2π h f (x)e−ikxi∞ −∞+ ik √ 2π Z ∞ −∞

Korollarium 1.3.1. Om f är en kontinuerlig funktion, n gånger styckvis deriverbar och f, f0_{, . . . , f}(n)_{∈ L}1_{(R), och}

lim|x|→∞f(m)(x) = 0 för m = 0, . . . , n − 1

då gäller att

Fnf(n) o

= (ik)nF {f } .

Denition 1.3.5. En faltning av två funktioner f(x), g(x) ges som (f ∗ g)(x) = √1 2π Z ∞ −∞ f (x − ξ)g(ξ)dξ = √1 2π Z ∞ −∞ f (ξ)g(x − ξ)dξ. Sats 1.3.3. För Fouriertransformen av en faltning mellan två funktioner f, g gäller

F (f ∗ g) = F (f )F (g) ⇔ f ∗ g = F−1(F (f )F (g)), eller, om man så vill

Z ∞ −∞ g(x − ξ)f (ξ)dξ = Z ∞ −∞ eikxg(k) ˆˆ f (k)dk. (1.25) Bevis. F (f ∗ g) = 1 2π Z ∞ −∞ e−ikx Z ∞ −∞ f (x − ξ)g(ξ)dξdx = 1 2π Z ∞ −∞ g(ξ) Z ∞ −∞ e−ikxf (x − ξ)dxdξ = 1 2π Z ∞ −∞ g(ξ) Z ∞ −∞ e−ik(x+ξ)f (x)dxdξ = √1 2π Z ∞ −∞ g(ξ)e−ikξdξ√1 2π Z ∞ −∞ e−ikxf (x)dx = F (g)F (f ).

För att kunna tillämpa dessa transformer på de dierentialekvationer vi är intresserade av, behöver vi redogöra för distributioner och Greenfunktio-ner.

Denition 1.3.6. En testfunktion är en oändligt deriverbar funktion på RN

som försvinner utanför något begränsat område. Rummet av alla testfunk-tioner skriver vi som D(RN_{) eller bara D.}

Denition 1.3.7. Med en distribution F på RN_{, menas en kontinuerlig}

linjär funktional på D(RN_{). Med andra ord, en avbildning F : D(R}N_{) → C}

1. F (aϕ + bψ) = aF (ϕ) + bF (ψ) för varje a, b ∈ C och ϕ, ψ ∈ D(RN₎

2. F (ϕn) → F (ϕ)(i C) då ϕn D

→ ϕ.

Rummet av alla distributioner skriver vi som D0_(RN_{) eller bara D}0_.

Denition 1.3.8. En Greenfunktion G(x − ξ, t − τ) av en linjär dierentia-loperator A = A(x, t) i 2D som verkar på distributioner, vid en punkt (ξ, τ), är en lösning till

AG(x − ξ, t − τ ) = δ(x − ξ)δ(t − τ )

där δ är Diracs deltafunktion. Denna egenskap ska visa sig kunna användas till att lösa partiella dierentialekvationer.

Diracs deltafunktion (se [9]) är en distribution, med följande egenskaper, 1. δ(x) = ( +∞, x = 0 0, x 6= 0 2. δ(x) = δ(−x) 3. R R_Ωf (x, t)δ(x − ξ)δ(t − τ )dxdt = f (ξ, τ ), om (ξ, τ) ∈ Ω 4. F(δ(x − ξ)) = _√1 2πe −ikξ 5. L(δ(t − τ)) = e−sτ_.

1.3.2 Lösningar med Greenfunktion, inhomogena fallet I det inhomogena fallet vill vi då lösa följande system

∂T ∂t − D · ∂2T ∂x2 = v(x, t), 0 < x < L, t > 0 T (x, 0) = f (x) T (0, t) = T (L, t) = 0 (1.26)

För att lösa detta delar vi upp (1.26) i två separata problem, ett med en homogen värmeledningsekvation med inhomogent begynnelsevillkor

∂u1 ∂t − D · ∂2u1 ∂x2 = 0, 0 < x < L, t > 0 u1(x, 0) = f (x) u1(0, t) = u1(L, t) = 0, (1.27)

och ett med en inhomogen värmeledningsekvation med homogent begynnel-sevillkor ∂u2 ∂t − D · ∂2u2 ∂x2 = v(x, t), 0 < x < L, t > 0 u2(x, 0) = 0 u2(0, t) = u2(L, t) = 0. (1.28)

Sätter vi nu T (x, t) = u1(x, t) + u2(x, t), kan vi visa att villkoren i (1.26)

stämmer. Vi betecknar den linjära dierentialoperatorn ∂ ∂t− D

∂2

∂x2 = Aoch

låter den verka på u1(x, t) + u2(x, t). Då får vi

AT = A(u1+ u2) = Au1+ Au2= 0 + v(x, t) = v(x, t)

T (0, t) = u1(0, t) + u2(0, t) = 0 = u1(L, t) + u2(L, t) = T (L, t)

T (x, 0) = u1(x, 0) + u2(x, 0) = f (x) + 0 = f (x).

(1.29)

Så vi kan vara säkra på att T (x, t) = u1(x, t) + u2(x, t). För u1(x, t) hämtar

vi nu en alternativ metod ur [9] och börjar med att applicera Fouriertrans-formen på problemet m.a.p. x, vilket ger oss en ordinär dierentialekvation

∂ ˆu1(k, t) ∂t = −Dk 2_uˆ 1(k, t) ˆ u1(k, 0) = ˆf (k), (1.30)

vars lösning ges som

ˆ

u1(k, t) = ˆf (k)e−Dk

2_t

. (1.31)

Sedan får vi applicera den inversa Fouriertransformen m.a.p. k för att få fram u1(x, t) = 1 √ 2π Z ∞ −∞ eikxf (k)eˆ −Dk2tdk. (1.32)

Enligt sats 1.3.3 så gäller u1(x, t) = Z ∞ −∞ g(x − ξ, t)f (ξ)dξ, (1.33) där g(x, t) = √1 2πF −1_(e−Dk2_t ) = 1 2π Z ∞ −∞ eikx−Dk2tdk = 1 2π Z ∞ −∞ exp " −Dt  k − ix 2Dt 2 − x 2 4Dt # dk = 1 2π  π Dt 1/2 exp  − x 2 4Dt  = √ 1 4πDtexp  − x 2 4Dt  . (1.34)

För att få fram u2(x, t), så betraktar vi först

AG(x − ξ, t − τ ) = δ(x − ξ)δ(t − τ ). (1.35) Vi multiplicerar med v(ξ, τ) och integrerar över vårt område m.a.p. ξ och τ.

Z L 0 Z t 0 AG(x − ξ, t − τ )v(ξ, τ )dξdτ = Z L 0 Z t 0 δ(x − ξ)δ(t − τ )v(ξ, τ )dξdτ. (1.36) Nu utnyttjar vi att A inte verkar m.a.p ξ och τ, samt egenskaper för delta-funktionen Z L 0 Z t 0 AG(x − ξ, t − τ )v(ξ, τ )dξdτ = A Z L 0 Z t 0 G(x − ξ, t − τ )v(ξ, τ )dξdτ = = Z L 0 Z t 0 δ(ξ − x)δ(τ − t)v(ξ, τ )dξdτ = v(x, t). (1.37) Och eftersom vi har att

Au2(x, t) = v(x, t) (1.38) får vi att Au2(x, t) = A Z L 0 Z t 0 G(x − ξ, t − τ )v(ξ, τ )dξdτ ⇔ u2(x, t) = Z L 0 Z t 0 G(x − ξ, t − τ )v(ξ, τ )dξdτ. (1.39)

Så att vi kan skriva vår lösning T (x, t) som

T (x, t) = Z L 0 g(x − ξ, t)f (ξ)dξ + Z L 0 Z t 0 G(x − ξ, t − τ )v(ξ, τ )dξdτ. (1.40)

Nu återstår att få fram Greenfunktionen G(x − ξ, t − τ). Om vi låter Gt

vara derivatan m.a.p. t och Gxx vara andraderivatan m.a.p. x, så vet vi från

denition 1.3.8 och (1.28) att G(x − ξ, t − τ) måste uppfylla Gt− DGxx= δ(x − ξ)δ(t − τ )

G(x, 0) = 0 (⇒ ˆG(k, 0) = 0). (1.41) Vi börjar med att applicera Fouriertransformen m.a.p. x,

F (G_t) − DF (Gxx) = δ(t − τ )F (δ(x − ξ)) F (Gt) + Dk2F (G) = δ(t − τ )e−ikξ ∂ ∂tF (G) + Dk 2_{F (G) = δ(t − τ )e}−ikξ ˆ G(k, t) + Dk2Gˆt(k, t) = δ(t − τ )e−ikξ. (1.42)

Sedan Laplacetransformen m.a.p. t,

L( ˆGt(k, t)) + Dk2L( ˆG(k, t)) = L(δ(t − τ ))e−ikξ sL( ˆG(k, t)) − ˆG(k, 0) + Dk2L( ˆG(k, t)) = e−sτe−ikξ sG(k, s) + Dkˆ˜ 2G(k, s) = eˆ˜ −(sτ +ikξ) ˆ ˜ G(k, s) = 1 s + Dk2e −(sτ +ikξ)_. (1.43)

Genom att ta inversen för Fouriertransformen och Laplacetransformen kom-mer man slutligen fram till

G(x − ξ, t − τ ) = 1 p4πD(t − τ )exp  − (x − ξ) 2 4D(t − τ )  . (1.44)

Så vi ser att g(x, t) = G(x, t) och lösningen kan till slut skrivas som

T (x, t) = Z L 0 G(x − ξ, t)f (ξ)dξ + Z L 0 Z t 0 G(x − ξ, t − τ )v(ξ, τ )dξdτ. (1.45)

Kapitel 2

Slumpvandring

2.1 Metod

2.1.1 Begreppet slumpvandring

För att förklara begreppet slumpvandring formulerar vi ett enkelt exempel: Exempel 2.1.1. Vi tänker oss att vi spelar ett spel, där vi börjar med W0kr.

Vi singlar en slant varpå vi vinner en kr om slanten visar krona och förlorar en kr om slanten visar klave. Vi antar rimligen att sannolikheten är 1/2 för vardera valör på myntet. Låt si vara +1 eller −1 vid försök i, beroende på

om vi får krona eller klave. Efter att vi har singlat slant i gånger kan vi beräkna antal kr vi har enligt

Wi= W0+ s1+ s2+ . . . + si = W0+ i

X

k=1

sk. (2.1)

Detta beskriver en slumpvandring med steglängd 1 och startpunkt W0.

Vi kommer hädanefter att kalla den punkt som vandrar för partikel.

2.1.2 Koppling till värmeledningsekvationen

För att ge en antydan om en förklaring till hur slumpvandringar kan an-vändas vid partiella dierentialekvationer (vilket behandlas utförligare i [6]), betraktar vi ett speciellt exempel. Låt P (x, t) beteckna sannolikheten att vi benner oss i punkt x efter t steg, med steglängder ∆x och ∆t. Eftersom vi räknar med att det är lika stor sannolikhet att vi rör oss till vänster som till höger, får vi att

P (x, t) = 1

2P (x + ∆x, t − ∆t) + 1

2P (x − ∆x, t − ∆t). (2.2) Genom att subtrahera P (x, t − ∆t) och dela med ∆t får vi

P (x, t) − P (x, t − ∆t)

∆t =

P (x + ∆x, t − ∆t) + P (x − ∆x, t − ∆t) − 2P (x, t − ∆t)

2∆t .

(2.3) Vi noterar att vänsterledet ser ut som derivatans denition i t om vi låter ∆t → 0. För att se vad som händer i högerledet, Taylorutvecklar vi P (x, t). Vi bara tar med de termer som vi behöver från utvecklingen nedan,

P (x, t − ∆t) ≈ P (x, t) −∂P ∂t · ∆t P (x ± ∆x, t − ∆t) ≈ P (x, t) ±∂P ∂x · ∆x + 1 2 ∂2P ∂x2 · (∆x) 2₋∂P ∂t · ∆t. (2.4)

Med detta kan vi skriva ekvation (2.3) som P (x, t) − P (x, t) +∂P_∂t · ∆t ∆t ≈ 1 2∆t·  P (x, t) + ∂P ∂x · ∆x + 1 2 ∂2P ∂x2 · (∆x) 2₋∂P ∂t · ∆t +P (x, t) −∂P ∂x · ∆x + 1 2 ∂2P ∂x2 · (∆x) 2₋∂P ∂t · ∆t −2P (x, t) + 2∂P ∂t · ∆t  . (2.5) Detta förenklas nu till

∂P ∂t ≈ (∆x)2 2∆t ∂2P ∂x2. (2.6)

Vid gränsen ∆x → 0, ∆t → 0 men där kvoten (∆x)2_{/∆t är begränsad, blir}

detta en exakt relation, detta ser nu precis ut som värmeledningsekvationen. Om vi sätter (∆x)2_{/2∆t = D, där D är diusionskonstanten, får vi}

∂P ∂t = D

∂2P

∂x2. (2.7)

Vi minns att ∆x och ∆t är våra valda steglängder i x och t. Denna relation ger oss ett samband som behöver uppfyllas när vi väljer våra steglängder,

för att ekvationen (2.7) ska vara uppfylld. Notera att vi nu tagit fram detta samband gällande sannolikhet för en partikel. När vi sedan använder detta samband i praktiken så har vi era partiklar som rör sig, varpå det snarare är en täthet vi kommer beräkna. Så vi låter framöver i rapporten P (x, t) vara tätheten av partiklar snarare än en sannolikhet.

2.2 Numerisk lösning, homogena randvillkor

Nu när vi har sett att vi kan lösa värmeledningsekvationen genom att lösa ett ekvivalent problem gällande sannolikhet, ska vi titta på hur vi lösa det homogena fallet med hjälp av slumpvandringar. Vi påminner om att vi vill beräkna tätheten av partiklar i punkten x vid tiden t. För att detta ska lösa värmeledningsekvationen, vill vi ha samma begynnelse-och randvillkor som i ekvation (1.6).

P (x, 0) = 1, 0 < x < L P (0, t) = 0, t > 0 P (L, t) = 0, t > 0

(2.8)

Villkoret P (x, 0) = 1 tolkas som att vi har en partikel i varje punkt initialt. Villkoren P (0, t) = P (L, t) = 0 tolkas som att en partikel försvinner vid kontakt med randen, vi säger att vi har en absorberande rand. Nedan ses ett exempel med x0= 0.4, som försvinner när den når randen vid x = 1.

Figur 2.2: Exempel på slumpvandring. Låt oss nu diskretisera vårt x-t-område enligt

∆x = xi+1− xi= 10−2 i = 0, 1, . . . , N − 1 x0 = 0, xN = 1

∆t = tj+1− tj = 5 · 10−5 j = 0, 1, . . . , M − 1 t0 = 0, tM = 1

(2.9) Här har vi för enkelhetens skull satt L = D = 1. Vi får följande villkor från (2.6) och (2.7), vilket ger har gett oss ovanstående samband mellan ∆x och ∆t. (∆x)2 2∆t = D ⇐⇒ ∆x = √ 2D∆t ⇐⇒ ∆t = (∆x) 2 2D . (2.10)

Så, för att beräkna P (xi, tj), initierar vi en slumpvandring i P (xi, t0), för

varje i och håller koll på hur varje partikel vandrar. För att få hög nog-grannhet gör vi detta många gånger och tar sedan ett medelvärde genom att dela med totala antalet vandringar vi har utfört. Detta kan beskrivas med en pseudokod.

Sätt m = antal inre punkter i x, n = antal slumpvandringar för varje partikel for l=1:m do Sätt startpunkt = xl, t = ∆t for k=1:n do while 0 < xi< 1 do Ta steg i = i ± 1 Lägg till 1 till P (i, j) Ta steg i t (tj = tj+1) end end Uppdatera startpunkt (xl= xl+1) end Sätt P = P/n

2.3 Resultat

Här visar vi ett antal gurer som jämför den analytiska lösningen till test-problemet (1.6) och den numeriska lösningen.

Figur 2.3: Jämförelse mellan analytisk lösning (med 100 termer) och nume-risk lösning till ekvation (1.6), för xerat x = 0.5, med 104 _{iterationer och}

∆x = 10−2.

Figur 2.4: Absoluta felet mellan analytisk lösning (med 100 termer) och numerisk lösning till ekvation (1.6), vid x = 0.5, med 104 _{iterationer och}

Figur 2.5: Jämförelse mellan analytisk lösning (med 100 termer) och nume-risk lösning till ekvation (1.6), med 104 _{iterationer och ∆x = 10}−2_.

Figur 2.6: Jämförelse mellan analytisk lösning (med 100 termer) och nume-risk lösning till ekvation (1.6), med 104 _{iterationer och ∆x = 10}−2_{, x-t-plan.}

Kapitel 3

Utvidgning till rörlig rand

3.1 Formulering

Vid tillämpningar av värmeledningsekvationen nns det många fall där ran-den inte är xerad, utan istället rör sig med tiran-den. Exempel på detta är när is bildas eller smälter. Vid isbildning till exempel, kommer isen bli tjockare allteftersom tiden går. Därför är det ofta av intresse att lösa värmelednings-ekvationen i dessa fall, då randen rör på sig med tiden.

Låt s(t) vara en funktion (dessa är problemberoende) som betecknar den övre gränsen i x. Då ges vårt problem som

T (x, 0) = 1, 0 < x < s(t), t > 0, s(0) = 1 T (0, t) = 0, t > 0

T (s(t), t) = 0, t > 0

(3.1)

och vårt ekvivalenta problem som

P (x, 0) = 1, 0 < x < s(t), t > 0, s(0) = 1 P (0, t) = 0, t > 0

P (s(t), t) = 0, t > 0

(3.2)

För att skapa ett enkelt testproblem börjar vi med något enkelt s(t) för att utvidga vår metod till rörlig rand, så vi sätter s(t) = 1 + t.

3.2 Slumpvandring

3.2.1 Metod för rörlig rand

Vi bygger vidare på vårat problem och minns från kapitel 2.2 att vi räknar med en absorberande rand, dvs att när vi träar randen med en partikel så försvinner den. Nu måste vi även ta i beaktning att randen föryttar sig med tiden. Detta gör att när vi tar ett steg i vår slumpvandring, måste vi ta reda på vart randen ligger för detta steg så att vi kan kontrollera om vi har nått randen eller ej. Således bör vi kunna använda vår metod för xerad rand, men vi modierar den så att den kontinuerligt uppdaterar positionen för randen i algoritmen.

Sätt m = antal inre punkter i x, n = antal iterationer for l=1:m do

Sätt startpunkt = xl, t = ∆t

for k=1:n do

while 0 < xi< s(tj) do

Ta steg i = i ± 1 Lägg till 1 till P (i, j) Ta steg i t (tj = tj+1) end end Uppdatera startpunkt (xl= xl+1) end Sätt P = P/n

3.2.2 Stefanproblemet, något att jämföra med

Ett av de mest välkända problemen med rörlig rand är Stefanproblemet, som handlar om att beräkna temperatur vid smältning av is. Detta avsnitt är hämtat från [7]. I detta fallet kommer vi endast betrakta temperaturen i vattnet som bildas vid smältningen, men det går även att beräkna hur temperaturen uppför sig i isen. Stefanproblemet är formulerat som

∂T ∂t = Dl· ∂2T ∂x2, 0 < x < s(t) T (0, t) = f (t), t > 0 T (x, 0) = 0 T (s(t), t) = 0 = Ts s(0) = 0 c1· ds dt = c2· ∂T ∂x (3.3)

där Dl är diusionskonstanten för vatten, c1, c2 är konstanter som bestäms

av vattnets fysikaliska egenskaper, s(t) är en funktion som beskriver hur randen (fasövergången mellan vatten och is) rör sig med tiden. Det sista villkoret kallas för Stefanvillkor (eng. Stefan condition) och beskriver ett samband i fasövergången mellan randen och temperaturen, ett samband som kan härledas ur den s.k. energiprincipen.

Figur 3.1: Illustration, randen (fasövergången) som rör sig med tiden Vi betraktar fallet då vi håller en konstant temperatur T0 vid x = 0, så

vi sätter f(t) = T0. Enligt [7] ges lösningen till detta system ges som

T (x, t) = T0− T0 erf(λ)·erf  x 2√Dlt  s(t) = 2λpDlt (3.4)

där λ är lösningen till ekvationen

λeλ2erf(λ) = ck∆T Lm

√

π. (3.5)

Här är ck värmekapaciteten för vatten, ∆T = T0− Ts, Lm latent smältvärme

för vatten och erf(x), kallat error function är en funktion denierad som erf(x) = _√2 π Rx 0 et 2 dt.

I den senare jämförelsen av den exakta och den numeriska lösningen har jag använt följande fysikaliska konstanter, angivna i SI-enheter.

Dl= 1.4868 · 10−7 [m2/s]

ck= 4.0897 · 103 [J/(kg · K)]

Lm= 334 · 103 [J/kg]

(3.6)

3.2.3 Numerisk beräkning av Stefanproblemet

När vi ska lösa Stefanproblemet med slumpvandringsmetoden måste vi mo-diera metoden vi använde oss av vid testproblemet. Minns att för det enkla testproblemet (se 3.1) hade vi begynnelsevillkoret T (x, 0) = 1 och rand-villkoret T (0, t) = 0. I termer av sannolikhet innebär begynnelserand-villkoret

P (x, 0) = 1att det benner sig en partikel i varje punkt x initialt. Nu kom-pliceras fallet genom att vi då har vårat diskretiserade problem med villkor givet som ∂P (xi, tj) ∂t = Dl· ∂2P (xi, tj) ∂(x)2 , 0 < xi< s(tj) P (x0, tj) = 20 = P0, t > 0 P (xi, t0) = 0 P (s(tj), tj) = 0 = Ps s(tj) = 2λpDltj. (3.7)

Här har vi att det benner sig 20 partiklar vid undre randen P (x0, tj) för

varje tj. Om vi nu räknar med inhomogen rand, får vi nu istället initiera

slumpvandringar i (x0, tj)för alla j istället. För att randvillkoret P (x0, tj) =

P0 ska stämma, får vi initiera P0· n slumpvandringar i (x0, tj)för alla j, där

när antalet iterationer. Nedan ses en pseudokod för algoritmen. Sätt m = antal inre punkter i t, n = antal iterationer for l=2:m do for k=1:(P0· n) do Sätt startpunkt x0 = 0, tj = tl while xi < s(tj) do Ta steg i = i ± 1 if 0 < xi < s(tj) then

Lägg till 1 till P (i, j) end if xi = 0 then Avbryt end Ta steg i t (tj = tj+1) end end Uppdatera startpunkt (tl = tl+1) end Sätt P = P/n

3.3 Felanalys för Stefanproblemet

Här återkopplar vi till (2.4) där vi Taylorutvecklade funktionen P (x, t), men då ignorerade vi för tillfället ordo-termerna. Vi ska nu titta mer formellt på utvecklingen. Enligt [8] kan vi skriva Taylorutvecklingen för en funktion f av två variabler som f (x + ∆x, t + ∆t) = f (x, t) + [fx(x, t)∆x + ft(x, t)∆t] + 1 2!fxx(x, t)(∆x) 2_{+ 2f} xt(x, t)∆x∆t + ftt(x, t)(∆t)2 + . . . (3.8)

Så låt oss nu, med anledning av (2.4), utveckla P (x, t − ∆t). Vi får då P (x, t−∆t) = P (x, t)−Pt(x, t)∆t+

1

2Ptt(x, t)(∆t)

2_{+. . . = P (x, t)−P}

t(x, t)∆t+O((∆t)2).

Och sedan utvecklar vi P (x ± ∆x, t − ∆t), P (x ± ∆x, t − ∆t) = P (x, t) ± Px(x, t)∆x + 1 2Pxx(x, t)(∆x) 2_{− P} t(x, t)∆t± 2Pxt(x, t)∆x∆t + 1 2Ptt(x, t)(∆t) 2_{± . . . =} P (x, t) ± Px(x, t)∆x + 1 2Pxx(x, t)(∆x) 2_{− P} t(x, t)∆t + O((∆t)2),

vilket tyder på att felet bör vara av ordning (∆t)2_{. Att den sista likheten}

ger O((∆t)2_{) är p.g.a. att ∆x < ∆t genom sambandet (∆x)}2_{/2∆t = D när}

vi väljer D som diusionskonstanten för vatten. För att analysera felet av den numeriska lösningen till Stefanproblemet, konstruerar vi först en matris som vi kan kalla för Texakt, med den exakta lösningen (se (3.4)) i de diskreta

punkterna. Sedan tar vi fram en matris E med felet i varje punkt enligt E(i, j) = |Texakt(i, j) − P (i, j)| . (3.9)

Denna vill vi nu beräkna normen för, för att få ett mått på felet. Man kan visa (se [11]) att den exakta lösningen till Stefanproblemet (3.4) ligger i L2(R),

därför använder vi `2-normen, denierad som

kEk_`2 =   X i X j |E(i, j)|2∆x∆t   1/2 . (3.10)

Nedan har jag samlat normen av felet i en tabell för olika val av ∆x och olika antal iterationer för slumpvandringarna för Stefanproblemet.

∆x ∆t Iterationer kEk`2 0.0005 0.8407 103 _0.8119 0.0005 0.8407 104 _0.8470 0.0005 0.8407 105 _0.8176 0.0003 0.3027 103 _0.4588 0.0003 0.3027 104 _0.5055 0.0003 0.3027 105 _0.4768 0.0001 0.0336 103 _0.1959 0.0001 0.0336 104 _0.1669 0.0001 0.0336 105 _0.1583

Från tabellen kan vi se att diskretiseringen av området har större påverkan i att minska felet, än antalet iterationer för slumpvandringarna på den nivån. Dock noterade jag att för Stefanproblemet så var felet stort i början av området, troligtvis för att området är väldigt litet i början vilket ger oss dålig statistisk mätdata där. Slutligen gjorde jag en gur med normen av felet som funktion av √∆t, vilket bör ge linjärt samband om felet följer (∆t)2.

Figur 3.2: Norm av felet mellan exakt lösning och numerisk lösning för Ste-fanproblemet, som funktion av √∆t.

3.4 Resultat, testproblem

Här visar vi gurer för numerisk lösning till testproblemet (3.2).

Figur 3.3: Numerisk lösning för ekvation (3.1), med 104_{iterationer och ∆x =}

10−2.

Figur 3.4: Numerisk lösning för ekvation (3.1), med 104_{iterationer och ∆x =}

3.5 Resultat, Stefanproblemet

Här visar vi ett antal gurer som jämför den analytiska lösningen till Ste-fanproblemet (3.3) och den numeriska lösningen.

Figur 3.5: Jämförelse, exakt lösning (3.4) och numerisk lösning med T0 =

20Co för t ∈ [0, 10] minuter.

Figur 3.6: Jämförelse, exakt lösning (3.4) och numerisk lösning med T0 =

Figur 3.7: Numerisk lösning till Stefanproblemet (3.7) med T0 = 20Co, 104

iterationer och ∆x = 10−4_{m för xed t ≈ 300s.}

Figur 3.8: Relativa felet mellan den numeriska och analytiska lösningen till Stefanproblemet med T0 = 20Co, 104 iterationer och ∆x = 10−4m för t ≈

Kapitel 4

Appendix

Här samlar vi matlab-kod som har använts till de numeriska lösningarna.

1 %% Numerisk beräkning av testproblemet med fixerad rand

2 syms x t

3 D=1; % Diffusionskonstant

4 L=1; % Ändpunkt i intervall

5 dx = 0.01; % Steglängd i x

6 dt = ((dx)^2) / (2*D); % Steglängd i t

7 nof_traj_t = 1e4; % Antal iterationer för varje partikel

8 nof_traj_x = (L/dx)-1; % Övre gräns för startpunkter

9 t_vector = 0:dt:L;

10 x_vector = 0:dx:L;

11 P = zeros(length(x_vector),length(t_vector)); % Preallokera ... minne till matris med slumpvandringar

12 x=dx; % Sätt startpunkt för slumpvandring

13 for l=2:nof_traj_x % Loop för startpunkter

14 for m=1:nof_traj_t % Loop för antal iterationer

15 x_traj=x; % Starta slumpvandring i x

16 i=(x_traj/dx)+1; % Sätt plats för x i matris P

17 j=2; % Sätt plats för t i matris P

18 while x_traj>0 && x_traj<L && j≤length(P) % Starta ... splumpvandring, avbryt vid kontakt med rand

19 pm = (-1)^round(rand); % pm är plus eller minus 1

20 x_traj = x_traj + pm*dx; % Ta steg i x för partikel

21 i=round(i+pm); % Uppdatera punkt för x i P

22 if i>1 && i<length(x_vector) % Registrera om partikel ... ej är på rand 23 P(i,j) = P(i,j) + 1; 24 end 25 j = j + 1; % Ta steg i t 26 end 27 end 28 x=x+dx; % Uppdatera startpunkt i x 29 end

30 P=P/nof_traj_t; % Beräkna täthet

1 %% Numerisk beräkning av testproblemet med rörlig rand 2 syms x t 3 D=1; % Diffusionskonstant 4 L=1; % Ändpunkt i intervall 5 dx = 0.01; % Steglängd i x 6 dt = ((dx)^2) / (2*D); % Steglängd i t

7 f = L+t; % Definiera övre randen som funktion av t

8 f_max = double(subs(f,L)); % Beräkna max av övre rand

9 rows = f_max/dx + 1; % Antal rader som krävs för matris

10 row_vector = 0:dx:f_max;

11 nof_traj_t = 1e4; % Antal iterationer för varje partikel

12 nof_traj_x = (L/dx)-1; % Övre gräns för startpunkter

13 t_vector = 0:dt:L;

14 x_vector = 0:dx:L;

15 P = zeros(rows,length(t_vector)); % Preallokera minne till ... matris med slumpvandringar

16 x=dx; % Sätt startpunkt för slumpvandring

17 t=0;

18 F=zeros(1,length(t_vector));

19

20 for i=1:length(t_vector) % Vektor med funktionsvärdet av ... f för varje diskret t

21 F(i)=double(subs(f,t));

22 t=t+dt;

23 end

24

25 for l=1:nof_traj_x % Loop för startpunkter

26 for m=1:nof_traj_t % Loop för antal iterationer

27 x_traj=x; % Starta slumpvandring i x

28 i=(x_traj/dx)+1; % Sätt plats för x i matris P

29 j=2; % Sätt plats för t i matris P

30 t=dt; % Ta steg i t

31 while x_traj>0 && x_traj<F(j-1) && j≤length(P) % Starta ... splumpvandring, avbryt vid kontakt med rand

32 pm = (-1)^round(rand); % pm är plus eller minus 1

33 x_traj = x_traj + pm*dx; % Ta steg i x för partikel

34 i=round(i+pm); % Uppdatera punkt för x i P

35 if i>1 && i<F(j)/dx % Registrera om partikel ej är på rand

36 P(i,j) = P(i,j) + 1; 37 end 38 j = j + 1; 39 end 40 end 41 x=x+dx; % Uppdatera startpunkt i x 42 end

43 P=P/nof_traj_t; % Beräkna täthet

1 %% Numerisk beräkning av Stefanproblemet

2 run trans_eq.m % Kör skript som numeriskt löser den ... transcendentala ekvationen

3 syms x t

4 D=1.4868E-7; % Diffusionskonstant för vatten

5 dx = 0.0001; % Steglängd i x

6 dt = ((dx)^2) / (2*D); % Steglängd i t

7 s = 2*lambda*sqrt(D*t); % Den rörliga randen

8 T0=20; % Temperatur vid undre rand T(0,t), Celsius

9 time=600; % Övre gräns för t, sekunder

10 s_max = double(subs(s,time)); % Beräkna max av övre rand

11 t_vector = 0:dt:time;

12 x_vector = 0:dx:s_max;

13 rows = length(x_vector); % Antal rader som krävs för matris

14 iterations = 1e4; % Antal iterationer för varje partikel

15 it=T0*iterations; % Antal iterationer för algoritm

16 P = zeros(rows,length(t_vector)); % Preallokera minne till ... matris med slumpvandringar

17 S=zeros(1,length(t_vector));

18 t=0;

19 for i=1:length(t_vector) % Vektor med funktionsvärdet av s ... för varje diskret t

20 S(i) = double(subs(s,t));

21 t=t+dt;

22 end

23 t=dt; % Sätt startpunkt i t

24 for k=2:length(t_vector) % Loop för startpunkter

25 for m=1:iterations % Loop för antal slumpvandringar

26 j=round(t/dt); % Sätt plats för t i matris P

27 x_traj=0; % Sätt startpunkt i x

28 i=1; % Sätt plats för x i matris P

29 while x_traj<S(j) && j<length(P) % Starta ...

splumpvandring, avbryt vid kontakt med övre rand

30 pm = (-1)^round(rand); % pm är plus eller minus 1.

31 x_traj = x_traj + pm*dx; % Ta steg i x för ... partikel

32 i=round(i+pm); % Uppdatera punkt för x i P

33 if i>1 && i<S(j)/dx % Registrera om ... partikel ej är på rand

34 P(i,j) = P(i,j) + 1;

35 end

36 if i==1 % Avbryt vid kontakt med undre rand

37 break 38 end 39 j=j+1; % Uppdatera punkt för t i P 40 end 41 end 42 t=t+dt; % Uppdatera startpunkt i t 43 end

44 P=P/iterations; % Beräkna täthet

1 %% Löser den transcendentala ekvationen med ... Newton-Raphson's metod

2 l=334*10^3; % Latent värme [J/kg]

3 C_l=4.0897E+3; _{% Värmekapacitet vid T0 [J/kg*K]}

4 beta=l/C_l;

5

6 syms lambda

7 f=sqrt(pi)*beta*lambda*exp((lambda)^2)*erf(lambda) - 20;

8 % Vi vill hitta den positiva roten till ... f med T0=20

9

10 tol=1e-04; % Tolerans

11 lambda0=1; % Initial gissning

12 g=diff(f,lambda); % Derivatan v f

13 f0=double(subs(f,lambda,lambda0)); % Initial evaluering ... av f i lambda0

14 g0=double(subs(g,lambda,lambda0)); % Initial evaluering ... av g i lambda0

15

16 while f0 > tol

17

18 lambda = lambda0 -(f0/g0); % Ta steg

19 lambda0 = lambda; % Uppdatera lambda

20 f0=double(subs(f)); % Uppdatera f i nytt lambda

21 g0=double(subs(g)); % Uppdatera g i nytt lambda

22 23 end

Litteraturförteckning

[1] Gui-Qiang Chen, Henrik Shahgholian, Juan-Luis Vazquez, Free boundary problems: the forefront of current and future developments, 2015.

[2] Huyen Pham, Feynman-Kac representation of fully nonlinear PDEs and applications, 2014, arXiv:1409.0625v1 [math.PR].

[3] Steven G. Krantz, Dierential Equations - theory, technique and practice, 2nd edition, CRC Press, 2015.

[4] Lokenath Debnath & Piotr Mikusinski, Introduction to Hilbert Spa-ces, 3rd edition, Elsevier Academic Press, 2005.

[5] Timothy Sauer, Numerical Analysis, 2nd edition, Pearson, 2014. [6] Mark Kac, Random Walk and the Theory of Brownian Motion, The

American Mathematical Monthly: Vol.54, No 7, P369:391, 1947. [7] A. C. Boucígueza, R. F. Lozanoa, M. A. Lara, About the exact solution

in two phase-Stefan problem, Thermal Engineering, Vol. 6, December 2007, p. 70-75.

[8] Persson & Böiers, Analys i era variabler, 3rd edition, Studentlitte-ratur AB, 2005.

[9] Lokenath Debnath & Dambaru Bhatta, Integral Transforms And Their Applications, 3rd edition, CRC Press, 2015.

[10] Annika Sparr & Gunnar Sparr, Kontinuerliga System, 2nd edition, Studentlitteratur AB, 2000.

[11] Antonio Fasano & Mario Primicerio, General Free-Boundary Pro-blems for the Heat Equation, III, Journal of mathematical analysis and applications 59, 1977.