Libor Kadrnka Schémata typu Moving Mesh pro řešení

Univerzita Karlova v Praze Matematicko-fyzikáln´ı fakulta

´ RSK ˇ ´ PRACE ´ BAKALA A

Libor Kadrnka Sch´ emata typu Moving Mesh pro ˇ reˇ sen´ı nestacion´ arn´ıch u ´ loh Katedra numerické matematiky

Vedouc´ı bakaláˇrské práce: Doc. RNDr. Jiˇr´ı Felcman, CSc. Studijn´ı program: Matematika, obecná matematika

2007

Chtˇel bych podˇekovat vedouc´ımu své bakaláˇrské práce panu doc. RNDr. Jiˇr´ımu Felcmanovi, CSc., za cenné rady a pˇripom´ınky bˇehem veden´ı mé bakaláˇrské práce, panu Mgr. Petru Kuberovi za pomoc pˇri programován´ı, pan´ı RNDr. Ivetˇe Hnˇetynkové, Ph.D., za zap˚ ujˇcen´ı poˇc´ıtaˇce a zejména sv´ ym rodiˇc˚ um za podporu bˇehem celého mého studia.

Prohlaˇsuji, ˇze jsem svou bakaláˇrskou práci napsal samostatnˇe a v´ yhradnˇe s pouˇzit´ım citovan´ ych pramen˚ u. Souhlas´ım se zap˚ ujˇcován´ım práce a jej´ım zveˇrejˇ nován´ım. V Praze dne 1. srpna 2007

Libor Kadrnka

2

Obsah ´ 1 Uvod

5

2 Popis metody 2.1 Vytvoˇren´ı s´ıtˇe pomoc´ı minimalizace funkcionálu . . . . . . . . . . . . . 2.1.1 Jednorozmˇern´ y pˇr´ıpad . . . 2.2 Popis algoritmu . . . . . . . . . . . 2.2.1 Metoda s´ıt´ı . . . . . . . . . 2.2.2 Metoda umˇelého ˇcasu . . . .

7 . . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

7 8 11 11 12

3 Numerick´ e experimenty 3.1 Funkce s pohybuj´ıc´ı se nespojitost´ı . . 3.2 Aplikace adaptivn´ı s´ıtˇe pˇri ˇreˇsen´ı PDR 3.2.1 Burgersova rovnice . . . . . . . 3.2.2 Riemann˚ uv problém . . . . . . 3.3 Srovnán´ı r˚ uzn´ ych adaptivn´ıch metod .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

18 18 25 25 26 30

Literatura

. . . . .

33

3

Název práce: Schémata typu Moving Mesh pro ˇreˇsen´ı nestacionárn´ıch u ´loh Autor: Libor Kadrnka Katedra (´ ustav): Katedra numerické matematiky Vedouc´ı bakaláˇrské práce: Doc. RNDr. Jiˇr´ı Felcman, CSc. e-mail vedouc´ıho: [email protected] Abstrakt: V pˇredloˇzené práci se zab´ yv´ ame vytvoˇren´ım adaptivn´ı v´ ypoˇcetn´ı s´ıtˇe, která se pouˇz´ıvá pˇri ˇreˇsen´ı parciáln´ıch diferenciáln´ıch rovnic (PDR). S´ıt’ budeme adaptovat v závislosti na chov´ an´ı pˇribliˇzného ˇreˇsen´ı této rovnice v ˇcase. Problém adaptace v´ ypoˇcetn´ı s´ıtˇe pˇrevedeme na minimalizaci jistého funkcionálu. Nutnou podm´ınkou k tomu, aby funkcion´ al nab´ yval minima, je splnˇen´ı jisté diferenciáln´ı rovnice. Hledané uzly adaptované v´ ypoˇcetn´ı s´ıtˇe jsou tvoˇreny hodnotami diskrétn´ıho ˇreˇsen´ı této rovnice. Na nespojit´ ych funkc´ıch simuluj´ıc´ıch v´ yvoj pˇribliˇzného ˇreˇsen´ı PDR v ˇcase ukáˇzeme, jak adaptace s´ıtˇe funguje. Poté s pomoc´ı adaptivn´ı v´ ypoˇcetn´ı s´ıtˇe najdeme pˇribliˇzné ˇreˇsen´ı dvou PDR a toto ˇreˇsen´ı srovn´ ame s ˇreˇsen´ım pˇresn´ ym a s ˇreˇsen´ım z´ıskan´ ym bez vyuˇzit´ı adaptivn´ı s´ıtˇe. Na závˇer srovn´ ame metodu popsanou v této práci s anizotropn´ı adaptivn´ı metodou vyv´ıjenou na katedˇre numerické matematiky MFF UK. Kl´ıˇcová slova: adaptivn´ı s´ıt’, geometrick´ y zákon zachov´ an´ı, metoda koneˇcn´ ych objem˚ u, pohybuj´ıc´ı se nespojitost

Title: Moving Mesh schemes for nonstationary problems Author: Libor Kadrnka Department: Department of Numerical Mathematics Supervisor: Doc. RNDr. Jiˇr´ı Felcman, CSc. Supervisor’s e-mail address: [email protected] Abstract: The presented work focuses on creating an adaptive computational mesh which is to be used for solving partial differential equations (PDE). The mesh is adapted according to the numerical solution behaviour of this equation in time. The problem of mesh adaptation is reduced to minimization of a certain functional. A necessary condition for the functional to have a minimum is that a certain differential equation holds. The grid points of the moving mesh on the next time level are given by the discrete solution of this equation. On discontinuous functions simulating the evolution of the PDE numerical solution in time we demonstrate how the mesh adaptation works. Then, using a moving mesh, we find a numerical solution to two PDEs and we compare this solution with the exact solution and with the solution obtained without using any adaptive mesh. Finally, we compare the adaptive method described in this work with an anisotropic mesh adaptation method developed at the Department of Numerical Mathematics MFF UK. Keywords: adaptive mesh, geometric conservation law, finite volume method, moving discontinuity

4

Kapitola 1 ´ Uvod Funkce, které jsou ˇreˇsen´ım parciáln´ıch diferenciáln´ıch rovnic, obecnˇe mohou obsahovat nespojitosti nebo mohou m´ıt v nˇekter´ ych bodech velkou hodnotu gradientu. Pˇri hledán´ı ˇreˇsen´ı numerickou metodou, jako je napˇr´ıklad metoda koneˇcn´ ych objem˚ u, m˚ uˇze pˇri pouˇzit´ı rovnomˇerné triangulace oblasti, na n´ıˇz rovnici ˇreˇs´ıme, doj´ıt ke zkreslen´ı tohoto ˇreˇsen´ı v m´ıstech, kde se tato nespojitost ˇci vysok´ y gradient vyskytuje. Proto se pouˇz´ıvaj´ı tzv. adaptivn´ı metody, kdy se v pr˚ ubˇehu v´ ypoˇctu triangulace oblasti v kaˇzdém kroku pˇrizp˚ usobuje pˇribliˇznému ˇreˇsen´ı z´ıskanému v kroku pˇredchoz´ım na základˇe daného kritéria adaptivn´ı metody. V naˇsem pˇr´ıpadˇe se v m´ıstech s vysok´ ym gradientem nebo v m´ıstech nespojitosti zv´ yˇs´ı poˇcet prvk˚ u triangulace, zat´ımco na zbytku oblasti se jejich poˇcet sn´ıˇz´ı, ˇza´dné nové prvky triangulace tedy nevytváˇr´ıme, pouze se na základˇe chován´ı z´ıskaného pˇribliˇzného ˇreˇsen´ı na v´ ypoˇcetn´ı oblasti snaˇz´ıme prvky triangulace rozloˇzit na v´ ypoˇcetn´ı oblasti v´ yˇse uveden´ ym zp˚ usobem. Pˇredpokládáme, ˇze pˇribliˇzné ˇreˇsen´ı parciáln´ı diferenciáln´ı rovnice se hledá metodou koneˇcných objem˚ u. Pˇri pouˇzit´ı této metody je pˇribliˇzné ˇreˇsen´ı z´ıskáno ve tvaru po ˇca´stech konstantn´ı funkce, pˇriˇcemˇz tato funkce je konstantn´ı ve vnitˇrku kaˇzdého prvku triangulace. V této práci se budeme zab´ yvat adaptivn´ı metodou popsanou v ˇclánku [3], kter´ y se t´ yká ˇreˇsen´ı jedno- a dvourozmˇern´ ych hyperbolick´ ych rovnic, my se vˇsak omez´ıme pouze na jednorozmˇern´ y pˇr´ıpad parciáln´ı diferenciáln´ı rovnice. Nejprve si tuto metodu podrobnˇe pˇredstav´ıme a následnˇe si pomoc´ı numerick´ ych experiment˚ u ukáˇzeme, jak funguje. Bˇehem experiment˚ u nejdˇr´ıve vyzkouˇs´ıme adaptivitu s´ıtˇe na funkc´ıch, u nichˇz známe jejich pr˚ ubˇeh, konkrétnˇe budeme u funkc´ı, které jsou po ˇca´stech lineárn´ı a jejichˇz nespojitost se v pr˚ ubˇehu ˇcasu pohybuje, zkoumat, zda se uzly s´ıtˇe taktéˇz pohybuj´ı spolu s nespojitost´ı a zda se souˇcasnˇe v daném ˇcasovém okamˇziku v m´ıstˇe nespojitosti nalézá v´ıce uzl˚ u neˇz v m´ıstech, kde takovou vlastnost funkce nemá a vˇse si ilustrujeme na obrázc´ıch. Poté pouˇzijeme adaptivn´ı s´ıt’ pˇri ˇreˇsen´ı parciáln´ıch diferenciáln´ıch rovnic, je-

