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

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

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 

/* 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;
}