第15回(2011年7月27日)


今日の演習



今回は各自実習をしてください.

期末テスト解答例

以下では、期末テスト問題の簡単な解説と解答例の紹介をします。どちらも基本的には 第12回、第13回 で扱った課題をベースにした問題でした。



第13回 課題6で扱った問題にやや変更が加わった問題です。変更点は

  方程式で、先の課題で -a u(x,t) と成っていた所が u(x,t) - u(x,t)3 となっている。
  境界条件、初期条件がやや異なる

ですので基本的には課題6の上にある(1変数反応拡散系に関する)解説通りに差分化して、そのアルゴリズムに従ってプログラムを作成(もしくは課題6のプログラムを変更)すれば、正解になります。境界条件については 第12回 課題4等で扱いました。

問題1 解答例






2つの変数 u(x,t), v(x,t) について計算して、同時に表示する問題なので、やや難しくなっていますが、差分化とアルゴリズムは基本的には、第13回 課題8とその上にある(2変数反応拡散系に関する)解説通りに進めていけば出来ます。あとは境界条件が左右で異なるところが注意点です。境界を振動的に時間変化させるのは、第12回 課題5や課題7 等で扱いました。

問題2 解答例




ところで、この期末テストで扱った偏微分方程式は「ギンツブルグ・ランダウ方程式」と呼ばれる、有名な方程式です。各自これがどのような方程式なのか(何を表し、どのような性質を持つのか)調べてみてください。


期末レポート (7/13 出題)

期末レポートを提出してもらいます。提出〆切は8/5の予定です。

課題は以下の通り。

(1) 反応拡散系、(第13回 課題8で扱っている)チューリングパターン (またはアラン・チューリング本人) について調べ, A4 用紙 2〜4枚程度に、手書きでまとめよ. (数学事務室のレポートボックスに提出すること)
 (自分で資料を探し、頭の中で何か考えてまとめる事. 参考文献を記す事. 必ず手書きの事.)

(2) 以下で紹介する波動方程式をに関し、以下の課題(9, 10, 11, 12) の何れかを解け。作成したプログラムはファイル名をrep0727.c として 下の提出用ページより提出する事. プログラムの先頭にコメント文で名前、学籍番号と選んだ課題番号. またプログラムの説明もコメント文として適宜記入する事.

提出用ページ

〆切はどちらも8/5中とします。




波動方程式の数値解法

(7/13 に紹介した内容です)

今回は1次元波動方程式(wave equation)を扱います。ここでは,1次元波動方程式をコンピュータを用いて解くために差分化を行い,コンピュータを用いて波動方程式を解くアルゴリズムを示します.今日の内容に関する詳しい情報は講義で説明があります(ました).演習では細かい事は気にせずに一気にアルゴリズムの紹介まで行きたいと思います.

さて,天下り的ですが1次元波動方程式は、下図のような弦や(深さのある水の)水面の時空間的な振動を表現する方程式系です(導出等は講義で説明されます. ).下記方程式 (W) における u(x,t) は弦(もしくは水面)の変位をあらわします.演習では1次元波動方程式を数値計算を用いて解き,その解 u(x,t) を画面上にアニメーションとして表示する事を目標とします.つまり,コンピュータ内部に弦の振動を再現し,それをアニメーションとして表示するわけです.




(W) 

時間ー空間2階偏微分方程式の差分化

さて,次の1次元波動方程式の初期値境界値問題、

(W) 

但し、c >0, f(0)=0, f(L )=0, g(0)=0, g(L )=0.

を考えます。境界条件,初期条件が適切に与えられれば(W)の解が弦の振動を表現していると言えるでしょう.これまで扱ってきた拡散方程式との違いは、
i) u(x,t) の時間 ( t ) に関する2階微分が含まれている事,
ii) それに伴って, 初期条件に u(x,t) の時間の1階微分に関する条件がつく事,
です

