1 / 58

８．３ 常微分方程式の数値解法

８．３ 常微分方程式の数値解法. 数値解法では，解析的に積分できないような問題を扱うことが多い. オイラーってどんな人？. 数学者・物理学者であり、天体物理学者。 スイスのバーゼル生まれ 現ロシアのサンクトペテルブルクで死んだ。. 1707/4/15-1783/9/18 Leonhard Euler. ８．３．１　オイラー法 （１）前進差分によるオイラー法. 微分方程式. ［例］前進差分近似. 次のような繰返し計算で近似する. あるいは，. （本来は，. は限りなく０に近い値）. 差分近似が微分の近似であることから. Excel でオイラー法（１）.

taurus
Download Presentation

８．３ 常微分方程式の数値解法

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

  1. ８．３ 常微分方程式の数値解法 数値解法では，解析的に積分できないような問題を扱うことが多い オイラーってどんな人？ 数学者・物理学者であり、天体物理学者。 スイスのバーゼル生まれ 現ロシアのサンクトペテルブルクで死んだ。 1707/4/15-1783/9/18 Leonhard Euler

  2. ８．３．１　オイラー法（１）前進差分によるオイラー法８．３．１　オイラー法（１）前進差分によるオイラー法 微分方程式 ［例］前進差分近似 次のような繰返し計算で近似する

  3. あるいは， （本来は， は限りなく０に近い値） 差分近似が微分の近似であることから

  4. Excelでオイラー法（１） もちろん，この式は次のようにして厳密解を求めることができるが，誤差評価のためにこの例を求める。 （例） ［Excelの定義例］ tとして現時刻をとる定義

  5. Excelでオイラー法（２） グラフを描いてみると （例） tとして現時刻をとる定義のほうが誤差が少ない

  6. （２）修正オイラー法（修正法１） データ点における微分係数で増分を計算し， 増分から得られた座標点から更に増分を計算し， 両増分の平均値を最終的な増分とする。

  7. Excelによる修正法１ 式定義 結果グラフ

  8. （３）修正オイラー法（修正法２） データ点における微分係数で増分を計算し， 増分の１／２から得られた座標点から増分を計算し， この増分を最終的な増分とする。

  9. Excelによる修正法２ 式定義 結果グラフ

  10. VBAによる定義　①データ宣言，データ設定および結果設定用シートへのX値設定VBAによる定義　①データ宣言，データ設定および結果設定用シートへのX値設定 Private DX As Double Private Y0 As Double Private XST As Double Private Y(11) As Double Private N As Integer Sub データ設定() N = 11: Y0 = 1: DX = 0.1: XST = 0 End Sub Sub X値設定() With Worksheets("Sheet4") X = XST For i = 1 To N .Cells(i + 1, 1) = X: X = X + DX Next End With End Sub

  11. VBAによる定義　②結果設定とボタンのClickイベントハンドラVBAによる定義　②結果設定とボタンのClickイベントハンドラ Sub 結果設定(ID) With Worksheets("Sheet4") For i = 1 To N .Cells(i + 1, ID) = Y(i) Next End With End Sub Sub ボタン1_Click() データ設定 X値設定 オイラー法前進差分 Y, Y0, XST, DX, N 結果設定 2 修正オイラー法1 Y, Y0, XST, DX, N 結果設定 3 修正オイラー法2 Y, Y0, XST, DX, N 結果設定 4 厳密解 Y, Y0, XST, DX, N 結果設定 5 End Sub

  12. VBAによる定義　③関数値，オイラー法，修正法１VBAによる定義　③関数値，オイラー法，修正法１ Function g(X, Y) As Double g = (X + 1) * Y End Function Sub オイラー法前進差分(Y, Y0, XST, DX, N) Y(1) = Y0: X = XST For i = 2 To N X = X + DX Y(i) = Y(i - 1) + g(X, Y(i - 1)) * DX Next End Sub Sub 修正オイラー法1(Y, Y0, XST, DX, N) Y(1) = Y0: X = XST For i = 2 To N K1 = DX * g(X, Y(i - 1)) K2 = DX * g(X + DX, Y(i - 1) + K1) X = X + DX Y(i) = Y(i - 1) + (K1 + K2) / 2 Next End Sub

  13. VBAによる定義　④修正法２，厳密解設定 Sub 修正オイラー法2(Y, Y0, XST, DX, N) Y(1) = Y0: X = XST For i = 2 To N K1 = DX * g(X, Y(i - 1)) K2 = DX * g(X + DX / 2, Y(i - 1) + K1 / 2) Y(i) = Y(i - 1) + K2 X = X + DX Next End Sub Sub 厳密解(Y, Y0, XST, DX, N) X = XST For i = 1 To N Y(i) = Exp(X * X / 2 + X) X = X + DX Next End Sub

  14. VBAによる定義による結果 結果（Excelによる結果より誤差が少ない）

  15. （４）修正オイラー法（一般化） 以下の係数に適当な値を代入する。修正法１，２はその代表的な例。 ［注］以下，式の展開による導出を示すが， 式の展開自体を覚える必要はない。 　　　どのような流れで導出されているかを理解し， 　　　最終的な結果を確認すること。 　　　式の展開自体が必要なときは，資料を参照すれば良い。

  16. 確認 をテーラ展開すると 第３式に代入すると テーラ展開の式と比較すると 修正法２ 修正法１

  17. ８．３．２　ルンゲ・クッタ法（Runge-Kutta method）（１）２次のルンゲ・クッタ法(その１） とおくと

  18. （１）２次のルンゲ・クッタ法（その２） であるからテーラ展開すると 定義から 係数を比較して これらを満足する値として を選択して この式はホイン（Heun）の２次公式とも呼ばれる。

  19. （２）３次のルンゲ・クッタ法 において とおくと，２次の場合と同様，係数比較によって を決めることができる。 この式はクッタ（Kutta）の公式と呼ばれる。

  20. （３）４次のルンゲ・クッタ法 において とおくと，２次の場合と同様，係数比較によって を決めることができる。 この式はルンゲ・クッタの式と呼ばれる。

  21. （４）ルンゲ・クッタ・ギルの公式（Runge-Kutta-Gill method） ４次のルンゲ・クッタ型の公式の一つとして，次の公式がある。

  22. ルンゲとクッタ • カール・ルンゲ(Runge,Carl David Tolme) ドイツの数学者で物理学者。 1856/08/30-1927/01/03 • ウイルヘルム・クッタ(Wilhelm Kutta)ドイツの数学者でかつ物理学者。1867~1944。渦による揚力の発生クッタ-ジェコフスキーの揚力理論で有名。

  23. VBAによる定義　①データ宣言，データ設定および結果設定用シートへのX値設定VBAによる定義　①データ宣言，データ設定および結果設定用シートへのX値設定 Private DX As Double Private Y0 As Double Private XST As Double Private Y(11) As Double Private N As Integer Sub データ設定() N = 11: Y0 = 1: DX = 0.1: XST = 0 End Sub Sub X値設定() With Worksheets("Sheet4") X = XST For i = 1 To N .Cells(i + 1, 1) = X: X = X + DX Next End With End Sub

  24. VBAによる定義　②結果設定とボタンのClickイベントハンドラVBAによる定義　②結果設定とボタンのClickイベントハンドラ Sub 結果設定(ID) With Worksheets("Sheet4") For i = 1 To N .Cells(i + 1, ID) = Y(i) Next End With End Sub Sub ボタン1_Click() データ設定 X値設定 ルンゲクッタ Y, Y0, XST, DX, N 結果設定 2 ルンゲクッタギル Y, Y0, XST, DX, N 結果設定 3 厳密解 Y, Y0, XST, DX, N 結果設定 4 End Sub

  25. VBAによる定義　③関数値，４次のルンゲ・クッタ法VBAによる定義　③関数値，４次のルンゲ・クッタ法 Function g(X, Y) As Double g = (X + 1) * Y End Function Sub ルンゲクッタ(Y, Y0, XST, DX, N) Y(1) = Y0: X = XST For i = 2 To N k1 = g(X, Y(i - 1)) k2 = g(X + DX / 2, Y(i - 1) + DX * k1 / 2) k3 = g(X + DX / 2, Y(i - 1) + DX * k2 / 2) k4 = g(X + DX, Y(i - 1) + DX * k3) Y(i) = Y(i - 1) + (k1 + 2 * k2 + 2 * k3 + k4) * DX / 6 X = X + DX Next End Sub

  26. VBAによる定義　④ルンゲ・クッタ・ギル法 Sub ルンゲクッタギル(Y, Y0, XST, DX, N) Y(1) = Y0: X = XST A1 = -(1 / 2 - 1 / Sqr(2)) A2 = (1 - 1 / Sqr(2)) B1 = -1 / Sqr(2) B2 = 1 + 1 / Sqr(2) C1 = 2 * (1 - 1 / Sqr(2)) C2 = 2 * (1 + 1 / Sqr(2)) For i = 2 To N k1 = g(X, Y(i - 1)) k2 = g(X + DX / 2, Y(i - 1) + DX * k1 / 2) k3 = g(X + DX / 2, Y(i - 1) + A1 * DX * k1 + A2 * DX * k2) k4 = g(X + DX, Y(i - 1) + B1 * DX * k2 + B2 * DX * k3) Y(i) = Y(i - 1) + (k1 + C1 * k2 + C2 * k3 + k4) * DX / 6 X = X + DX Next End Sub

  27. VBAによる定義　⑤厳密解の設定 Sub 厳密解(Y, Y0, XST, DX, N) X = XST For i = 1 To N Y(i) = Exp(X * X / 2 + X) X = X + DX Next End Sub

  28. VBAによる定義　④修正法２，厳密解設定 Sub 修正オイラー法2(Y, Y0, XST, DX, N) Y(1) = Y0: X = XST For i = 2 To N K1 = DX * g(X, Y(i - 1)) K2 = DX * g(X + DX / 2, Y(i - 1) + K1 / 2) Y(i) = Y(i - 1) + K2 X = X + DX Next End Sub Sub 厳密解(Y, Y0, XST, DX, N) X = XST For i = 1 To N Y(i) = Exp(X * X / 2 + X) X = X + DX Next End Sub

  29. ルンゲクッタ法による結果 非常に良い結果が得られていることに注意

  30. ８．３．３　連立微分方程式（１）手順 を求める。 ②増分 ～ 次の連立１次微分方程式を考える。 を求める。 ③ を求める。 ①増分ベクトル

  31. （２）プログラム例 次の連立１次微分方程式を解くプログラムを例とする。

  32. ①データ宣言，データ設定および結果設定用シートへのX値設定①データ宣言，データ設定および結果設定用シートへのX値設定 Private DX As Double Private Y0() As Double Private XST As Double Private Y(11, 2) As Double Private N As Integer Private M As Integer Sub データ設定() N = 11: DX = 0.1: XST = 0 M = 2: ReDim Y0(M) For k = 1 To M Y0(k) = 0 Next End Sub Sub X値設定() With Worksheets("Sheet2") X = XST For i = 1 To N .Cells(i + 1, 1) = X: X = X + DX Next End With End Sub

  33. ②結果設定およびボタンのClickイベントハンドラ②結果設定およびボタンのClickイベントハンドラ Sub 結果設定() With Worksheets("Sheet2") For i = 1 To N For k = 1 To M .Cells(i + 1, k + 1) = Y(i, k) Next Next End With End Sub Sub Sheet2_ボタン1_Click() データ設定 X値設定 連立ルンゲクッタ Y, Y0, XST, DX, N, M 結果設定 End Sub

  34. ③関数設定および処理本体（その１） Function g(ID, X, Y) As Double If ID = 1 Then g = X + Y(2) Else g = X End If End Function Sub 連立ルンゲクッタ(Y, Y0, XST, DX, N, M) Dim K1() As Double: ReDim K1(M) Dim K2() As Double: ReDim K2(M) Dim K3() As Double: ReDim K3(M) Dim K4() As Double: ReDim K4(M) Dim YY() As Double: ReDim YY(M) For k = 1 To M Y(1, k) = Y0(k) Next

  35. ⑤処理本体（その３） X = XST For i = 2 To N For k = 1 To M YY(k) = Y(i - 1, k) Next For k = 1 To M K1(k) = g(k, X, YY) Next For k = 1 To M YY(k) = Y(i - 1, k) + DX * K1(k) / 2 Next For k = 1 To M K2(k) = g(k, X + DX / 2, YY) Next For k = 1 To M YY(k) = Y(i - 1, k) + DX * K2(k) / 2 Next For k = 1 To M K3(k) = g(k, X + DX / 2, YY) Next

  36. ⑤処理本体（その４） For k = 1 To M YY(k) = Y(i - 1, k) + DX * K3(k) Next For k = 1 To M K4(k) = g(k, X + DX, YY) Next For k = 1 To M Y(i, k) = Y(i - 1, k) + _ (K1(k) + 2 * K2(k) + 2 * K3(k) + K4(k)) * DX / 6 Next X = X + DX Next End Sub

  37. ８．３．４　高階の常微分方程式（１）考え方８．３．４　高階の常微分方程式（１）考え方 次の微分方程式を考える。 各次数の時間微分を別々の変数と考えると という連立１次微分方程式の形式であるため， ルンゲ・クッタ法を適用することができる。

  38. ［例題］ 解析的に解けるが， 厳密解が分かっているので， 誤差検証のために あえてこの例題を取り上げる。

  39. （連立ルンゲ・クッタ法の変更部分） ・ 　・(中略） 　・ Function g(ID, X, Y) As Double If ID = 1 Then g = Y(2) Else g = -Y(1) End If End Function 　・ 　・(後略） 　・ （下記，赤線部分の変更） Private DX As Double Private Y0() As Double Private XST As Double Private Y(60, 2) As Double Private N As Integer Private M As Integer Sub データ設定() N = 60: DX = 0.1: XST = 0 M = 2: ReDim Y0(M) Y0(1) = 0 Y0(2) = 4 End Sub

  40. 実行結果 誤差が少ないことを確認する。

  41. （２）ルンゲ・クッタ法を使わない高次微分方程式（オイラー法による例題）（２）ルンゲ・クッタ法を使わない高次微分方程式（オイラー法による例題） ここで，s = nhとなる有限の刻み幅hを導入し， 次のような数式モデルを設定 とおくと とすることができるので， と変形しておく。

  42. 処理の手順 ① DT = 刻み幅 ② V0 = 初速度，X0 = 0，T = 0とする． ③ T時刻の速度，位置をV0，X0とする． ④ T = T + DTとする． ⑤ T ≦ 最終時刻の間，以下を繰り返す． ・X = X0 ＋ V0 × DT ・V = V0 － X0 × DT ・T時刻の速度，位置をV，Xとする． ・X0 = X，V0=V ・T = T + DT

  43. 単純な方法による式定義（厳密解との比較を含む）単純な方法による式定義（厳密解との比較を含む）

  44. 単純な方法による結果（誤差が蓄積している）単純な方法による結果（誤差が蓄積している）

  45. この理由（エネルギー保存則） ＝一定， したがって ＝一定 ＝一定 一方，数式モデルから 刻み幅が0であれば一定となるが，今刻み幅は有限の値にしているので 保存則が成立していない。

  46. 式の改良 現在時刻における速度と位置の値は，お互いに現在時刻における値と前時刻との平均値を使用して計算するものとする。 現在時刻における値がお互いに入っているので，計算上はループしてしまうので，お互いに代入して整理すると， 両式を二乗して加えると

  47. 確認（エネルギー保存則が満足されている）

  48. 考慮した結果

  49. （３）かえる飛び法 次のようにお互いに相手の中間値を使用する方法 論理的には，中心差分であることに注意しよう。

  50. かえるとび法（元の式は変わらないが誤差が少ない）かえるとび法（元の式は変わらないが誤差が少ない）

More Related