矩形格子データへの双３次スプラインフィット

ホーム > 製品 > nAG数値計算ライブラリ > nAG Fortranライブラリ > サンプルソースコード集 > 矩形格子データへの双３次スプラインフィット

Keyword: 双３次, スプライン, フィット

概要

本サンプルは矩形格子データへの双３次スプラインフィットを行うFortranによるサンプルプログラムです。本サンプルは以下に示されるデータについて双３次スプラインフィットを行いノットとBスプライン係数を求めて出力します。

双３次スプラインフィットのデータ　

※本サンプルはnAG Fortranライブラリに含まれるルーチン e02dcf() のExampleコードです。本サンプル及びルーチンの詳細情報は e02dcf のマニュアルページをご参照ください。
ご相談やお問い合わせはこちらまで

入力データ

このデータをダウンロード

E02DCF Example Program Data
 11    9      MX, MY, number of grid points on the X and Y axes
 0.0000E+00  5.0000E-01  1.0000E+00  1.5000E+00  2.0000E+00
 2.5000E+00  3.0000E+00  3.5000E+00  4.0000E+00  4.5000E+00
 5.0000E+00   End of MX grid points
 0.0000E+00  5.0000E-01  1.0000E+00  1.5000E+00  2.0000E+00
 2.5000E+00  3.0000E+00  3.5000E+00  4.0000E+00   End of MY grid points
 1.0000E+00  8.8758E-01  5.4030E-01  7.0737E-02 -4.1515E-01
-8.0114E-01 -9.7999E-01 -9.3446E-01 -6.5664E-01  1.5000E+00
 1.3564E+00  8.2045E-01  1.0611E-01 -6.2422E-01 -1.2317E+00
-1.4850E+00 -1.3047E+00 -9.8547E-01  2.0600E+00  1.7552E+00
 1.0806E+00  1.5147E-01 -8.3229E-01 -1.6023E+00 -1.9700E+00
-1.8729E+00 -1.4073E+00  2.5700E+00  2.1240E+00  1.3508E+00
 1.7684E-01 -1.0404E+00 -2.0029E+00 -2.4750E+00 -2.3511E+00
-1.6741E+00  3.0000E+00  2.6427E+00  1.6309E+00  2.1221E-01
-1.2484E+00 -2.2034E+00 -2.9700E+00 -2.8094E+00 -1.9809E+00
 3.5000E+00  3.1715E+00  1.8611E+00  2.4458E-01 -1.4565E+00
-2.8640E+00 -3.2650E+00 -3.2776E+00 -2.2878E+00  4.0400E+00
 3.5103E+00  2.0612E+00  2.8595E-01 -1.6946E+00 -3.2046E+00
-3.9600E+00 -3.7958E+00 -2.6146E+00  4.5000E+00  3.9391E+00
 2.4314E+00  3.1632E-01 -1.8627E+00 -3.6351E+00 -4.4550E+00
-4.2141E+00 -2.9314E+00  5.0400E+00  4.3879E+00  2.7515E+00
 3.5369E-01 -2.0707E+00 -4.0057E+00 -4.9700E+00 -4.6823E+00
-3.2382E+00  5.5050E+00  4.8367E+00  2.9717E+00  3.8505E-01
-2.2888E+00 -4.4033E+00 -5.4450E+00 -5.1405E+00 -3.5950E+00
 6.0000E+00  5.2755E+00  3.2418E+00  4.2442E-01 -2.4769E+00
-4.8169E+00 -5.9300E+00 -5.6387E+00 -3.9319E+00   End of data values
 0.1         S, smoothing factor
 6    0.0   5.0
 5    0.0   4.0

１行目はタイトル行で読み飛ばされます。
２行目にx軸の格子点の数(mx)、y軸の格子点の数(my)を指定しています。
３～５行目にはx座標を指定しています。
６～７行目にはy座標を指定しています。
８～２７行目にはデータ値(f)を指定しています。
２８行目に平滑化因子(s)を指定しています。
２９行目に矩形格子のx軸に沿った点の数(npx)、最小値(xlo)、最大値(xhi)を指定しています。
３０行目に矩形格子のy軸に沿った点の数(npy)、最小値(ylo)、最大値(yhi)を指定しています。

出力結果

（本ルーチンの詳細はe02dcf のマニュアルページを参照）