5

jichˇz pˇresné ˇreˇsen´ı známe a ukáˇzeme, ˇze adaptivita s´ıtˇe v´ yraznˇe zvyˇsuje pˇresnost pˇribliˇzného ˇreˇsen´ı oproti pˇr´ıpadu, kdy se s´ıt’ neadaptuje. Vˇse bude opˇet ilustrováno na obrázc´ıch. Na závˇer srovnáme zkoumanou adaptivn´ı metodu s jinou adaptivn´ı metodou, a to anizotropn´ı metodou zjemˇ nován´ı v´ ypoˇcetn´ı s´ıtˇe popsanou v ˇclánku [2] a na internetové adrese http://www.karlin.mff.cuni.cz/ e dolejsi/angen/angen3.1.htm a na obrázc´ıch srovnáme pohyb uzl˚ u v´ ypoˇcetn´ı s´ıtˇe pˇri pouˇzit´ı obou adaptivn´ıch metod.

6

Kapitola 2 Popis metody V této kapitole si podrobnˇe pˇredstav´ıme adaptivn´ı metodu popsanou v ˇclánku [3], veˇskeré vzorce v této kapitole jsou citovány z [3].

2.1

Vytvoˇ ren´ı s´ıtˇ e pomoc´ı minimalizace funkcion´ alu

Pˇredpokládejme, ˇze N ≥ 1 oznaˇcuje dimenzi prostoru, v nˇemˇz se nacház´ı fyzik´ aln´ı oblast Ωp ⊂ IRN , tj. oblast, v n´ıˇz máme definovanou parciáln´ı diferenciáln´ı rovnici, pˇri jej´ımˇz ˇreˇsen´ı budeme adaptovat v´ ypoˇcetn´ı s´ıt’. Necht’ souˇradnice fyzikáln´ı oblasti jsou dány jako x = (x1 , x2 , . . . , xN ), pˇredpokládejme dále, ˇze fyzikáln´ı oblast je obrazem výpoˇcetn´ı oblasti Ωc ⊂ IRN pˇri vzájemnˇe jednoznaˇcném zobrazen´ı z Ωc na Ωp , oznaˇcme ξ = (ξ1 , ξ2 , . . . , ξN ) souˇradnice této v´ ypoˇcetn´ı oblasti. Transformaci souˇradnic z oblasti Ωc do oblasti Ωp oznaˇc´ıme x = x(ξ),

ξ ∈ Ωc ,

inverzn´ı zobrazen´ı z oblasti Ωp do oblasti Ωc oznaˇc´ıme ξ = ξ(x),

x ∈ Ωp .

Pˇri hledán´ı ˇreˇsen´ı zmiˇ nované parciáln´ı diferenciáln´ı rovnice provedeme triangulaci fyzikáln´ı oblasti Ωp a následnˇe hledáme pˇribliˇzné ˇreˇsen´ı. Pˇresné ˇreˇsen´ı m˚ uˇze m´ıt v nˇekter´ ych m´ıstech velk´ y gradient ˇci dokonce nespojitost a nen´ı-li triangulace oblasti dostateˇcnˇe jemná, nemus´ı m´ıt numerické ˇreˇsen´ı námi poˇzadovanou pˇresnost. Na druhou stranu, je-li triangulace oblasti rovnomˇerná a má-li b´ yt dostateˇcnˇe jemná, m˚ uˇzeme t´ımto poˇzadavkem na jemnost s´ıtˇe znaˇcnˇe zv´ yˇsit v´ ypoˇcetn´ı nároˇcnost, coˇz je neˇzádouc´ı. Proto je vhodné triangulaci zjemˇ novat pouze v tˇech m´ıstech, ve kter´ ych docház´ı k velk´ ym zmˇenám v ˇreˇsen´ı. Naˇs´ım c´ılem je tedy nalézt takové zobrazen´ı z pomocné v´ ypoˇcetn´ı oblasti Ωc do oblasti 7

Ωp , tj. takové rozloˇzen´ı prvk˚ u triangulace, aby se v takto problematick´ ych m´ıstech nacházel vˇetˇs´ı poˇcet prvk˚ u triangulace neˇzli v m´ıstech, kde v´ yˇse zm´ınˇené problémy nenastávaj´ı. Toho lze dosáhnout napˇr´ıklad minimalizac´ı následuj´ıc´ıho funkcionálu, jak je uvedeno v ˇclánku [3]. Minimalizujme tedy funkcionál (2.1)

E(ξ) =

N Z 1X ∇ξkT G−1 k ∇ξk dx, 2 k=1 Ωp

kde ∇ := (∂x1 , ∂x2 , . . . , ∂xN )T a Gk , 1 ≤ k ≤ N jsou dané symetrické pozitivnˇe definitn´ı matice zvané monitorovac´ı funkce. Nutnou podm´ınkou, aby funkcionál E nab´ yval minima, je Euler-Lagrangeova rovnice funkcionálu (2.1), a to (2.2)

∇.(G−1 k ∇ξk ) = 0,

1 ≤ k ≤ N.

Za monitorovac´ı funkce Gk zvol´ıme Gk = ωI, 1 ≤ k ≤ q N , kde I je jednotková matice a ω je kladná váhová funkce, napˇr. ω = 1 + |∇u| 2 , kde u je ˇreˇsen´ım parciáln´ı diferenciáln´ı rovnice, které hledáme. Rovnost (2.2) pak lze zapsat ve tvaru µ

(2.3)

2.1.1

∇.

¶

1 ∇ξk = 0, ω

1 ≤ k ≤ N.

Jednorozmˇ ern´ y pˇ r´ıpad

Zab´ yvejme se nyn´ı pˇr´ıpadem, kdy N = 1. Necht’ x oznaˇcuje fyzikáln´ı souˇradnice definované v odstavci 2.1. Necht’ x ∈ Ωp = [a, b], a, b ∈ IR, a < b. Triangulace intervalu [a, b] je v tomto jednoduchém pˇr´ıpadˇe dˇelen´ı intervalu. V této práci se budeme zab´ yvat adaptivn´ım zjemˇ nován´ım v´ ypoˇcetn´ı s´ıtˇe na intervalu [a, b] pˇri hledán´ı ˇreˇsen´ı parciáln´ı diferenciáln´ı rovnice tvaru (2.4)

∂ (f (u)) ∂u (x, t) + (x, t) = 0, ∂t ∂x

t ∈ (0, T ],

x ∈ [a, b],

s poˇca´teˇcn´ı podm´ınkou (2.5)

u(x, 0) = u0 (x),

x ∈ [a, b],

a vhodnou okrajovou podm´ınkou. Pˇredpokládáme, ˇze f ∈ C 1 (IR) a ˇze funkce u0 , která má kompaktn´ı nosiˇc 8

v IR, splˇ nuje u0 ∈ L∞ (IR) ∩ BV (IR). O pˇresném ˇreˇsen´ı u této rovnice budeme pˇredpokládat, ˇze u ∈ L∞ (IR) a ˇze u je slab´ ym po ˇca´stech hladk´ ym ˇreˇsen´ım rovnice (2.4) (definice pojmu slabé po ˇc´ astech hladké ˇreˇsen´ı viz kapitola 2.3 v knize [1]). Definujme si pomocnou v´ ypoˇcetn´ı oblast Ωc = [0, 1], oznaˇcme ξ ∈ [0, 1] souˇradnice na [0, 1] a definujme také transformace souˇradnic (2.6)

x = x(ξ), x(0) = a,

ξ ∈ [0, 1], x(1) = b.

ξ = ξ(x), ξ(a) = 0,

x ∈ [a, b], ξ(b) = 1.

a (2.7)

Nyn´ı najdeme zobrazen´ı z [0, 1] na [a, b] takové, abychom pomoc´ı nˇej zkonstruovali adaptivn´ı v´ ypoˇcetn´ı s´ıt’ na [a, b]. To provedeme pomoc´ı minimalizace funkcionálu (2.1), konkrétnˇe podm´ınky (2.3), která má v pˇr´ıpadˇe N = 1 tvar (ω −1 ξx )x = 0 na [a, b],

(2.8)

kde monitorovac´ı funkci ω definujeme jako (2.9)

ω(x) =

v u u t

Ã

!2

∂u 1+ (x) ∂x

,

x ∈ [a, b],

kde u je pˇresné ˇreˇsen´ı rovnice (2.4). V´ yraz na levé stranˇe rovnice (2.8) je funkce promˇenné x. Kdyby se nám podaˇrilo tuto rovnici transformovat do takového tvaru, aby levá strana transformované rovnice byla funkc´ı promˇenné ξ, diskrétn´ı pˇribliˇzné ˇreˇsen´ı x transformované rovnice nám následnˇe umoˇzn´ı, jak dále uvid´ıme, zkonstruovat uzly adaptivn´ıho dˇelen´ı intervalu [a, b]. Takto z´ıskané diskrétn´ı ˇreˇsen´ı bude totiˇz minimalizovat funkcionál (2.1), protoˇze rovnice (2.8) je pouze jednorozmˇernou podm´ınkou (2.2). Transformujme tedy rovnici (2.8). K tomu vyuˇzijeme následuj´ıc´ı postup. Necht’ f : X → IR, g : Y → X, g −1 : X → Y , f ∈ C 1 (X), g ∈ C 1 (Y ), g −1 ∈ C 1 (X), kde X, Y jsou uzavˇrené intervaly v IR a g je transformace souˇradnic z Y na X. Plat´ı tedy x = g(ξ), ξ = g −1 (x) ∀x ∈ −1 6= 0 v Y a dg 6= 0 v X. Rovnˇeˇz plat´ı ∀x ∈ X, ∀ξ ∈ Y X, ∀ξ ∈ Y , dg dξ dx df dg df dg −1 d (f ◦ g)(ξ) = (g(ξ)) . (ξ) = (x) . (g (x)) = dξ dx dξ dx dξ 9

