当前位置：首页 > news >正文

C#，数值计算——插值和外推，曲线插值（Curve_interp）的计算方法与源程序

news 来源：原创 2024/5/20 13:56:24

1 文本格式

using System;

namespace Legalsoft.Truffer
{
/// <summary>
/// Object for interpolating a curve specified by n points in dim dimensions.
/// </summary>
public class Curve_interp
{
private int dim { get; set; }
private int n { get; set; }
private int xin { get; set; }
private bool cls { get; set; }
private double[,] pts { get; set; }
private double[] s { get; set; }
private double[] ans { get; set; }
private Spline_interp[] srp { get; set; }// = new Spline_interp[];

/// <summary>
/// The n x dim matrix ptsin inputs the data points.Input close as 0 for an
/// open curve, 1 for a closed curve. (For a closed curve, the last data point
/// should not duplicate the first 鈥?the algorithm will connect them.)
/// </summary>
/// <param name="ptsin"></param>
/// <param name="close"></param>
/// <exception cref="Exception"></exception>
public Curve_interp(double[,] ptsin, bool close = false)
{
this.n = ptsin.GetLength(0);
this.dim = ptsin.GetLength(1);
this.xin = close ? 2 * n : n;
this.cls = close;
this.pts = new double[dim, xin];
this.s = new double[xin];
this.ans = new double[dim];
this.srp = new Spline_interp[dim];

//int i;
//int ii;
//int im;
//int j;
//int ofs;
//double ss;
//double soff;
//double db;
//double de;
int ofs = close ? n / 2 : 0;
s[0] = 0.0;
for (int i = 0; i < xin; i++)
{
int ii = (i - ofs + n) % n;
int im = (ii - 1 + n) % n;
for (int j = 0; j < dim; j++)
{
pts[j, i] = ptsin[ii, j];
}
if (i > 0)
{
//s[i] = s[i - 1] + rad(ptsin.GetRow(ii).ToArray(), ptsin.GetRow(im).ToArray());
s[i] = s[i - 1] + Globals.dist(Globals.CopyFrom(ii, ptsin), Globals.CopyFrom(im, ptsin));
if (s[i] == s[i - 1])
{
throw new Exception("error in Curve_interp");
}
}
}
double ss = close ? s[ofs + n] - s[ofs] : s[n - 1] - s[0];
double soff = s[ofs];
for (int i = 0; i < xin; i++)
{
s[i] = (s[i] - soff) / ss;
}
for (int j = 0; j < dim; j++)
{
double db = xin < 4 ? 1.0e99 : fprime(s, 0, Globals.CopyFrom(j, pts), 0, 1);
double de = xin < 4 ? 1.0e99 : fprime(s, xin - 1, Globals.CopyFrom(j, pts), xin - 1, -1);
srp[j] = new Spline_interp(s, pts[j, 0], db, de);
}
}

/// <summary>
/// Interpolate a point on the stored curve.The point is parameterized by t,
/// in the range[0, 1]. For open curves, values of t outside this range will
/// return extrapolations(dangerous!). For closed curves, t is periodic with
/// period 1.
/// </summary>
/// <param name="t"></param>
/// <returns></returns>
public double[] interp(double t)
{
if (cls)
{
t = t - Math.Floor(t);
}
for (int j = 0; j < dim; j++)
{
ans[j] = (srp[j]).interp(t);
}
return ans;
}

/// <summary>
/// Utility for estimating the derivatives at the endpoints.x and y point to
/// the abscissa and ordinate of the endpoint.If pm is C1, points to the right
/// will be used (left endpoint); if it is 1, points to the left will be used
/// (right endpoint).
/// </summary>
/// <param name="x"></param>
/// <param name="y"></param>
/// <param name="pm"></param>
/// <returns></returns>
public double fprime(double[] x, double[] y, int pm)
{
double s1 = x[0] - x[pm * 1];
double s2 = x[0] - x[pm * 2];
double s3 = x[0] - x[pm * 3];
double s12 = s1 - s2;
double s13 = s1 - s3;
double s23 = s2 - s3;
return -(s1 * s2 / (s13 * s23 * s3)) * y[pm * 3] + (s1 * s3 / (s12 * s2 * s23)) * y[pm * 2] - (s2 * s3 / (s1 * s12 * s13)) * y[pm * 1] + (1.0 / s1 + 1.0 / s2 + 1.0 / s3) * y[0];
}

private double fprime(double[] x, int off_x, double[] y, int off_y, int pm)
{
double s1 = x[off_x + 0] - x[off_x + pm * 1];
double s2 = x[off_x + 0] - x[off_x + pm * 2];
double s3 = x[off_x + 0] - x[off_x + pm * 3];
double s12 = s1 - s2;
double s13 = s1 - s3;
double s23 = s2 - s3;
return -(s1 * s2 / (s13 * s23 * s3)) * y[off_y + pm * 3] + (s1 * s3 / (s12 * s2 * s23)) * y[off_y + pm * 2] - (s2 * s3 / (s1 * s12 * s13)) * y[off_y + pm * 1] + (1.0 / s1 + 1.0 / s2 + 1.0 / s3) * y[off_y + 0];
}
/*
public double rad(double[] p1, double[] p2)
{
double sum = 0.0;
for (int i = 0; i < dim; i++)
{
sum += Globals.SQR(p1[i] - p2[i]);
}
if (sum <= float.Epsilon)
{
return 0.0;
}
return Math.Sqrt(sum);
}
*/
}
}

