2004/12/07 解析学とコンピュータ 先週は、オイラー法、今週は、ルンゲ=クッタ法 オイラー法は、テイラー展開の一次元の項を取ったもの [テイラー展開] f(x+h) = f(x) + h f'(x) + \frac{h^2}{2!} f''(x) + .. + [オイラー法] f(x+h) = f(x) + h f'(x) で近似 ( h^2 以降の項を無視してる.. ) テイラー展開は、(無限次元連続という仮定の元に..) 特定の点による全ての 情報を含んでいる => 無限の計算が可能なら、誤差は入らない 計算機で計算するためには、有限に落すために、どこかで切る必要がある 一次の項目で切れば、オイラー法 一般に、n 次の項目で切るという立場を考えることができる [Text から..] k_1 = h f( x_i, y_i ) k_2 = h f( x_i + \thita h, y_i + \thita h ) where 0 <= \thita <= 1 # (x_0, f(x_0)) から、( x_0 + h, f(x_0+h) ) を計算するのに、 # x_0 と x_0 + h の間のどこの点をとるか ? を \thita で表す としたときに、 y_{i+1} = y_i + a_1 k_1 + a_2 k_2 # where a_1 + a_2 = 1, a_2 \thita = 1/2 # a_1 = 1, a_2 = 0 とすれば、オイラー法と同じ結果 ## つまり、ルンゲ=クッタは、オイラー法の一般化になっている とする。 # ルンゲ=クッタは、4 次元までは、4 回の関数計算で求められる # 5 次の場合は、6 回、計算が必要 !! y_{i+1} = y_i + h(a_1+a_2)f(x_i,y_i) + h^2 a_2 \thita { \frac{\partial}{\partial x}.. # テーラ展開と比較し、係数を色々と揃えて、特定の次数まで係数が # 合うように、a_i を選ぶ ホインの方法 ( 6.17 式 ) a_1 + a_2 = 1 # 1 を積分したら 1 になる必要がある.. a_2 \thita = 1/2 を満すもので、特に a_1 = a_2 = 1/2 \thita = 1 を利用すると.. # これは区間 (a,b) の傾きを、a の傾き # と b の傾きの平均をとっている。 ## 次のステップを計算する時に、一度計 ## 算した値をもう一度利用できる。 # 他のパターンは創りにくいが、作ろうとすれば作れる ## 一次 ( オイラー ) では、工夫のしようがない ### 二次だと工夫は可能だが、まともなのは、一つ ( ホイン ) だけ #### 三次以降だと、色々とありうる 4 次元だと、色々と工夫が可能 標準的な制約としては.. 0 <= \thita_i <= 1 \thita_1 <= \thita_2 <= \thita_3 <= \thita_4 a_i > 0 ( 単純な値が望ましいが... ) \Sum_{i=1}^{4} a_i = 1 位が普通 ( だが.. ) # 次元数が上ると、\Sum_{i=1}^{4} a_i = 1 が守れなくなる # a_i の値は複雑になる.. ## 高次の公式の場合は、係数を色々と選べば、公 ## 式はいくらでも作れる。ただし、精度が保てない ルンゲ=クッタの場合は.. \thita_1 = 0 \thita_2 = \thita_3 = 1/2 \thita_4 = 1 で計算する。 # 宇宙ロケットの起動計算では、宇宙空間に出てしまえば、高次のルンゲ=クッタで Okey # 大気の下の方では、低次の解法を.. # # 実は、中間段階の所が大変で、大気にむらがあるらしい.. ## スペースシャトルは、形が非対称なので、さらに大変 ルンゲ=クッタの一般形式 k_1 = h f(x_i,y_i) k_2 = h f(x_i + \alpha_1 h, y_i + \beta_1 k_1 ) k_3 = h f(x_i + \alpha_2 h, y_i + \beta_21 k_1 + \beta_22 k_2 ) k_4 = h f(x_i + \alpha_3 h, y_i + \beta_31 k_1 + \beta_32 k_2 + \beta_33 k_3 ) y_{i+1} = y_i + \mu_1 k_1 + \mu_2 k_2 + \mu_3 k_3 + \mu_4 k_4 # 次の方針で、係数を選ぶ ( かなりの任意性がある.. ) # 係数を簡単にする、 # 関数値の計算回数をできるだけ減らす よく利用されるルンゲ=クッタ法のパラメータ \alpha_1 = \alpha_2 = 1/2 \alpha_3 = 1 \beta_1 = 1/2 \beta_21 = 0 \beta_22 = 1/2 \beta_31 = 0 \beta_32 = 0 \beta_33 = 1 \mu_1 = \mu_4 = 1/6 \mu_2 = \mu_3 = 2/6 よって、 k_1 = h f(x_i,y_i) k_2 = h f(x_i + h/2, y_i + k_1/2 ) k_3 = h f(x_i + h/2, y_i + k_2/2 ) k_4 = h f(x_i + h, y_i + k_3 ) y_{i+1} = y_i + k_1/6 + k_2/3 + k_3/3 + k_4/6 となる。 真中の点を利用しているので、比較的精度も良い # 対称性がある !! ## 対称性があると、精度もよく安定している # 計算機ができる前に考えられた手法なので、単純さが重要 ## 項目が増えると手計算だと大変 !! # 6 次の公式は、ルンゲ=クッタ=マーチンの公式 ルンゲ=クッタ=ジル法 メモリの消費量が少い ( 係数は厄介だが.. ) # 半分位の量で済む == これまでの解法の応用範囲 \frac{dy}{dx} = f(x,y) # 一階の微分方程式だけ ニ階の微分方程式はどうするか ? [例] \frac{d^2x}{dt^2} + \frac{dx}{dt} + x = a この場合は、補助変数 v を次のように導入する v = \frac{dx}{dt} すると、元の式は、 \frac{dv}{dt} = a - v - x \frac{dx}{dt} = v という連立方程式となる。 一般には、高階の微分方程式も、n 個の連立一次微分方程式にできる 個々の一次元の微分方程式は、これまでの方法で解くことができる # 実は、偏微分方程式も同じような考えかたで、解けてしまう。 ## t (時間) の方向は、ルンゲ=クッタ、x,y (空間) の方向は、差分化すればよい。 ルンゲ＝クッタは、4 次なので、h を 1/2 にすると、誤差は 1/16 となる # 連立方程式の方程式間の関係は、ルンゲ=クッタの中では、関数 f の計算の方に # に吸収されるため、アルゴリズムとしてのルンゲ=クッタには影響しない ## => この為に、汎用なアルゴリズムを設計できる == 講義でやった内容の中で、常微分方程式の解法は比較的、よく出てくる。 連立方程式の係数行列の固有値を考えると、h の取り方との関係が見えてくる ?? 実際の計算の時には「 h をどうするか」という問題がある 計算時間を小くするには、h を大きくすることが望ましいが.. h を大きくしすぎると、関数の値を近似しない ex) y = \sin{x} で、 h = 2 \pi とすると、単なる直線に.. h = \pi/2 とすると、三角波 h = \pi/4 で、なんとか三角関数風.. 一般に、周期の 1/10 位にする 振動数がわかれば、h を决めることができる # しかし、良く考えると、近似する関数の中身がわ # からないので、その振動数がわかるわけがない.. ## 地震波のように、観測値があれば、話は別だが.. h を决めるということは、分解能を决めている h を大きくしすぎると、見たい現象をみることができなくなる 現象の大きさに対応した適当な大きさを决める必要がある # カメラの画素数と同じ 一つの方程式の場合は、h を决めるのであればそれほど問題はない 少しやってみて、だめなら、調整すればよい 複数の現象 ( 速い現象と遲い現象 ) が混ざっている場合は、こまる 周波数の大きい ( 遲い ) 問題と周波数の小さい問題 スティッフ問題 h の最小値は、周波数の小さい方に併せる必要がある ( 波長の 1/3 位 ) 区間の最大値は、周波数の大きい方に併せる必要がある ( 波長の 3 倍位 ) h が小く、区間が大きいの、計算量が増えてしまう。 == [演習] ルンゲ=クッタ法 p.135 にサンプルがある # 連立方程式の形になっていることに注意する 微分方程式は、前回のオイラー法と同じ 精度の違いを比較する == 数値計算の近似系は、テーラ展開が基本 # 解の存在は仮定する必要があるが.. == Q. オイラー法と同様に、Excel でできるか ? A. f の計算が Excel で簡単にできるならば、Okey == # 数値計算をしてわかったこと.. リニアーモータの計算を ルンゲ=クッタ でといたら ? 元々、不安定ものなので、制御系をいれた 地震がおきたら ? # ふっ飛ぶとおもったら.. 安定のための制御系が働き、かえって安定に.. 電車のブレーキ 電気を回収する仕組になっている 不安定な時は ? ラッシュの時かとおもったら .. 始発や終電の時 台数が少いので、発電きばかりで電気を消費する部分がなくなり不安定に.. ラッシュ時は、台数が多いので、どこかの電車が電気を消費している コンセントを抜くと、放電がおきる もし放電がなければ.. ? => 無限の電圧が.. => 回路がこわれてしまう.. 放電があると、わずかながら電流が流れるので、回路が壊れない