df (x) . dx

= kde posledn´ı ˇclen

dg (g −1 (x)) dξ

1 dg −1 dx

1

=

6= 0 ∀x ∈ X, a proto plat´ı

dg −1 (x) dx

následuj´ıc´ı ekvivalence: ∀ξ ∈ Y :

,

(x)

d (f ◦ g)(ξ) = 0 dξ

⇔

∀x ∈ X :

df (x) = 0. dx

Aplikac´ı pˇredchoz´ıho postupu na funkce g(ξ) := x(ξ), g −1 (x) := ξ(x),

f (x) :=

ξ ∈ [0, 1], dg −1 (x), dx

1 . ω(x)

x ∈ [a, b]

a s vyuˇzit´ım následuj´ıc´ıho faktu, pˇri jehoˇz odvozen´ı se pouˇzije vˇeta o derivaci inverzn´ı funkce: dξ (x) 1 ξx = dx = . ω ω(x) ω(x(ξ))

1 dx (ξ) dξ

=

1 . e ω(ξ)

1 dx (ξ) dξ

=

1 , e ξ ωx

∀ξ ∈ [0, 1] ∀x ∈ [a, b], kde e ω(ξ) = ω(x(ξ)),

ξ ∈ [0, 1].

Dostáváme pro x = x(ξ), ξ = ξ(x), x ∈ [a, b], ξ ∈ [0, 1]: Ã −1

0 = (ω ξx )x ⇔

⇔

0=

1 e ξ ωx

e ξ xξ + ωx e ξξ = 0 ω

!

= ξ

⇔

