第10回,第11回(2009年7月1日,8日)


今回の演習 


連結バネ系, 仕上げ(7/1)

前回紹介した連結バネ系について, 以下の事柄に留意し, 仕上げをしてください.

連結バネ系のレポートが完成した人は, 今日の演習後に一度回収しますので, レポートフォルダに rep0624.c というファイル名で置いてください. 出来なかった人も, 来週火曜 (7/7) の夕方4時までに提出してください.

連結バネ系は, 「オイラー法」でやってみてください. (ルンゲ・クッタ法で出来た人はそれで提出で構いませんが, 現時点で無理してルンゲ・クッタ法にこだわる必要はありません.)

オイラー法の場合, h = 0.0001 (T = 300 程度) とし, if 文等を使って, 1000 ステップに一回おもりの位置を表示するようにしてください. h は 0.0001 程度以下でないと, 計算精度が悪く, あっというまに発散してしまいます. またこの程度 h を小さくすると, 計算回数が多くなるため, 毎回おもりの位置を GLSC で表示させると非常に動きが遅くなります. 動きを滑らかにするために, 「 if( j%1000 == 0 ) ... 」などとして, 描写する, しない, の条件分岐をさせてください. あと, g_sleep も使って調整してください. この辺りはアニメーションの基本で, 結果の見やすさはものすごく重要です. (そのアニメーションが, (時に解析的には得られない) 現象や方程式系の解の情報をもれなく示しているのだから...) )

第9回


波動方程式の数値解法(7/1, 7/8)

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

なお,1次元波動方程式については今週と来週の2回にわたって演習を行う予定です.来週,1次元波動方程式に関するレポート課題を出す予定です.

さて,天下り的ですが1次元波動方程式は弦の振動をあらわします(導出等は講義で説明されます. また前回の「連結バネ系」からも導出できます. 時間があれば紹介します.).下記方程式 (W) における u(x,t) は弦の変位をあらわします.演習では1次元波動方程式を数値計算を用いて解き,その解 u(x,t) を画面上にアニメーションとして表示する事を目標とします.つまり,コンピュータ内部に弦の振動を再現し,それをアニメーションとして表示するわけです.


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

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

(W) 

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

を考えます。境界条件,初期条件が適切に与えられれば(W)の解が弦の振動を表現していると言えるでしょう.さて,数値計算する為にまず(W)を離散化する必要があります。分割方法は以下の2種類が考えられますが,

ここでは(講義の説明に従い,)時間方向の離散化は(1)を,空間方向の離散化は(2)を採用します.(これまでの演習で常微分方程式の数値解法を扱ったが,そこでは時間方向の離散化として(1)を採用してきた)

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

と書くことにすると、(W)の第一式は次のように離散化されます。

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

という式が得られます。2階微分の差分化については,このページ最後のスライドを参照してください(右辺最後の λ2(...) という項は, どこかで見覚えがあると思います.(連結バネ系と同じです.)).


初期条件

常備分方程式同様, 初期条件を与える必要があります. これは空間全域に対して与える必要があります. また上記の式はこの式は、時刻 tk-1, tkでの U から次の時刻 tk+1での U が計算できるという形をしています。したがって、初期条件として、 が必要となります。まず は(W)の第三式より、

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

となります。(W)より

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

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

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

と求まります.


境界条件

空間に関しても端を適度に扱う必要があり, そのために境界条件が与えられます. 境界条件は(W)の第二式ですが,まずは境界条件について触れる必要があります.詳しくは,講義で説明がありますが,ここでは,境界の処理の為に, という架空のセルを用意することにします.(以下講義のプリント)

様々な境界条件がありますが,今回の(W)における境界条件は Dirichlet 条件と呼ばれるものです.空間離散化の方法として(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) は課題1,2を参照 )

(境界条件)


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

を繰り返す ( g(x) は課題1,2を参照 )

(境界条件)


(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 の順に

を繰り返す


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

を繰り返す

(境界条件)


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

を繰り返す

(境界条件)


を繰り返す


課題1

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 を適度に選び, 描写頻度も適度に選ぶ事.)


課題2

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))

としたり

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))

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

レポート課題

課題1〜2のについてのプログラムを rep0708.c というファイル名で report フォルダにおいてください. (プログラムの先頭にコメント文で名前と学籍番号と課題の番号を書く事. またプログラムの説明もコメント文として適宜記入する事.)

自動回収日時:2009年7月21日(火)午後4時



課題3

他の境界条件,Neumann 境界条件,周期境界条件について復習し(講義を受けていない人は図書館なので調べる),それぞれの境界条件を適応したプログラムを作成し,実際に数値計算してみよ.


課題4

弦の振動ではなく,膜の振動を考えたい.すなわち,1次元波動方程式ではなく2次元波動方程式を考える.どのような方程式になるであろうか?


課題5

2次元波動方程式を差分化し、プログラムを作成し、数値計算せよ.


差分近似について(講義のスライドより).

導関数の差分近似について

関数 u(x) を x のまわりでTaylor展開すると

      ・・・(S1)

      ・・・(S2)

となり、両辺を加えると、

      ・・・(S3)

となります。O(d4)を無視すると、

      ・・・(S4)

この式の右辺の誤差は d2 程度となります。



期末レポート予告

期末レポートを提出してもらいます。提出〆切は7月末の予定です。(詳細は次回)

課題は以下の通り。

(1) 反応拡散系, チューリングパターンに関することを, A4 2枚程度にまとめよ.(自分の頭の中で何か考えてまとめる事. 必ず手書きの事.)

(2) チューリングパターンを生成する偏微分方程式系を数値的に解き, チューリングパターンを生成するプログラムを作成せよ

   (2) については、ヒントとなる方程式系を後の演習(13回目くらい)で扱いますので、それを解くプログラムを作成し、チューリングパターンが生成される条件を見つけてください.

   (1) は数学事務室に提出、(2) はいつも通りプログラムを自動回収します.

   成績は [これまでのレポートの点数]×[期末レポートの点数] でつけます. (つまり未提出では単位はありませんのでご注意ください.)


この後は, 第12回(7/15)に拡散方程式、13回(7/17 <- 金曜が水曜の授業に振替になってます. ご注意を.)以降で、反応拡散系、時空パターン形成の初歩に触れる予定です.

戻る