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

C#，数值计算——有理函数插值和外推（Rational_interp）的计算方法与源程序

news 来源：原创 2024/5/9 7:36:53

1 文本格式

using System;

namespace Legalsoft.Truffer
{
/// <summary>
/// 有理函数插值和外推
/// Rational Function Interpolation and Extrapolation
/// Given a value x, and using pointers to data xx and yy, this routine returns
/// an interpolated value y, and stores an error estimate dy. The returned value
/// is obtained by mm-point polynomial interpolation on the subrange
/// xx[jl..jl + mm - 1].
/// </summary>
public class Rational_interp : Base_interp
{
private double dy { get; set; }

public Rational_interp(double[] xv, double[] yv, int m) : base(xv, yv[0], m)
{
this.dy = 0.0;
}

/// <summary>
/// Given a value x, and using pointers to data xx and yy, this routine returns
/// an interpolated value y, and stores an error estimate dy. The returned
/// value is obtained by mm-point diagonal rational function interpolation on
/// the subrange xx[jl..jl + mm - 1].
/// </summary>
/// <param name="jl"></param>
/// <param name="x"></param>
/// <returns></returns>
/// <exception cref="Exception"></exception>
public override double rawinterp(int jl, double x)
{
const double TINY = 1.0e-99;
int ns = 0;
double[] c = new double[mm];
double[] d = new double[mm];
double hh = Math.Abs(x - xx[jl + 0]);
for (int i = 0; i < mm; i++)
{
double h = Math.Abs(x - xx[jl + i]);
//if (h == 0.0)
if (Math.Abs(h) <= float.Epsilon)
{
dy = 0.0;
return yy[jl + i];
}
else if (h < hh)
{
ns = i;
hh = h;
}
c[i] = yy[jl + i];
d[i] = yy[jl + i] + TINY;
}
double y = yy[jl + ns--];
for (int m = 1; m < mm; m++)
{
for (int i = 0; i < mm - m; i++)
{
double w = c[i + 1] - d[i];
double h = xx[jl + i + m] - x;
double t = (xx[jl + i] - x) * d[i] / h;
double dd = t - c[i + 1];
//if (dd == 0.0)
if (Math.Abs(dd) <= float.Epsilon)
{
throw new Exception("Error in routine ratint");
}
dd = w / dd;
d[i] = c[i + 1] * dd;
c[i] = t * dd;
}
y += (dy = (2 * (ns + 1) < (mm - m) ? c[ns + 1] : d[ns--]));
}
return y;
}
}
}

2 代码格式

using System;namespace Legalsoft.Truffer
{/// <summary>/// 有理函数插值和外推/// Rational Function Interpolation and Extrapolation/// Given a value x, and using pointers to data xx and yy, this routine returns/// an interpolated value y, and stores an error estimate dy. The returned value/// is obtained by mm-point polynomial interpolation on the subrange/// xx[jl..jl + mm - 1]./// </summary>public class Rational_interp : Base_interp{private double dy { get; set; }public Rational_interp(double[] xv, double[] yv, int m) : base(xv, yv[0], m){this.dy = 0.0;}/// <summary>/// Given a value x, and using pointers to data xx and yy, this routine returns/// an interpolated value y, and stores an error estimate dy. The returned/// value is obtained by mm-point diagonal rational function interpolation on/// the subrange xx[jl..jl + mm - 1]./// </summary>/// <param name="jl"></param>/// <param name="x"></param>/// <returns></returns>/// <exception cref="Exception"></exception>public override double rawinterp(int jl, double x){const double TINY = 1.0e-99;int ns = 0;double[] c = new double[mm];double[] d = new double[mm];double hh = Math.Abs(x - xx[jl + 0]);for (int i = 0; i < mm; i++){double h = Math.Abs(x - xx[jl + i]);//if (h == 0.0)if (Math.Abs(h) <= float.Epsilon){dy = 0.0;return yy[jl + i];}else if (h < hh){ns = i;hh = h;}c[i] = yy[jl + i];d[i] = yy[jl + i] + TINY;}double y = yy[jl + ns--];for (int m = 1; m < mm; m++){for (int i = 0; i < mm - m; i++){double w = c[i + 1] - d[i];double h = xx[jl + i + m] - x;double t = (xx[jl + i] - x) * d[i] / h;double dd = t - c[i + 1];//if (dd == 0.0)if (Math.Abs(dd) <= float.Epsilon){throw new Exception("Error in routine ratint");}dd = w / dd;d[i] = c[i + 1] * dd;c[i] = t * dd;}y += (dy = (2 * (ns + 1) < (mm - m) ? c[ns + 1] : d[ns--]));}return y;}}
}