2005/11/29 数値解析学と演習 福井 先生 # 先週から、常微分方程式を始めて、オイラー法を学んだ 常微分方程式 初期値問題 正規形を利用すると、多階の微分方程式は一階の微分方程式の連立系になる => 一般化しやすい ( 陽的解法 ) 単独のライブラリーで多くの問題が解ける # 偏微分方程式は、個々にプログラムを作成する必要がある 微分方程式の解釈 dy/dt = f(x,y) ^^^^^^ <= 問題を記述している ^^^^^ <= 積分の必要性 この二つの問題は分離して解決できる # 陽的解法の特徴 例 (オイラー法 : text p.478) y(x+h) = y(x) + h f(x,y(x)) ^^^^^^^^^ <= 問題の記述 ^ <= オイラー法の本質的な部分 # Routin として、積分の部分(オイラー / ルンゲクッタ)と問題の計算 ( f(x,y(x)) の計算 ) を分離できるので、一般的なものが作成できる オイラー法 ( 復習 ) y(x+h) = y(x) + h f(x,y(x)) 初期値 (f(0,y(0))) から、順に計算 => O(h) の精度しかでない 一次近似だから.. # 精度は h を小くすれば良いが、そうすると、今度は丸め誤差の問題があるので、h はあまり小くできない 一次 オイラー法 二次 修正オイラー法 ( ホインの方法 : text p.423 ) オイラー法の計算を台形を用いて計算 三次 四次 ルンゲ＝クッタ法 原理的には、更に高階のものも作れる シャンクスの公式 ( text p.494 / p.498 ) 滑らかな関数の近似に向いている # cf ハヤブサの軌道計算 低次 高次 -------------------> 滑らか ( h を大きくできる ) <------------------- タフ ( どんな場合でも利用できる ) 一般に、m 段 n 次の公式 m は、関数の計算回数 n は、近似の次数 オイラー法は、 1 段 1 次の公式 ( m=1, n=1 ) シャンクスの公式 10 段 8 次の公式 ( m=10, n=8 ) # m < n が望ましいのだが.. ありえない # 一般的には # m >= n ( 等号は、 n <= 4 の時のみ ) # となることが解っている 常微分方程式の公式を作成する時の方針 ( 必須ではないが望ましい ) => 係数をシンプルにする 例 : オイラーは、係数が「1」 整数、あるいは、分数 無限小数だと倍精度と単精度で係数をかえる必要がある # 昔は手計算だから、間違いが少いように ## プログラムを作る場合も、同様 => 係数を正の値にする 負の係数があると、桁落 ( 情報落 ) がおきる可能性がある # ルンゲ=クッタ法でも、この原則が守られている [ルンゲ=クッタ法] y(x+h) = y(x) + \frac{ d_1 + 2d_2 + 2d_3 + d_4){6} 1 2 2 1 <= 係数は自然数 # 合計は 6 # 理論的には、1 6 -2 1 でも *あり* # ただ、負の数がでると桁落が置きる ## 全ての係数を正数にしても、関数値が負になる場合は、結局、桁落がおきるのだが.. 関数の計算回数 ( 4 ) に対して、精度が高い ( 4 ) し、公式も簡単なので使い安い公式 # 微係数が解っていれば、直接「テーラー展開法」ができるのだが、一般には、No なので.. ルンゲ=クッタ法の原理 # まず、テーラ展開を行う f(x+h)+y(x)+\frac{h}{1!}f'(x)+\frac{h^2}{2!}f''(x)+\frac{h^3}{3!}f^{(3)}(x)+\frac{h^4}{4!}f^{(4)}(x)+\frac{h^5}{5!}f^{(5)}(x)+.. 区間 [x,x+h] での計算をするのに、その中点、更にその中点(1/4点) を考える d_0 = f(x,y(x)) d_1 = h f(x,y(x)) # この式は、オイラー法と同じ d_2 = h f(x+\frac{h}{2},y(x)+\frac{d_1}{2}) # ここが、数値積分と違う点 # y の値が必要なので d_1 を転用している d_3 = h f(x+\frac{h}{2},y(x)+\frac{d_2}{2}) # d_2 での計算の補正を d_3 で補う d_4 = h f(x+h,y(x)+d_3) # x (独立変数) の方は、数値積分と同様 # y (従属変数) の方は、段々と精度を良くした結果を利用する # 最後に、平均を求める ## 非線形問題として、精度が高くなるように係数を定める ### 係数を决めるのは凄く大変 ( 人手だったらもっと.. 1900 頃 ) #### 係数が正数になるようにすると、「制約」になるので、人手で計算する時には、楽になる要因 #### あとは「対称性」があるとよさそうだということは解る # ルンゲ=クッタ : 百年も利用されている # (数値計算の世界では..) 単純なものが生き残る m,n に関しては、色々なパターンが作れる ( 単段法 ) 一般に、段数をふやせば、精度は上げられるが、段数をふやした計算量に比較して、それほど精度は良くならない。 [単段法/多段法] 単段法では、一ステップでの計算で用いた d_1 .. d_4 は、そのステップだけでしか利用しない。 例: ルンゲー=クッタ # 微分方程式の計算では、出来るだけ、関数計算の回数を減らして、精度を上たい # できるだけ少い m で、多くの n !! 多段法 単段法での計算した関数計算を、再度利用する => 予測子修正子法 [予測子修正子法] まず、新しいステップの値を計算 ( 予測子 ) これまでの計算した値を使って修正する ( 修正子 ) [例] ミルン法 # 計算回数を減らすことができる # 2 回計算して、4 次の公式が作れる # 確かに、計算量が減るが、良いことばかりではない # => 不安定性が生れる ( ミルン自身が指摘 ) [常微分方程式の解法への解釈] 公式を漸化式と考える オイラー法は、一次関数 ルンゲ＝クッタは、三次関数 => 一般に、多項式で近似 常微分方程式の解法 微分方程式を離散化して、差分方程式にしている # 実空間から、離散空間への変換を行っている 微分方程式の根が微分方程式に一致するという保証はない # h を 0 に近付けると一致するということ証明する必要がある 解法で、作られる近似した多項式(特性多項式)の根に 1 がないと、h を 0 にした時に、一致しないことが解っている。 # 関数論の成果 ## 良く知られている解法は、大丈夫 ## 自分で離散化する場合に注意 ## 素直に離散化すれば大丈夫だが.. ## 変な離散化をするとそこで変な事に.. ## cf. 宇野 利雄 先生の「数値計算」浅倉書店 オイラー法 一次式なので、根は一つで 1 しか答がない => (誤差は多いが..) 必ず、正い近似値になる 他の場合 多次式なので、根が複数ある ( 1 でないものもある ) 1 以外の根が、単位円内にないと幻影解が出る 単段法の場合は、区間の内側の点の値しか使わない => 多項式の根が、全て、単位円内にあることが保証できる 多段法の場合は、区間の外側の点の値にも利用する => 多項式の根の中で、単位円の外の根が現れることがある => 幻影解が出る可能性がある ( 不安定 !! ) # 一般には、安定な単段法が望ましい。慣れてきたらい、速い、多段法を用いる == [演習] 11/29 (text. p.448) y'' + \sin{y} = 0 0 <= x <= 16 y(0)=3 y'(0)=3 [解法] \frac{d y_0}{d x} = y_1 \frac{d^2 y_0}{d x^2} = \frac{d y_1}{d x} を利用し、 \frac{d y_1}{x} = -\sin{y_0} \frac{d y_0}{x} = y_1 の形で、ルンゲ=クッタで解く [刻み幅] h = 1/8 == 演習は、今日まで 次回以降は、偏微分方程式をやるが、偏微分方程式は、実習が大変。 ただ、熱伝導程度だと、Excel でできるのですが、Excel が使える人がいなければ演習はなしに.. == 防災研究所の Web ( 微分方程式の応用例 ) 地震の直後にアクセスすると、様々な情報が得られる == 講議はあと四回数 偏微分方程式 2 回 トピック 1 回 まとめ 1 回