さて,数値計算する為にまず(W)を離散化する必要があります。拡散方程式のときにも紹介しましたが、分割方法は以下の2種類が考えられますが,

ここでは時間方向の離散化は「京都方式」を,空間方向の離散化は「札幌方式」を採用します.(拡散方程式と同様です.)

[0,L ] 区間を N 等分するとし,上記離散化方法に従い, (δx = ) h = L /N, τ >0 とし,それを用いて xj = ( j - 1/2)h, tk = として,

と書くことにします。(W)の第一式は、左辺・右辺それぞれ u(x,t) の t に関する2階微分、x に関する2階微分を含んでいます。 これらは拡散方程式における, 空間に関する2階微分の差分化(以下のスライドを参照)と同様、次のように離散化されます。

よって λ=/h として、まとめると

という式が得られます。


境界条件

拡散方程式と同様に、空間に関しても端を適度に扱う必要があり, そのために境界条件が与えられます. 境界の処理の仕方も拡散方程式と同様です,

上記の Dirichlet 条件は、両端を u(x,t) = 0 に固定した場合の物と成っています、拡散方程式の時と同様に
  すなわち 

  すなわち 

で与えると、両端の u の値を好きな値 A , B で固定した境界条件を与えられます。


初期条件

拡散方程式同様, 初期条件を空間全域に対して与える必要があります. ここで拡散方程式と波動方程式の大きな違いが現れます。今回の波動方程式を差分化した上記の式では、 時刻 tk-1 , tk での U ( つまり Uik-1 , Uik ) から次の時刻 tk+1 での U ( つまり Uik+1 ) が計算できるという形をしています。したがって計算を開始するためには、初期値として、 が必要となります。

まず は(W)の第三式より、

となります。 を決める方法はいくつか考えられますが、例えば次のようにします。まず、テイラーの公式を用いると、

となります。(W)より

であり,さらに t=0 でも(W)の第一式が成り立つとすると、

となります。この式の右辺は2階差分商を用いて

と近似できます。つまり,

と求まります.以下ではこのように求まった を用います。


解くべき(差分)方程式

以上の考察から解くべき式をまとめると、以下のようになります.

(2)

但し、 ( λ=/h ) が安定性のための必要十分条件(安定性条件とは,数値計算がうまく行く為の条件です.なぜにこういった条件が必要であるかは,講義で説明があります).

めでたく差分方程式が得られました。(2)をコンピューターを用いて解くのが今週の目標です。


アルゴリズム

1次元波動方程式を解く為のアルゴリズムは次のようになります。安定性条件を満たすため λ ≦ 1 となるように h, τ, を与えます.

アルゴリズム中の 「 xj 」, 「 uj 」, 「 old_uj 」, 「 new_uj 」 はそれぞれ実数型の配列として定義すると良い.境界条件処理用に配列要素はN個よりも余分に2つ必要となり,計 N+2 個必要となる事に注意(new_uj については,実際には境界処理用の要素は用いない).例えば uj について,C 言語風に配列を定義すると,

float u[N+2];

となる.このようにすれば,u[0], u[1], ... , u[N+1] が用意される.実際そのように用意すれば良い.

アルゴリズム

(1)初期パラメータ設定

T(時間上限), M(時間刻み数), L(空間幅), N(空間刻み数), c を設定する

h = L/N, τ=T/M, λ = τ*c/h (λ ≦ 1 となるように. M, N を調整, プログラム中では, λは「ramda」, τは「tau」などとするとよい.)

(1.5) GLSC の準備

(2) 初期値設定

j := 1,....,N の順に

xj = j*h - 0.5*h

を繰り返す ( f(x) は課題9以降を参照 )

(境界条件)


j := 1,....,N の順に

を繰り返す ( g(x) は課題9以降を参照 )

(境界条件)


(3)計算, 動画描画

k := 0,1,.....,M-1 の順に

(x[1], u[1]) 〜 (x[N], u[N]) を画面に表示 (g_move + g_plot, x[0], u[0], x[N+1], u[N+1]

は境界処理用なので表示しない)

