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)
      logical small
      parameter (small = .true.)
      real*8 a(lda,maxn), eval(maxn), c, h, h4inv
      external devasb
      COMMON /WORKSP/  RWKSP
      REAL RWKSP(MAXW)
*
      CALL IWKIN(MAXW)
c
      write(*,*) ' NSTART,NEND,NSTEP='
      read(*,*) N1,N2,NSTEP
      if (N1 .lt. 6) then
        N1 = 6
      endif
c
      do n = N1,N2,NSTEP
        m = n - 3
        h = 1.0d0 / n
        h4inv = 1.0d0 / h**4
c
        a(1,1)=0.0d0
        a(1,2)=0.0d0
        c = 1.0d0 * h4inv
        do i=3,m
          a(1,i) = c
        end do
c
        a(2,1) = 0.0d0
        c = -4.0d0 * h4inv
        do i=2,m-1
          a(2,i)= c
        end do
        a(2,m)=-2.0 * h4inv
c
        c = 6.0d0 * h4inv
        do i=1,m-2
          a(3,i)= c
        end do
        a(3,m-1)=5.0d0 * h4inv
        a(3,m)=1.0d0 * h4inv
c
        call devasb(m, 5, a, lda, ncoda, small, eval)
c
        write(*,*) ' N=', N
        write(*,*) (eval(i),i=1,5)
      end do
      end

さて、結果はどうなるかというと、、、

oyabun%  test-FDM2
 NSTART,NEND,NSTEP=
1000 10000 1000
 N=  1000
    12.462191808892    489.42847340945    3837.1448241929    14734.664365594
    40264.212562460
 N=  2000
    12.395044092660    487.34666340961    3821.7483886199    14675.816251591
    40103.824563559
 N=  3000
    13.134802409501    486.94196351073    3816.7670363691    14656.253512514
    40050.098696306
 N=  4000
    11.946737822504    486.33007299197    3813.4893156216    14645.588853085
    40023.221212918
 N=  5000
    19.102980867991    494.31233271986    3818.8659709859    14646.849156585
    40013.678099232
 N=  6000
    28.244970754340    503.45967693861    3832.3625652045    14659.027784007
    40018.703056702
 N=  7000
   -3.6834530007044    492.44054690281    3820.1818452542    14648.801176336
    39988.834453875
 N=  8000
   -63.982736504047    408.39409769491    3740.3961247383    14541.477295313
    39945.408837434
 N=  9000
    84.591125935711    508.35975169377    3837.5120895015    14655.358069608
    39998.978912310
 N=  10000
   -52.363946332173    418.17647792973    3803.9786566901    14652.305272351
    40010.437030957
oyabun%

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