この出力例をダウンロード

 E02DCF Example Program Results

 Calling with smoothing factor S =   1.0000E-01: NX =   10, NY =   13.

                I   Knot LAMDA(I)      J     Knot MU(J)

                4      0.0000          4      0.0000
                5      1.5000          5      1.0000
                6      2.5000          6      2.0000
                7      5.0000          7      2.5000
                                       8      3.0000
                                       9      3.5000
                                      10      4.0000

 The B-spline coefficients:

    0.9918   1.5381   2.3913   3.9845   5.2138   5.9965
    1.0546   1.5270   2.2441   4.2217   5.0860   6.0821
    0.6098   0.9557   1.5587   2.3458   3.3860   3.7716
   -0.2915  -0.4199  -0.7399  -1.1763  -1.5527  -1.7775
   -0.8476  -1.3296  -1.8521  -3.3468  -4.3628  -5.0085
   -1.0168  -1.5952  -2.4022  -3.9390  -5.4680  -6.1656
   -0.9529  -1.3381  -2.2844  -3.9559  -5.0032  -5.8709
   -0.7711  -1.0914  -1.8488  -3.2549  -3.9444  -4.7297
   -0.6476  -1.0373  -1.5936  -2.5887  -3.3485  -3.9330

 Sum of squared residuals FP =   1.0004E-01

 Values of computed spline:

           X    0.00    1.00    2.00    3.00    4.00    5.00
      Y
     4.00      -0.65   -1.36   -1.99   -2.61   -3.25   -3.93
     3.00      -0.98   -1.97   -2.91   -3.91   -4.97   -5.92
     2.00      -0.42   -0.83   -1.24   -1.66   -2.08   -2.48
     1.00       0.54    1.09    1.61    2.14    2.71    3.24
     0.00       0.99    2.04    3.03    4.01    5.02    6.00

３行目に呼び出される際の平滑化因子、NX（変数xに関するスプラインのノットの数）、NY（変数yに関するスプラインのノットの数）が出力されています。
５～１３行目にx方向のノットの位置（LAMDA）とy方向のノットの位置（MU）が出力されています。
１７～２５行目にBスプライン係数が出力されています。
２７行目に残差平方和が出力されています。
３１～３７行目に計算されたスプラインの値が出力されています。

ソースコード

（本ルーチンの詳細はe02dcf のマニュアルページを参照）

※本サンプルソースコードは科学技術・統計計算ライブラリである「nAG Fortranライブラリ」のルーチンを呼び出します。
サンプルのコンパイル及び実行方法

このソースコードをダウンロード

!   E02DCF Example Program Text
!   Mark 23 Release. nAG Copyright 2011.
    MODULE e02dcfe_mod

!      E02DCF Example Program Module:
!             Parameters and User-defined Routines

!      .. Use Statements ..
       USE nag_library, ONLY : nag_wp
!      .. Implicit None Statement ..
       IMPLICIT NONE
!      .. Parameters ..
       INTEGER, PARAMETER                  :: nin = 5, nout = 6
    CONTAINS
       SUBROUTINE cprint(c,ny,nx,nout)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. Scalar Arguments ..
          INTEGER, INTENT (IN)                :: nout, nx, ny
!         .. Array Arguments ..
          REAL (KIND=nag_wp), INTENT (IN)     :: c(ny-4,nx-4)
!         .. Local Scalars ..
          INTEGER                             :: i
!         .. Executable Statements ..
          WRITE (nout,*)
          WRITE (nout,*) 'The B-spline coefficients:'
          WRITE (nout,*)

          DO i = 1, ny - 4
             WRITE (nout,99999) c(i,1:(nx-4))
          END DO

          RETURN

99999     FORMAT (1X,8F9.4)
       END SUBROUTINE cprint
    END MODULE e02dcfe_mod
    PROGRAM e02dcfe

!      E02DCF Example Main Program

!      .. Use Statements ..
       USE nag_library, ONLY : e02dcf, e02dff, nag_wp
       USE e02dcfe_mod, ONLY : cprint, nin, nout
!      .. Implicit None Statement ..
       IMPLICIT NONE
!      .. Local Scalars ..
       REAL (KIND=nag_wp)                  :: delta, fp, s, xhi, xlo, yhi, ylo
       INTEGER                             :: i, ifail, j, liwrk, lwrk, mx,    &
                                              my, npx, npy, nx, nxest, ny, nyest
       CHARACTER (1)                       :: start
!      .. Local Arrays ..
       REAL (KIND=nag_wp), ALLOCATABLE     :: c(:), f(:), fg(:), lamda(:),     &
                                              mu(:), px(:), py(:), wrk(:),     &
                                              x(:), y(:)
       INTEGER, ALLOCATABLE                :: iwrk(:)
!      .. Intrinsic Functions ..
       INTRINSIC                              max, min, real
