Keyword: 双３次, スプライン, フィット
概要
本サンプルは矩形格子データへの双３次スプラインフィットを行うFortranによるサンプルプログラムです。 本サンプルは以下に示されるデータについて双３次スプラインフィットを行いノットとBスプライン係数を求めて出力します。
※本サンプルはnAG Fortranライブラリに含まれるルーチン e02dcf() のExampleコードです。本サンプル及びルーチンの詳細情報は e02dcf のマニュアルページをご参照ください。
ご相談やお問い合わせはこちらまで
入力データ
（本ルーチンの詳細は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