j := 1,....,N の順に

を繰り返す


( uj , old_uj 両方の値を更新します。以下の順序を守らないと計算がうまく行きません。(なぜだか分かりますか?))

j := 1,....,N の順に

を繰り返す

(境界条件)


j := 1,....,N の順に

を繰り返す

(境界条件)


を繰り返す


課題9

L=1,  c=1,  f(x)=2*x*(1-x),  g(x)=0,  u(0,t)=u(1,t)=0  として (2) を解くプログラムを作成してください。(N は100程度, T ~ 100 で λ << 1 となるように M を適度に選び, 描写頻度も適度に選ぶ事.)

実行サンプル (mpg)


課題10

L=1, c=1, u(0,t)=u(1,t)=0 (N = 100 程度) はそのままで、f(x) および g(x) を他のもにかえてみよ。例えば関数 g(x)=0、関数 f(x) を

f(x)=exp(-100*(x-0.4)*(x-0.4)) ... (i)

としたり

f(x)=exp(-100*(x-0.4)*(x-0.4)) - exp(-50*(x-0.7)*(x-0.7)) + exp(-10*(x-0.5)*(x-0.5)) ... (ii)

とする等。
(局在化した初期値からはじめると、波の伝播が観察される。)

実行サンプル (i)  実行サンプル (ii) (mpg)


課題11

他の境界条件,(i) Neumann 境界条件もしくは (ii) 周期境界条件、を適応したプログラムを作成し,実際に数値計算してみよ.

実行サンプル (i) 実行サンプル (ii) (mpg)


課題12

境界条件として (端に適当な値を与える) ディレクレ条件で, A = sin(wt), B = 0 とした場合について、様々な c, w について数値的に解き ( 初期条件 : f(x)=0 , g(x)=0 )、、c, w を適度に選んで、(バネ振動のときに見たような)「共鳴(共振)」を観察せよ. (以前の「バネ振動子」に比べ、より「吊り橋の共振」に近い状況がシミュレーションできます.)

実行サンプル(非共鳴的) 実行サンプル(共鳴的) (mpg)


おまけ 1

境界条件として Neumann 条件を用いて,初期条件として
g(x)=0、
f(x)
  = 1 ( 0.4 ≦ x ≦ 0.6 ) ,
  = 0 ( Otherwise )
等として計算しその様子を観察せよ.

実行サンプル (N = 5000 としています. ) (mpg)


先日大きな被害を出した津波は、地震によって生じた海底の隆起(沈降)によって生じる急激な海面の変形(今回の初期条件のような形状)が波として伝搬したものです。その様子は波動方程式で観察されます。



浅い所では波と海底が擦れて海底付近の海水の流れが遅くなるため、波の上部が多いかぶさるようになります.このような現象は、今回扱っている波動方程式の適用外になっています。




おまけ 2(フラクタル図形を描く)

(拡散方程式が習得できて余裕のある方、トライしてみてください.)

(1) Ukj が以下の差分方程式に従うとする.

   Ujk+1 = F( Ukj-1 ,  Ukj ,  Ukj+1 )
ただし k = 0, 1, 2, ...M (正整数), j = 1, 2, ..., N (正整数) 


ここで F( Ukj-1 ,  Ukj ,  Ukj+1 ) = Ukj + α( Ukj-1 -2Ukj + Ukj+1 )
と与えられれば、先に扱った拡散方程式の差分化したものになります.


