I 君の固有値問題
桂田 祐史

1992 年 11 月 20 日 (火)古い記録なのだが、懐かしいのと、現在でも参考になるところが多いので…

発端

1

それは K 先生のところの学生 I 君からの相談で始まった。
数学概論応用編、C1 第 6 章 p440 にある固有値問題∂ 4u= λu in [0, 1],∂x4
u(0) = u′ (0) = 0,u′′ (1) = u′′′ (1) = 0の数値計算が K 先生から彼に与えられた課題。本に書いてある計算法は次の通り。数列 {aij },{bij } をi(i + 1)j(j + 1)
1
aij =, bij =
,
i+j−1
i+j+3
で定めて、各 N = 1, 2, · · · に対して N 次正方行列 A = A(N ) , B = B (N ) を

A = (aij ),B = (bij )で定めた時の一般化固有値問題

det(A − λB) = 0
(N )

(N )

(N )

の解を小さい順に並べたものを λ1 , λ2 , · · · , λN(N )λi = lim λiN →+∞とする。これから(i ∈ N)によって得られる {λi ; i ∈ N} が最初の固有値問題の解になる。
(N )行列の次数を N として、N が小さい時の λi の計算値は載っている。まず N = 1 の時、(1)

λ1 = 20,N = 2 の時
(2)

λ1 = 12.480,1(2)

λ2 = 1211.520

加藤敏夫先生の変分法だ！
！

1

となるそうだ。
一般化固有値問題というと、どういう問題であるか聞いたことはあるものの、具体的にどうやって解くのがカシコイというような知識はほとんどなかった。典型的な “教えながら勉強する” パターンになりそうだ、と思った。

困ったもんだ、見つからない2

まずは自室にある邦書で解説を探したが、ほとんど記述なし。どこかで読んだ覚えがあったのだけど（後で分かったが一松先生の本だった）。
Numerical Recipes には少し載っていたが、慌てていたのであまり読まなかった（本当はちゃんと読むべきだった）。とりあえず I 君には Jacobi 法のプログラムの載った本2 を渡して、勉強しておきなさい、次は出張から帰ってきたらね、と逃げておいた。
B −1 をかけたらどうかなぁ、とか思いつつ。いわきへと旅立つ。

いわきにて

3

質問出来る人を探す。偉い人3 を紹介されもしたが、ちょっとおそれ多いので、適当なとこ
ろで K さんに SOS。K さんによれば、逆行列をかけてでもいいでしょ、でも直接やるのでも色々ある、文献としては Cleve Moler が 1970 年代後半に SIAM Review に書いたもの、さらに
Van Loan and Gene H. Golub, Matrix Computations, John Hopkins University Press (1988)
でしょう、これは良く売れている本ですよ、とのこと。ははあ。聞いて良かった。さらにすぐ使えるサブルーチンが、EISPACK や IMSL の中に見つかるでしょう、とおっしゃいました（急に敬語になる）。

4

再び明治にて —作戦を立てる—

4.1

数学科資料室

早速 K さん推薦の文献を探す、、、ない！書籍の方も有隣社のリストにあがっていない。困ったものだ。書籍はそのうち何とか買うとして、文献（雑誌）の方は暇を見つけて東大数学科にコピーしに行くしかないな。

4.2

一松先生の本

本棚にある本を端から順番に邦書をチェックしていって、ようやく見覚えのある記述に遭遇した。一松信著、数値解析、朝倉書店の第 3 章 (P61) に曰く。

2
3

戸川隼人, 数値計算, 岩波書店。今回は直接役には立たなかったが、いい本である。
森正武、名取亮、両先生の名前があがった。

2

§11. 数学では普通にはあまりうるさく区別しないが、実は行列の固有値問題にはいろいろの意味がある。定数係数連立線型常微分方程式系 y ′ = Ay の係数行列 A
は混合テンソルであって、座標系の変換の標準化による Jordan の標準系が自然な意味をもつ。これに対して、2 次形式の係数行列 A の固有値は、標準的な正値対
称行列 B を定めて初めて意味がある。普通には B を単位行列 I としているのだが、その意味では、“一般化された固有値問題” とよばれる det(A − λB) = 0 の方がかえって本質的である。§16 の m × n 次の長方形行列 A の特異値も、m 次, n
次の標準正値対称行列 B, C を定めて初めて意味がある。これも普通には B, C をそれぞれの次元の単位行列に定めて扱う。以上は数値解析自体の問題ではないかも知れないが、実用上には重要な背景である。なおこの点を明確にご注意下さった
東京大学伊理正夫教授に感謝する。
とか。やはり伊理先生は偉いんだ、と再認識したが、この問題に関しては直接得るところは
なし。