e ξ xξ − ωx e ξξ −ω e ξ )2 (ωx

⇔

e ξ )ξ = 0. (ωx

Vid´ıme, ˇze rovnost (2.8) je ekvivalentn´ı rovnosti (2.10)

e ξ )ξ = 0 na [0, 1], (ωx

ˇc´ımˇz dostáváme poˇzadovanou transformaci rovnice (2.8), kde v´ yraz na levé stranˇe závisel na promˇenné x, na rovnici, kde je jiˇz v´ yraz na levé stranˇe funkc´ı promˇenné ξ. Konstrukci uzl˚ u adaptivn´ı s´ıtˇe na intervalu [a, b] provedeme v následuj´ıc´ım odstavci.

10

2.2

Popis algoritmu

Na základˇe odstavce 2.1.1 jiˇz nyn´ı m˚ uˇzeme zkonstruovat uzly adaptivn´ı v´ ypoˇcetn´ı s´ıtˇe v jednorozmˇerném pˇr´ıpadˇe parciáln´ı diferenciáln´ı rovnice (2.4). Toho dosáhneme nalezen´ım pˇribliˇzného ˇreˇsen´ı rovnice (2.11)

e ξ )ξ = 0 na [0, 1] (ωx

s okrajov´ ymi podm´ınkami (2.12)

x(0) = a,

x(1) = b.

Ze z´ıskaného diskrétn´ıho ˇreˇsen´ı uˇz m˚ uˇzeme jednoduˇse zkonstruovat uzly adaptivn´ı s´ıtˇe, jeˇz budou splˇ novat naˇse poˇzadavky, to znamená ˇze v m´ıstech, kde bude m´ıt funkce u velk´ y gradient nebo nespojitost, bude vyˇsˇs´ı poˇcet uzl˚ u dˇelen´ı intervalu [a, b], na nˇemˇz je funkce definována, neˇz v m´ıstech, kde funkce u takovou vlastnost nemá. Zaved’me oznaˇcen´ı. Uvaˇzujme interval [0, 1], necht’ ξ ∈ [0, 1] a bud’ 1 n ∈ IN, n ≥ 2. Necht’ h = n−1 je krok rovnomˇerného dˇelen´ı intervalu [0, 1] a ξj = (j − 1)h, 1 ≤ j ≤ n jsou uzly s´ıtˇe. Oznaˇcme 0 = t0 < t1 < . . . < tl = T , l ∈ IN dˇelen´ı intervalu [0, T ]. V kaˇzdém ˇcasovém kroku tk , k ∈ {0, . . . , l}, hledán´ı pˇribliˇzného ˇreˇsen´ı rovnice (2.4) oznaˇcme uk (x) := u(x, tk ), x ∈ [a, b]. Protoˇze ke zjemˇ nován´ı v´ ypoˇcetn´ı s´ıtˇe docház´ı pouze na základˇe pˇribliˇzn´ ych hodnot funkce u z´ıskan´ ych v daném ˇcasovém kroku tk , a zjemˇ nován´ı tedy nezávis´ı na pˇredchoz´ıch hodnotách funkce u v ˇcasech t1 , . . . , tk−1 , bez nebezpeˇc´ı nedorozumˇen´ı budeme dále index k u funkce uk vynechávat.

2.2.1

Metoda s´ıt´ı

Nejdˇr´ıve zkus´ıme nalézt ˇreˇsen´ı rovnice (2.11) metodou s´ıt´ı. Pˇribliˇzné ˇreˇsen´ı budeme hledat ve vnitˇrn´ıch bodech dˇelen´ı intervalu [0,1] takov´ ych, aby vˇsechny aproximace byly dobˇre definovány. Pro 1 ≤ j ≤ n oznaˇcme x(ξj ) ≈ xj aproximaci hodnot funkce x v bodˇe ξj a pro 1 ≤ j ≤ n − 1 x(ξj+ 1 ) ≈ xj+ 1 aproximaci hodnot funkce x 2

2

v bodech ξj+ 1 = ξj +ξ2j+1 , 1 ≤ j ≤ n − 1. Dále pro 1 ≤ j ≤ n − 1 oznaˇcme 2 u(xbj+ 1 ) ≈ uj+ 1 aproximaci hodnot funkce u v bodˇe xbj+ 1 = xj +x2 j+1 , 2 2 2 1 ≤ j ≤ n − 1. Aproximujme derivace rovnice (2.11) pomoc´ı centráln´ıch diferenc´ı: xξ (ξj ) =

xj+ 1 − xj− 1 dx 2 2 (ξj ) ≈ , dξ h

2 ≤ j ≤ n − 1.

Dále aproximujme druhou derivaci centráln´ımi diferencemi za pomoci 11

pˇredchoz´ı aproximace: e ξ )ξ (ξj ) ≈ (ωx

e j+ 1 )(xj+1 − xj ) − ω(ξ e j− 1 )(xj − xj−1 ) ω(ξ 2

2

h2

,

2 ≤ j ≤ n − 1. e j+ 1 ) = ω(xj+ 1 ), 2 ≤ j ≤ n − 1. Zb´ yvá aproximovat hodnoty ω(ξ 2 2 Zkusme to opˇet pomoc´ı centráln´ıch diferenc´ı takto, udˇelejme to pouze pro 3 ≤ j ≤ n − 2: v  2 u u uj+ 3 − uj− 1 u 2 2  e j+ 1 ) = ω(xj+ 1 ) ≈ t1 +  ω(ξ = 2 2 b b 3 1 xj+ − xj− 2

=

v u u t

Ã

1+4

!2

uj+ 3 − uj− 1 2

2

2

xj+2 + xj+1 − xj − xj−1

,

3 ≤ j ≤ n − 2.

Je vidˇet, ˇze abychom mohli aproximovat ω(xj+ 1 ), 3 ≤ j ≤ n − 2, 2 mus´ıme znát hodnoty uj+ 1 , 3 ≤ j ≤ n − 2, pˇriˇcemˇz uj+ 1 ≈ u(xbj+ 1 ) = 2 2 2 u( xj +x2 j+1 ), 3 ≤ j ≤ n − 2, tedy koeficienty uj+ 1 , 3 ≤ j ≤ n − 2 s´ıt’ové 2 rovnice s

µ 1+4

u −u j+ 3 j− 1 2 2 xj+2 +xj+1 −xj −xj−1

s

¶2 (xj+1 − xj ) −

µ 1+4

h2

u −u j+ 1 j− 3 2 2 xj+1 +xj −xj−1 −xj−2

¶2 (xj − xj−1 ) =0 ,

3 ≤ j ≤ n − 2,

závis´ı na s´ıt’ov´ ych uzlech xj , xj+1 , které vˇsak poˇc´ıtáme, tedy hodnoty uj+ 1 , 3 ≤ j ≤ n−2 pˇredem neznáme. Protoˇze z rovnice (2.9) dostáváme, 2 ˇze funkce ω závis´ı nejen na hodnotˇe promˇenné x, ale také na hodnotˇe derivace funkce u promˇenné x, a dále z konstrukce aproximace ω(xj+ 1 ), 2 3 ≤ j ≤ n − 2, je vidˇet, ˇze s´ıt’ová rovnice daná touto aproximac´ı bude vˇzdy obsahovat koeficienty, které budou záviset na uzlech s´ıtˇe, tj. na ˇreˇsen´ı, které hledáme. To znaˇcnˇe komplikuje v´ ypoˇcet, nebot’ dostáváme soustavu nelineárn´ıch rovnic a nalezen´ı ˇreˇsen´ı takové soustavy je v´ ypoˇcetnˇe mnohem nároˇcnˇejˇs´ı neˇz u soustavy lineárn´ı. Uˇzit´ı metody s´ıt´ı na problém (2.11) tedy nen´ı vhodné, a proto pouˇzijeme následuj´ıc´ı, i kdyˇz ménˇe efektivn´ı metodu, a to metodu umˇelého ˇcasu.

2.2.2

Metoda umˇ el´ eho ˇ casu

Metoda umˇelého ˇcasu je zaloˇzena na hledán´ı stacionárn´ıho ˇreˇsen´ı u ´lohy, kterou z´ıskáme modifikac´ı rovnice (2.11). Zaved’me funkci x(ξ, τ ) 12

promˇenn´ ych ξ a τ , kde τ je tzv. umˇelý ˇcas, bez nebezpeˇc´ı nedorozumˇen´ı ji budeme znaˇcit stejnˇe jako funkci p˚ uvodn´ı, a ˇreˇsme rovnici (2.13)

e ξ )ξ , xτ = (ωx

ξ ∈ (0, 1), τ > 0

s okrajov´ ymi podm´ınkami x(0, τ ) = a, x(1, τ ) = b,

τ >0

a poˇca´teˇcn´ı podm´ınkou x(ξ, 0) = (b − a)ξ + a,

ξ ∈ [0, 1],

která odpov´ıdá poˇcáteˇcn´ımu rovnomˇernému rozloˇzen´ı uzl˚ u v´ ypoˇcetn´ı s´ıtˇe na [a, b]. Bude-li ˇcasová derivace funkce x(ξ, τ ) nulová, bude schéma zaloˇzené na metodˇe umˇelého ˇcasu ˇreˇsit rovnici (2.11). Oznaˇcme ∆τ > 0 velikost kroku diskretizace v promˇenné τ , bud’ [ν] τν = ν ∆τ , ν ∈ IN0 , krok této diskretizace, definujme x(ξj , τν ) ≈ xj , 1 ≤ j ≤ n, ν ∈ IN0 , aproximaci funkce x(ξ, τ ) v bodˇe (ξj , τν ), 1 ≤ j ≤ n, [ν] [ν] ν ∈ IN0 . Pro 1 ≤ j ≤ n − 1 a ν ∈ IN0 oznaˇcme u(xbj+ 1 ) ≈ uj+ 1 apro2 [ν]

2

[ν]

xj +xj+1 , 2

[ν] xbj+ 1 2

[ν] xbj+ 1 , 2

kde = 1 ≤ j ≤ ximaci hodnot funkce u v bodˇe e j+ 1 ) z odstavce 2.2.1 nyn´ı budeme znaˇ n − 1, ν ∈ IN0 . Funkci ω(ξ cit jako 2 ω(uj+ 1 ), 1 ≤ j ≤ n − 1, aby byla patrná závislost funkce ω na funkci u 2 (ve skuteˇcnosti ovˇsem funkce ω závis´ı na ∇u). Diskretizujme rovnici (2.13), dostaneme: [ν]

(2.14)

[ν]

[ν]

[ν]

[ν]

[ν]

[ν] ω(uj+ 1 )(xj+1 − xj ) − ω(uj− 1 )(xj − xj−1 ) − xj 2 2 = , ∆τ h2

[ν+1]

xj

2 ≤ j ≤ n − 1,

ν ∈ IN0

z ˇcehoˇz dostáváme diferenˇcn´ı schéma (2.15)

[ν+1]

xj

[ν]

= xj +

∆τ [ν] [ν] [ν] [ν] [ν] [ν] [ω(uj+ 1 )(xj+1 − xj ) − ω(uj− 1 )(xj − xj−1 )] , 2 h 2 2 2 ≤ j ≤ n − 1,

s okrajov´ ymi podm´ınkami Funkci ω(uj+ 1 ) = 2

v u u t

[ν] x0

ν ∈ IN0

= a a x[ν] n = b ∀ν ∈ IN0 .

Ã

!2

∂u (xb 1 ) 1+ ∂x j+ 2

,

1≤j ≤n−1

aproximujeme napˇr´ıklad pomoc´ı centráln´ıch diferenc´ı.

13

Nyn´ı ukáˇzeme, ˇze schéma (2.15) ˇreˇs´ı rovnici (2.11). Skuteˇcnˇe, jestliˇze [ν+1] [ν] | xj −xj | ≈ 0, je pravá strana rovnice (2.14) také pˇribliˇznˇe rovna nule [ν+1] a uzly {xj } splˇ nuj´ı difereˇcn´ı schéma pro pˇribliˇzné ˇreˇsen´ı rovnice (2.11).

V praxi se tato podm´ınka aplikuje tak, ˇze pˇri hledán´ı stacionárn´ıho ˇreˇsen´ı rovnice (2.13) pˇredep´ıˇseme toleranci ε a provád´ıme iterace schématu [ν+1] [ν] (2.15) pro ν = 0, 1, . . . , dokud nen´ı splnˇeno | xj − xj | < ε. Schéma (2.15) m˚ uˇzeme pˇrepsat do tvaru (2.16)

[ν+1]

xj

[ν]

[ν]

[ν]

= αj+ 1 xj+1 + (1 − αj+ 1 − αj− 1 )xj + αj− 1 xj−1 , 2

2

2

2

2 ≤ j ≤ n − 1, kde ν ∈ IN0 je krok iterace a αj+ 1 = 2

∆τ [ν] ω(uj+ 1 ), 2 2 h

2 ≤ j ≤ n − 1.

Uvaˇzujeme následuj´ıc´ı podm´ınku stability max αj+ 1 ≤

(2.17)

1≤j≤n−1

2

1 . 2

V praxi je ovˇsem schéma (2.16) pˇr´ıliˇs pomalé a nedává dobré v´ ysledky. M´ısto schématu (2.16) proto budeme uvaˇzovat schéma [ν]

(2.18)

[ν+1] xj

=

[ν]

[ν]

[ν+1]

ω(uj+ 1 )xj+1 + ω(uj− 1 )xj−1 2

2

[ν] ω(uj+ 1 ) 2

+

[ν] ω(uj− 1 ) 2

,

2 ≤ j ≤ n − 1.

V praktick´ ych v´ ypoˇctech je tato metoda mnohem rychlejˇs´ı. Schéma (2.18) lze odvodit z pˇredeˇslého schématu. Vˇsimnˇeme si, ˇze ve vzorci (2.16) se [ν+1] [ν] [ν] [ν] hodnota uzlu xj z´ıská pomoc´ı jisté kombinace uzl˚ u xj+1 , xj a xj−1 , tedy pro v´ ypoˇcet hodnoty uzlu v iteraci ν + 1 pouˇz´ıváme hodnoty pouze z pˇredchoz´ıho kroku ν. Pokusme se v´ ypoˇcet urychlit t´ım, ˇze k v´ ypoˇctu [ν+1] xj pouˇzijeme jiˇz známé hodnoty z kroku ν + 1. Dostaneme rovnici [ν+1]

xj

[ν]

[ν+1]

= αj+ 1 xj+1 + (1 − αj+ 1 − αj− 1 )xj 2

2

2

[ν+1]

+ αj− 1 xj−1 , 2

2 ≤ j ≤ n − 1, [ν+1]

z n´ıˇz si vyjádˇr´ıme xj a z´ıskáme rovnici (2.18). V ˇclánku [3] se pouˇz´ıvá jeˇstˇe tzv. zhlazován´ı“ monitorovac´ı funkce ” ω ve tvaru 1 ω(uj+ 1 ) ← [ω(uj+ 3 ) + 2ω(uj+ 1 ) + ω(uj− 1 )], 2 2 2 2 4 14

2 ≤ j ≤ n − 2.

Funkce ω má totiˇz takov´ y pr˚ ubˇeh, ˇze nab´ yvá velk´ ych hodnot v m´ıstech, ve kter´ ych nab´ yvá velk´ ych hodnot gradient u a je-li gradient u mal´ y, je funkce ω malá. V závislosti na hodnotách funkce ω se následnˇe adaptuje s´ıt’, ˇc´ım vyˇsˇs´ı jsou tyto hodnoty, t´ım v´ıce uzl˚ u se v daném m´ıstˇe nashromáˇzd´ı. Problém nastává, je-li gradient u v jistém smyslu velk´ y. V takovém pˇr´ıpadˇe nab´ yvá funkce ω v daném m´ıstˇe extrémnˇe vysok´ ych hodnot a pro splnˇen´ı podm´ınky stability (2.17) je nutné pouˇz´ıt velmi mal´ y krok umˇelého ˇcasu ∆τ , coˇz zpomaluje v´ ypoˇcet. Vhodné pouˇzit´ı v´ yˇse zm´ınˇeného zhlazován´ı má za následek, ˇze k nab´ yván´ı takto extrémn´ıch hodnot u funkce ω nedocház´ı a podm´ınka stability je splnˇena i pro vyˇsˇs´ı hodnoty ∆τ , a tedy v´ ypoˇcet prob´ıhá rychleji. Pˇri v´ ypoˇctu nov´ ych hodnot uzl˚ u potˇrebujeme obvykle vykonat nˇekolik iterac´ı schématu (2.18), abychom z´ıskali dostateˇcnˇe jemnou s´ıt’. V takovém pˇr´ıpadˇe vˇsak vzniká problém. Zmˇen´ıme-li v´ ypoˇcetn´ı s´ıt’, neznáme hodnoty pˇribliˇzného ˇreˇsen´ı, které jsou definované na pˇredchoz´ı s´ıti, na s´ıti nové. Jedn´ım z ˇreˇsen´ı je hodnoty na p˚ uvodn´ı s´ıti pˇrepoˇc´ıtat na hodnoty pro novou s´ıt’ v kaˇzdé iteraci zjemˇ nován´ı s´ıtˇe tak, aby nové hodnoty splˇ novaly urˇcité podm´ınky. My pouˇzijeme následuj´ıc´ı postup. Oznaˇcme [ν] uj+ 1 hodnotu uj+ 1 v ν-té iteraci adaptivn´ı metody, 1 ≤ j ≤ n − 1, 2

2

[ν+1]

ν ∈ IN0 . Hodnotu uj+ 1 spoˇcteme jako 2

[ν+1]

(2.19)

[ν] [ν]

[ν]

[ν]

[ν]

uj+ 1 = βj uj+ 1 − γj ((ccu)j+1 − (ccu)j ), 2

2

1 ≤ j ≤ n − 1, kde

[ν]

[ν+1]

[ν+1] −1

γj = (xj+1 − xj

)

[ν]

,

[ν]

[ν]

[ν]

βj = γj (xj+1 − xj ),

1 ≤ j ≤ n − 1, a numerický tok (ccu)j poˇc´ıtáme jako [ν]

[ν]

[ν]

(ccu)j =

|cj | + cj − (u+ (uj − u− j + uj ) − j ), 2 2

2 ≤ j ≤ n − 2.

Protoˇze ∀ν ∈ IN0 plat´ı x1 = a a xn = b, poloˇzme [ν]

[ν]

(ccu)1 := 0,

(ccu)n−1 := 0.

V definici numerického toku dále máme [ν]

[ν]

[ν+1]

cj = xj − xj a

,

1 ≤ j ≤ n,

1 e u± j = uj± 12 + (xj − xj±1 )Sj± 12 , 2

15

3 ≤ j ≤ n − 2,

pˇriˇcemˇz hodnoty Sej+ 1 definujeme 2

j+ 12

=

¶ + − e e sgn(Sj+ 1 ) + sgn(Sj+ 1 )

+ e− |Sej+ | 1S j+ 1

µ

Se

2

2

2

2

+ e− | |Sej+ 1 | + |S j+ 1 2

+ Sej+ 1 = 2

uj+ 3 − uj+ 1 2

2

xj+ 3 − xj+ 1 2

− Sej+ 1 = 2

,

2 ≤ j ≤ n − 2,

2

,

1 ≤ j ≤ n − 2,

,

2 ≤ j ≤ n − 1.

2

uj+ 1 − uj− 1 2

2

xj+ 1 − xj− 1 2

2

Hodnoty u± akladˇe j pro j ∈ {1, 2, n − 1, n} definujeme analogicky na z´ okrajov´ ych podm´ınek ˇreˇsené parciáln´ı diferenciáln´ı rovnice. Jedn´ım z d˚ uvod˚ u pouˇzit´ı této metody pˇrepoˇc´ıtáván´ı hodnot uj+ 1 je, ˇze 2 splˇ nuje geometrický zákon zachován´ı, tj. plat´ı n−1 X

[ν+1]

[ν+1]

(xj+1 − xj

[ν+1]

)uj+ 1 = 2

j=1

n−1 X j=1

[ν]

[ν]

[ν]

(xj+1 − xj )uj+ 1

2

∀ν ∈ IN0 .

Algoritmus ˇ reˇ sen´ı PDR s vyuˇ zit´ım adaptivn´ı s´ıtˇ e: 0. Inicializace Zvolme n ∈ IN, n ≥ 2, poˇcet uzl˚ u fyzikáln´ı oblasti Ωp = [a, b]. Definujme rovnomˇerné rozloˇzen´ı bod˚ u ξj = (j − 1)h, 1 ≤ j ≤ n, 1 v oblasti Ωc = [0, 1], h = n−1 . Definujme dˇelen´ı intervalu [0, T ] jako 0 = t0 < t1 < t2 < . . . < tl = T pro nˇejaké l ∈ IN . Bud’ ν ∈ IN0 a k ∈ IN0 . Necht’ {xkj }nj=1 oznaˇcuje uzly v´ ypoˇcetn´ı s´ıtˇe tvoˇr´ıc´ı v reálném ˇcase k k k tk koneˇcné objemy Dj = [xj , xj+1 ] na [a,b] pro 1 ≤ j ≤ n − 1. Oznaˇcme aproximaci u( . , tk )|Djk ≈ ukj+ 1 , 1 ≤ j ≤ n − 1, pˇresného 2

ˇreˇsen´ı u rovnice (2.4) na koneˇcném objemu Djk = [xkj , xkj+1 ] v reálném ˇcase tk . k[ν] Necht’ xj , 1 ≤ j ≤ n, jsou uzly v´ ypoˇcetn´ı s´ıtˇe v ν-tém kroku adaptace s´ıtˇe (2.18) - (2.19), kterou provád´ıme v reálném ˇcase tk , k[0] pˇr´ısluˇsné poˇcáteˇcn´ı rozloˇzen´ı uzl˚ u je xj := xkj . Pro jednoduchost k[ν] [ν] budeme index k vynechávat, tj. m´ısto xj budeme psát xj . k[ν]

Symbolem uj+ 1 (a opˇet pro jednoduchost vynecháváme v horn´ım 2

indexu k) budeme rozumˇet aproximaci u( . , tk )|[x[ν] ,x[ν] j

j+1 ]

[ν]

≈ uj+ 1 , 2

1 ≤ j ≤ n − 1, pˇresného ˇreˇsen´ı u rovnice (2.4) na koneˇcném objemu [ν] [ν] [xj , xj+1 ] pˇr´ısluˇsnému reálnému ˇcasu tk a ν-tému kroku adaptace 16

[ν]

s´ıtˇe (2.18) - (2.19) v ˇcase tk , tedy uj+ 1 znaˇc´ı ν-tou iteraci schématu (2.19) s poˇca´teˇcn´ı podm´ınkou

[0] uj+ 1 2

2

:= ukj+ 1 . 2

Definujme poˇca´teˇcn´ı rovnomˇerné rozloˇzen´ı uzl˚ u x0j , 1 ≤ j ≤ n, v´ ypoˇcetn´ı s´ıtˇe ve fyzikáln´ı oblasti Ωp = [a, b]. Vypoˇc´ıtejme hodnoty u0j+ 1 jako integráln´ı pr˚ umˇer poˇcáteˇcn´ı pod2 m´ınky u(x, 0) = u0 (x) na pˇr´ısluˇsném intervalu, tj. u0j+ 1 2

xZj+1 1 = u0 (x)dx , xj+1 − xj x

1 ≤ j ≤ n − 1.

j

[0]

[0]

Poloˇzme ν = 0 a k = 0 a dále xj := x0j pro 1 ≤ j ≤ n a uj+ 1 := 2

u0j+ 1 pro 1 ≤ j ≤ n − 1. 2

1. Adaptace s´ıtˇ e [ν] [ν+1] Adaptujme s´ıt’ {xj } na novou s´ıt’ {xj } podle (2.18) a pˇrepoˇc´ı[ν] [ν+1] tejme hodnoty {uj+ 1 } na hodnoty {uj+ 1 } pomoc´ı (2.19). Poloˇzme 2 2 ν := ν + 1 a opakujme tento krok pro pˇredem dan´ y poˇcet iterac´ı [m] m. V praxi vol´ıme poˇcet iterac´ı m = 1 - 5. Dostaneme {xj } a [m] {uj+ 1 }. 2

2. V´ ypoˇ cet pˇ ribliˇ zn´ eho ˇ reˇ sen´ı PDR (2.4) [m] [m] k Poloˇzme uj+ 1 := uj+ 1 pro 1 ≤ j ≤ n − 1 a xkj := xj pro 1 ≤ 2 2 j ≤ n. Uˇzit´ım metody koneˇcn´ ych objem˚ u vypoˇc´ıtejme na základˇe k k v ˇcase tk+1 . hodnot uj+ 1 a s´ıtˇe {xj } pˇribliˇzné ˇreˇsen´ı PDR uk+1 j+ 1 2

2

3. Ukonˇ covac´ı podm´ınka [0] [0] Jestliˇze tk+1 ≤ T , definujme uj+ 1 := uk+1 a xj := xkj , poloˇzme j+ 12 2 ν = 0, k = k + 1 a jdˇeme na krok 1.

17

Kapitola 3 Numerick´ e experimenty V pˇredchoz´ıch kapitolách jsme adaptivn´ı metodu zkoumali teoreticky. V této kapitole se zamˇeˇr´ıme na praktickou ˇcást. Pˇredvedeme si pouze adaptivn´ı metodu (2.18), protoˇze metoda umˇelého ˇcasu (2.16) je znaˇcnˇe neefektivn´ı a nebudeme se j´ı tedy zab´ yvat. Pomoc´ı vhodnˇe definovan´ ych ’ funkc´ı si ukáˇzeme, ˇze se s´ıt adaptuje v závislosti na poloze nespojitosti, tj. ˇze se uzly shlukuj´ı“ v m´ıstˇe nespojitosti, poté adaptivn´ı s´ıt’ apli” kujeme pˇri ˇreˇsen´ı nˇekter´ ych parciáln´ıch diferenciáln´ıch rovnic a nakonec srovnáme pohyb uzl˚ u adaptivn´ı s´ıtˇe pˇri pouˇzit´ı algoritmu uˇz´ıvaj´ıc´ıho schéma (2.18) a algoritmu anizotropn´ıho zjemˇ nován´ı v´ ypoˇcetn´ı s´ıtˇe popsaného v ˇclánku [2].

3.1

Funkce s pohybuj´ıc´ı se nespojitost´ı

V tomto odstavci budeme vhodnou definic´ı funkc´ı, u nichˇz známe jejich pr˚ ubˇeh, modelovat situace, kdy se nespojitost funkce v pr˚ ubˇehu ˇcasu pohybuje nebo se velikost skoku v m´ıstˇe nespojitosti zvˇetˇsuje ˇci naopak zmenˇsuje. Na tyto funkce následnˇe aplikujeme metodu adaptivn´ıho zjemˇ nován´ı s´ıtˇe a na obrázc´ıch si ukáˇzeme, jak funguje. Pˇ r´ıklad 3.1. (Pohybuj´ıc´ı se nespojitost) Definujme si funkci

f1 (x, t) =

 (  x ,     x − 1,   (       

x , x−1 ,

0 ≤ x ≤ 1+t 1+t < x ≤ 3

,

0 ≤ t ≤ 1,

0 ≤ x ≤ 3−t 3−t < x ≤ 3

,

1 < t ≤ 3.

Pro pevné t ∈ [0, 3] je funkce f1t (x) = f1 (x, t), x ∈ [0, 3] po ˇcástech lineárn´ı, pˇriˇcemˇz graf funkce f1t obsahuje pouze jeden bod nespojitosti, jehoˇz poloha závis´ı na hodnotˇe parametru t, parametr t budeme naz´ yvat 18

ˇcas. Nespojitost v ˇcase t = 0 se nacház´ı v bodˇe x = 1 (viz obrázek 3.1), s rostouc´ım t se posouvá do bodu x = 2 pro t = 1 a následnˇe se vrac´ı do bodu x = 1 v ˇcase t = 3. V analogii ke kapitole 2 funkce f1 (x, t) odpov´ıdá funkci u(x, t) ze zmin ˇované kapitoly, funkce f1t tedy simuluje v´ yvoj pˇribliˇzného ˇreˇsen´ı nˇejaké parciáln´ı diferenciáln´ı rovnice v reálném ˇcase t, a to následuj´ıc´ım zp˚ usobem. Funkci f1t budeme zkoumat v pr˚ ubˇehu ˇcasu pro t = 0 aˇz t = 3 s ˇcasov´ ym krokem 0.1 (coˇz odpov´ıdá ˇcasovému kroku pˇri ˇreˇsen´ı PDR). Máme-li D0 poˇca´teˇcn´ı rovnomˇernou s´ıt’ v reálném ˇcase t = 0 tvoˇrenou body xj , 1 ≤ j ≤ 50, pak po proveden´ı pˇeti iterac´ı schématu (2.18) (2.19) dostaneme novou s´ıt’ D1 , která vypadá tak, ˇze se uzly s´ıtˇe posouvaj´ı k bodu nespojitosti, jehoˇz poloha odpov´ıdá ˇcasu t = 0 (na obrázku 3.2 ovˇsem tˇemto uzl˚ um odpov´ıdá zobrazen´ı uzl˚ u v ˇcase t = 0.1, nebot’ na obrázku znázorˇ nujeme polohu uzl˚ u odpov´ıdaj´ıc´ı adaptaci s´ıtˇe t na základˇe hodnot f1 z pˇredchoz´ıho ˇcasového kroku, tedy rozloˇzen´ı uzl˚ u s´ıtˇe v ˇcase t odpov´ıdá tomu, jak´ ych hodnot funkce f1t nab´ yvala v kroku pˇredchoz´ım). Tato s´ıt’ se následnˇe pouˇzije k nalezen´ı pˇribliˇzného ˇreˇsen´ı dané PDR, v naˇsem pˇr´ıpadˇe jsou ovˇsem t´ımto tzv. ˇreˇsen´ım“ hodnoty ” f1 (xj , t) pro xj ∈ D1 a t = 0.1. Nyn´ı opˇet na s´ıt’ D1 aplikujeme pˇet iterac´ı schématu (2.18) - (2.19) a dostaneme s´ıt’ D2 , pro kterou zjist´ıme hodnoty f1 (xj , 0.2). Dalˇs´ı v´ ypoˇcet prob´ıhá analogicky aˇz do t = 3. V´ ysledn´ y pohyb uzl˚ u je znázornˇen na obrázku 3.2. Z nˇej je patrné, ˇze uzly se shlukuj´ı v m´ıstˇe nespojitosti, která se pohybuje v ˇcase, coˇz je pˇresnˇe to, co od adaptivn´ı metody vyˇzadujeme. Pˇ r´ıklad 3.2. (Zvˇetˇsuj´ıc´ı se skok) Pro 0 ≤ t ≤ 1.5 definujme funkci (

f2 (x, t) =

(1 + t)x , 0 ≤ x ≤ 0.5 (1 + t)x − (1 + t) , 0.5 < x ≤ 1 .

Pro pevné t ∈ [0, 2] definujme f2t (x) = f2 (x, t), x ∈ [0, 1] Vid´ıme, ˇze pro 0 ≤ t ≤ 1.5 se bod nespojitosti funkce f2t nacház´ı neustále v bodˇe x = 0.5, pˇriˇcemˇz s rostouc´ım t se velikost skoku v bodˇe x = 0.5 (tj. pro pevné t ∈ [0, 1.5] | lim f2t (x) − lim f2t (x)|) zvˇetˇsuje (viz obrázek 3.3). x→0.5−

x→0.5+

S t´ım souvis´ı i pohyb uzl˚ u s´ıtˇe na obrázku 3.4, které se se vzr˚ ustaj´ıc´ım ˇcasem shlukuj´ı u bodu x = 0.5. Celkov´ y poˇcet uzl˚ u s´ıtˇe je 30, ˇcasov´ y krok vol´ıme opˇet 0.1, provád´ıme 5 iterac´ı adaptivn´ı metody (2.18) - (2.19).

19

Pˇ r´ıklad 3.3. (Zmenˇsuj´ıc´ı se skok) Pro 0 ≤ t ≤ 2 definujme funkci (

f3 (x, t) =

(2 − t)x + 0.5t (2 − t)x − 2 + 1.5t

, 0 ≤ x ≤ 1 , 1 < x ≤ 2.

Funkce f3t (x) = f3 (x, t), x ∈ [0, 2] pro pevné t ∈ [0, 2] se narozd´ıl od funkce f2t (x) chová tak, ˇze se nyn´ı skok mezi hodnotami v bodˇe x = 1 s rostouc´ım ˇcasem nezvˇetˇsuje, n´ ybrˇz naopak zmenˇsuje, dokonce v ˇcase t t = 2 je funkce f3 (x) ≡ 1 pro x ∈ [0, 2] (viz obrázek 3.5). Uzly se tedy zpoˇcátku pohybuj´ı smˇerem k nespojitosti,v d˚ usledku zmenˇsován´ı velikosti skoku se vˇsak jejich pohyb zpomaluje, dokonce pro t → 2 se uzly pohybuj´ı smˇerem od nespojitosti, velikost skoku uˇz je totiˇz tak malá, ˇze nen´ı nutné v tomto m´ıstˇe s´ıt’ adaptovat, pohyb uzl˚ u je znázornˇen na obrázku 3.6. Poˇcet uzl˚ u s´ıtˇe je 30, ˇcasov´ y krok je opˇet 0.1, provád´ıme 4 iterace adaptivn´ı metody (2.18) - (2.19). Pˇ r´ıklad 3.4. (Pohybuj´ıc´ı se nespojitost) Pro 0 ≤ t ≤ 2 definujme funkci g(t) = 0.5 + (exp(−1) − exp(−1 − t0.7 )) a dále ( x , 0 ≤ x ≤ g(t) f4 (x, t) = x−1 , g(t) < x ≤ 2 . Nespojitost funkce f4t (x) = f4 (x, t), x ∈ [0, 2], pro pevné t ∈ [0, 2], se nacház´ı v bodˇe x = g(t) pro t ∈ [0, 2], z ˇcehoˇz je vidˇet, ˇze se poloha nespojitosti s rostouc´ım ˇcasem mˇen´ı, avˇsak rychlost pohybu se pro t → 2 zmenˇsuje (viz obrázek 3.7), coˇz dobˇre modeluje situaci, s n´ıˇz se setkáme napˇr´ıklad pˇri ˇreˇsen´ı Burgersovy rovnice, a to, ˇze se poloha nespojitosti od urˇcitého ˇcasového okamˇziku mˇen´ı uˇz velmi pomalu. Poˇcet uzl˚ u s´ıtˇe je 50, ˇcasov´ y krok je volen 0.1, provád´ıme 5 iterac´ı adaptivn´ı metody (2.18) - (2.19). Pohyb uzl˚ u s´ıtˇe viz obrázek 3.8.

20

2

1.8

1.6

1.4

1.2

1

0.8

0.6

0.4

0.2

0

0

0.5

1

1.5

2

2.5

3

Obrázek 3.1: Pˇr´ıklad (3.1): Graf funkce f1tk (x) v ˇcase tk = 0.1k, k = 0, 1, 2, . . . , 10, osa x zn´ azorˇ nuje hodnoty promˇenné x, osa y hodnoty funkce f1tk (x). 3

2.5

2

1.5

1

0.5

0

0

0.5

1

1.5

2

2.5

3

Obrázek 3.2: Pˇr´ıklad (3.1): Pohyb uzl˚ u s´ıtˇe v z´ avislosti na pohybuj´ıc´ı se nespojitosti, osa x zn´ azorˇ nuje souˇradnice uzl˚ u, osa y ˇcas t.

21

1

0.5

0

−0.5

−1

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Obrázek 3.3: Pˇr´ıklad (3.2): Graf funkce f2t (x) v ˇcase t = 0 (ˇcárkovaná ˇcára) a v ˇcase t = 1.5 (pln´ a ˇc´ ara), osa x zn´ azorˇ nuje hodnoty promˇenné x, osa y hodnoty funkce f2t (x). 1.5

1

0.5

0

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Obrázek 3.4: Pˇr´ıklad (3.2): Pohyb uzl˚ u s´ıtˇe v z´ avislosti na zvˇetˇsuj´ıc´ı se velikosti skoku, osa x zn´ azorˇ nuje souˇradnice uzl˚ u, osa y ˇcas t.

22

2

1.8

1.6

1.4

1.2

1

0.8

0.6

0.4

0.2

0

0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

Obrázek 3.5: Pˇr´ıklad (3.3): Graf funkce f3t (x) v ˇcase t = 0 (ˇcárkovaná ˇcára) a v ˇcase t = 1.9 (pln´ a ˇc´ ara), osa x zn´ azorˇ nuje hodnoty promˇenné x, osa y hodnoty funkce f3t (x). 2

1.8

1.6

1.4

1.2

1

0.8

0.6

0.4

0.2

0

0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

Obrázek 3.6: Pˇr´ıklad (3.3): Pohyb uzl˚ u s´ıtˇe v z´ avislosti na zmenˇsuj´ıc´ı se velikosti skoku, osa x zn´ azorˇ nuje souˇradnice uzl˚ u, osa y ˇcas t.

23

1

0.8

0.6

0.4

0.2

0

−0.2

−0.4

−0.6

−0.8

−1

0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

Obrázek 3.7: Pˇr´ıklad (3.4): Graf funkce f4tk (x) v ˇcase tk = 0.1k, k = 0, 1, 2, . . . , 20, osa x zn´ azorˇ nuje hodnoty promˇenné x, osa y hodnoty funkce f4tk (x). 2

1.8

1.6

1.4

1.2

1

0.8

0.6

0.4

0.2

0

0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

Obrázek 3.8: Pˇr´ıklad (3.4): Pohyb uzl˚ u s´ıtˇe v z´ avislosti na pohybuj´ıc´ı se nespojitosti, osa x zn´ azorˇ nuje souˇradnice uzl˚ u, osa y ˇcas t.

24

3.2

Aplikace adaptivn´ı s´ıtˇ e pˇ ri ˇ reˇ sen´ı PDR

V této ˇca´sti se zamˇeˇr´ıme na aplikaci adaptivn´ı metody pˇri ˇreˇsen´ı parciáln´ıch diferenciáln´ıch rovnic a ukáˇzeme si na konkrétn´ıch pˇr´ıkladech, jak adaptivita v´ ypoˇcetn´ı s´ıtˇe zvyˇsuje pˇresnost pˇribliˇzného ˇreˇsen´ı ve srovnán´ı s pouˇzit´ım rovnomˇerné s´ıtˇe bˇehem celého v´ ypoˇctu. Budeme hledat pˇribliˇzné ˇreˇsen´ı u ´lohy (2.4) s poˇcáteˇcn´ı podm´ınkou (2.5). Pro nalezen´ı pˇribliˇzného ˇreˇsen´ı rovnice (2.4) pouˇzijeme tzv. MUSCL (monotone upstream-centered scheme for conservation laws) metodu koneˇcn´ ych objem˚ u, jej´ıˇz definici a schéma algoritmu lze nalézt na str. 494 ˇclánku [3]. Pod´ıvejme se jiˇz na konkrétn´ı pˇr´ıklady.

3.2.1

Burgersova rovnice

Pˇredpis Burgersovy rovnice je následuj´ıc´ı: 2

(3.1)

∂ ( u2 ) ∂u (x, t) + (x, t) = 0, ∂t ∂x

t > 0,

0 ≤ x ≤ 2π,

se spojitou poˇca´teˇcn´ı podm´ınkou (3.2)

u(x, 0) = 0.5 + sin(x),

x ∈ [0, 2π],

a periodickou okrajovou podm´ınkou. Uˇzit´ım MUSCL algoritmu s adaptivn´ım zjemˇ nován´ım s´ıtˇe metodou (2.18) - (2.19) splˇ nuj´ıc´ı geometrick´ y zákon zachován´ı jsme z´ıskali pˇribliˇzné ˇreˇsen´ı rovnice (3.1) v ˇcase t = 2, které je na obrázku 3.9 znázornˇeno pomoc´ı bod˚ u, pˇresné ˇreˇsen´ı je znázornˇeno plnou ˇcarou. Graf pˇresného ˇreˇsen´ı byl z´ıskán metodou charakteristik. Na obrázku 3.10 je zobrazeno pˇribliˇzné ˇreˇsen´ı, které jsme z´ıskali stejn´ ym MUSCL algoritmem za uˇzit´ı stejn´ ych parametr˚ u, avˇsak bez zjemˇ nován´ı s´ıtˇe, tzn. ˇze bˇehem celého v´ ypoˇctu jsme pracovali na rovnomˇerné s´ıti na intervalu [0, 2π]. V m´ıstech, na jejichˇz okol´ı je pˇresné ˇreˇsen´ı hladké, dostáváme pˇri pouˇzit´ı obou postup˚ u obdobné hodnoty pˇribliˇzného ˇreˇsen´ı, avˇsak v m´ıstˇe nespojitosti je pˇribliˇzné ˇreˇsen´ı z´ıskané pomoc´ı zjemˇ nován´ı s´ıtˇe mnohem pˇresnˇejˇs´ı (ve smyslu, kdy pˇresnost urˇcujeme jako absolutn´ı hodnotu rozd´ılu funkˇcn´ıch hodnot pˇresného a pˇribliˇzného ˇreˇsen´ı), narozd´ıl od rovnomˇerné s´ıtˇe, kdy na okol´ı bodu nespojitosti docház´ı k v´ yraznému zkreslen´ı pˇribliˇzného ˇreˇsen´ı. Pohyb uzl˚ u adaptivn´ı s´ıtˇe je zachycen na obrázku 3.11. V pr˚ ubˇehu v´ ypoˇctu docház´ı k jejich postupnému shromaˇzd’ován´ı kolem bodu nespojitosti a pˇribliˇznˇe od ˇcasu t = 1.3 uˇz jen nejbliˇzˇs´ı skupina uzl˚ u sleduje pohyb nespojitosti, uzly, které jsou dostateˇcnˇe daleko od bodu nespojitosti, se témˇeˇr nepohybuj´ı a zachovávaj´ı takˇrka rovnomˇerné rozloˇzen´ı, nebot’ 25

gradient pˇribliˇzného ˇreˇsen´ı je vzhledem k bodu nespojitosti v ostatn´ıch bodech znaˇcnˇe mal´ y. ’ V´ ypoˇcetn´ı s´ıt byla v tomto pˇr´ıpadˇe tvoˇrena 30-ti uzly, v kaˇzdém ˇcasovém kroku MUSCL algoritmu jsme pouˇzili 5 iterac´ı schématu (2.18) - (2.19) pro zjemnˇen´ı s´ıtˇe, ˇcasov´ y krok MUSCL metody byl volen 0.0005.

3.2.2

Riemann˚ uv probl´ em

Rovnice pro tzv. Riemann˚ uv problém pro hyperbolickou rovnici zákona zachován´ı je tvaru (3.3)

∂u ∂ (f (u)) (x, t) + (x, t) = 0, ∂t ∂x

t > 0,

−1 ≤ x ≤ 1,

kde

1 f (u) = (u2 − 1)(u2 − 4), 4 s nespojitou poˇca´teˇcn´ı podm´ınkou (3.4)

u(x, 0) = −2 sgn(x),

−1 ≤ x ≤ 1.

a okrajovou podm´ınkou

(3.5)

u(−1, t) = −2, u(1, t) = 2,

t > 0, t > 0.

Tato parciáln´ı diferenciáln´ı rovnice a graf jej´ıho pˇresného ˇreˇsen´ı jsou pˇrevzaty z ˇclánku [3]. Pˇresné ˇreˇsen´ı v pˇr´ıpadˇe této rovnice nen´ı spojité jiˇz od ˇcasu t = 0, jak je vidˇet z definice poˇca´teˇcn´ı podm´ınky a v ˇcase t = 1.2 obsahuje dva body nespojitosti. Opˇet jsme pouˇzili MUSCL algoritmus k nalezen´ı pˇribliˇzného ˇreˇsen´ı rovnice (3.3). Na v´ ypoˇcetn´ı s´ıti, která je tvoˇrena 50-ti uzly, jsme opˇet pomoc´ı adaptivn´ıho zjemˇ nován´ı z´ıskali mnohem pˇresnˇejˇs´ı ˇreˇsen´ı neˇz v pˇr´ıpadˇe pouˇzit´ı pouze rovnomˇerné s´ıtˇe po celou dobu v´ ypoˇctu, jak je vidˇet z obrázk˚ u 3.13 a 3.14. Ke zkreslen´ı pˇribliˇzného ˇreˇsen´ı v pˇr´ıpadˇe uˇzit´ı rovnomˇerné s´ıtˇe docház´ı opˇet pouze v m´ıstech nespojitosti pˇresného ˇreˇsen´ı. V pˇr´ıpadˇe adaptivn´ı s´ıtˇe se uzly akumuluj´ı na okol´ı m´ısta nespojitosti a absolutn´ı hodnota rozd´ılu hodnot pˇribliˇzného a pˇresného ˇreˇsen´ı je menˇs´ı neˇz v pˇr´ıpadˇe uˇzit´ı pouze rovnomˇerné s´ıtˇe. Pohyb uzl˚ u v´ ypoˇcetn´ı s´ıtˇe je znázornˇen na obrázku 3.12. V kaˇzdém ˇcasovém kroku MUSCL algoritmu jsme adaptovali v´ ypoˇcet’ n´ı s´ıt pouze pomoc´ı jedné iterace schématu (2.18) sak pouˇzili q- (2.19), avˇ θ 2 jsme modifikovanou monitorovac´ı funkci ω = 1 + θ|ux | s parametˇ rem θ > 0, v naˇsem pˇr´ıpadˇe jsme volili θ = 0.6. Casov´ y krok MUSCL algoritmu byl volen 0.0001. 26

1.6 1.4 1.2 1 0.8 0.6 0.4 0.2 0 −0.2 −0.4

0

Obrázek 3.9:

1

2

3

4

5

6

Graf pˇ resn´ eho (pln´ a ˇ c´ ara) a pˇ ribliˇzn´ eho (body) ˇ reˇsen´ı Burgersovy rovnice (3.1)

na intervalu [0, 2π] v ˇ case t = 2 z´ıskan´ eho pomoc´ı adaptivn´ıho zjemˇ nov´ an´ı v´ ypoˇ cetn´ı s´ıtˇ e.

1.6 1.4 1.2 1 0.8 0.6 0.4 0.2 0 −0.2 −0.4

0

Obrázek 3.10:

1

2

3

4

5

6

Graf pˇ resn´ eho (pln´ aˇ ca ´ra) a pˇ ribliˇzn´ eho (body) ˇreˇsen´ı Burgersovy rovnice (3.1)

na intervalu [0, 2π] v ˇ case t = 2 z´ıskan´ eho na rovnomˇ ern´ e v´ ypoˇ cetn´ı s´ıti.

27

2

1.8

1.6

1.4

1.2

1

0.8

0.6

0.4

0.2

0

0

1

2

3

4

5

6

Obrázek 3.11: Pohyb uzl˚u výpoˇcetn´ı s´ıtˇe v závislosti na ˇcase t bˇehem hledán´ı pˇribliˇzného ˇreˇsen´ı Burgersovy rovnice (3.1) na intervalu [0, 2π]. Osa x reprezentuje souˇ radnice uzl˚ u, osa y ˇ cas t.

1

0.8

0.6

0.4

0.2

0 −1

−0.8

−0.6

−0.4

−0.2

0

0.2

0.4

0.6

0.8

1

Obrázek 3.12: Pohyb uzl˚u výpoˇcetn´ı s´ıtˇe v závislosti na ˇcase t bˇehem hledán´ı pˇribliˇzného ˇreˇsen´ı rovnice (3.3) na intervalu [−1, 1]. Osa x reprezentuje souˇ radnice uzl˚ u, osa y ˇ cas t.

28

2

1.5

1

0.5

0

−0.5

−1

−1.5

−2 −1

−0.8

−0.6

−0.4

−0.2

0

0.2

0.4

0.6

0.8

1

Obrázek 3.13: Graf pˇresného (plná ˇca´ra) a pˇribliˇzného (body) ˇreˇsen´ı rovnice (3.3) na intervalu [−1, 1] v ˇ case t = 1.2 z´ıskan´ eho pomoc´ı adaptivn´ıho zjemˇ nov´ an´ı v´ ypoˇ cetn´ı s´ıtˇ e.

2

1.5

1

0.5

0

−0.5

−1

−1.5

−2 −1

−0.8

−0.6

−0.4

−0.2

0

0.2

0.4

0.6

0.8

1

Obrázek 3.14: Graf pˇresného (plná ˇca´ra) a pˇribliˇzného (body) ˇreˇsen´ı rovnice (3.3) na intervalu [−1, 1] v ˇ case t = 1.2 z´ıskan´ eho na rovnomˇ ern´ e v´ ypoˇ cetn´ı s´ıti.

29

3.3

Srovn´ an´ı r˚ uzn´ ych adaptivn´ıch metod

Na závˇer srovnáme adaptivn´ı metodu zkoumanou v této práci a anizotropn´ı adaptivn´ı metodu z ˇclánku [2] stejnˇe jako v odstavci 3.1 pro funkce s nespojitost´ı pohybuj´ıc´ı se v ˇcase. Z obrázk˚ u 3.15 - 3.18 je patrné, ˇze u obou metod docház´ı ke zhuˇstˇen´ı s´ıtˇe v m´ıstˇe nespojitosti a toto zhuˇstˇen´ı“ prvk˚ u s´ıtˇe se pohybuje spolu ” s nespojitost´ı. Na zbytku v´ ypoˇcetn´ı s´ıtˇe obˇe adaptivn´ı metody zachovávaj´ı rovnomˇerné dˇelen´ı, nedocház´ı k pohybu uzl˚ u v m´ıstech, kde nen´ı potˇreba s´ıt’ zjemnit. V tomto ohledu tedy obˇe metody vykazuj´ı obdobné vlastnosti a jsou vhodné k pouˇzit´ı napˇr´ıklad pˇri ˇreˇsen´ı parciáln´ıch diferenciáln´ıch rovnic, jak jsme vidˇeli na pˇredcházej´ıc´ıch pˇr´ıkladech. Ze srovnán´ı v´ ypoˇcetn´ı nároˇcnosti vycház´ı jako rychlejˇs´ı, zhruba o jeden ˇra´d, ˇ anizotropn´ı adaptivn´ı metoda, viz tabulka 3.1. Casovou nároˇcnost jsme zkoumali v následuj´ıc´ım smyslu. Pro kaˇzd´ y algoritmus zvláˇst’ jsme bˇehem pˇr´ısluˇsného ˇcasového intervalu adaptovali v´ ypoˇcetn´ı s´ıt’ pro funkce definované v pˇr´ıkladech (3.1) (3.4), v pˇr´ıpadˇe adaptivn´ı metody (2.18) - (2.19) jsme volili stejné parametry jako v pˇr´ıkladech (3.1) - (3.4), v pˇr´ıpadˇe anizotropn´ı adaptivn´ı metody jsme na stejném v´ ypoˇcetn´ım intervalu bˇehem stejného ˇcasového intervalu volili stejn´ y poˇcet uzl˚ u jako u pˇredcházej´ıc´ı metody, parametry byly voleny takto: MAXITER = 1, e1 = 100, e2 = 1 (v´ yznam parametr˚ u viz pˇriloˇzené CD a internetové stránky http://www.karlin.mff.cuni.cz/ e dolejsi/angen/angen3.1.htm). Tento krok v´ ypoˇctu s´ıtˇe jsme opakovali desetkrát a v´ ysledn´ y ˇcas, bˇehem kterého probˇehl v´ ypoˇcet, jsme vydˇelili deseti a tak z´ıskali hodnoty uvedené v tabulce 3.1. Pro urˇcen´ı v´ ypoˇcetn´ıho ˇcasu jsme uˇzili standardn´ı matlabovské pˇr´ıkazy tic a toc, v´ ysledné ˇcasy jsou uvedené v sekundách. Zkoumán´ı v´ ypoˇcetn´ı nároˇcnosti jsme provádˇeli na poˇc´ıtaˇci Intel(R) Pentium(R) M, 1.73GHz / 504 MB RAM v programu Matlab verze 7.2.0.232 (R2006a). Funkce f1 f2 f3 f4

Schéma (2.18) - (2.19)

0.2238 0.0791 0.0938 0.1080

Anizotropie 0.0227 0.0077 0.0097 0.0105

Tabulka 3.1: Srovnán´ı ˇcasové nároˇcnosti jednotlivých algoritm˚ u pˇri aplikaci pro funkce definované v pˇr´ıkladech (3.1) - (3.4). V levém sloupci je uvedena daná funkce, v prostˇredn´ım sloupci ˇcas schématu (2.18) - (2.19) a v pravém sloupci ˇcas anizotropn´ı adaptivn´ı metody, ˇcasy jsou uvedeny v sekundách.

30

3

2.5

2

1.5

1

0.5

0

0

0.5

Obrázek 3.15:

1

1.5

2

2.5

3

Funkce z pˇ r. (3.1): Pln´ aˇ ca ´ra zn´ azorˇ nuje pohyb uzl˚ u metody (2.18), teˇ ckovan´ a

pohyb uzl˚ u anizotropn´ı adaptivn´ı metody, osa x reprezentuje souˇ radnice uzl˚ u, osa y ˇ cas t.

1.5

1

0.5

0

0

0.1

Obrázek 3.16:

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1



31

2

1.8

1.6

1.4

1.2

1

0.8

0.6

0.4

0.2

0

0

0.2

Obrázek 3.17:

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2



2

1.8

1.6

1.4

1.2

1

0.8

0.6

0.4

0.2

0

0

0.2

Obrázek 3.18:

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2



32

Literatura [1] M. Feistauer, J. Felcman, and I. Straˇskraba. Mathematical and Computational Methods for Compressible Flow. Oxford University Press, Oxford, 2003. [2] J. Felcman and P. Kubera. Computational aspects of the mesh adaptation for the time marching procedure. In A. Berm´ udez de Castro, D. Gómez, P. Quintela, and P. Salgado, editors, Numerical Mathematics and Advanced Applications, ENUMATH 2005, pages 225– 232, Santiago de Compostela, Spain, July 2005 2006. Springer Verlag Berlin Heidelberg. ISBN-10 3-54034287-7 Springer Berlin Heidelberg New York, ISBN-13 978-3-540-34287-8 Springer Berlin Heidelberg New York. [3] H. Tang and T. Tang. Adaptive mesh methods for one- and twodimensional hyperbolic conservation laws. SIAM J. Appl. Math., pages 487–515, 2003.

33

Libor Kadrnka Schémata typu Moving Mesh pro řešení

Recommend Documents