大気環境シミュレーション

（Ｑ1） 各自excelを用いて、次の漸化式

x

_n+1

= x

_n2

– 2

の解の初期値依存性を調べよ．ｎは50まで。 (1) x₀=2.10 (2) x₀=2.15 (3) x₀=2.20 締切 2014年1月6日（月）夕方まで 提出先 3417室 オーバーフロー 失敗 ゴメンなさい

（Ｑ1）‘ 各自excelを用いて、次の漸化式

x

_n+1

= x

_n2

– 2

の解の初期値依存性を調べよ．ｎは50まで。 (1) x₀=1.330 (2) x₀=1.331 (3) x₀=1.332 (4) x₀=1.333 締切 2014年1月14日（火）第1時限まで 提出先 3417室

0

2 2



















x

f

a

x

f

c

t

f

を差分式にしてみる。 時間項 ＝ 移流項 ＋ 拡散項 2 2

x

f

a

x

f

c

t

f



















今回と次回にわたって、次の式を

計算してもらう

  0 2 2 2 2 1 1 1 1 1 1 1 1 1                     x f f f a x f f c t f f _jn _jn _jn _jn _jn _jn _jn 005 . 0 2 . 0 1 . 0 , 1 . 0       a c x t





_{ }



1



1 1 1 1 2 1 1 1 1 2 2 _         _{ }         n j n j n j n j n j n j n j f f f x t a f f x t c f f どうしてこのような数値を選ぶのか理解する 時間項 ： リープフロッグ 空間差分： _中央差分