!      .. Executable Statements ..
       WRITE (nout,*) 'E02DCF Example Program Results'

!      Skip heading in data file
       READ (nin,*)

!      Input the number of X, Y co-ordinates MX, MY.

       READ (nin,*) mx, my
       nxest = mx + 4
       nyest = my + 4
       lwrk = 4*(mx+my) + 11*(nxest+nyest) + nxest*my + max(my,nxest) + 54
       liwrk = 3 + mx + my + nxest + nyest
       ALLOCATE (x(mx),y(my),f(mx*my),lamda(nxest),mu(nyest),wrk(lwrk), &
          iwrk(liwrk),c((nxest-4)*(nyest-4)))

       READ (nin,*) x(1:mx)
       READ (nin,*) y(1:my)

!      Input the MX*MY function values F at the grid points.

       READ (nin,*) f(1:mx*my)

       start = 'C'

       READ (nin,*) s

!      Determine the spline approximation.

       ifail = 0
       CALL e02dcf(start,mx,x,my,y,f,s,nxest,nyest,nx,lamda,ny,mu,c,fp,wrk, &
          lwrk,iwrk,liwrk,ifail)

       DEALLOCATE (wrk,iwrk)

       WRITE (nout,*)
       WRITE (nout,99999) 'Calling with smoothing factor S =', s, ': NX =', &
          nx, ', NY =', ny, '.'
       WRITE (nout,*)
       WRITE (nout,*) '               I   Knot LAMDA(I)      J     Knot MU(J)'
       WRITE (nout,*)

       DO j = 4, max(nx,ny) - 3

          IF (j<=nx-3 .AND. j<=ny-3) THEN
             WRITE (nout,99997) j, lamda(j), j, mu(j)
          ELSE IF (j<=nx-3) THEN
             WRITE (nout,99997) j, lamda(j)
          ELSE IF (j<=ny-3) THEN
             WRITE (nout,99996) j, mu(j)
          END IF

       END DO

       CALL cprint(c,ny,nx,nout)

       WRITE (nout,*)
       WRITE (nout,99998) 'Sum of squared residuals FP =', fp

       IF (fp==0.0E+0_nag_wp) THEN
          WRITE (nout,*) '(The spline is an interpolating spline)'
       ELSE IF (nx==8 .AND. ny==8) THEN
          WRITE (nout,*) &
             '(The spline is the least-squares bi-cubic polynomial)'
       END IF

!      Evaluate the spline on a rectangular grid at NPX*NPY points
!      over the domain (XLO to XHI) x (YLO to YHI).

       READ (nin,*) npx, xlo, xhi
       READ (nin,*) npy, ylo, yhi

       lwrk = min(4*npx+nx,4*npy+ny)

       IF (4*npx+nx>4*npy+ny) THEN
          liwrk = npy + ny - 4
       ELSE
          liwrk = npx + nx - 4
       END IF

       ALLOCATE (px(npx),py(npy),fg(npx*npy),wrk(lwrk),iwrk(liwrk))

       delta = (xhi-xlo)/real(npx-1,kind=nag_wp)

       DO i = 1, npx
          px(i) = min(xlo+real(i-1,kind=nag_wp)*delta,xhi)
       END DO

       delta = (yhi-ylo)/real(npy-1,kind=nag_wp)

       DO i = 1, npy
          py(i) = min(ylo+real(i-1,kind=nag_wp)*delta,yhi)
       END DO

       ifail = 0
       CALL e02dff(npx,npy,nx,ny,px,py,lamda,mu,c,fg,wrk,lwrk,iwrk,liwrk, &
          ifail)

       WRITE (nout,*)
       WRITE (nout,*) 'Values of computed spline:'
       WRITE (nout,*)
       WRITE (nout,99995) '          X', (px(i),i=1,npx)
       WRITE (nout,*) '     Y'

       DO i = npy, 1, -1
          WRITE (nout,99994) py(i), (fg(npy*(j-1)+i),j=1,npx)
       END DO

99999  FORMAT (1X,A,1P,E13.4,A,I5,A,I5,A)
99998  FORMAT (1X,A,1P,E13.4)
99997  FORMAT (1X,I16,F12.4,I11,F12.4)
99996  FORMAT (1X,I39,F12.4)
99995  FORMAT (1X,A,7F8.2)
99994  FORMAT (1X,F8.2,3X,7F8.2)
    END PROGRAM e02dcfe

矩形格子データへの双３次スプラインフィット

Fortranによるサンプルソースコード : 使用ルーチン名：e02dcf

概要

入力データ

出力結果

ソースコード