ちょっと暇ができたので、scheme分を補充しておこう……。Durand-Kerner法(DK法)というのはニュートン法の改良版で、代数方程式の全ての解を一度に求めることができる。

数式の表示にはMathMLを使っています。また僕は数値計算にあまり詳しくないため、間違いなどありましたら指摘してください。

Durand-Kerner法で多項式の近似解を求めるプログラム

(require (lib "1.ss" "srfi"));SRFI-1のロード

(define PI 3.1415)
(define (square x) (* x x))
  
(define (dka input)     ;Durand-Kerner法
  (define formula (if (proper-list? input) 
                      (drop-while zero? input) ;最高次係数が0なら削除
                      (error "INVALID INPUT")))
  (define degree (- (length formula) 1))
  (define degree_list (iota degree 1 1))　
                        ;次数nに対して、'(1 2 ... n)というリストを作成
  (define zcenter 
    (if (>= 0 degree) (error "INVALID DEGREE")
        (- (/ (cadr formula) (* (car formula) degree)))))
  
  (define (get-R list)  ;解を囲む円の半径の候補を求める
    (let ((denom (car list)))
      (map (lambda (x k) 
             (expt (/ (abs x) denom) (/ 1 k))) 
           (cdr list) degree_list)))

  (define trial-list　　　　;Aberthの初期値
    (let* ((angle_base 
            (/ (* 2 PI (make-rectangular 0 1)) degree))
           (R (* degree 
                 (fold max 0　(get-R (kumitate-exp formula zcenter))))))
      (map (lambda (x) 
             (+ zcenter (* R (exp (* angle_base (- x 0.75)))))) 
           degree_list)))

  (define (dka-iter list)　　;反復
    (let ((deriv-list　
           (map (lambda (x)
                  (subst　(get-form (kumitate formula x)) x))
                list)))　　　　;それぞれの近似値における導関数値を求める
      (cond ((< (abs-list 
                 (map (lambda (x) (subst formula x)) list))
                0.0001) list);それぞれの近似値を与式に代入し、誤差を調べる
            (#t (dka-iter 
                 (map (lambda (x y) 
                        (if (= y 0) x (- x (/ (subst formula x) y)))) 
                      list deriv-list))))))
  (cond ((= degree 1) zcenter)
        (#t (dka-iter trial-list))))
             
(define (kumitate-exp list x) ;組み立て除法を用いて式を変形
  (define (kumitate-t-iter list ans)
  (cond ((null? list) ans)
        (#t (let ((form (kumitate list x)))
              (kumitate-t-iter (get-form form) (cons (last form) ans))))))
  (kumitate-t-iter list '()))

(define (kumitate list x);組み立て除法
  (define (kumitate-iter list v)
    (cond ((null? list) '())
          (#t (cons (+ (car list) (* v x)) 
                    (kumitate-iter (cdr list) (+ (car list) (* v x)) )))))
  (cons (car list) (kumitate-iter (cdr list) (car list))))

(define (get-form list)	;組み立て除法後のリストlistから、式の部分を取得
  (drop-right list 1))

(define (abs-list list)
  (fold + 0 (map (lambda (x) (+ (magnitude x))) list)))

(define (subst form v)	;式formに値vを代入した結果を返す
  (last (kumitate form v)))

リスト操作ライブラリ

SRFI-1を使う。proper-list?、 last、drop-while、drop-right、fold、iotaあたりがそう。便利すぎ。

式の表現

多項式は最高次係数を先頭とするリストであらわすことにした。つまり$x^2-5x+6$は'(1 -5 6)。

式にたいして値が代入できなきゃどうしようもないので、それを行う関数substを定義。計算には組み立て除法を用いる。

> (subst '(1 -5 6) 2)
0
> (subst '(1 -5 6) 4)
2

組み立て除法

組み立て除法というのはある多項式を (x-a) といった一次式で割るのに適した方法。それを使って商関数Q(x)と剰余Rを求めると、 $f(x)=(x-a)Q(x)+R$ より、 $f(a)=R$

> (kumitate '(1 -4 -2 34 2 4 3) 2)
(1 -2 -6 22 46 96 195)
> (kumitate '(1 -2 -6 22 46 96) 2)
(1 0 -6 10 66 228)

組み立て除法によって導関数値を求める

組み立て除法を繰り返すと、 $x^6-4x^5-2x^4+34x^3+2x^2+4x+3=(x-2)(x^5-2x^4-6x^3+22x^2+46x+96)+195={(x-2)}^2(x^4-6x^2+10x+66)+(x-2)228+195$ ...と更なる変形が得られる。これを使うとaに対するm次の導関数値が比較的簡単に求められる。

いま、f(x)が $$f(x)={(x-a)}^{m}Q(x)+{(x-a)}^{(m-1)}{R}_m...$$ と変形できたとして、点 x = a における m 次導関数値を求めたい。m 回微分すると次数がm-1以下の項はすべて消えることから、 $$f^{(m)}(x)={({(x-a)}^{m}Q(x))}^{(m)}$$ としてよい。

$${(g(x)h(x))}^{(m)}=\sum _{k=1}^m\left(\genfrac{}{}{0ex}{}{m}{k}\right)g^{(m-k)}(x)h^{(m)}(x)$$ を考える。これを展開したとき、( x-a )が含まれる項はxにaを代入したとき全て0になるので、結局 $$f^{(m)}(a)=m!Q_m(a)$$ が言える。まーここでは一次導関数値しか使わないけど。

Aberthの初期値

ある代数方程式の解が含まれる領域を推定する方法 。直感的には、解の分布の中心zcenterから、その円内にすべての解が含まれるように半径Rの円を描いて、その円周上に初期値を配置する。こうすることで近似値それぞれが真値に近づくスピードをそろえることができる。

解の中心は簡単に求められる。なぜなら、$a(x-b)(x-c)(x-d)=a{x}^3-a(b+c+d)x^2+...$

半径を求めるには、与式を解の中心zcenterで組み立て除法を用いて変形し、r = ( z - zcenter ) と置換することで円の半径 r が解となるような方程式を立てればよい。

結局、n個の初期値を円周上に等間隔に配置するには $$z_j=-\frac{a_1}{n{a}_0}+Rexp(\frac{2\pi i}{n}(j-\frac{3}{4}))$$ を j = 1～n について計算すればよい。

こうして求めた初期値はリストとしてtrial-listに格納される。

  (define (get-R list)  ;解を囲む円の半径の候補を求める
    (let ((denom (car list)))
      (map (lambda (x k) 
             (expt (/ (abs x) denom) (/ 1 k))) 
           (cdr list) degree_list)))

  (define trial-list　　　　;Aberthの初期値
    (let* ((angle_base 
            (/ (* 2 PI (make-rectangular 0 1)) degree))
           (R (* degree 
                 (fold max 0　(get-R (kumitate-exp formula zcenter))))))
      (map (lambda (x) 
             (+ zcenter (* R (exp (* angle_base (- x 0.75)))))) 
           degree_list)))

Durand-Kerner法

上の方法で求めた初期値それぞれについて、ニュートン法による近似を行う。本来は導関数値をさらに近似する。

$P_N(x)=a_0(x-a_1)(x-a_2)...(x-a_n)$とすると、その一次導関数は $$P_N^{(1)}(x)=a_0(x-a_2)\cdots (x-a_n)+a_0(x-a_1)(x-a_3)\cdots (x-a_n)+a_0(x-a_1)(x-a_2)(x-a_4)\cdots (x-a_n)+\cdots =a_0\sum _{i=1}^n\prod _{j=1,j\ne i}^n(x-a_j)$$ ここで、$z_i\approx a_i$ とすると、 $a_i$ の含まれる項は全て消えるので、 $$P_N^{(1)}(z_i)==a_0\prod _{j=1,j\ne i}^n(z_i-a_j)$$ このとき、$z_j\approx a_j$もまた言えるので、 $$P_N^{(1)}(z_i)==a_0\prod _{j=1,j\ne i}^n(z_i-z_j)$$

反復処理

  (define (dka-iter list)　　;反復
    (let ((deriv-list　
           (map (lambda (x)
                  (subst　(get-form (kumitate formula x)) x))
                list)))　　　　;それぞれの近似値における導関数値を求める
      (cond ((< (abs-list 
                 (map (lambda (x) (subst formula x)) list))
                0.0001) list);それぞれの近似値を与式に代入し、誤差を調べる
            (#t (dka-iter 
                 (map (lambda (x y) 
                        (if (= y 0) x (- x (/ (subst formula x) y)))) 
                      list deriv-list))))))

まず組み立て除法を用いて各近似値に対応する導関数値をそれぞれ計算し、deriv-listとする。それぞれの近似値を与式に代入してみて、全てが0に近ければ十分収束しているとして解を出力、終了。まだダメそうならmapを用いて、近似値全てをニュートン法で再近似。

導関数値リストはもともと別の場所でも使うはずだったので一次変数にしてあるが、その部分を消したためdka-iterの反復呼び出しの際に直接計算してもいい。でもそうするとすごいわかりにくくなると思うのでそのまま。上で示したように、点aにおける一階導関数値を求めるにはまずkumitateを用いて一回変形し、変形後のリストからget-formを用いて式の部分のみを抜き出し、それにaを代入すればよい。

結果

$$x^2-2x+1=0$$ $$x^3-3x^2-10x+24=0$$ $$x^4-5x^3-21x^2-22x-16=0$$ $$8x^5+24x^4-25x^3+3x^2+98=0$$ について、

> (dka '(1 -2 1))
(1 1)

> (dka '(1 -3 -10 24))
(3.9999998318967864+2.8387151406426836e-008i
 -3.000000000000001+1.441359404959568e-017i
 2.0000000001780975+1.4040173748418967e-010i)

> (dka '(1 -5 -21 -22 -16))
(8.0+0.0i
 -0.4999999999911601+0.8660254037577851i
 -2.0000000003879568+5.235939426694963e-010i
 -0.500000298422228-0.8660254066797882i)

> (dka '(8 24 -25 3 0 98))
(1.151630579359255+0.817203444821391i
 -0.7099165021391798+1.038229298587666i
 -3.883428154439853-9.183549615799121e-041i
 -0.7099165035189984-1.0382292988192343i
 1.151630579277084-0.8172034448086737i)

それぞれ誤差は、

0
2.3890323896364474e-006
7.669790640016116e-006
4.0979890195062324e-007

もっとシンプルにできそうな予感はするんだけど思いつかないなー。

集合への30講

集合への30講を読み終えました。整列集合あたりからは僕の想像力の限界を突破してるのでほとんど眼で追っただけだけれども、その難しさはわかった気がする。本書のコラムの中に、2つの異なる実数をとるとは、現実に何を意味しているのだろうか。という言葉があるんですが、うーん、何を意味しているんでしょうね。位相への30講が終わったらもう一回読んでみようと思います。

枝葉の部分でちょっと気になったのは、この本は連続体仮説に連続体「仮設」って文字をあててるんですよね。英語がCH、Continuum Hypothesisだから訳としては仮説なんだろうけど(Googleなんかでも仮説が圧倒的に多い)、証明も反証もできないという意味を考えて仮設にしてあるのかな。

数学を勉強しています。集合・位相に手をつけようと思い、とりあえず図書館へ行って志賀浩二先生の「集合への30講」、「位相への30講」二冊を借りてきた。アレフってどこかでやったな、と思って昔の教科書を調べると、計算機科学入門(アービブ)という本に8Pくらい載ってた(同名でゴールドシュレーガーという人が書いてるのもある。こっちはハードウェアや人工知能にも触れる内容)。どうやら情報基礎学という講義で習ったことがあるようだ。読み返してみるとこの本けっこう盛りだくさんにいろいろ載ってる。まーそんなことはどうでもよくて、集合への30講を読み始めましたよ、という話。

空でない集合Γから、任意の集合Aに対して写像 $\gamma (\in \mathrm{\Gamma })\to A_{\gamma }(\in P(A))$が与えられたとき、Γを添数（てんすう）とする集合族が与えられたといい、これを $$\{A_{\gamma }{\}}_{\gamma \in \mathrm{\Gamma }}$$ と書く。

ある中学校の生徒全員の集合をAとして、学年{1, 2, 3}をΓとすれば、 $\gamma \in \mathrm{\Gamma }$ によってある1学年の生徒の集合が取り出せる。

このとき、その直積集合は $$\prod _{\gamma \in \mathrm{\Gamma }}{A}_{\gamma }$$と書き、上の例では1学年からひとりずつ生徒を選ぶ組み合わせを意味する。