） （ ） （ （○） ） （ ） （ （◯）                                         01 . 0 , , 1 . 0 , 1 . 0 ) 6 ( 01 . 0 , 5 . 1 , 1 . 0 , 1 . 0 ) 5 ( 01 . 0 , 8 . 0 , 1 . 0 , 1 . 0 ) 4 ( 01 . 0 , 2 . 0 , 1 . 0 , 0 . 1 ) 3 ( 01 . 0 , 2 . 0 , 0 . 1 , 1 . 0 ) 2 ( 01 . 0 , 2 . 0 , 1 . 0 , 1 . 0 ) 1 ( a f c x t a c x t a c x t a c x t a c x t a c x t n j





_{ }



1



1 1 1 1 2 1 1 1 1 2 2 _         _{ }         n j n j n j n j n j n j n j f f f x t a f f x t c f f 宿題: (1)から(6)までの（ｘ、ｔ、振幅）の図を描け 2014年1月20日（月） 午後5時まで 3417室

時間発展

 時間微分を含む次の方程式 について差分化する。時間間隔を Δt とし、n ステップ目の値を f n _{と書き表すと、以下の} ような差分式で表せる ) ( f g t f _   ) 3 ( ) ( ) 2 /( ) ( ) 2 ( ) ( / ) ( ) 1 ( ) ( / ) ( 1 1 1 1 1 n n n n n n n n n f g t f f f g t f f f g t f f               t_n t_n+1 t_n-1 t_n+2 Δt Δt Δt t

時間発展

(1) 陽解法（explicit scheme） nステップの値からn+1ステップの値が直ちに求まる。 解は必ずしも安定ではない (2) 陰解法（implicit scheme） n+1ステップの計算は複雑。解は安定 (3) リープ・フロッグ法（leap-frog scheme） n-1、nステップの値が必要。解は中立 ) 3 ( ) ( ) 2 /( ) ( ) 2 ( ) ( / ) ( ) 1 ( ) ( / ) ( 1 1 1 1 1 n n n n n n n n n f g t f f f g t f f f g t f f               t_n t_n+1 t_n-1 t_n+2 Δt Δt Δt t

線形移流方程式

 _{f についての移流方程式} 移流項 u(∂f/∂x) は非線形のため複雑なので、 ここでは簡単のため、移流速度は c (≥ 0) で一定であるとする。 移流項は線形化され、解析解が求まる。 0       x f u t f 0       x f c t f

線形移流方程式の解

において変数 X, T を のように置くと、  これを利用して 線形移流方程式 を書き直すと、 → t = 0 における f の分布を維持し たまま速さ c で 移動 0       x f c t f       t T ct x X                                       X X x X T x T x X c T X t X T t T t 0    T f

線形移流方程式の解

c

拡散方程式

 拡散方程式  拡散方程式の解は、上 に凸（_∂2_{f /∂x}2 _{< 0）のとこ} ろで減少（_{∂f /∂t < 0）し、} 下に凸（_∂2_{f /∂x}2 _{> 0）の} ところで増加（_{∂f /∂t > 0）} するので、分布を平滑化 させる解となる 2 2 x f a t f      a（≥ 0）は拡散係数で定数とする f（ｘ）

拡散方程式の厳密解

0 1 ) 0 ( sin ) ( ) , ( 2 2         t at H x k t H t x f x f a t f a（≥ 0）は拡散係数で定数とする x k e t x f e H e A H k a b H k a t d H d t b t b t b sin ) , ( 2 2            問題： ｘ方向には周期的と仮定する 解：

線形

移流拡散

方程式の厳密解

において変数 X, T を のように置くと、  これを利用して線 形移流拡散方程 式_{を書き直すと、} →速さ c で移動しな がら拡散するよう な解 0 2 2          x f a x f c t f       t T ct x X                                       X X x X T x T x X c T X t X T t T t 2 2 X f a T f     

線形移流拡散方程式の解

c

線形移流方程式の差分近似

（復習）  時間微分については、陽解法を採用  空間微分に中心差分を用いると  _fj±1/2 の補間を とすると、 t f f t f jn n j      1 x f f x f jn n j        ₂1 2 1 2 1 2 1 n j n j n j f f f _    x f f x f jn n j         2 1 1

線形移流方程式の差分近似

 線形移流方程式 を中心差分_{で差分近似すると、} 整理すると、 0       x f c t f 0 2 1 1 1        _{ }  x f f c t f f_jn _jn _jn _jn n j n j n j n j f x t c f f x t c f 1 _{1 1} 2 2           n n+1 j–1 j j+1

x

t

c

差分解の例（

中心差分

）

c = 0.3 m/s

t = 0

Δt = 1 s

Δx = 1 m

差分解の例（

中心差分

）

c = 0.3 m/s

t = 10

Δt = 1 s

Δx = 1 m

差分解の例（

中心差分

）

c = 0.3 m/s

t = 20

Δt = 1 s

差分解の例（

中心差分

）

c = 0.3 m/s

t = 30

Δt = 1 s

差分解の例（

中心差分

）

c = 0.3 m/s

t = 40

Δt = 1 s

差分解の例（

中心差分

）

c = 0.3 m/s

t = 50

Δt = 1 s

拡散方程式の差分近似

（復習）  拡散方程式 を陽解法と中心差分_で 差分近似すると、 整理すると、 2 2 x f a t f       2 1 1 1 2 x f f f a t f f_jn _jn _jn _jn _jn       _{ }        jn n j n j n j f x t a f x t a f x t a f 1 _{2 1} 1 2 _{2 2 1}                 n n+1 j–1 j j+1

x

t

拡散方程式の差分解

t = 0

a = 0.3

Δt = 1 s

拡散方程式の差分解

t = 10

a = 0.3

Δt = 1 s

拡散方程式の差分解

t = 20

a = 0.3

Δt = 1 s

拡散方程式の差分解

t = 30

a = 0.3

Δt = 1 s

拡散方程式の差分解

t = 40

a = 0.3

Δt = 1 s

拡散方程式の差分解

t = 50

a = 0.3

Δt = 1 s

拡散方程式の差分解

a = 0.3

Δt = 1 s

安定性－ノイマン法

（復習）  変数 f (x, t) を、空間方向についてフーリエ 変換すると  ノイマン（von Neumann）の安定条件： 時間発展とともにあらゆる波数 k の振幅 A_k(t) が増幅しないこと



 k k t ikx A t x f ( , ) ( )exp( )

安定性－ノイマン法

 差分点（x_{j = jΔx, tn = nΔt）における f の値} f_jn を とおき、_{λ = A}n+1 _{/ A}n _{の大きさを求めて以下の} ように判別する 全ての k について安定なら、無条件安定 ) exp(i k j x A f _jn  n        不安定 1 安定　 1 | |   

 λ = λ_r＋ ｉ λ_i （ λ_r，λ_iは実数）  |λ|:  これまで講義したことを復習する： ・移流スキーム：陽解法＋中央差分 ・拡散スキーム：陽解法＋中央差分 2 2

|

|







_r





_i

(1)移流方程式の安定性

 陽解法、中心差分での移流方程式の差分形式 についての波数 k の振幅の式は x f f c t f f _jn _jn _jn _jn        _{ }  2 1 1 1 n n n

A

x

x

k

ic

t

A

A















)

sin(

1

x

f

c

t

f













)

