C#，数值计算——插值和外推，二维三次样条插值（Spline2D

1 文本格式

using System;

namespace Legalsoft.Truffer
{
/// <summary>
/// 二维三次样条插值
/// Object for two-dimensional cubic spline interpolation on a matrix.Construct
/// with a vector of x1 values, a vector of x2 values, and a matrix of tabulated
/// function values yij.Then call interp for interpolated values.
/// </summary>
public class Spline2D_interp
{
private int[,] bcucof_wt = new int[,] {
{ 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 },
{ 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0 },
{ -3, 0, 0, 3, 0, 0, 0, 0, -2, 0, 0, -1, 0, 0, 0, 0 },
{ 2, 0, 0, -2, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0 },
{ 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 },
{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0 },
{ 0, 0, 0, 0, -3, 0, 0, 3, 0, 0, 0, 0, -2, 0, 0, -1 },
{ 0, 0, 0, 0, 2, 0, 0, -2, 0, 0, 0, 0, 1, 0, 0, 1 },
{ -3, 3, 0, 0, -2, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 },
{ 0, 0, 0, 0, 0, 0, 0, 0, -3, 3, 0, 0, -2, -1, 0, 0 },
{ 9, -9, 9, -9, 6, 3, -3, -6, 6, -6, -3, 3, 4, 2, 1, 2 },
{ -6, 6, -6, 6, -4, -2, 2, 4, -3, 3, 3, -3, -2, -1, -1, -2 },
{ 2, -2, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 },
{ 0, 0, 0, 0, 0, 0, 0, 0, 2, -2, 0, 0, 1, 1, 0, 0 },
{ -6, 6, -6, 6, -3, -3, 3, 3, -4, 4, 2, -2, -2, -2, -1, -1 },
{ 4, -4, 4, -4, 2, 2, -2, -2, 2, -2, -2, 2, 1, 1, 1, 1 }
};

private int m { get; set; }
private int n { get; set; }
private double[,] y { get; set; }
private double[] x1 { get; set; }
private double[] yv { get; set; }
private Spline_interp[] srp { get; set; }

public Spline2D_interp(double[] x1v, double[] x2v, double[,] ym)
{
this.m = x1v.Length;
this.n = x2v.Length;
this.y = ym;
this.yv = new double[m];
this.x1 = x1v;
this.srp = new Spline_interp[m];
for (int i = 0; i < m; i++)
{
srp[i] = new Spline_interp(x2v, y[i, 0]);
}
}

public double interp(double x1p, double x2p)
{
for (int i = 0; i < m; i++)
{
yv[i] = (srp[i]).interp(x2p);
}
Spline_interp scol = new Spline_interp(x1, yv);
return scol.interp(x1p);
}

/// <summary>
/// Given arrays y[0..3], y1[0..3], y2[0..3], and y12[0..3], containing the
/// function, gradients, and cross-derivative at the four grid points of a
/// rectangular grid cell(numbered counterclockwise from the lower left), and
/// given d1 and d2, the length of the grid cell in the 1 and 2 directions,
/// this routine returns the table c[0..3][0..3] that is used by routine bcuint
/// for bicubic interpolation.
/// </summary>
/// <param name="y"></param>
/// <param name="y1"></param>
/// <param name="y2"></param>
/// <param name="y12"></param>
/// <param name="d1"></param>
/// <param name="d2"></param>
/// <param name="c"></param>
private void bcucof(double[] y, double[] y1, double[] y2, double[] y12, double d1, double d2, double[,] c)
{
//int[,] bcucof_wt = new int[16, 16]{ bcucof_wt_d };
double d1d2 = d1 * d2;
double[] cl = new double[16];
double[] x = new double[16];

for (int i = 0; i < 4; i++)
{
x[i] = y[i];
x[i + 4] = y1[i] * d1;
x[i + 8] = y2[i] * d2;
x[i + 12] = y12[i] * d1d2;
}
for (int i = 0; i < 16; i++)
{
double xx = 0.0;
for (int k = 0; k < 16; k++)
{
xx += bcucof_wt[i, k] * x[k];
}
cl[i] = xx;
}
int l = 0;
for (int i = 0; i < 4; i++)
{
for (int j = 0; j < 4; j++)
{
c[i, j] = cl[l++];
}
}
}

/// <summary>
/// Bicubic interpolation within a grid square.Input quantities are
/// y, y1, y2, y12 (as described in bcucof); x1l and x1u, the lower and upper
/// coordinates of the grid square in the 1 direction; x2l and x2u likewise for
/// the 2 direction; and x1, x2, the coordinates of the desired point for the
/// interpolation.The interpolated function value is returned as ansy, and the
/// interpolated gradient values as ansy1 and ansy2.This routine calls bcucof.
/// </summary>
/// <param name="y"></param>
/// <param name="y1"></param>
/// <param name="y2"></param>
/// <param name="y12"></param>
/// <param name="x1l"></param>
/// <param name="x1u"></param>
/// <param name="x2l"></param>
/// <param name="x2u"></param>
/// <param name="x1"></param>
/// <param name="x2"></param>
/// <param name="ansy"></param>
/// <param name="ansy1"></param>
/// <param name="ansy2"></param>
/// <exception cref="Exception"></exception>
public void bcuint(double[] y, double[] y1, double[] y2, double[] y12, double x1l, double x1u, double x2l, double x2u, double x1, double x2, ref double ansy, ref double ansy1, ref double ansy2)
{
double d1 = x1u - x1l;
double d2 = x2u - x2l;
double[,] c = new double[4, 4];
bcucof(y, y1, y2, y12, d1, d2, c);
//if (x1u == x1l || x2u == x2l)
if (Math.Abs(x1u - x1l) <= float.Epsilon || Math.Abs(x2u - x2l) <= float.Epsilon)
{
throw new Exception("Bad input in routine bcuint");
}
double t = (x1 - x1l) / d1;
double u = (x2 - x2l) / d2;
ansy = ansy2 = ansy1 = 0.0;
for (int i = 3; i >= 0; i--)
{
ansy = t * ansy + ((c[i, 3] * u + c[i, 2]) * u + c[i, 1]) * u + c[i, 0];
ansy2 = t * ansy2 + (3.0 * c[i, 3] * u + 2.0 * c[i, 2]) * u + c[i, 1];
ansy1 = u * ansy1 + (3.0 * c[3, i] * t + 2.0 * c[2, i]) * t + c[1, i];
}
ansy1 /= d1;
ansy2 /= d2;
}

}
}

