1 文本格式
using System;
namespace Legalsoft.Truffer
{
/// <summary>
/// Routine implementing the extended trapezoidal rule.
/// </summary>
public class Trapzd : Quadrature
{
/// <summary>
/// Limits of integration and current value of integral.
/// </summary>
private double a { get; set; } = 0.0;
private double b { get; set; } = 0.0;
private double s { get; set; } = 0.0;
private UniVarRealValueFun funx;
public Trapzd()
{
}
/// <summary>
/// The constructor takes as inputs func, the function or functor to be
/// integrated between limits a and b, also input.
/// </summary>
/// <param name="funcc"></param>
/// <param name="aa"></param>
/// <param name="bb"></param>
public Trapzd(UniVarRealValueFun funcc, double aa, double bb)
{
this.funx = funcc;
this.a = aa;
this.b = bb;
n = 0;
}
/// <summary>
/// Returns the nth stage of refinement of the extended trapezoidal rule.On
/// the first call(n= 1),the routine returns the crudest estimate of S(a, b)f(x)dx.
/// Subsequent calls set n=2,3,... and improve the accuracy by adding 2n-2
/// additional interior points.
/// </summary>
/// <returns></returns>
public override double next()
{
n++;
if (n == 1)
{
return (s = 0.5 * (b - a) * (funx.funk(a) + funx.funk(b)));
}
else
{
int it = 1;
for (int j = 1; j < n - 1; j++)
{
it <<= 1;
}
double tnm = it;
double del = (b - a) / tnm;
double x = a + 0.5 * del;
double sum = 0.0;
for (int j = 0; j < it; j++, x += del)
{
sum += funx.funk(x);
}
s = 0.5 * (s + (b - a) * sum / tnm);
return s;
}
}
/// <summary>
/// Returns the integral of the function or functor func from a to b.The
/// constants EPS can be set to the desired fractional accuracy and JMAX so
/// that 2 to the power JMAX-1 is the maximum allowed number of steps.
/// Integration is performed by the trapezoidal rule.
/// </summary>
/// <param name="funcc"></param>
/// <param name="a"></param>
/// <param name="b"></param>
/// <param name="eps"></param>
/// <returns></returns>
/// <exception cref="Exception"></exception>
public static double qtrap(UniVarRealValueFun funcc, double a, double b, double eps = 1.0e-10)
{
const int JMAX = 20;
double s;
double olds = 0.0;
Trapzd t = new Trapzd(funcc, a, b);
for (int j = 0; j < JMAX; j++)
{
s = t.next();
if (j > 5)
{
//if (Math.Abs(s - olds) < eps * Math.Abs(olds) || (s == 0.0 && olds == 0.0))
if (Math.Abs(s - olds) < eps * Math.Abs(olds) ||
(Math.Abs(s) <= float.Epsilon && Math.Abs(olds) <= float.Epsilon)
)
{
return s;
}
}
olds = s;
}
throw new Exception("Too many steps in routine qtrap");
}
/// <summary>
/// Returns the integral of the function or functor func from a to b.The
/// constants EPS can be set to the desired fractional accuracy and JMAX so
/// that 2 to the power JMAX-1 is the maximum allowed number of steps.
/// Integration is performed by Simpson's rule.
/// </summary>
/// <param name="funcc"></param>
/// <param name="a"></param>
/// <param name="b"></param>
/// <param name="eps"></param>
/// <returns></returns>
/// <exception cref="Exception"></exception>
public static double qsimp(UniVarRealValueFun funcc, double a, double b, double eps = 1.0e-10)
{
const int JMAX = 20;
double ost = 0.0;
double os = 0.0;
Trapzd t = new Trapzd(funcc, a, b);
for (int j = 0; j < JMAX; j++)
{
double st = t.next();
double s = (4.0 * st - ost) / 3.0;
if (j > 5)
{
//if (Math.Abs(s - os) < eps * Math.Abs(os) || (s == 0.0 && os == 0.0))
if (Math.Abs(s - os) < eps * Math.Abs(os) || (Math.Abs(s) <= float.Epsilon && Math.Abs(os) <= float.Epsilon))
{
return s;
}
}
os = s;
ost = st;
}
throw new Exception("Too many steps in routine qsimp");
}
public static double qromb(UniVarRealValueFun func, double a, double b)
{
return qromb(func, a, b, 1.0e-10);
}
/// <summary>
/// Returns the integral of the function or functor func from a to b.
/// Integration is performed by Romberg's method of order 2K, where, e.g.,
/// K=2 is Simpson's rule.
/// </summary>
/// <param name="funcc"></param>
/// <param name="a"></param>
/// <param name="b"></param>
/// <param name="eps"></param>
/// <returns></returns>
/// <exception cref="Exception"></exception>
public static double qromb(UniVarRealValueFun funcc, double a, double b, double eps)
{
int JMAX = 20;
int JMAXP = JMAX + 1;
int K = 5;
double[] s = new double[JMAX];
double[] h = new double[JMAXP];
Poly_interp polint = new Poly_interp(h, s, K);
h[0] = 1.0;
Trapzd t = new Trapzd(funcc, a, b);
for (int j = 1; j <= JMAX; j++)
{
s[j - 1] = t.next();
if (j >= K)
{
double ss = polint.rawinterp(j - K, 0.0);
if (Math.Abs(polint.dy) <= eps * Math.Abs(ss)) return ss;
}
h[j] = 0.25 * h[j - 1];
}
throw new Exception("Too many steps in routine qromb");
}
}
}
2 代码格式
using System;
namespace Legalsoft.Truffer
{
/// <summary>
/// Routine implementing the extended trapezoidal rule.
/// </summary>
public class Trapzd : Quadrature
{
/// <summary>
/// Limits of integration and current value of integral.
/// </summary>
private double a { get; set; } = 0.0;
private double b { get; set; } = 0.0;
private double s { get; set; } = 0.0;
private UniVarRealValueFun funx;
public Trapzd()
{
}
/// <summary>
/// The constructor takes as inputs func, the function or functor to be
/// integrated between limits a and b, also input.
/// </summary>
/// <param name="funcc"></param>
/// <param name="aa"></param>
/// <param name="bb"></param>
public Trapzd(UniVarRealValueFun funcc, double aa, double bb)
{
this.funx = funcc;
this.a = aa;
this.b = bb;
n = 0;
}
/// <summary>
/// Returns the nth stage of refinement of the extended trapezoidal rule.On
/// the first call(n= 1),the routine returns the crudest estimate of S(a, b)f(x)dx.
/// Subsequent calls set n=2,3,... and improve the accuracy by adding 2n-2
/// additional interior points.
/// </summary>
/// <returns></returns>
public override double next()
{
n++;
if (n == 1)
{
return (s = 0.5 * (b - a) * (funx.funk(a) + funx.funk(b)));
}
else
{
int it = 1;
for (int j = 1; j < n - 1; j++)
{
it <<= 1;
}
double tnm = it;
double del = (b - a) / tnm;
double x = a + 0.5 * del;
double sum = 0.0;
for (int j = 0; j < it; j++, x += del)
{
sum += funx.funk(x);
}
s = 0.5 * (s + (b - a) * sum / tnm);
return s;
}
}
/// <summary>
/// Returns the integral of the function or functor func from a to b.The
/// constants EPS can be set to the desired fractional accuracy and JMAX so
/// that 2 to the power JMAX-1 is the maximum allowed number of steps.
/// Integration is performed by the trapezoidal rule.
/// </summary>
/// <param name="funcc"></param>
/// <param name="a"></param>
/// <param name="b"></param>
/// <param name="eps"></param>
/// <returns></returns>
/// <exception cref="Exception"></exception>
public static double qtrap(UniVarRealValueFun funcc, double a, double b, double eps = 1.0e-10)
{
const int JMAX = 20;
double s;
double olds = 0.0;
Trapzd t = new Trapzd(funcc, a, b);
for (int j = 0; j < JMAX; j++)
{
s = t.next();
if (j > 5)
{
//if (Math.Abs(s - olds) < eps * Math.Abs(olds) || (s == 0.0 && olds == 0.0))
if (Math.Abs(s - olds) < eps * Math.Abs(olds) ||
(Math.Abs(s) <= float.Epsilon && Math.Abs(olds) <= float.Epsilon)
)
{
return s;
}
}
olds = s;
}
throw new Exception("Too many steps in routine qtrap");
}
/// <summary>
/// Returns the integral of the function or functor func from a to b.The
/// constants EPS can be set to the desired fractional accuracy and JMAX so
/// that 2 to the power JMAX-1 is the maximum allowed number of steps.
/// Integration is performed by Simpson's rule.
/// </summary>
/// <param name="funcc"></param>
/// <param name="a"></param>
/// <param name="b"></param>
/// <param name="eps"></param>
/// <returns></returns>
/// <exception cref="Exception"></exception>
public static double qsimp(UniVarRealValueFun funcc, double a, double b, double eps = 1.0e-10)
{
const int JMAX = 20;
double ost = 0.0;
double os = 0.0;
Trapzd t = new Trapzd(funcc, a, b);
for (int j = 0; j < JMAX; j++)
{
double st = t.next();
double s = (4.0 * st - ost) / 3.0;
if (j > 5)
{
//if (Math.Abs(s - os) < eps * Math.Abs(os) || (s == 0.0 && os == 0.0))
if (Math.Abs(s - os) < eps * Math.Abs(os) || (Math.Abs(s) <= float.Epsilon && Math.Abs(os) <= float.Epsilon))
{
return s;
}
}
os = s;
ost = st;
}
throw new Exception("Too many steps in routine qsimp");
}
public static double qromb(UniVarRealValueFun func, double a, double b)
{
return qromb(func, a, b, 1.0e-10);
}
/// <summary>
/// Returns the integral of the function or functor func from a to b.
/// Integration is performed by Romberg's method of order 2K, where, e.g.,
/// K=2 is Simpson's rule.
/// </summary>
/// <param name="funcc"></param>
/// <param name="a"></param>
/// <param name="b"></param>
/// <param name="eps"></param>
/// <returns></returns>
/// <exception cref="Exception"></exception>
public static double qromb(UniVarRealValueFun funcc, double a, double b, double eps)
{
int JMAX = 20;
int JMAXP = JMAX + 1;
int K = 5;
double[] s = new double[JMAX];
double[] h = new double[JMAXP];
Poly_interp polint = new Poly_interp(h, s, K);
h[0] = 1.0;
Trapzd t = new Trapzd(funcc, a, b);
for (int j = 1; j <= JMAX; j++)
{
s[j - 1] = t.next();
if (j >= K)
{
double ss = polint.rawinterp(j - K, 0.0);
if (Math.Abs(polint.dy) <= eps * Math.Abs(ss)) return ss;
}
h[j] = 0.25 * h[j - 1];
}
throw new Exception("Too many steps in routine qromb");
}
}
}
