2004/06/08 先週 : ニュートン法 今週 : 代数方程式 これまでの解法は、「非線形方程式の解法」 範囲が広い ( = 条件が少い.. ) ので、解法も限られる 代数方程式も、非線形方程式の一種 代数方程式 ( f(x) が x の多項式の場合 ) f(x)=0 f(x) = \Sum_{i=0}^n{a_i x^i} # 代数方程式は、非線形方程式より条件 ( = 情報 ) が多いので、特別な解法がある # 後期に行う、「有限要素法」も利用できる n < 5 の時には、根の公式がある n = 1 a_1 x + a_0 = 0 n = 2 a_2 x^2 + a_1 x + a_0 = 0 n = 3 a_3 x^3 + a_2 x^2 + a_1 x + a_0 = 0 (カルダン法) n = 4 a_4 x^4 + a_3 x^3 + a_2 x^2 + a_1 x + a_0 = 0 (フェラリ法) n = 5 a_5 x^5 + a_4 x^4 + a_3 x^3 + a_2 x^2 + a_1 x + a_0 = 0 (根の公式は存在しない) # 結局、n > 4 の場合は、反復法を適用する必要がある # 代数方程式も、非線形方程式なので、普通にニュートン法などが利用できるが.. # 代数方程式の性質を利用して、更に、「旨い」解法を考えることができる n 次代数方程式の性質 ( 代数の基本方程式 ) n 個の ( 複素 ) 根を持つ ( 重根も含めて数える )。 # 代数方程式を実数だけで、解こう (ベアスト法) とすると、破綻することが多い # 代数方程式の複素数で閉じた性質を持つ # 複素係数の複素根を求めるニュートン法が一番素直 ベアスト法 実数係数の代数方程式を「ほとんど実数だけ」で解く 実数係数の代数方程式の根の一つが複素数ならばその複素数の共役複素数も根 つまり元の方程式は、その複素根を根とする実数係数の二次方程式の因子を持つ 二次方程式の係数 ( 実数 ) を求める方法 # 計算機で複素数を扱うのは不便なので.. # 実数係数だけで求めることができるので便利そうだが.. # あまり上手くゆかない ニュートン法は根を一つ求める方式 ベアスト法は、根を同時に二つもとめようとしている # そこに無理が.. ( 複素空間を実数空間に押し込めようとすると無理が生じている [問題を捻じまげている..] ) # 解こうとすると觧けないことが多い # 物事は、1 は簡単、2 は難しい、3 は、2 と同じようなもの.. 普通に、ニュートン法を複素数に適用すれば問題はない ( 素直に扱うことができる ) # 実数部と虚数部を別に扱う必要があるので、計算量は増える 問題を素直に考える 引き算 : 整数 割り算 : 有理数 平方根 : 複素空間 代数方程式を解くには、複素数が「必須」と考えてよい !! # これが「素直」なアプローチ # 問題に、特別な性質があれば話はべつだが.. # 元の方程式が、「実根しかないと判っている」とか.. == 解法 根の推定 、修正 # 修正の方法はある ( 二分法やニュートン法は、その具体的な方法 ) # 推定の一般的な方法はない # そもそも、推定ができる ( 根が判っている ) なら計算は不要 !! # 一般的な方法もあるが、当らないこともある.. ## 代数方程式の場合は、根の推定が可能.. f(x) の計算 能率よく計算するには、どうすればよいか ? という問題 f(x) の計算の効率化 代数方程式の計算個数 n * (かけ算) + (足し算) 1 1 1 2 3 2 3 6 3 4 10 4 5 15 5 i i(i+1)/2 i なんとか工夫できないか ? n = 1 の時 a_1 x + a_0 どうにもならない n = 2 の時 a_2 x^2 + a_1 x + a_0 ( a_2 x + a_1 ) x + a_0 足し算の数は変わらない かけ算の数が減っている n = 3 の時 a_3 x^3 + a_2 x^2 + a_1 x + a_0 ( ( a_3 x + a_2 ) x + a_1 ) x + a_0 かけ算の回数は n 次式で n 回数 ( 足し算も n で変らない ) となる => ホーナー法 ( 組立除法 ) n(n+1)/2 が n になるので (n+1)/2 倍も速くなる ホーナー法の性質 計算量が少い プログラムも簡単 ( すなおにかける ) (実は..) 計算精度も良い # 足し算をするときに、同じ位の大きさの数を加えている # 普通は、計算量を減らす ( 手間を惜しむ.. ) と、精度が悪くなる # ホーナー法は、計算量を減らしつつ、精度もよくなる (珍しい) 良い方法 ホーナー法の精度が良い理由 x^n は、直に大きくなる ( オーバフローする ) x が大きい場合は、x^n が計算できない でも f(x) が 0 に近いはずなので.. a_n x - a_{n-1} も 0 に近いはず !! 直接計算すると x^n がオーバフローする場合も ホーナー法だと途中の計算で値がキャンセルして計算可能 == 複素数ニュートン法 # 色々工夫するより、これが一番、素直 代数方程式 f(x) のニュートン法も基本は同じ x_{k+1} = x_k - \frac{f(x_k)}{f'{x_k)} ただ、f(x) が多項式なので、微分が容易であるため、計算がやりやすい # f(x) の計算と f'(x) の計算を並列して計算する # Fortran は、Complex 型が利用できるので、実数 ( double ) と同 # じように考えればよい。 ## 昔の Fortran だと、絶対値 ( abs ) の計算を ( cabs など ) 別 ## の関数で行う必要があった今の Fortran なら、自動的に選択され ## るので、問題なし。 == 収束判定 どこで、計算をストップすべきか ? # 二分法 : 一次収束なので、有効桁数の bit 数できまる 一般には、十分に根に近ければ終りにしたい # 無駄に回しても時間が無駄になるだけ !! f(x_k) が 0 に近い ? x_k が根に近い ( これは本来的な意味 ) x_{k+1} - x_k が 0 に近い ? x_k が収束しつつある ( これ以上計算してもあまり改良されない.. ) y = f(x) と考えると.. 前者 : y で考えている 後者 : x で考えている どちらが望ましいか ? f'(x) の振舞いによる ? y = f(x) と x 軸の交わり方による 45 度位ならばどっちでも同じ x と y がほぼ直交 交わりが浅い場合 x と y が直交していない x の変化に対して y の変化が小さい場合が問題 !! ニュートン法の場合 !! 結果的に同じ ? !! f'(x) で割っているだけ ??? 二次方程式の場合 根が近接していない 交わりが急 どっちも可 根が近接していいる ( 重根に近い ) 交わりが斜め 方程式の判別式 D = b^2-4ac が小さい 根は、+- \sqrt{D} の差をもっている 平方根の振舞い 1 より大きい数は、小くなる 1 より小さい数は、大きくなる 結果的に 1 に近くなる 誤差を拡大する結果になる n bit の平方根は n/2 bit の精度しかない # 有限桁の計算であることが「ポイント」 重根ということは、(誤差が拡大されている) \sqrt{D} の値が根の値を決め手いる 精度が悪い !! x で比較すると、精度が悪い場合に、収束しない !! f(x) で比較すると、精度が悪い場合でも、収束する # 代数方程式では、重根があると、半分しか精度がでない !!!! # 三重根だと、1/3 の精度しかでない.. # n 重根だと、1/n の精度しかでない.. # 精度がでないので、「空回り」してしまう.. ## 有限桁での計算では問題となる ( 無限桁が半分でも無限なので気にしない.. ) ### 後で、 \frac{f(x+h)-f(x)}{h} の精度が 1/2 になることも示す == 初期値の問題 n 次のニュートン法 # 日大の生産工学部 ( 平野菅保 先生 発案 ) どの初期値から初めても、必ず根に収束する # 数学的な証明もあり !! テーラ展開 f(x+h) = \Sum_{n=0}^{\infty} \frac{1}{n!}\frac{df^n}{dx^n} これをどう解釈するか ? 一点での情報が全て ( 無限回数の微分係数 ) が判れば、全ての値が解る # 一般には、無限回数の微分係数はわからないので無意味 ## k 次以後は無視している ### x の近傍 ( k による ) だけがわかる #### その廻りなら根が収束する n 次の代数方程式は、n 次微分係数が 0 なので、全ての情報が得られる 全ての情報が得られているので、全ての場合が Okey n 次の代数方程式は、n 次のニュートン法で必ず全ての空間で、根が収束する # 普通のニュートン法は、一次なので、根が収束しない場所が生じてしまう.. == [まとめ] 代数方程式は、複素ニュートン法で解くのが素直 代数方程式の関数値の計算は、ホーナー法を利用する 初期値の定め方 旨い方法はない ( 平野の方法ならば、Don't care.. ) == 課題 ホーナー法のプログラム 多項式の係数 a_i と、独立変数 x_0 の値を与えて、関数値 f(x_0) = \Sum_{i=0}^n a_i x_0^i を計算する [ヒント] S = a_n として、 S = S * x + a_{n-1} を繰り返すだけ [締切] 6/14