第9回 微分方程式と差分法:Euler法(11月29日)



*出席確認

<出席確認メール> 添付ファイルなし 受付:本日14:30〜16:15

本日のキーワードは授業中に連絡します。


宛先:miwamoto[at]riko.shimane-u.ac.jp

件名:MMM2 20161129

本文: 学生番号 名前 「本日のキーワード」


*グループ分け


常微分方程式の数値解法(オイラー法)


次のような1階常微分方程式の初期値問題を数値的に解き、結果をグラフにしよう。

 

   


求めたいのは、 の範囲の関数 であり、関数 および定数   は与えられているとする。また、 の値は が与えられればいつでも計算できるものとする。


コンピューターでは数値を離散的に(とびとびに)しか扱うことができない。そこで、上記の微分方程式式をコンピューターで扱えるようにする必要がある。今回は、その方法として最も代表的なオイラー法を学ぶ。



差分商


   の関数としたとき、 の   における微分係数は


             ......(1)


で定義される。ここで、定義中の極限操作を取り払い、h を有限にとどめた

         ......(2)


を考えると、  を十分小さな正の実数にとれば、(2)式は(1)式の近似となっていると考えられる。この  のように、関数のいくつかの点における値の差を用いてその関数の微分係数を近似することを差分近似といい、(2)式の右辺のような量を差分商という。この場合は1階微分係数を近似する1階差分商である。


差分方程式


 では、このような差分商を用いて1階常微分方程式離散化してみよう。

      ......(3)

(3)の微分係数を上記の差分商で置き換えると、

       .......(4)

ただし、h は正の数とする。

(3)式と(4)式は異なる方程式なので、(3)式の解 は一般に(4)式を満たさない。

そこで、混乱を防ぐために(4)式の を   と書き換えておく。


      .......(5)

(5)式のように差分を含む方程式を差分方程式という。


次に、1階常微分方程式の  を   に置き換えた初期条件

      .......(6)

の下で(5)式を解くことを考える。



(5)式は、

     ......(7)


と変形できるので、  を代入すると、

 

       

となり、  から直ちに   が求まる。


注意: は初期条件として与えられているので既知


同様にして(7)式を繰り返し用いると、

 

 

 

............

というように、j=1,2,3,.... とすると Y((j+1)h) を Y(jh) から計算できる。

h をいったん決めると、 t = jh 以外の時刻の  Y の値は(7)式からは求めることができない。

このように、差分方程式を解くと従属変数はとびとびの時刻で値が定まる。そのようなとびとびの時刻を格子点と呼ぶ。Y は格子点でのみ意味があるので、そのことを明示するために  Y(jh) を Yj と書き換え、 Tj = jh とすると、(5)式と初期条件は、


 

 

となり、常微分方程式の初期値問題が、Yjに関する漸化式の問題に置き換えられる。この方法をオイラー法と呼ぶ。





オイラー法のアルゴリズム


必要な変数の宣言(T, h, a, Y, Y_new, tは実数型(float)、Nは整数型(int))


void setup()

{

ウィンドウサイズの設定

背景色の設定

T,N,aを設定

刻み幅の設定 h = T/(N-1)

初期値の設定 Y = a

}


void draw()

{

t = frameCount*h;

Y_new = Y + h*f(t, Y)

Y = Y_new


座標変換(下記参照)

色指定

描画


 if(frameCount*h >= T) noLoop();

}



 写像(座標変換)


写像(座標変換)

map(x, a,b,c,d)

Processingnにおいて、座標はピクセル数なので、実際に方程式を解いた求まる値をそのまま描画すると小さすぎたり、また縦方向は正の方向が逆向きになってしまう。そこで、map()関数を使うと、変数xの取りうる範囲が[a,b]のところを[c,d]の範囲へ変換してくれる。

float value = 100;

float m = map(value, 0, 100, -1, 1);

println(m);


*tmapやymapなど変換後の値を格納する変数を用意しておくとよい。



練習1 常微分方程式の初期値問題

 

 

をオイラー法を用いて 0≦ t ≦20 の範囲での解を求めるプログラムを作成し、

0≦ t ≦20 の範囲を2000等分割して、グラフ化せよ。


練習2 分割数を変えてみよ。


練習3 手計算にて解を求め、数値計算結果が正しいか確認せよ。

練習4 f(t,y)について独立した関数を用意せよ。



練習5

次の常微分方程式の初期値問題をオイラー法を用いて解け。いろいろな異なる時間刻み幅 h について計算し、結果について考察せよ。

 

 


練習6

次の常微分方程式の初期値問題をオイラー法を用いて解け。いろいろな異なる時間刻み幅 h について計算し、結果について考察せよ。

 

 


 


オイラー法による2階常微分方程式解法


2階常微分方程式をオイラー法を用いて解く方法を考える。

 

 

 


 これは単振動を記述する2階の常微分方程式である。2階の常微分方程式は次のように連立の1階の常微分方程式に帰着できる。

 ,    とおくと、


 

 

 

 


と書き換えることができる。

これにより、1階常微分方程式と同様、オイラー法を用いて解を求めることができる。




2階常微分方程式のオイラー法



<変数の宣言>

 y1, y2, y1_new, y2_newをfloatで用意する。


<setup関数で>

初期値の設定。y1 = a1, y2 = a2


<draw関数で>

    y1_new = y1 + h*f1(t, y1, y2)

    y2_new = y2 + h*f2(t, y1, y2)

    y1 = y1_new

    y2 = y2_new

 

練習7

次の常微分方程式の初期値問題をオイラー法をオイラー法で解き、グラフを描け。グラフは横軸が t であり、縦軸を y(t) とせよ。いろいろな異なる時間刻み幅 h について計算し、結果について考察せよ。



 

 






本日の課題1

下の図のように壁に一端を固定されたバネの他端におもりがつけられている状況を考える。

このおもりを少しだけ右に引っ張ってから手をはなす.おもりの中心の時刻 t におけるつり合いの位置からのずれを y(t) とする。また、時刻 0 でのおもりの位置を y(0)=1、おもりの速度を y'(0)=0 とすると、練習7と全く同じ、つぎの単振動の問題となる。

 

 

 

   

これではおもりに働く力はバネの復元力のみで、ある意味非現実的なので、今回は復元力以外におもりの速度 y' に比例する抵抗が働く場合を考えましょう(おもりと床との間の摩擦力のようなものをイメージ)。

改良されたモデルは次のようになる(2ky' の項が摩擦をあらわしている)。

 

 

 

ここで k,ω はそれぞれ抵抗力、復元力の強さを表す正の定数である。

この問題を連立1階常微分方程式に変換し、オイラー法を用いて解き、グラフにせよ。

k < ω でのグラフと k > ω でのグラフを見比べて違いに関して考察してみて下さい。


本日の課題2

下の図のように壁に一端を固定されたバネの他端におもりがつけられている状況を考える。

課題1の微分方程式を用いて、おもりが動くアニメーションを作成せよ。

このブラウザは埋め込まれたオーディオファイルを再生できません




このブラウザは埋め込まれたオーディオファイルを再生できません



このブラウザは埋め込まれたオーディオファイルを再生できません