Keyword: 双３次, スプライン, フィット
概要
本サンプルは矩形格子データへの双３次スプラインフィットを行うC言語によるサンプルプログラムです。 本サンプルは以下に示されるデータについて双３次スプラインフィットを行いノットとBスプライン係数を求めて出力します。
※本サンプルはnAG Cライブラリに含まれる関数 nag_2d_spline_fit_grid() のExampleコードです。本サンプル及び関数の詳細情報は nag_2d_spline_fit_grid のマニュアルページをご参照ください。
ご相談やお問い合わせはこちらまで
入力データ
（本関数の詳細はnag_2d_spline_fit_grid のマニュアルページを参照）| このデータをダウンロード |
nag_2d_spline_fit_grid (e02dcc) Example Program Data 11 9 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 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 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 0.1 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)を指定しています。
出力結果
（本関数の詳細はnag_2d_spline_fit_grid のマニュアルページを参照）| この出力例をダウンロード |
nag_2d_spline_fit_grid (e02dcc) Example Program Results
Calling with smoothing factor s = 1.0000e-01: spline.nx = 10, spline.ny = 13.
Distinct knots in x direction located at
0.0000 1.5000 2.5000 5.0000
Distinct knots in y direction located at
0.0000 1.0000 2.0000 2.5000 3.0000
3.5000 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
- ３行目に呼び出される際の平滑化因子、spline.nx（変数xに関するスプラインのノットの数）、spline.ny（変数yに関するスプラインのノットの数）が出力されています。
- ５〜６行目にx方向のノットの位置が出力されています。
- ８〜１０行目にy方向のノットの位置が出力されています。
- １２〜２２行目にBスプライン係数が出力されています。
- ２４行目に残差平方和が出力されています。
- ２６〜３４行目に計算されたスプラインの値が出力されています。
ソースコード
（本関数の詳細はnag_2d_spline_fit_grid のマニュアルページを参照）
※本サンプルソースコードはnAG数値計算ライブラリ（Windows, Linux, MAC等に対応）の関数を呼び出します。
サンプルのコンパイル及び実行方法
| このソースコードをダウンロード |
/* nag_2d_spline_fit_grid (e02dcc) Example Program.
*
* CLL6I261D/CLL6I261DL Version.
*
* Copyright 2017 Numerical Algorithms Group.
*
* Mark 26.1, 2017.
*
*/
#include <nag.h>
#include <stdio.h>
#include <nag_stdlib.h>
#include <nage02.h>
int main(void)
{
Integer exit_status = 0, i, j, mx, my, npx, npy, nx, ny;
Nag_2dSpline spline;
Nag_Comm warmstartinf;
Nag_Start start;
double delta, *f = 0, *fg = 0, fp, *px = 0, *py = 0, s, *x = 0, xhi,
xlo, *y = 0, yhi;
double ylo;
NagError fail;
INIT_FAIL(fail);
/* Initialize spline */
spline.lamda = 0;
spline.mu = 0;
spline.c = 0;
warmstartinf.nag_w = 0;
warmstartinf.nag_iw = 0;
printf("nag_2d_spline_fit_grid (e02dcc) Example Program Results\n");
scanf("%*[^\n]"); /* Skip heading in data file */
/* Input the number of x, y co-ordinates mx, my. */
scanf("%ld%ld", &mx, &my);
if (mx >= 4 && my >= 4) {
if (!(f = nAG_ALLOC(mx * my, double)) ||
!(x = nAG_ALLOC(mx, double)) || !(y = nAG_ALLOC(my, double)))
{
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
}
else {
printf("Invalid mx or my.\n");
exit_status = 1;
return exit_status;
}
/* Input the x co-ordinates followed by the y co-ordinates. */
for (i = 0; i < mx; i++)
scanf("%lf", &x[i]);
for (i = 0; i < my; i++)
scanf("%lf", &y[i]);
/* Input the mx*my function values f at the grid points. */
for (i = 0; i < mx * my; i++)
scanf("%lf", &f[i]);
start = Nag_Cold;
scanf("%lf", &s);
/* Determine the spline approximation. */
/* nag_2d_spline_fit_grid (e02dcc).
* Least squares bicubic spline fit with automatic knot
* placement, two variables (rectangular grid)
*/
nag_2d_spline_fit_grid(start, mx, x, my, y, f, s, mx + 4, my + 4,
&fp, &warmstartinf, &spline, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_2d_spline_fit_grid (e02dcc).\n%s\n", fail.message);
exit_status = 1;
goto END;
}
nx = spline.nx;
ny = spline.ny;
printf("\nCalling with smoothing factor s = %13.4e:"
" spline.nx = %2ld, spline.ny = %2ld.\n\n",
s, nx, ny);
/* Print the knot sets, lamda and mu. */
printf("Distinct knots in x direction located at\n");
for (j = 3; j < spline.nx - 3; j++)
printf("%12.4f%s", spline.lamda[j],
((j - 3) % 5 == 4 || j == spline.nx - 4) ? "\n" : " ");
printf("\nDistinct knots in y direction located at\n");
for (j = 3; j < spline.ny - 3; j++)
printf("%12.4f%s", spline.mu[j], ((j - 3) % 5 == 4 || j == spline.ny - 4)
? "\n" : " ");
printf("\nThe B-spline coefficients:\n\n");
for (i = 0; i < ny - 4; i++) {
for (j = 0; j < nx - 4; j++)
printf("%9.4f", spline.c[i + j * (ny - 4)]);
printf("\n");
}
printf("\nSum of squared residuals fp = %13.4e\n", fp);
if (fp == 0.0)
printf("\nThe spline is an interpolating spline\n");
else if (nx == 8 && ny == 8)
printf("\nThe spline is the least squares bi-cubic polynomial\n");
/* Evaluate the spline on a rectangular grid at npx*npy points
* over the domain (xlo to xhi) x (ylo to yhi).
*/
scanf("%ld%lf%lf", &npx, &xlo, &xhi);
scanf("%ld%lf%lf", &npy, &ylo, &yhi);
if (npx >= 1 && npy >= 1) {
if (!(fg = nAG_ALLOC(npx * npy, double)) ||
!(px = nAG_ALLOC(npx, double)) || !(py = nAG_ALLOC(npy, double)))
{
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
}
else {
printf("Invalid npx or npy.\n");
exit_status = 1;
return exit_status;
}
delta = (xhi - xlo) / (npx - 1);
for (i = 0; i < npx; i++)
px[i] = MIN(xlo + i * delta, xhi);
for (i = 0; i < npy; i++)
py[i] = MIN(ylo + i * delta, yhi);
/* nag_2d_spline_eval_rect (e02dfc).
* Evaluation of bicubic spline, at a mesh of points
*/
nag_2d_spline_eval_rect(npx, npy, px, py, fg, &spline, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_2d_spline_eval_rect (e02dfc).\n%s\n",
fail.message);
exit_status = 1;
goto END;
}
printf("\nValues of computed spline:\n");
printf("\n x");
for (i = 0; i < npx; i++)
printf("%7.2f ", px[i]);
printf("\n y\n");
for (i = npy - 1; i >= 0; i--) {
printf("%7.2f ", py[i]);
for (j = 0; j < npx; ++j)
printf("%8.2f ", fg[npy * j + i]);
printf("\n");
}
END:
nAG_FREE(spline.lamda);
nAG_FREE(spline.mu);
nAG_FREE(spline.c);
nAG_FREE(warmstartinf.nag_w);
nAG_FREE(warmstartinf.nag_iw);
nAG_FREE(f);
nAG_FREE(x);
nAG_FREE(y);
nAG_FREE(fg);
nAG_FREE(px);
nAG_FREE(py);
return exit_status;
}