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,

>    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 ,

                                                6
>    7.85416773748417801457297775449574357188 10 ,

                                                 7
>    1.104159009168505078752128124330832849583 10 ,

                                                 7
>    2.761431598064332497648305858001351817508 10 ,

                                                7
>    4.34651133206433406410291335571975185642 10 ,

                                                 8
>    2.080831345259133798313121604931666684158 10 }

In[2]:= Quit

となります。まあ、ここらへんが満足のしどころかな。