sin(

1

1

x

k

x

t

c

i

A

A

n n

















(1) 移流方程式の安定性

 振幅の大きさ |λ| は となり、従って、陽解法、中心差分による移流方 程式の差分解法は、無条件不安定_である 1 ) ( sin 1 | | 2 2             k x x t c  c = 0.3 m/s Δt = 1 s Δx = 1 m

(2)拡散方程式の安定性

 拡散方程式の差分形式:  陽解法、  現在値の拡散 についてフーリエ変換を適用し、安定性を調 べる。波数 k の振幅の式は

 

2 1 1 1 2 x f f f a t f f _jn _jn _jn _jn _jn       _{ }    n n n A x x k a t A A 2 1 1 ) cos( 2          sin 2 4 1 ₂ 2 1 x k x t a A A n n         2 cos 1 2 sin2     2 2 x f a t f     

(2) 拡散方程式の安定性

 安定条件は |λ| ≤ 1 であ るから  結局、あらゆる波数 k に ついて安定となる条件は  条件付き安定   sin 2 1 4 1 1 ₂ 2        k x x t a   a x t 2 2    a = 0.3 a = 0.5 a = 0.7 Δt = 1 (Δx)2_{/2a = 1.67} Δt = 1 (Δx)2_{/2a = 1.00} Δt = 1 (Δx)2_{/2a = 0.71}

どうしてマイナス×マイナスはプラスなのか

最初に負の数を習ったときに、負の数同士を掛け合わせ ると正の数になると教わった。例えば、（-2）×（-3） = 6。 不思議に思う人がずいぶん多い。

(3)移流方程式の安定性

 リープフロッグ法、中心差分での移流方程 式の差分形式 についての波数 k の振幅の式は x f f c t f f _jn _jn _jn _jn        _{ }  2 2 1 1 1 1 n n n

A

x

x

k

ic

t

A

A













 

)

sin(

2

1 1

x

f

c

t

f













0

1

)

sin(

2

2

















k

x

x

t

c

i

(3) 移流方程式の安定性

 振幅の大きさ λ は となり、従って、リープフロッグ法、中心差分に よる移流方程式の差分解法は、条件付安 定_である unstable neutral neutral x k x t c x k x t c x k x t c i 1 1 1 ) ( sin ) ( sin 1 ) sin( 2 2                             

(3) 移流方程式の安定性

リープフロッグ法、中心差分による移流方程 式の差分解法の安定条件： クーラン数： C = cΔt / Δx → 情報伝達距離（cΔt）と格子幅（Δx）の比。 一般に、C ≤ 1 は計算安定性の必要条件。こ れは Δx / Δt ≥ c であり、「情報が伝播する速 さ」が「実際の現象の進む速さ」以上でなけれ ばならないことを示す。この条件を_CFL （Courant-Friedrichs-Lewy）条件という

1







x

t

c

(4) 拡散方程式の安定性

 拡散方程式の差分形式:  リープフロッグ法、  現在値の拡散 についてフーリエ変換を適用し、安定性を調べる。 波数 k の振幅の式は

 

2 1 1 1 1 2 2 x f f f a t f f _jn _jn _jn _jn _jn        _{ }    n n n A x x k a t A A 2 1 1 1 ) cos( 2 2           sin 2 1 0 8 2 2 2         k x x t a 2 cos 1 2 sin2      2 2 x f a t f     

(4) 拡散方程式の安定性

 |λ| ＞ 1であるので、いつも不安定である。 リープフロッグ法で、拡散項に現在値を用い るのは無条件不安定_である。 拡散項には過去値_{をもちいるべき。}  2 1 1 1 1 ) cos( 2 2          n _n n A x x k a t A A     sin 2 8 1 0 2 sin 8 1 ₂ 2 ₂ 2 2 k x x t a x k x t a               

 

 

01 . 0 0125 . 0 8 1 . 0 1 . 0 , 1 . 0 1 8 1 2 sin 8 1 2 2 2                a a x t x t a x k x t a の値であれば拡散項に関しては 安定である

0 2 2          x f a x f c t f   0 2 2 2 2 1 1 1 1 1 1 1 1 1                     x f f f a x f f c t f f _jn _jn _jn _jn _jn _jn _jn





_{ }



1



1 1 1 1 2 1 1 1 1 2 2 _         _{ }         n j n j n j n j n j n j n j f f f x t a f f x t c f f ・時間： リープフロッグ法がよい （２次の精度、 移流スキームはCFL条件を満たせばよい） ・移流スキーム：中央差分 ・拡散項：過去値を使う

01 . 0 2 . 0 1 . 0 , 1 . 0       a c x t





_{ }



1



1 1 1 1 2 1 1 1 1 2 2 _         _{ }         n j n j n j n j n j n j n j f f f x t a f f x t c f f １次元線形移流拡散式_プログラム作成 (2) のエクセル参照

） （ ） （ （○） ） （ ） （ （◯）                                         01 . 0 , , 1 . 0 , 1 . 0 ) 6 ( 01 . 0 , 5 . 1 , 1 . 0 , 1 . 0 ) 5 ( 01 . 0 , 8 . 0 , 1 . 0 , 1 . 0 ) 4 ( 01 . 0 , 2 . 0 , 1 . 0 , 0 . 1 ) 3 ( 01 . 0 , 2 . 0 , 0 . 1 , 1 . 0 ) 2 ( 01 . 0 , 2 . 0 , 1 . 0 , 1 . 0 ) 1 ( a f c x t a c x t a c x t a c x t a c x t a c x t n j





_{ }



1



1 1 1 1 2 1 1 1 1 2 2 _         _{ }         n j n j n j n j n j n j n j f f f x t a f f x t c f f 宿題: (1)から(6)までの（ｘ、ｔ、振幅）の図を描け 2014年1月20日（月） 午後5時まで 3417室