5.1 初日

とりあえずI君にやらせて、途中から真打ち (我輩のことだ) 登場。 $a_{ij}$ , $b_{ij}$ の計算を紙の上で計算しておったのを、 ちゃんと do loop で計算するようにした。

$N$ の小さいところは快調に求まるも $N=7$ でエラーが出た。

oyabun% test2
 n=
7

 *** FATAL    ERROR 2 from GVLSP.  Matrix B is not positive definite.
oyabun% 

これはどういうことか？正定値であることは間違いないはずなのに。 ピンと来るものがあって、試しに倍精度にしてみたら、 大丈夫になった。さては、 ということで GVLSP を call する前に $B$ の固有値を計算させてみたら

oyabun% test2d
 n=
7
    1.3158636382684D-11    2.8279747505052D-09    2.8588867231650D-07
    1.8236197276094D-05    8.3397167029818D-04    2.9314481226287D-02
   0.71712364005519
    12.362363368353    485.51884143459    3808.7559124218    14669.662848705
    49899.723376200    126886.20894806    2765041.2146235

となった。なるほど、$B$ は対称で、 固有値も計算してみたら確かにみな正なのだけれど、絶対値がとても小さく、 それで正定値でないと判断されてしまったわけだ。 倍精度にした場合は限界まで少し余裕が出たというわけだね。 とはいえ、かなり条件の厳しい問題らしく、 倍精度にしても $N=11$ 程度が限界で、 $N=12$ ではやはり ``Matrix B is not positive definite.'' と言われてしまった。

それでも、絶対値の小さい方の固有値の収束ぶりはまあまあのものだ。 試しに $N=10$ と $N=11$ の場合の結果を並べて見よう。

oyabun% test2d
 n=
10
    12.362363368326    485.51881851566    3806.5462784675    14617.793474394
    39954.444527442    91907.607463908    184992.73985194    598922.97555703
    1177146.5931475    30675173.939123
oyabun% test2d
 n=
11
    12.362363368326    485.51881851341    3806.5462775007    14617.275831296
    39954.439347729    89277.089283761    184984.26168742    352114.83160918
    1177043.8311788    2683602.5548179    72724488.719649
oyabun%

ちなみにプログラムの中身は

        integer i,j,n, maxn
        parameter (maxn = 100)
        real*8 a(maxn,maxn), b(maxn,maxn), eval(maxn)
        external dgvlsp
c
        write(*,*) ' n='
        read(*,*) n
c
        do 100 i=1,n
           do 110 j=1,n
              a(i,j)=i*(i+1)*j*(j+1.0d0)/(i+j-1.0d0)
              b(i,j)=1.0d0/(i+j+3.0d0)
 110       continue
 100    continue
c
c        call dwrrrn('b=', n, n, b, maxn, 0)
c
c        call devlsf(n, b, maxn, eval)
c        write(*,*) (eval(i),i=1,n)
c
        call dgvlsp (n, a, maxn, b, maxn, eval)
c
        write(*,*) (eval(i),i=1,n)
        end

という感じである。

ここでI君の指摘で加藤先生の本の記述に帰る(P.440 下〜 p.441)。 何でも $\lambda_i$ は方程式

$$\cos\theta\cosh\theta + 1 = 0$$

の解の $4$ 乗であって、例えば $\lambda_1= 12.36236$ , $\lambda_2=485.5217$ であって、 別の方法⁵で得られた $\lambda^{(2)}_2=515.86$ という値は $6\%$ 程度の誤差がある、と。

私は $\lambda_2=485.5217$ という値のシッポはちょっと信じられない、 と思った。$N=11$ までの計算を見ると、 収束している雰囲気があって、 $\lambda^{(2)}_2 = 485.5188185$ 程度はかなりの精度であると思われた。 何といっても古い本だし、非線形方程式の解を疑ってみよう、 と考えた。そこで Mathematica の登場。 最初は Solve を用いたが、駄目だった。無限個の解があるから無理か。 でもそれならば範囲を指定すれば解けるのでは？ と考えて探していたら FindRoot というのを見つけた。

oyabun% math
Mathematica 2.0 for SPARC
Copyright 1988-91 Wolfram Research, Inc.
 -- X11 windows graphics initialized --

In[1]:= FindRoot[Cos[x]Cosh[x]+1==0,{x,4}]

Out[1]= {x -> 4.69409}

In[2]:= FindRoot[Cos[x]Cosh[x]+1==0,{x,4.69},WorkingPrecision->100]

Out[2]= {x -> 4.6940911329741745764363917780198120493898967375457668289728032\

>      77849077936801125617639769164807883401}

In[3]:= 4.6940911329741745764363917780^4

Out[3]= 485.5188185133710378116913838

In[4]:= Quit
oyabun%

やはり、と思った (こんな細かい計算でツッパルのは大人げないな)。