第13章 微分方程式その2 13.1 予測子・修正子法 オイラー法・ルンゲクッタ法 既知の𝑥𝑗 、𝑥𝑗+1 と𝑦𝑗 から 𝑦𝑗+1 を直接求める方法 予測子・修正子法 2段階で計算 予測子（陽的公式）で近似値の予測 修正子（陰的公式）で修正 𝑥𝑗−𝑘 , 𝑥𝑗+1 で積分 一般的に予測子・修正子法では積分で求める (13.1) 左辺を積分すると (13.2) 式(13.1)の左辺、すなわち式(13.2)の右辺はyが含まれているので積分できな い。式(13.1)の右辺の積分は台形公式を使えば求められる。 ●k=0のとき、 𝑥𝑗 , 𝑥𝑗+1 で積分 式(13.2)と式(13.3)より、 k=0のとき、式(13.1)は、 𝑦𝑗+1 を例えばオイラー法で予測する 予測値を𝑦 ∗𝑗+1 予測子（陽的公式） 既知𝑥𝑗 と𝑥𝑗+1 から未知𝑦 ∗𝑗+1 を求める 修正子（陰的公式） 予測した𝑦 ∗𝑗+1 から𝑦𝑗+1 を求める ●ホイン法に似ている。二段階で求める 𝑦𝑗+1 をオイラー法で予測してその予測値を𝑦 ∗𝑗+1 xj+1の傾きを 𝑓 𝑥𝑗+1 , 𝑦𝑗 + ℎ𝑓(𝑥𝑗 , 𝑦𝑗 ) ≅ 𝑦𝑗+1 y 𝑦(𝑥) 𝑦𝑗+1 yj とおく。これは既知 xjの傾きだと この分が誤差 xjの傾き 𝑓 𝑥𝑗 , 𝑦𝑗 この高さ 𝑦𝑗 + ℎ𝑓(𝑥𝑗 , 𝑦𝑗 ) をyj+1と仮定している。 ℎ 𝑥𝑗 𝑥𝑗+1 x 式 13.5 の第1式を第2式に代入して、𝑦 ∗𝑗+1 を消去して、直接求めると、 →ホイン法と一致 既知の𝑥𝑗 、𝑥𝑗+1 と𝑦𝑗 から 𝑦𝑗+1 を直接求めている!! ●ミルン法 既知の(xj-3,yj-3)から(xj-2,yj-2),(xj-1,yj-1)、および、(xj,yj)をルンゲクッタ法で 求め、予測子として高精度の公式(13.7)を用いy*j+1を予測。 ●修正子としてシンプソンの公式を用いる k=1の場合、式(13.2)の右辺は、3点𝑥𝑗−1 、 𝑥𝑗 、 𝑥𝑗+1 を用いて近似 高精度の予測子 二階微分まで考慮 詳細は専門書参照 シンプソンの公式による修正子 (13.6) (13.7) 13.2 アダムス・バッシュフォース法 前述の式(13.1) (13.1) 積分区間 𝑥𝑗 , 𝑥𝑗+1 において、 (13.8) f(x,y)を2点 𝑥𝑗−1 , 𝑓𝑗−1 と 𝑥𝑗 , 𝑓𝑗 を通る直線P(x)で近似 j-1とjのみを使ってこの点 𝑓𝑗 − 𝑓𝑗−1 𝑃 𝑥 = 𝑓𝑗−1 + 𝑥 − 𝑥𝑗−1 を外挿!! ℎ 𝑓 ∗𝑗+1 傾きf j+ 1 y 𝑦𝑗+1 P(x) 𝑦(𝑥) yj+1－yj 𝑓𝑗+1 𝑓𝑗 𝑦𝑗 傾きfj f(x,y) 𝑦𝑗−1 𝑥 𝑓 𝑥, 𝑦(𝑥) 𝑑𝑥 傾きfj-1 𝑥 𝑓𝑗−1 𝑗+1 𝑗 𝑥𝑗−1 ℎ 𝑥𝑗 ℎ 𝑥𝑗+1 x 𝑥𝑗−1 𝑥𝑗 𝑥𝑗+1 P(x)を式(13.8)に代入して、積分すると、 𝑦𝑗+1 ℎ = 𝑦𝑗 + 3𝑓(𝑥𝑗 , 𝑦𝑗 ) − 𝑓(𝑥𝑗−1 , 𝑦𝑗−1 ) 2 → 2次精度 (13.9) 𝑥𝑗−1 , 𝑦𝑗−1 , 𝑥𝑗 , 𝑦𝑗 のみを使って𝑦𝑗+1 を計算 𝑥𝑗+1 , 𝑓𝑗+1 を使っていない!! 同様に3点 𝑥𝑗−2 , 𝑓𝑗−2 、 𝑥𝑗−1 , 𝑓𝑗−1 、 𝑥𝑗 , 𝑓𝑗 を通る放物線で置き換え て、積分すると、 𝑦𝑗+1 = 𝑦𝑗 ℎ + 23𝑓(𝑥𝑗 , 𝑦𝑗 ) − 16𝑓 𝑥𝑗−1 , 𝑦𝑗−1 + 5𝑓(𝑥𝑗−2 , 𝑦𝑗−2 ) 12 → 3次精度 (13.10) 𝑥𝑗−2 , 𝑦𝑗−2 , 𝑥𝑗−1 , 𝑦𝑗−1 , 𝑥𝑗 , 𝑦𝑗 のみを使って𝑦𝑗+1 を計算 6 𝑥𝑗+1 , 𝑓𝑗+1 を使っていない!! 𝑥𝑗+1 , 𝑓𝑗+1 を使う方法 ●アダムス・ムルトン法 式(13.8) の𝑓 𝑥, 𝑦 を3点 𝑥𝑗−1 , 𝑓𝑗−1 、 𝑥𝑗 , 𝑓𝑗 、 𝑥𝑗+1 , 𝑓𝑗+1 を通る放物線 𝑃(𝑥) = 𝑓𝑗−1 (𝑥 − 𝑥𝑗 )(𝑥 − 𝑥𝑗+1 ) (𝑥 − 𝑥𝑗−1 )(𝑥 − 𝑥𝑗+1 ) (𝑥 − 𝑥𝑗−1 )(𝑥 − 𝑥𝑗 ) + 𝑓𝑗 +𝑓𝑗+1 (𝑥𝑗−1 − 𝑥𝑗 )(𝑥𝑗−1 − 𝑥𝑗+1 ) (𝑥𝑗 − 𝑥𝑗−1 )(𝑥𝑗 − 𝑥𝑗+1 ) (𝑥𝑗+1 − 𝑥𝑗−1 )(𝑥𝑗+1 − 𝑥𝑗 ) で置き換えて区間 𝑥𝑗 , 𝑥𝑗+1 で積分。 修正子𝑓𝑗+1 = 𝑓 𝑥𝑗+1 , 𝑦 ∗𝑗+1 𝑓 ∗𝑗+1 未知の𝑦𝑗+1 が右辺にも存在 傾きfj+1 y 𝑦𝑗+1 𝑦𝑗 𝑦𝑗−1 𝑦(𝑥) 傾きfj yj+1－yj 𝑓𝑗+1 𝑓𝑗 傾きfj-1 𝑥𝑗−1 ℎ ※𝑓𝑗+1 = 𝑓(𝑥𝑗+1 , 𝑦𝑗+1 ) この点は不明!! この放物 線f(x,y) 𝑥𝑗+1 𝑓𝑗−1 𝑥𝑗 ℎ x 𝑓 𝑥, 𝑦(𝑥) 𝑑𝑥 𝑥𝑗 h 𝑥 (13.11) を用いる h 𝑥𝑗−1 𝑥𝑗 𝑥𝑗+1 修正子公式の列挙 𝑓𝑗+1 = 𝑓 𝑥𝑗+1 , 𝑦 ∗𝑗+1 ⇐アダムス・バッシュフォース法 2次精度予測子 ⇐台形公式修正子 ⇐アダムス・バッシュ フォース法3次精度予 測子 12 ⇐アダムス・ムルトン法 修正子 13.3 連立微分方程式 ●オイラー法を連立微分方程式に適用する 微分を前進差分で近似 (13.18) 式(13.16)は、 yj+1, zj+1 xj, yj, zj, xj+1 ●3元以上の連立微分方程式への適用 (13.20) yj+1, zj+1 , uj+1 xj, yj, zj, uj, xj+1 (12.18)式の4次のルンゲ・クッタ法を適用すると、 ⇐ホイン法のs1と同じ。 ‥‥(13.22) s1～s4はxj、yjおよびzjで書け る!! r1～r4も同じ!! ⇐s1とr1使ってs2を表す。 ⇐s1とr1使ってr2を表す。 ⇐s2とr2を使ってs3を表す。 ⇐s2とr2を使ってr3を表す。 ⇐ホイン法のs2と似ている。 ただしs3とr3を用いる。 ⇐ホイン法のs2と似ている。 ただしs3とr3を用いる。 ●例1 連立微分方程式の初期値問題 式(13.19)の前進差分式より、 初期値より、 後は陽的に 解法 この例題の場合、𝑦𝑗+1 ± 𝑧𝑗+1 は、 なので、𝑦𝑗 ± 𝑧𝑗 は、 したがって、𝑦𝑗 と𝑧𝑗 の一般式は、 ●例1の解析解双曲線関数 𝑒 𝑥 + 𝑒 −𝑥 𝑦 = cosh𝑥 = 2 𝑒 𝑥 − 𝑒 −𝑥 𝑧 = sinh𝑥 = 2 双曲線関数の性質 cosh𝑥 + sinh𝑥 = 𝑒 𝑥 cosh𝑥 − sinh𝑥 = 𝑒 −𝑥 cosh2 𝑥 − sinh2 𝑥 = 1 cosh𝑥 tanh𝑥 (cosh𝑥, sinh𝑥) sinh𝑥 𝑥 2 𝑥 2 ここの面積が 13.4 高階微分方程式 ●オイラー法を高階微分方程式の初期値問題にも適用 1階微分をzとおくと、 式(13.23)と式(13.24)は、 y 例1 ファン・デル・ポールの方程式 減衰力負の1自由度ばね質量ダンパ系 の運動方程式 𝑚𝑦 − 𝑐𝑦 + 𝑘𝑦 = 0 振幅が指数関数的に成長⇒振幅 がある値を超えると減衰力の符号 を+にする。⇒リミットサイクル m c k ‥‥(13.28) 自励振動の例 タコマ橋の崩壊 https://ja.wikipedia.org とおくと ‥‥(A) (A)式を式(12.6)のオイラー法で解くと、 𝑥𝑗+1 = 𝑥𝑗 + ℎ 𝑦𝑗+1 = 𝑦𝑗 + ℎ 𝑧𝑗 右辺は 左辺はxj+1, 𝑧𝑗+1 = 𝑧𝑗 + ℎ 0.25 1 − 𝑦𝑗2 𝑧𝑗 − 𝑦𝑗 + sin 𝑥𝑗 xj, yj, zj yj+1, zj+1 初期値y(0)=0, z(0)=y’(0)=0, h=0.1の場合 解(x,y) A 解(y,z) C →振幅 と思う →時間と思う B B 図13.2 →振幅の 微分と思う →振幅と思う C 図13.3 A
© Copyright 2024 ExpyDoc