2 代码格式

using System;namespace Legalsoft.Truffer
{/// <summary>/// 二维三次样条插值/// Object for two-dimensional cubic spline interpolation on a matrix.Construct/// with a vector of x1 values, a vector of x2 values, and a matrix of tabulated/// function values yij.Then call interp for interpolated values./// </summary>public class Spline2D_interp{private int[,] bcucof_wt = new int[,] {{  1,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0, 0 },{  0,  0,  0,  0,  0,  0,  0,  0,  1,  0,  0,  0,  0,  0,  0,  0 },{ -3,  0,  0,  3,  0,  0,  0,  0, -2,  0,  0, -1,  0,  0,  0,  0 },{  2,  0,  0, -2,  0,  0,  0,  0,  1,  0,  0,  1,  0,  0,  0,  0 },{  0,  0,  0,  0,  1,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0 },{  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  1,  0,  0,  0 },{  0,  0,  0,  0, -3,  0,  0,  3,  0,  0,  0,  0, -2,  0,  0, -1 },{  0,  0,  0,  0,  2,  0,  0, -2,  0,  0,  0,  0,  1,  0,  0,  1 },{ -3,  3,  0,  0, -2, -1,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0 },{  0,  0,  0,  0,  0,  0,  0,  0, -3,  3,  0,  0, -2, -1,  0,  0 },{  9, -9,  9, -9,  6,  3, -3, -6,  6, -6, -3,  3,  4,  2,  1,  2 },{ -6,  6, -6,  6, -4, -2,  2,  4, -3,  3,  3, -3, -2, -1, -1, -2 },{  2, -2,  0,  0,  1,  1,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0 },{  0,  0,  0,  0,  0,  0,  0,  0,  2, -2,  0,  0,  1,  1,  0,  0 },{ -6,  6, -6,  6, -3, -3,  3,  3, -4,  4,  2, -2, -2, -2, -1, -1 },{  4, -4,  4, -4,  2,  2, -2, -2,  2, -2, -2,  2,  1,  1,  1,  1 }};private int m { get; set; }private int n { get; set; }private double[,] y { get; set; }private double[] x1 { get; set; }private double[] yv { get; set; }private Spline_interp[] srp { get; set; }public Spline2D_interp(double[] x1v, double[] x2v, double[,] ym){this.m = x1v.Length;this.n = x2v.Length;this.y = ym;this.yv = new double[m];this.x1 = x1v;this.srp = new Spline_interp[m];for (int i = 0; i < m; i++){srp[i] = new Spline_interp(x2v, y[i, 0]);}}public double interp(double x1p, double x2p){for (int i = 0; i < m; i++){yv[i] = (srp[i]).interp(x2p);}Spline_interp scol = new Spline_interp(x1, yv);return scol.interp(x1p);}/// <summary>/// Given arrays y[0..3], y1[0..3], y2[0..3], and y12[0..3], containing the/// function, gradients, and cross-derivative at the four grid points of a/// rectangular grid cell(numbered counterclockwise from the lower left), and/// given d1 and d2, the length of the grid cell in the 1 and 2 directions,/// this routine returns the table c[0..3][0..3] that is used by routine bcuint/// for bicubic interpolation./// </summary>/// <param name="y"></param>/// <param name="y1"></param>/// <param name="y2"></param>/// <param name="y12"></param>/// <param name="d1"></param>/// <param name="d2"></param>/// <param name="c"></param>private void bcucof(double[] y, double[] y1, double[] y2, double[] y12, double d1, double d2, double[,] c){ //int[,] bcucof_wt = new int[16, 16]{ bcucof_wt_d };double d1d2 = d1 * d2;double[] cl = new double[16];double[] x = new double[16];for (int i = 0; i < 4; i++){x[i] = y[i];x[i + 4] = y1[i] * d1;x[i + 8] = y2[i] * d2;x[i + 12] = y12[i] * d1d2;}for (int i = 0; i < 16; i++){double xx = 0.0;for (int k = 0; k < 16; k++){xx += bcucof_wt[i, k] * x[k];}cl[i] = xx;}int l = 0;for (int i = 0; i < 4; i++){for (int j = 0; j < 4; j++){c[i, j] = cl[l++];}}}/// <summary>/// Bicubic interpolation within a grid square.Input quantities are/// y, y1, y2, y12 (as described in bcucof); x1l and x1u, the lower and upper/// coordinates of the grid square in the 1 direction; x2l and x2u likewise for/// the 2 direction; and x1, x2, the coordinates of the desired point for the/// interpolation.The interpolated function value is returned as ansy, and the/// interpolated gradient values as ansy1 and ansy2.This routine calls bcucof./// </summary>/// <param name="y"></param>/// <param name="y1"></param>/// <param name="y2"></param>/// <param name="y12"></param>/// <param name="x1l"></param>/// <param name="x1u"></param>/// <param name="x2l"></param>/// <param name="x2u"></param>/// <param name="x1"></param>/// <param name="x2"></param>/// <param name="ansy"></param>/// <param name="ansy1"></param>/// <param name="ansy2"></param>/// <exception cref="Exception"></exception>public void bcuint(double[] y, double[] y1, double[] y2, double[] y12, double x1l, double x1u, double x2l, double x2u, double x1, double x2, ref double ansy, ref double ansy1, ref double ansy2){double d1 = x1u - x1l;double d2 = x2u - x2l;double[,] c = new double[4, 4];bcucof(y, y1, y2, y12, d1, d2, c);//if (x1u == x1l || x2u == x2l)if (Math.Abs(x1u - x1l) <= float.Epsilon || Math.Abs(x2u - x2l) <= float.Epsilon){throw new Exception("Bad input in routine bcuint");}double t = (x1 - x1l) / d1;double u = (x2 - x2l) / d2;ansy = ansy2 = ansy1 = 0.0;for (int i = 3; i >= 0; i--){ansy = t * ansy + ((c[i, 3] * u + c[i, 2]) * u + c[i, 1]) * u + c[i, 0];ansy2 = t * ansy2 + (3.0 * c[i, 3] * u + 2.0 * c[i, 2]) * u + c[i, 1];ansy1 = u * ansy1 + (3.0 * c[3, i] * t + 2.0 * c[2, i]) * t + c[1, i];}ansy1 /= d1;ansy2 /= d2;}}
}