4.3

Numerical Recipes

出張前に Numerical Recipes をきちんと読まなかったのは失敗だった。Cleve Moler の文献
まではなかったが、Van Loan & Gene H. Golub の本 (正確には K さんが勧めてくれた本の前の版に相当するもの）は References に載っていたし、何より固有値問題に対する戦略について色々指示があって、それに従っていれば K 大明神の御神託を聞かなくても実用的な観点からは何とかなっていた。
我々に関係ある部分の要旨を大雑把に述べると、(1) 固有値問題は難しく、既成のパッケージの利用を避けることを勧めない数少ない例である（英国人て回りくどくて、やあね。つまり既成のパッケージを使いましょうということだな。）。この本では、その種のパッケージの中で何をやっているかの理解の助けになるようなプログラムを載せる。(2) 一般化固有値問題には、
B の平方根を求めて普通の固有値問題に帰着させるという解法もあるが、既存のパッケージを利用するのが簡単で効率的だ。大抵のパッケージには一般化固有値問題のルーチンがある。
具体的なパッケージの名前の中に IMSL があった。この時点で私の頭の中は IMSL の利用
で固まっていった。

4.4

I 君来る

ここで突如 I 君登場 (作ったような話だが本当)。彼は K 先生 (K さんにあらず…同じイニ
シャルだけど) から渡されたコピーを持ってきた。そこには Jacobi 法で一般化固有値問題を解く方法の説明が載っていた。これは本質的には Numerical Recipes に説明のあった平方根を用いる方法であった。なるほど、解き方が載っている本があったか、でも、一からプログラムを書き起こすのは面倒だな、どうしようか、と少し迷った。なお、この解説に書いてあったのだが、A, B が対称であっても、B −1 A は対称とは限らないので、B −1 A の固有値問題に帰着させるのはうまくない、とのこと。なるほど、と感じていたところ、I 君から質問があった。
B −1 A が非対称になっても、この間渡された本に載っていた非対称行列の固有値を求めるサブルーチンを使えば解けるのではないか、と。おお、鋭い (人の心中が読めるのだろうか)。自3

力で解く方法を見出したのは偉い、その線でやらせようか、と一瞬思ったが、やはり非対称行
列用の固有値ルーチンは (使用を避ける道がある限り) 使うべきでないという正論を説明。先日の研究集会での某先生の言葉が甦る。曰く、非対称行列の固有値問題を解くのはアートである4 、と。で、二人で IMSL を調べましょうということになった。

4.5

IMSL

マニュアルを引いたら、大した苦労もなしに GVLSP なるルーチンが見つかった。能書きに曰く、Purpose: Compute all of the eigenvalues of the generalized real symmetric eigenvalue
problem A*z = w*B*z, with B symetric positive definite. これですね。サンプル・プログラムも短いので、IMSL の使い方の復習を兼ねて打ち込んで実行してみた。
コンパイルは

f77o sample393.f -limslib
で OK でした。実行して即 OK となりました。
うん、これで何とかなるでしょう。

5

戦い

5.1

初日

とりあえず I 君にやらせて、途中から真打ち (我輩のことだ) 登場。aij , bij の計算を紙の上で計算しておったのを、ちゃんと do loop で計算するようにした。
N の小さいところは快調に求まるも N = 7 でエラーが出た。

oyabun% test2n=7*** FATAL
oyabun%ERROR 2 from GVLSP.

Matrix B is not positive definite.

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

oyabun% test2dn=71.3158636382684D-11
1.8236197276094D-05
42.8279747505052D-09
8.3397167029818D-04

本当は英語だった。

4

2.8588867231650D-07
2.9314481226287D-02

0.71712364005519
12.362363368353
49899.723376200

485.51884143459
126886.20894806

3808.7559124218
2765041.2146235

14669.662848705となった。なるほど、B は対称で、固有値も計算してみたら確かにみな正なのだけれど、絶対値がとても小さく、それで正定値でないと判断されてしまったわけだ。倍精度にした場合は限界まで少し余裕が出たというわけだね。とはいえ、かなり条件の厳しい問題らしく、倍精度にしても N = 11 程度が限界で、N = 12 ではやはり “Matrix B is not positive definite.” と言われてしまった。
それでも、絶対値の小さい方の固有値の収束ぶりはまあまあのものだ。試しに N = 10 とN = 11 の場合の結果を並べて見よう。

oyabun% test2dn=1012.362363368326
39954.444527442
1177146.5931475
oyabun% test2dn=1112.362363368326
39954.439347729
1177043.8311788
oyabun%485.51881851566
91907.607463908
30675173.939123

3806.5462784675
184992.73985194

14617.793474394
598922.97555703

485.51881851341
89277.089283761
2683602.5548179

3806.5462775007
184984.26168742
72724488.719649

14617.275831296
352114.83160918

ちなみにプログラムの中身はinteger i,j,n, maxn
parameter (maxn = 100)
real*8 a(maxn,maxn), b(maxn,maxn), eval(maxn)
external dgvlsp
cwrite(*,*) ’ n=’read(*,*) n
c

110
100
c
c

do 100 i=1,ndo 110 j=1,na(i,j)=i*(i+1)*j*(j+1.0d0)/(i+j-1.0d0)b(i,j)=1.0d0/(i+j+3.0d0)
continue
continue
call dwrrrn(’b=’, n, n, b, maxn, 0)
5

c
c
c
c

call devlsf(n, b, maxn, eval)
write(*,*) (eval(i),i=1,n)call dgvlsp (n, a, maxn, b, maxn, eval)cwrite(*,*) (eval(i),i=1,n)endという感じである。
ここで I 君の指摘で加藤先生の本の記述に帰る (P.440 下〜 p.441)。何でも λi は方程式

cos θ cosh θ + 1 = 0の解の 4 乗であって、例えば λ1 = 12.36236, λ2 = 485.5217 であって、別の方法5 で得られた(2)
λ2 = 515.86 という値は 6% 程度の誤差がある、と。
私は λ2 = 485.5217 という値のシッポはちょっと信じられない、と思った。N = 11 までの(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

5改良された基底関数を取った場合。

6

In[4]:= Quitoyabun%
やはり、と思った (こんな細かい計算でツッパルのは大人げないな)。

5.2

二日目

これを書いているのは実は三日目なのだが、大変な問題に関わってしまったと思い始めている。
初日ので大体話は済んだと思っていたのだが、I 君は別の計算を始めていた。そもそもこの
計算は Galerkin 法に基づいているのだが、これまで使ってきた基底関数は境界条件を満足していないので、あまりうまくないのだという。境界条件を満たす基底関数を使うと収束が良くなる、と本に書いてあった。ところがそこには基底関数が 2 つしか書かれていなかった:

w1 (x) = x2 (6 − 4x + x2 ),w2 (x) = x3 (10 − 10x + 3x2 )そこで基底関数の一般形を求めることが問題。はっきりとは書いていないが、加藤先生の頭にあったのは

wn (x) = xn × (x の 2 次式で 1 における 2 階、3 階の微分係数が 0 になるもの)ということらしい。こういうものを具体的に求めるのは高校数学の問題だが、面倒だね。I 君に紙の上で計算させたが、時間がかかりそうだったので、Mathematica を持ち出して計算を始めてみた。実際にどういう計算をしたのかは、ログに残っているが、エッセンスをプログラムにまとめると、

u=x^(n+1)(x^2+a x+b)
DDu=D[u,{x,2}]
u2=Simplify[DDu /. x->1]
DDDu=D[u,{x,3}]
u3=Simplify[DDDu /. x->1]
sol=Simplify[Solve[{u2==0,u3==0},{a,b}]]ということである。その結果はoyabun% math
Mathematica 2.0 for SPARC
Copyright 1988-91 Wolfram Research, Inc.
-- X11 windows graphics initialized -In[1]:= << ishibashi1.m
2-2 (3 + n)
6 + 5 n + nOut[1]= {{a -> ----------, b -> ------------}}
7

1 + n2n + nIn[2]:=となる。そこでwn (x) = xn+1 [n(n + 1)x2 − 2n(n + 3)x + (n + 2)(n + 3)]という基底関数の一般形が求まる。これは n = 1, 2 の時に加藤先生のあげたものと一致する(n = 1 の時は加藤先生の関数を 2 倍したもの)。
僕がこの計算を終える頃には I 君も大体計算を終えていて、同じ結果を得た (まあ間違っていない限り同じになるのは当たり前だけど)。
これで話が済んだわけではない。この基底関数を用いて行列 A, B を計算しないといけない。実はaij = (wi′′ , wj′′ ), bij = (wi , wj )ということだそうだ。ここで (, ) は区間 (0, 1) 上の内積を表す。
あまり手でやりたくない計算だ、ということで Mathematica との格闘が始まった。これも一部始終がログに残っているが、エッセンスをまとめると、%
u=x^(i+1)((i+1)i x^2-2i(i+3)x+(i+2)(i+3))v=x^(j+1)((j+1)j x^2-2j(j+3)x+(j+2)(j+3))
intuv=Integrate[u v, {x,0,1}]
uv=Simplify[intuv]
Table[uv, {i,5}, {j,5}]
FortranForm[uv]Print[Table[uv, {i,15}, {j,15}]]
Du=Simplify[D[u,{x,2}]]Print["Du=",Du]
Dv=Simplify[D[v,{x,2}]]Print["Dv=",Dv]
DuDv = Simplify[Du Dv]
Print["DuDv=",DuDv]
tmp=Integrate[DuDv /. i+j-2->m, {x,0,1}]
tmp=Simplify[tmp]u2v2 = tmp/. m->i+j-2
Print["u2v2=", u2v2]
FortranForm[u2v2]
Print[Table[u2v2, {i,15}, {j,15}]]ということになる。

oyabun% math8

Mathematica 2.0 for SPARC
Copyright 1988-91 Wolfram Research, Inc.
-- X11 windows graphics initialized -In[1]:= << foo.m
1 + i
2
u=x
((2 + i) (3 + i) - 2 i (3 + i) x + i (1 + i) x )
1 + j
2
v=x
((2 + j) (3 + j) - 2 j (3 + j) x + j (1 + j) x )2
34
Bij=(8 (3780 + 3564 i + 1209 i + 174 i + 9 i + 3564 j + 2505 i j +
2
>

565 i

3
j + 40 i
3>

40 i j2
j + 1209 j

2+ 565 i j

2 2
3+ 65 i j + 174 j +4+ 9 j )) />((3 + i + j) (4 + i + j) (5 + i + j) (6 + i + j) (7 + i + j))
TeXForm={{8\,\left( 3780 + 3564\,i + 1209\,{i^2} + 174\,{i^3} + 9\,{i^4} +>3564\,j + 2505\,i\,j + 565\,{i^2}\,j + 40\,{i^3}\,j + 1209\,{j^2} +>565\,i\,{j^2} + 65\,{i^2}\,{j^2} + 174\,{j^3} + 40\,i\,{j^3} +>9\,{j^4} \right) }\over>{\left( 3 + i + j \right) \,\left( 4 + i + j \right) \,>\left( 5 + i + j \right) \,\left( 6 + i + j \right) \,>\left( 7 + i + j \right) }}FortranForm=8*(3780 + 3564*i + 1209*i**2 + 174*i**3 + 9*i**4 + 3564*j +>2505*i*j + 565*i**2*j + 40*i**3*j + 1209*j**2 + 565*i*j**2 +>65*i**2*j**2 + 174*j**3 + 40*i*j**3 + 9*j**4)/>((3 + i + j)*(4 + i + j)*(5 + i + j)*(6 + i + j)*(7 + i + j))2
3
2 -1 + i
Du=i (6 + 11 i + 6 i + i ) (-1 + x) x
2
3
2 -1 + j
9Dv=j (6 + 11 j + 6 j+ j ) (-1 + x) x
2
32
3
4 -2 + i + j
DuDv=i (6 + 11 i + 6 i + i ) j (6 + 11 j + 6 j + j ) (-1 + x) x
23
23
Aij=(24 i (6 + 11 i + 6 i + i ) j (6 + 11 j + 6 j + j )) /2>

(120 + 274 (-2 + i + j) + 225 (-2 + i + j)

3
+ 85 (-2 + i + j)

+

4

5
>
15 (-2 + i + j) + (-2 + i + j) )TeXForm={{24\,i\,\left( 6 + 11\,i + 6\,{i^2} + {i^3} \right) \,j\,>>\left( 6 + 11\,j + 6\,{j^2} + {j^3} \right) }\over{120 + 274\,\left( -2 + i + j \right)+>

225\,{{\left( -2 + i + j \right) }^2} +>85\,{{\left( -2 + i + j \right) }^3} +>15\,{{\left( -2 + i + j \right) }^4} +>{{\left( -2 + i + j \right) }^5}}}
FortranForm=24*i*(6 + 11*i + 6*i**2 + i**3)*j*(6 + 11*j + 6*j**2 + j**3)/>(120 + 274*(-2 + i + j) + 225*(-2 + i + j)**2 + 85*(-2 + i + j)**3 +>15*(-2 + i + j)**4 + (-2 + i + j)**5)
416 2608 8636 1123 13384
2608 5216 227 904 17540{{---, ----, ----, ----, -----}, {----, ----, ---, ---, -----},
45
315
1155 165
2145
315
693
33
143 3003>

8636 227 908 5890 5480
1123 904 5890 2356 401{----, ---, ---, ----, ----}, {----, ---, ----, ----, ---},
1155 33
143 1001 1001
165
143 1001 429
7813384 17540 5480 401 3208>{-----, -----, ----, ---, ----}}
2145
3003
1001 78663
576
576
960 1080{{---, 96, ---, 72, 64}, {96, ---, ----, 160, 160},
10

5

7

7

7

>

576 1080 1440
2880
3360 3920{---, ----, ----, 240, ----}, {72, 160, 240, ----, ----},
7
7
7
11
11
11

>

2880 3920 62720{64, 160, ----, ----, -----}}
11
11
143

In[2]:= uv
後は I 君に任せよう。Mathematica 万歳。

5.3

二日目の晩

以前からの懸案だった EISPACK のインストールにやっと成功して、IMSL の代わりに EISPACK でも計算できた。結果は IMSL 内のルーチンを利用した場合と少し食い違う。少なく
とも単精度計算の結果は IMSL のそれの方が良いみたい (倍精度計算の結果と照合しての判断)。ふうむ、さすが商品、か。中で反復法でもやっているのか。それとも、何か本質的な改良してあるのか。もっとも、どの N まで計算できるかという限界は、両者に差がなかった。
ところで、僕は次の問について考える必要があると気がついた。
すべての固有値が必要なのか、それとも小さい方の数個があれば十分なのか？
これは加藤先生の本を真面目に読む必要があるかなぁ？、ということになった。

5.4

三日目 – 差分法で挑戦 –昨日ゼミ卒研があって、久しぶりに菊地＆山本の本を読んでいて、常微分方程式の固有値問
題の章があるのを思い出した。そうだ、これをやってみようということで、実行。
まずちょっと走らせてみて、あまり精度が低いのでコーディング・ミスを疑ったが、普通のルーチンと、帯行列用のルーチンで同じ答が出る。単純なコーディング・ミスではなさそう。
それではアルゴリズムのミスか？。ところが、次数を上げていくと段々精度が上がって来てそれらしい結果を出すようになる。正しいのかな？この時点で、次数を大きくするためにも帯行列用のルーチンでないと駄目だと分かった。さらに固有値を全部計算するのはしんどいので、小さい方から数個求めるようにしないと話にならないことも分かった。というわけで、次のよ
うなプログラムを走らせることになった。

program fdm
integer i,n,m,maxn,N1,N2,NSTEP,MAXW
parameter (maxn = 10000, ncoda=2, lda=ncoda+1)
parameter (MAXW = 2*maxn*(ncoda+4)+maxn)
11

logical small
parameter (small = .true.)
real*8 a(lda,maxn), eval(maxn), c, h, h4inv
external devasb
COMMON /WORKSP/ RWKSPREAL RWKSP(MAXW)*CALL IWKIN(MAXW)cwrite(*,*) ’ NSTART,NEND,NSTEP=’
read(*,*) N1,N2,NSTEP
if (N1 .lt. 6) thenN1 = 6endifcdo n = N1,N2,NSTEPm = n - 3
h = 1.0d0 / n
h4inv = 1.0d0 / h**4ca(1,1)=0.0d0a(1,2)=0.0d0
c = 1.0d0 * h4inv
do i=3,ma(1,i) = c
end doca(2,1) = 0.0d0
c = -4.0d0 * h4inv
do i=2,m-1
a(2,i)= c
end doa(2,m)=-2.0 * h4inv
cc = 6.0d0 * h4inv
do i=1,m-2
a(3,i)= c
end doa(3,m-1)=5.0d0 * h4inv
a(3,m)=1.0d0 * h4inv
ccall devasb(m, 5, a, lda, ncoda, small, eval)
12cwrite(*,*) ’ N=’, Nwrite(*,*) (eval(i),i=1,5)
end doend
さて、結果はどうなるかというと、、、

oyabun% test-FDM2
NSTART,NEND,NSTEP=
1000 10000 1000
N= 1000
12.462191808892
40264.212562460
N= 2000
12.395044092660
40103.824563559
N= 3000
13.134802409501
40050.098696306
N= 4000
11.946737822504
40023.221212918
N= 5000
19.102980867991
40013.678099232
N= 6000
28.244970754340
40018.703056702
N= 7000-3.6834530007044
39988.834453875
N= 8000-63.982736504047
39945.408837434
N= 9000
84.591125935711
39998.978912310
N= 10000-52.363946332173
40010.437030957
oyabun%489.42847340945

3837.1448241929

14734.664365594

487.34666340961

3821.7483886199

14675.816251591

486.94196351073

3816.7670363691

14656.253512514

486.33007299197

3813.4893156216

14645.588853085

494.31233271986

3818.8659709859

14646.849156585

503.45967693861

3832.3625652045

14659.027784007

492.44054690281

3820.1818452542

14648.801176336

408.39409769491

3740.3961247383

14541.477295313

508.35975169377

3837.5120895015

14655.358069608

418.17647792973

3803.9786566901

14652.305272351

うーん。Galerkin 法に惨敗している。N = 2000 程度で頭打ちみたいで、その際の結果が悪い。

13

5.5

再び Mathematica

この時点の状況を見ると Galerkin 法の優秀性が際だっている。差分法で大きな計算をしても Galerkin 法の小さな計算にかなわないのならば、Galerkin 法で、N = 10, 20 程度の計算をして満足するべきなのではないか？(最初のうちは暗に、Galerkin 法で N = 100, 200 というような計算をする気でいたわけだ。だから Jacobi 法の採用や、非対称行列の固有値問題への帰着を蹴って来たわけだね。)
そういうわけで、全部 Mathematica で計算する方法を考えることになった。Mathematicaに一般化固有値問題を解く組み込み手続きはないようだが、問題が小さければ非対称の問題になっても構わないのでないか？そこで一度は捨てた方法、すなわち B −1 A の固有値の計算をしてみることにした。プログラムは% ishibashi3.m
u=x^(i+1)((i+1)i x^2-2i(i+3)x+(i+2)(i+3))v=x^(j+1)((j+1)j x^2-2j(j+3)x+(j+2)(j+3))
Print["u=",u]Print["v=",v]
Bij=Simplify[Integrate[u v, {x,0,1}]]
Print["Bij=",Bij]
DDu=Simplify[D[u,{x,2}]]
DDv=Simplify[D[v,{x,2}]]
DDuDDv = Simplify[DDu DDv]
Print["DDuDDv=",DDuDDv]
Aij=Simplify[Integrate[DDuDDv /. i+j-2->m, {x,0,1}] /. m->i+j-2]
Print["Aij=", Aij]
B5=Table[Bij, {i,5}, {j,5}]A5=Table[Aij, {i,5}, {j,5}]X5=N[Inverse[B5] . A5, 40]
Print["(N=5) Eigenvalues of B^(-1) * A =", Sort[Eigenvalues[X5]]]
B20=Table[Bij, {i,20}, {j,20}]
A20=Table[Aij, {i,20}, {j,20}]
X20=N[Inverse[B20] . A20, 40]
Print["(N=20) Eigenvalues of B^(-1) * A =", Sort[Eigenvalues[X20]]]
Remove[u,v,DDu,DDv,DDuDDv,A5,B5,X5,A20,B20,X20]というものである。固有値は(N=5) Eigenvalues of inverse B^(-1) * A =>{12.36236336840952940231787585067154160033,>

485.5189744083006116625558974720579893284,>

3811.768990034377066087752823718141246277,
14

>14784.4917006380991631897111399552268624,>

57102.4650674393490172220221233488604026}(N=20) Eigenvalues of B^(-1) * A =>{12.36236336832619021871926166118869013603,>

485.518818513371037811691383834744388645,>

3806.546266391451058088482228842380210229,>

14617.27330511878066373978433933587907848,>

39943.8317785094667476023381845908568265,>

89135.405071423787031412816497573328912,>

173881.315656156365530902969367235780677,>

308208.4524575434976297314325665570192543,>

508481.5496567107232917183383377889822887,>

793407.140343260324224285231205468035843,
6>

1.184039968145607930407722269718423310352 10 ,6>

1.70619097696839393429262456204757437922 10 ,6>

2.386493083493206213818601932268233908251 10 ,6>

3.392589029623894784136767942540982255932 10 ,6>

4.63729964903979336008569892106668076025 10 ,

15

6
>

7.85416773748417801457297775449574357188 10 ,7>

1.104159009168505078752128124330832849583 10 ,7>

2.761431598064332497648305858001351817508 10 ,7>4.34651133206433406410291335571975185642 10 ,8>2.080831345259133798313121604931666684158 10 }In[2]:= Quitとなります。まあ、ここらへんが満足のしどころかな。

6

戦い済んで
やはり不勉強を痛感する。とりあえず、いくつか身につけることの出来たものがある。

• Mathematica の使い方。使えるという自信。
• IMSL の使用経験。
• 使える状態の EISPACK。EISPACK の使用経験。
• 常微分方程式の固有値問題の数値解法についての知識。
• 素朴な Galerkin 法の経験。
• 固有値問題への知見。必要な文献等。
IMSL は使い続けることにしよう (バージョンアップを続けるようにお願いしよう)。そして
EISPACK や NUMPAC 等も整備しておくことにする (NUMPAC 入手については動き出さな
いといけないな)。加藤先生の変分法は一度きちんと読む (必要—価値) があるだろう。数値計

16

算に関する基本的な文献は早めに揃えておくことにしよう。また、学生用の参考書も複数準備したい。
最後に自室でなく、計算機室備えつけの書籍には、一般化固有値問題の解説もちゃんと載っていることが明らかになった。
（ああ、恥ずかし。不明を恥じる。）特に渡辺＆名取＆小国の本
はプログラム付きだった。これもちゃんとインストールしておく必要があるな。

A

差分方程式

∆ を中心差分演算子とする。vj := ∆2 uj とする。すなわち
vj = uj−1 − 2uj + uj+1 .
∆4 uj = ∆2 vj = vj−1 − 2vj + vj+1
= uj−2 − 2uj−1 + uj − 2(uj−1 − 2uj + uj+1 ) + (uj − 2uj+1 + uj+2 )= uj−2 − 4uj−1 + 6uj − 4uj+1 + uj+2

(2 ≤ j ≤ N − 2).

ゆえにuj−2 − 4uj−1 + 6uj − 4uj+1 + uj+2 = λuj(1)(2 ≤ j ≤ N − 2).

一方で u(0) = u′ (0) = 0 より u0 = 0, u1 − u0 = 0 として、

(2)

u0 = u1 = 0.

また u′′ (0) = u′′′ (0) = 0 より

uN −2 − 2uN −1 + uN = 0,

(uN −3 − 2uN −2 + uN −1 ) − (uN −2 − 2uN −1 + uN ) = 0.

これから

uN −2 − 2uN −1 + uN = 0,

uN −3 − 2uN −2 + uN −1 = 0.

整理して

(3) uN −1 = 2uN −2 −uN −3 ,uN = 2uN −1 −uN −2 = 2(2uN −2 −uN −3 )−uN −2 = 3uN −2 −2uN −3 .

よって最後から 2 行目はuN −5 − 4uN −4 + 6uN −3 − 4uN −2 + uN −1 = uN −5 − 4uN −4 + 6uN −3 − 4uN −2 + (2uN −2 − uN −3 )= uN −5 − 4uN −4 + 5uN −3 − 2uN −2 .
最後の行はuN −4 − 4uN −3 + 6uN −2 − 4uN −1 + uN = uN −4 − 4uN −3 + 6uN −2 − 4(2uN −2 − uN −3 ) + (3uN −2 − 2uN −3 )= uN −4 − 2uN −3 + uN −2 .

17

よって





u2
6u2 −4u3 +u4 u3 

−4u2 +6u3 −4u4 +u5
 u4 

 u2
−4u
+6u
−4u
+u
3
4
5
6




 u5 
u3−4u4 +6u5 −4u6
+u7
1 




= λ . .


.
.
.
.
.
4
..
..
..
..
..
 .. 

h 





u

uN −6 −4uN −5 +6uN −4 −4uN −3 +uN −2 
 N −4 






uN −5
−4uN −4 +5uN −3 −2uN −2 uN −3 

uN −2
uN −4
−2uN −3 +uN −2

これは N − 3 次実対称行列の固有値問題である。

18