1 文本格式
using System;
namespace Legalsoft.Truffer
{
/// <summary>
/// 构造n点等间隔求积的权重
/// Constructs weights for the n-point equal-interval quadrature
/// from O to(n-1)h of a function f(x) times an arbitrary
/// (possibly singular) weight function w(x). The indefinite-integral
/// moments Fn(y) of w(x) are provided by the user-supplied function
/// kermom in the quad object.
/// </summary>
public class Wwghts
{
private double h { get; set; }
private int n { get; set; }
//private Quad_matrix quad { get; set; }
UniVarRealMultiValueFun quad;
private double[] wghts { get; set; }
public Wwghts(double hh, int nn, UniVarRealMultiValueFun q)
{
this.h = hh;
this.n = nn;
this.quad = q;
this.wghts = new double[n];
}
public double[] weights()
{
double hi = 1.0 / h;
for (int j = 0; j < n; j++)
{
wghts[j] = 0.0;
}
if (n >= 4)
{
double[] w = new double[4];
double[] wold = quad.funk(0.0);
double b = 0.0;
for (int j = 0; j < n - 3; j++)
{
double c = j;
double a = b;
b = a + h;
if (j == n - 4)
{
b = (n - 1) * h;
}
double[] wnew = quad.funk(b);
double fac = 1.0;
for (int k = 0; k < 4; k++, fac *= hi)
{
w[k] = (wnew[k] - wold[k]) * fac;
}
wghts[j] += (((c + 1.0) * (c + 2.0) * (c + 3.0) * w[0] - (11.0 + c * (12.0 + c * 3.0)) * w[1] + 3.0 * (c + 2.0) * w[2] - w[3]) / 6.0);
wghts[j + 1] += ((-c * (c + 2.0) * (c + 3.0) * w[0] + (6.0 + c * (10.0 + c * 3.0)) * w[1] - (3.0 * c + 5.0) * w[2] + w[3]) * 0.5);
wghts[j + 2] += ((c * (c + 1.0) * (c + 3.0) * w[0] - (3.0 + c * (8.0 + c * 3.0)) * w[1] + (3.0 * c + 4.0) * w[2] - w[3]) * 0.5);
wghts[j + 3] += ((-c * (c + 1.0) * (c + 2.0) * w[0] + (2.0 + c * (6.0 + c * 3.0)) * w[1] - 3.0 * (c + 1.0) * w[2] + w[3]) / 6.0);
for (int k = 0; k < 4; k++)
{
wold[k] = wnew[k];
}
}
}
else if (n == 3)
{
double[] w = new double[3];
double[] wold = quad.funk(0.0);
double[] wnew = quad.funk(h + h);
w[0] = wnew[0] - wold[0];
w[1] = hi * (wnew[1] - wold[1]);
w[2] = hi * hi * (wnew[2] - wold[2]);
wghts[0] = w[0] - 1.5 * w[1] + 0.5 * w[2];
wghts[1] = 2.0 * w[1] - w[2];
wghts[2] = 0.5 * (w[2] - w[1]);
}
else if (n == 2)
{
double[] wold = quad.funk(0.0);
double[] wnew = quad.funk(h);
wghts[0] = wnew[0] - wold[0] - (wghts[1] = hi * (wnew[1] - wold[1]));
}
return wghts;
}
}
}
2 代码格式
using System;
namespace Legalsoft.Truffer
{
/// <summary>
/// 构造n点等间隔求积的权重
/// Constructs weights for the n-point equal-interval quadrature
/// from O to(n-1)h of a function f(x) times an arbitrary
/// (possibly singular) weight function w(x). The indefinite-integral
/// moments Fn(y) of w(x) are provided by the user-supplied function
/// kermom in the quad object.
/// </summary>
public class Wwghts
{
private double h { get; set; }
private int n { get; set; }
//private Quad_matrix quad { get; set; }
UniVarRealMultiValueFun quad;
private double[] wghts { get; set; }
public Wwghts(double hh, int nn, UniVarRealMultiValueFun q)
{
this.h = hh;
this.n = nn;
this.quad = q;
this.wghts = new double[n];
}
public double[] weights()
{
double hi = 1.0 / h;
for (int j = 0; j < n; j++)
{
wghts[j] = 0.0;
}
if (n >= 4)
{
double[] w = new double[4];
double[] wold = quad.funk(0.0);
double b = 0.0;
for (int j = 0; j < n - 3; j++)
{
double c = j;
double a = b;
b = a + h;
if (j == n - 4)
{
b = (n - 1) * h;
}
double[] wnew = quad.funk(b);
double fac = 1.0;
for (int k = 0; k < 4; k++, fac *= hi)
{
w[k] = (wnew[k] - wold[k]) * fac;
}
wghts[j] += (((c + 1.0) * (c + 2.0) * (c + 3.0) * w[0] - (11.0 + c * (12.0 + c * 3.0)) * w[1] + 3.0 * (c + 2.0) * w[2] - w[3]) / 6.0);
wghts[j + 1] += ((-c * (c + 2.0) * (c + 3.0) * w[0] + (6.0 + c * (10.0 + c * 3.0)) * w[1] - (3.0 * c + 5.0) * w[2] + w[3]) * 0.5);
wghts[j + 2] += ((c * (c + 1.0) * (c + 3.0) * w[0] - (3.0 + c * (8.0 + c * 3.0)) * w[1] + (3.0 * c + 4.0) * w[2] - w[3]) * 0.5);
wghts[j + 3] += ((-c * (c + 1.0) * (c + 2.0) * w[0] + (2.0 + c * (6.0 + c * 3.0)) * w[1] - 3.0 * (c + 1.0) * w[2] + w[3]) / 6.0);
for (int k = 0; k < 4; k++)
{
wold[k] = wnew[k];
}
}
}
else if (n == 3)
{
double[] w = new double[3];
double[] wold = quad.funk(0.0);
double[] wnew = quad.funk(h + h);
w[0] = wnew[0] - wold[0];
w[1] = hi * (wnew[1] - wold[1]);
w[2] = hi * hi * (wnew[2] - wold[2]);
wghts[0] = w[0] - 1.5 * w[1] + 0.5 * w[2];
wghts[1] = 2.0 * w[1] - w[2];
wghts[2] = 0.5 * (w[2] - w[1]);
}
else if (n == 2)
{
double[] wold = quad.funk(0.0);
double[] wnew = quad.funk(h);
wghts[0] = wnew[0] - wold[0] - (wghts[1] = hi * (wnew[1] - wold[1]));
}
return wghts;
}
}
}
