第６回、サブルーチンの利用

目次

• 小テスト
• サブルーチンの利用
• 今日の課題
• レポート提出用form

サブルーチンの利用

最小２乗法

例題１

n組のデータ(x_i,y_i)(i=1,2,...,n)に対して、 直線 y=a+bx を最小2乗法であてはめる。

解説

n 組の観測値 (x_i,y_i) (i=1,2,...,n) があった ときに、それらに対して直線 y = a + bx をあてはめることを考える。 y_i は誤差を含んでいるが、 x_i は誤差を含んでいないものとする。最小2乗法は残差 e_i = y_i - ( a + b x_i ) の２乗和 Σe_i² が最小になるように a,b の値を 決める方法である。x,y の平均値、x の平方和、x,yの積和
< x > = 1/n Σx_i , < y > = 1/n Σy_i , S_x = Σ(x_i - < x >)² , S_xy = (x_i - < x >)(y_i - < y>)
を使うと
b = S_xy / S_x , a = < y > - b < x >
となる。

解答

      implicit real*8(a-h,o-z)
      real*8 a, b
      integer num, i, n
      parameter(num=100)
      real*8 datax(num), datay(num)
c
      read(5,'(i5)') n
      read(5,'(2f10.5)') (datax(i), datay(i), i=1,n)
c
c     subroutine call *** LSQ ***
c       (least square method)
c
      call lsq(datax, datay, n, a, b)
c
      write(6,'(a,f10.5,a,f10.5)') 'A = ',a, ' B = ',b
c
      end
c================================================================
      subroutine lsq(x, y, n, a, b)
      implicit real*8(a-h,o-z)
      real*8 x(*), y(*), a, b, sx, sy, xm, ym, sxx, sxy
      integer i
c
      sx = 0.0d0
      sy = 0.0d0
      do i = 1 , n
         sx = sx + x(i)
         sy = sy + y(i)
      end do
c
      xm = sx / n
      ym = sy / n
c
      sxx = 0.0d0
      sxy = 0.0d0
      do i = 1 , n
         sxx = sxx + (x(i)-xm)**2
         sxy = sxy + (x(i)-xm)*(y(i)-ym)
      end do
c
      b = sxy / sxx
      a = ym - b*xm
      return
      end

１６行目

主プログラムの中には特に難しいところはないと思います。しいてあげれば この行で、サブルーチン"lsq"を呼び出します。その後に"( )"で囲まれている ものはサブルーチンに必要な引数です。気をつけないといけないのは、変数の 型、順番をサブルーチンのものと一致させておかなければいけない、という ことです。そうでなければエラーとなります。変数名は変えてもかまいません。 このプログラムの場合では

• datax < -- > x 横軸のデータ
• datay < -- > y 縦軸のデータ
• n < -- > n データ総数
• a < -- > a 切片
• b < -- > b 傾き

となります。

サブルーチン中３行目

"real*8 x(*), y(*)" としているところがみなれないものだと思います。 サブルーチンの引数として配列をとった場合、その変数の数はサブルーチンを 呼び出す側で定義されているのが普通です。（定義されていないといけません） つまりサブルーチン中の"x, y"はどのくらいの大きさの配列かはサブルーチンの 中ではわかりません。ですから
real*8 x(*), y(*),....
として宣言します。ただし教科書によっては
real*8 x(n), y(n),....
とされているかもしれません。サブルーチンの引数に配列の大きさを引き渡して いるのですからこれでもよいかもしれません。ただ私はこういう書き方はきらい です。

サブルーチンの残り

必要なデータを計算していきます。

連立１次方程式（はき出し法）

例題２

はきだし法を用いて n元連立1次方程式を解く。

解説

n 行 (n+1) 列の係数行列を考える。はきだし法は係数行列の対角成分だけ を残し、他の要素を順次消去してゆく 方法で、以下の手順を k=1,2,...,n について繰り返すことによって、その解は 第 (n+1) 列に求められる。

1. a_kj < -- a_kj/a_kk (j=k+1,...,n+1)
2. a_ij < -- a_ij - a_ik a_kj (i=1,2,...,k-1,k+1,..,n ; j=k+1,...,n+1)

解答

      implicit real*8(a-h,o-z)
      integer m
      parameter (m=3)
      real*8 a(m,m+1), x(m)
c
      do i = 1 , m
         read(5,'(4f10.5)') (a(i,j),j=1,m+1)
      end do
c
      call gauss(a,x,m)
c
      write(6,'(a,i3,a,f10.5)') ('X(',i,') = ',x(i), i=1,m)
      end
c=============================================================
      subroutine gauss(a,x,n)
      implicit real*8(a-h,o-z)
      integer n
      real*8 a(n,*), x(*)
      integer i, j, k
c
      do k = 1 , n
         do j = k+1, n+1
            a(k,j) = a(k,j)/a(k,k)
         end do
c
         do i = 1 , n
            if( i .ne. k ) then
               do j = k+1 , n+1
                  a(i,j) = a(i,j) - a(i,k)*a(k,j)
               end do
            end if
         end do
      end do
c
      do k = 1 , n
         x(k) = a(k,n+1)
      end do
c
      return
      end

今日の課題

1. 小テスト
2. [6-1]安息香酸の水-ベンゼン系での分配係数を測定する実験において、 以下のような結果が得られた。これより会合度 n と分配係数 K を求めよ。 ただし K,n の間の関係式は、水中、ベンゼン中の安息香酸の濃度を それぞれC_a,C_bとすると
   C_aⁿ/C_b = K である。
   

   Ca/N	0.0121	0.0132	0.0179	0.0202	0.0215	0.0229
   Cb/N	0.233	0.288	0.364	0.525	0.588	0.639

3. [6-2]30 度におけるショ糖の転化速度を測定する実験において、以下のような 結果が得られた。これより速度定数 k(sec^-1) を求めよ。 ただし k,時間 t, 旋光度αの間の関係式は
   k t = ln(α₀-α_∞) - ln(α_t-α_∞) である。(α_∞=-0.678)
   

   t/min.	0	10	20	30	40	50	60
   αt/deg.	2.316	2.248	2.167	2.078	2.000	1.919	1.842

4. [6-3]上の例題２の連立1次方程式を解くためのサブルーチンには重要な欠陥が ある。それを指摘し、対処するサブルーチンに変更せよ。 (ヒント：”ピボットの選択”がキーワードである)
5. [6-4]例題１の発展問題である。 m 個のデータ(x_i,y_i)を x についての n 次式
   y = a + bx + ... + cxⁿ へあてはめるプログラムを作成せよ。(一般の n 次式に対応した プログラムがのぞましいが、難しいようなら 2 次式と限定してもよい)
6. [6-5]例題２を応用して、与えられた行列 A の逆行列 A^-1を求める プログラムを作成せよ。