しかし一般的に F( Ukj-1 ,  Ukj ,  Ukj+1 ) は、もっといろいろな形のものが考えられて良いはずです。
そこで今回、以下のような差分方程式系を考えてみます。

 Ukj は 0 もしくは 1 の2値のみをとり、次の差分式に従う.

   Ujk+1 = F( Ukj-1 ,  Ukj ,  Ukj+1 )

   ただし
       F(0, 0, 0) = 0, F(0, 0, 1) = 1, F(0, 1, 0) = 0, F(0, 1, 1) = 1,
       F(1, 0, 0) = 1, F(1, 0, 1) = 0, F(1, 1, 0) = 1, F(1, 1, 1) = 0
 もしくは
     F( Ukj-1 ,  Ukj ,  Ukj+1 ) = 1  ( Ukj-1 + Ukj+1 = 1 の時 )
     F( Ukj-1 ,  Ukj ,  Ukj+1 ) = 0  ( Ukj-1 + Ukj+1 = 1 以外の時 )
(どちらも同じです.)

境界条件: Uk0 = 0, UkN+1 = 0

初期条件: U0j = 1 ( j = N/2 ), U0j = 0 ( j = N/2 以外 )

とする. その解の様子を図示せよ.

表示方法は, 縦軸 k, 横軸 j, とし, Ukj = 1 のとき, g_circle, g_marker 等で ( j, k ) に丸(点)を描く等とすること.
( M = 200, N =400 程度でこんな絵が描けます.(点は半径 0.5 の円.)

 

(2)「フラクタル」について調べよ

(3) Cellular Automata (セルオートマトン) について調べよ




おまけ 3(交通渋滞のモデル)

(1) 一車線道路上の車の動きがシンプルなルールに従うとして、交通渋滞の様子をモデル化してみます。簡単のため、空間は j = 1, 2, ... N の格子点に区切られていて、各格子点上には車が1台だけ居られるとします。時刻も k = 0, 1, 2, ... と離散的な量として扱います。車は前方( 格子点番号 j が大きくなる方向 )に進むとします.

ある時刻 m で n 番目の格子に車が居る場合、その車の動きは以下の規則に従うとします。

 n+1 番目の格子(前方の格子点)に別の車が居なければ、次の時刻(m+1)に n+1 番目の格子に移動する。 (前が空いていたら進む。)

 n+1 番目の格子(前方の格子点)に別の車が居る場合、次の時刻(m+1)も n 番目の格子にとどまる。 (前が空い無ければ進まない。)

 


ここで時刻 k, 場所 j に車が存在したら Ukj = 1, いなければ Ukj = 0 とすると、Ukj の変化は以下のような差分方程式系で記述できることが分かります。

   Ujk+1 = F( Ukj-1 ,  Ukj ,  Ukj+1 )

       F(0, 0, 0) = 0, F(0, 0, 1) = 0, F(0, 1, 0) = 0, F(0, 1, 1) = 1,
       F(1, 0, 0) = 1, F(1, 0, 1) = 1, F(1, 1, 0) = 0, F(1, 1, 1) = 1
 

この方程式系について、境界条件を周期境界条件: Uk0 = UkN, UkN+1 = Uk1 とし、初期条件を

U0j = 1 ( sin( 0.5*j ) + P > 0 の場合 )

U0j = 0 ( それ以外 )

として. その様子を図示してみてください. (P : -1 〜 1 で大きいほど車の台数が増えます。)

表示方法は, 縦軸 k, 横軸 j, とし, Ukj = 1 のとき, g_circle, g_marker 等で ( j, k ) に丸(点)を描く等とすること.
( M = 50, N =50 程度でこんな絵が描けます.(点は半径 0.25 の円と長さ 0.5 の線. 時間の方向は下→上、車は左→右に進む。)

 

 P = -0.95(車少ない)各車はスムーズに前方に進む

 

 P = -0.5(車少し多い)各車は(まだ)スムーズに前方に進む (1格子点分の車間距離は空いているので、止まらない。)

 

 P = 0.5(車多い)渋滞ができて進めない車がいる。渋滞自身は後方に移動していく(車は前方にすすむ。)。

このモデルは、非常に単純化されていますが現実の交通流のもつ特徴を良く捉えています。この方法をもっと複雑な道路網やネットワークに拡張して、都市交通やインターネットのパケットの渋滞等の解析も行われています。



戻る