2 代码格式

using System;namespace Legalsoft.Truffer
{/// <summary>/// Object for interpolating a curve specified by n points in dim dimensions./// </summary>public class Curve_interp{private int dim { get; set; }private int n { get; set; }private int xin { get; set; }private bool cls { get; set; }private double[,] pts { get; set; }private double[] s { get; set; }private double[] ans { get; set; }private Spline_interp[] srp { get; set; }// = new Spline_interp[];/// <summary>/// The n x dim matrix ptsin inputs the data points.Input close as 0 for an/// open curve, 1 for a closed curve. (For a closed curve, the last data point/// should not duplicate the first 鈥?the algorithm will connect them.)/// </summary>/// <param name="ptsin"></param>/// <param name="close"></param>/// <exception cref="Exception"></exception>public Curve_interp(double[,] ptsin, bool close = false){this.n = ptsin.GetLength(0);this.dim = ptsin.GetLength(1);this.xin = close ? 2 * n : n;this.cls = close;this.pts = new double[dim, xin];this.s = new double[xin];this.ans = new double[dim];this.srp = new Spline_interp[dim];//int i;//int ii;//int im;//int j;//int ofs;//double ss;//double soff;//double db;//double de;int ofs = close ? n / 2 : 0;s[0] = 0.0;for (int i = 0; i < xin; i++){int ii = (i - ofs + n) % n;int im = (ii - 1 + n) % n;for (int j = 0; j < dim; j++){pts[j, i] = ptsin[ii, j];}if (i > 0){//s[i] = s[i - 1] + rad(ptsin.GetRow(ii).ToArray(), ptsin.GetRow(im).ToArray());s[i] = s[i - 1] + Globals.dist(Globals.CopyFrom(ii, ptsin), Globals.CopyFrom(im, ptsin));if (s[i] == s[i - 1]){throw new Exception("error in Curve_interp");}}}double ss = close ? s[ofs + n] - s[ofs] : s[n - 1] - s[0];double soff = s[ofs];for (int i = 0; i < xin; i++){s[i] = (s[i] - soff) / ss;}for (int j = 0; j < dim; j++){double db = xin < 4 ? 1.0e99 : fprime(s, 0, Globals.CopyFrom(j, pts), 0, 1);double de = xin < 4 ? 1.0e99 : fprime(s, xin - 1, Globals.CopyFrom(j, pts), xin - 1, -1);srp[j] = new Spline_interp(s, pts[j, 0], db, de);}}/// <summary>/// Interpolate a point on the stored curve.The point is parameterized by t,/// in the range[0, 1]. For open curves, values of t outside this range will/// return extrapolations(dangerous!). For closed curves, t is periodic with/// period 1./// </summary>/// <param name="t"></param>/// <returns></returns>public double[] interp(double t){if (cls){t = t - Math.Floor(t);}for (int j = 0; j < dim; j++){ans[j] = (srp[j]).interp(t);}return ans;}/// <summary>/// Utility for estimating the derivatives at the endpoints.x and y point to/// the abscissa and ordinate of the endpoint.If pm is C1, points to the right/// will be used (left endpoint); if it is 1, points to the left will be used/// (right endpoint)./// </summary>/// <param name="x"></param>/// <param name="y"></param>/// <param name="pm"></param>/// <returns></returns>public double fprime(double[] x, double[] y, int pm){double s1 = x[0] - x[pm * 1];double s2 = x[0] - x[pm * 2];double s3 = x[0] - x[pm * 3];double s12 = s1 - s2;double s13 = s1 - s3;double s23 = s2 - s3;return -(s1 * s2 / (s13 * s23 * s3)) * y[pm * 3] + (s1 * s3 / (s12 * s2 * s23)) * y[pm * 2] - (s2 * s3 / (s1 * s12 * s13)) * y[pm * 1] + (1.0 / s1 + 1.0 / s2 + 1.0 / s3) * y[0];}private double fprime(double[] x, int off_x, double[] y, int off_y, int pm){double s1 = x[off_x + 0] - x[off_x + pm * 1];double s2 = x[off_x + 0] - x[off_x + pm * 2];double s3 = x[off_x + 0] - x[off_x + pm * 3];double s12 = s1 - s2;double s13 = s1 - s3;double s23 = s2 - s3;return -(s1 * s2 / (s13 * s23 * s3)) * y[off_y + pm * 3] + (s1 * s3 / (s12 * s2 * s23)) * y[off_y + pm * 2] - (s2 * s3 / (s1 * s12 * s13)) * y[off_y + pm * 1] + (1.0 / s1 + 1.0 / s2 + 1.0 / s3) * y[off_y + 0];}/*public double rad(double[] p1, double[] p2){double sum = 0.0;for (int i = 0; i < dim; i++){sum += Globals.SQR(p1[i] - p2[i]);}if (sum <= float.Epsilon){return 0.0;}return Math.Sqrt(sum);}*/}
}