８．３ 常微分方程式の数値解法 数値解法では，解析的に積分できないような問題を扱うことが多い オイラーってどんな人？

Presentation on theme: "８．３ 常微分方程式の数値解法 数値解法では，解析的に積分できないような問題を扱うことが多い オイラーってどんな人？"— Presentation transcript:

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

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

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

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

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

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

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

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

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

10 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

11 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

12 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) K1 = DX * g(X, Y(i - 1)) K2 = DX * g(X + DX, Y(i - 1) + K1) Y(i) = Y(i - 1) + (K1 + K2) / 2

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)

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

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

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

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

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

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

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

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

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

23 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

24 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

25 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 / / Sqr(2)) A2 = (1 - 1 / Sqr(2)) B1 = -1 / Sqr(2) B2 = / 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)

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

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

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

32 ①データ宣言，データ設定および結果設定用シートへの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 End With

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

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 K1(k) = g(k, X, YY) YY(k) = Y(i - 1, k) + DX * K1(k) / 2 K2(k) = g(k, X + DX / 2, YY) YY(k) = Y(i - 1, k) + DX * K2(k) / 2 K3(k) = g(k, X + DX / 2, YY)

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

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

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

39 （連立ルンゲ・クッタ法の変更部分） ・ （下記，赤線部分の変更） ・(中略） Private DX As Double ・
　・(中略） 　・ 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 かえるとび法（元の式は変わらないが誤差が少ない）

51 （４）テーラ展開による方法 まず，微分方程式から高次の微分を求め， テーラ展開の式に代入する。 ［例］ とすると
これを次のようなテーラ展開の式に代入する。

52 式の変形 次のように近似できる。 更に次のように変形できる。

53 式の整理 すなわち， 。同様にして

54 単純化 一方，定義から 最終的な数式モデル

55 テーラ展開による方法

56 テーラ展開による方法（結果）

57 ［例題２］

58 ［結果］