1 文本格式
using System;
namespace Legalsoft.Truffer
{
public class GaussianWeights
{
private static double[] x = {
0.1488743389816312, 0.4333953941292472, 0.6794095682990244,
0.8650633666889845, 0.9739065285171717
};
private static double[] w = {
0.2955242247147529, 0.2692667193099963, 0.2190863625159821,
0.1494513491505806, 0.0666713443086881
};
public GaussianWeights()
{
}
/*
public static double qgaus(UniVarRealValueFun f1, double a, double b)
{
//func = f1;
return qgaus(a,b);
}
*/
/// <summary>
/// Returns the integral of the function or functor func between a and b, by
/// ten-point Gauss-Legendre integration: the function is evaluated exactly
/// ten times at interior points in the range of integration.
/// </summary>
/// <param name="a"></param>
/// <param name="b"></param>
/// <returns></returns>
public static double qgaus(UniVarRealValueFun func, double a, double b)
{
double xm = 0.5 * (b + a);
double xr = 0.5 * (b - a);
double s = 0;
for (int j = 0; j < 5; j++)
{
double dx = xr * x[j];
s += w[j] * (func.funk(xm + dx) + func.funk(xm - dx));
}
return s *= xr;
}
/// <summary>
/// Given the lower and upper limits of integration x1 and x2, this routine
/// returns arrays x[0..n - 1] and w[0..n - 1] of length n, containing the
/// abscissas and weights of the Gauss-Legendre n-point quadrature formula.
/// </summary>
/// <param name="x1"></param>
/// <param name="x2"></param>
/// <param name="x"></param>
/// <param name="w"></param>
public static void gauleg(double x1, double x2, double[] x, double[] w)
{
const double EPS = 1.0e-14;
double z1;
double pp;
int n = x.Length;
int m = (n + 1) / 2;
double xm = 0.5 * (x2 + x1);
double xl = 0.5 * (x2 - x1);
for (int i = 0; i < m; i++)
{
double z = Math.Cos(3.141592654 * (i + 0.75) / (n + 0.5));
do
{
double p1 = 1.0;
double p2 = 0.0;
for (int j = 0; j < n; j++)
{
double p3 = p2;
p2 = p1;
p1 = ((2.0 * j + 1.0) * z * p2 - j * p3) / (j + 1);
}
pp = n * (z * p1 - p2) / (z * z - 1.0);
z1 = z;
z = z1 - p1 / pp;
} while (Math.Abs(z - z1) > EPS);
x[i] = xm - xl * z;
x[n - 1 - i] = xm + xl * z;
w[i] = 2.0 * xl / ((1.0 - z * z) * pp * pp);
w[n - 1 - i] = w[i];
}
}
/// <summary>
/// Given alf, the parameter a of the Laguerre polynomials, this routine
/// returns arrays x[0..n - 1] and w[0..n - 1] containing the abscissas and weights
/// of the n-point Gauss-Laguerre quadrature formula.The smallest abscissa is
/// returned in x[0], the largest in x[n - 1].
/// </summary>
/// <param name="x"></param>
/// <param name="w"></param>
/// <param name="alf"></param>
/// <exception cref="Exception"></exception>
public static void gaulag(double[] x, double[] w, double alf)
{
const int MAXIT = 10;
const double EPS = 1.0e-14;
double p2 = 0.0;
double pp = 0.0;
double z = 0.0;
int n = x.Length;
for (int i = 0; i < n; i++)
{
if (i == 0)
{
z = (1.0 + alf) * (3.0 + 0.92 * alf) / (1.0 + 2.4 * n + 1.8 * alf);
}
else if (i == 1)
{
z += (15.0 + 6.25 * alf) / (1.0 + 0.9 * alf + 2.5 * n);
}
else
{
int ai = i - 1;
z += ((1.0 + 2.55 * ai) / (1.9 * ai) + 1.26 * ai * alf / (1.0 + 3.5 * ai)) * (z - x[i - 2]) / (1.0 + 0.3 * alf);
}
int its = 0;
for (; its < MAXIT; its++)
{
double p1 = 1.0;
p2 = 0.0;
for (int j = 0; j < n; j++)
{
double p3 = p2;
p2 = p1;
p1 = ((2 * j + 1 + alf - z) * p2 - (j + alf) * p3) / (j + 1);
}
pp = (n * p1 - (n + alf) * p2) / z;
double z1 = z;
z = z1 - p1 / pp;
if (Math.Abs(z - z1) <= EPS)
{
break;
}
}
if (its >= MAXIT)
{
throw new Exception("too many iterations in gaulag");
}
x[i] = z;
w[i] = -Math.Exp(Globals.gammln(alf + n) - Globals.gammln((double)n)) / (pp * n * p2);
}
}
/// <summary>
/// This routine returns arrays x[0..n - 1] and w[0..n - 1] containing the
/// abscissas and weights of the n-point Gauss-Hermite quadrature formula.
/// The largest abscissa is returned in x[0], the most negative in x[n - 1].
/// </summary>
/// <param name="x"></param>
/// <param name="w"></param>
/// <exception cref="Exception"></exception>
public static void gauher(double[] x, double[] w)
{
const double EPS = 1.0e-14;
const double PIM4 = 0.7511255444649425;
const int MAXIT = 10;
double p2 = 0.0;
double pp = 0.0;
double z = 0.0;
int n = x.Length;
int m = (n + 1) / 2;
for (int i = 0; i < m; i++)
{
if (i == 0)
{
z = Math.Sqrt((double)(2 * n + 1)) - 1.85575 * Math.Pow((double)(2 * n + 1), -0.16667);
}
else if (i == 1)
{
z -= 1.14 * Math.Pow((double)n, 0.426) / z;
}
else if (i == 2)
{
z = 1.86 * z - 0.86 * x[0];
}
else if (i == 3)
{
z = 1.91 * z - 0.91 * x[1];
}
else
{
z = 2.0 * z - x[i - 2];
}
int its = 0;
for (; its < MAXIT; its++)
{
double p1 = PIM4;
p2 = 0.0;
for (int j = 0; j < n; j++)
{
double p3 = p2;
p2 = p1;
p1 = z * Math.Sqrt(2.0 / (j + 1)) * p2 - Math.Sqrt((double)j / (j + 1)) * p3;
}
pp = Math.Sqrt((double)(2 * n)) * p2;
double z1 = z;
z = z1 - p1 / pp;
if (Math.Abs(z - z1) <= EPS)
{
break;
}
}
if (its >= MAXIT)
{
throw new Exception("too many iterations in gauher");
}
x[i] = z;
x[n - 1 - i] = -z;
w[i] = 2.0 / (pp * pp);
w[n - 1 - i] = w[i];
}
}
/// <summary>
/// Given alf and bet, the parameters a and b of the Jacobi polynomials, this
/// routine returns arrays x[0..n - 1] and w[0..n - 1] containing the abscissas and
/// weights of the n-point Gauss-Jacobi quadrature formula.The largest
/// abscissa is returned in x[0], the smallest in x[n - 1].
/// </summary>
/// <param name="x"></param>
/// <param name="w"></param>
/// <param name="alf"></param>
/// <param name="bet"></param>
/// <exception cref="Exception"></exception>
public static void gaujac(double[] x, double[] w, double alf, double bet)
{
const int MAXIT = 10;
const double EPS = 1.0e-14;
double p2 = 0.0;
double pp = 0.0;
double temp = 0.0;
double z = 0.0;
int n = x.Length;
for (int i = 0; i < n; i++)
{
if (i == 0)
{
double an = alf / n;
double bn = bet / n;
double r1 = (1.0 + alf) * (2.78 / (4.0 + n * n) + 0.768 * an / n);
double r2 = 1.0 + 1.48 * an + 0.96 * bn + 0.452 * an * an + 0.83 * an * bn;
z = 1.0 - r1 / r2;
}
else if (i == 1)
{
double r1 = (4.1 + alf) / ((1.0 + alf) * (1.0 + 0.156 * alf));
double r2 = 1.0 + 0.06 * (n - 8.0) * (1.0 + 0.12 * alf) / n;
double r3 = 1.0 + 0.012 * bet * (1.0 + 0.25 * Math.Abs(alf)) / n;
z -= (1.0 - z) * r1 * r2 * r3;
}
else if (i == 2)
{
double r1 = (1.67 + 0.28 * alf) / (1.0 + 0.37 * alf);
double r2 = 1.0 + 0.22 * (n - 8.0) / n;
double r3 = 1.0 + 8.0 * bet / ((6.28 + bet) * n * n);
z -= (x[0] - z) * r1 * r2 * r3;
}
else if (i == n - 2)
{
double r1 = (1.0 + 0.235 * bet) / (0.766 + 0.119 * bet);
double r2 = 1.0 / (1.0 + 0.639 * (n - 4.0) / (1.0 + 0.71 * (n - 4.0)));
double r3 = 1.0 / (1.0 + 20.0 * alf / ((7.5 + alf) * n * n));
z += (z - x[n - 4]) * r1 * r2 * r3;
}
else if (i == n - 1)
{
double r1 = (1.0 + 0.37 * bet) / (1.67 + 0.28 * bet);
double r2 = 1.0 / (1.0 + 0.22 * (n - 8.0) / n);
double r3 = 1.0 / (1.0 + 8.0 * alf / ((6.28 + alf) * n * n));
z += (z - x[n - 3]) * r1 * r2 * r3;
}
else
{
z = 3.0 * x[i - 1] - 3.0 * x[i - 2] + x[i - 3];
}
double alfbet = alf + bet;
int its = 1;
for (; its <= MAXIT; its++)
{
temp = 2.0 + alfbet;
double p1 = (alf - bet + temp * z) / 2.0;
p2 = 1.0;
for (int j = 2; j <= n; j++)
{
double p3 = p2;
p2 = p1;
temp = 2 * j + alfbet;
double a = 2 * j * (j + alfbet) * (temp - 2.0);
double b = (temp - 1.0) * (alf * alf - bet * bet + temp * (temp - 2.0) * z);
double c = 2.0 * (j - 1 + alf) * (j - 1 + bet) * temp;
p1 = (b * p2 - c * p3) / a;
}
pp = (n * (alf - bet - temp * z) * p1 + 2.0 * (n + alf) * (n + bet) * p2) / (temp * (1.0 - z * z));
double z1 = z;
z = z1 - p1 / pp;
if (Math.Abs(z - z1) <= EPS)
{
break;
}
}
if (its > MAXIT)
{
throw new Exception("too many iterations in gaujac");
}
x[i] = z;
w[i] = Math.Exp(Globals.gammln(alf + n) + Globals.gammln(bet + n) - Globals.gammln(n + 1.0) - Globals.gammln(n + alfbet + 1.0)) * temp * Math.Pow(2.0, alfbet) / (pp * p2);
}
}
/// <summary>
/// Computes the abscissas and weights for a Gaussian quadrature formula from
/// the Jacobi matrix.On input, a[0..n - 1] and b[0..n - 1] are the coefficients
/// of the recurrence relation for the set of monic orthogonal polynomials. The
/// quantity u0= S(a, b)W(x) dx is input as amu0.The abscissas x[0..n - 1] are
/// returned in descending order, with the corresponding weights in w[0..n - 1].
/// The arrays a and b are modified.Execution can be speeded up by modifying
/// tqli and eigsrt to compute only the zeroth component of each eigenvector.
/// </summary>
/// <param name="a"></param>
/// <param name="b"></param>
/// <param name="amu0"></param>
/// <param name="x"></param>
/// <param name="w"></param>
public static void gaucof(double[] a, double[] b, double amu0, double[] x, double[] w)
{
int n = a.Length;
for (int i = 0; i < n; i++)
{
if (i != 0)
{
b[i] = Math.Sqrt(b[i]);
}
}
Symmeig sym = new Symmeig(a, b);
for (int i = 0; i < n; i++)
{
x[i] = sym.d[i];
w[i] = amu0 * sym.z[0, i] * sym.z[0, i];
}
}
/// <summary>
/// Computes the abscissas and weights for a Gauss-Radau quadrature formula.On
/// input, a[0..n - 1] and b[0..n - 1] are the coefficients of the recurrence
/// relation for the set of monic orthogo- nal polynomials corresponding to the
/// weight function. (b[0] is not referenced.) The quantity u0=S(a, b)W(x)dx is
/// input as amu0.x1 is input as either endpoint of the interval.The
/// abscissas x[0..n - 1] are returned in descending order, with the
/// corresponding weights in w[0..n - 1]. The arrays a and b are modified.
/// </summary>
/// <param name="a"></param>
/// <param name="b"></param>
/// <param name="amu0"></param>
/// <param name="x1"></param>
/// <param name="x"></param>
/// <param name="w"></param>
public static void radau(double[] a, double[] b, double amu0, double x1, double[] x, double[] w)
{
int n = a.Length;
if (n == 1)
{
x[0] = x1;
w[0] = amu0;
}
else
{
double p = x1 - a[0];
double pm1 = 1.0;
double p1 = p;
for (int i = 1; i < n - 1; i++)
{
p = (x1 - a[i]) * p1 - b[i] * pm1;
pm1 = p1;
p1 = p;
}
a[n - 1] = x1 - b[n - 1] * pm1 / p;
gaucof(a, b, amu0, x, w);
}
}
/// <summary>
/// Computes the abscissas and weights for a Gauss-Lobatto quadrature formula.
/// On input, the vectors a[0..n - 1] and b[0..n - 1] are the coefficients of the
/// recurrence relation for the set of monic orthogonal polynomials
/// corresponding to the weight function. (b[0] is not referenced.) The
/// quantity u0=S(a, b)W(x)dx is input as amu0.x1 amd xn are input as the
/// endpoints of the interval.The abscissas x[0..n - 1] are returned in
/// descending order, with the corresponding weights in w[0..n - 1]. The arrays a
/// and b are modified.
/// </summary>
/// <param name="a"></param>
/// <param name="b"></param>
/// <param name="amu0"></param>
/// <param name="x1"></param>
/// <param name="xn"></param>
/// <param name="x"></param>
/// <param name="w"></param>
/// <exception cref="Exception"></exception>
public static void lobatto(double[] a, double[] b, double amu0, double x1, double xn, double[] x, double[] w)
{
int n = a.Length;
if (n <= 1)
{
throw new Exception("n must be bigger than 1 in lobatto");
}
double pl = x1 - a[0];
double pr = xn - a[0];
double pm1l = 1.0;
double pm1r = 1.0;
double p1l = pl;
double p1r = pr;
for (int i = 1; i < n - 1; i++)
{
pl = (x1 - a[i]) * p1l - b[i] * pm1l;
pr = (xn - a[i]) * p1r - b[i] * pm1r;
pm1l = p1l;
pm1r = p1r;
p1l = pl;
p1r = pr;
}
double det = pl * pm1r - pr * pm1l;
a[n - 1] = (x1 * pl * pm1r - xn * pr * pm1l) / det;
b[n - 1] = (xn - x1) * pl * pr / det;
gaucof(a, b, amu0, x, w);
}
}
}
2 代码格式
using System;
namespace Legalsoft.Truffer
{
public class GaussianWeights
{
private static double[] x = {
0.1488743389816312, 0.4333953941292472, 0.6794095682990244,
0.8650633666889845, 0.9739065285171717
};
private static double[] w = {
0.2955242247147529, 0.2692667193099963, 0.2190863625159821,
0.1494513491505806, 0.0666713443086881
};
public GaussianWeights()
{
}
/*
public static double qgaus(UniVarRealValueFun f1, double a, double b)
{
//func = f1;
return qgaus(a,b);
}
*/
/// <summary>
/// Returns the integral of the function or functor func between a and b, by
/// ten-point Gauss-Legendre integration: the function is evaluated exactly
/// ten times at interior points in the range of integration.
/// </summary>
/// <param name="a"></param>
/// <param name="b"></param>
/// <returns></returns>
public static double qgaus(UniVarRealValueFun func, double a, double b)
{
double xm = 0.5 * (b + a);
double xr = 0.5 * (b - a);
double s = 0;
for (int j = 0; j < 5; j++)
{
double dx = xr * x[j];
s += w[j] * (func.funk(xm + dx) + func.funk(xm - dx));
}
return s *= xr;
}
/// <summary>
/// Given the lower and upper limits of integration x1 and x2, this routine
/// returns arrays x[0..n - 1] and w[0..n - 1] of length n, containing the
/// abscissas and weights of the Gauss-Legendre n-point quadrature formula.
/// </summary>
/// <param name="x1"></param>
/// <param name="x2"></param>
/// <param name="x"></param>
/// <param name="w"></param>
public static void gauleg(double x1, double x2, double[] x, double[] w)
{
const double EPS = 1.0e-14;
double z1;
double pp;
int n = x.Length;
int m = (n + 1) / 2;
double xm = 0.5 * (x2 + x1);
double xl = 0.5 * (x2 - x1);
for (int i = 0; i < m; i++)
{
double z = Math.Cos(3.141592654 * (i + 0.75) / (n + 0.5));
do
{
double p1 = 1.0;
double p2 = 0.0;
for (int j = 0; j < n; j++)
{
double p3 = p2;
p2 = p1;
p1 = ((2.0 * j + 1.0) * z * p2 - j * p3) / (j + 1);
}
pp = n * (z * p1 - p2) / (z * z - 1.0);
z1 = z;
z = z1 - p1 / pp;
} while (Math.Abs(z - z1) > EPS);
x[i] = xm - xl * z;
x[n - 1 - i] = xm + xl * z;
w[i] = 2.0 * xl / ((1.0 - z * z) * pp * pp);
w[n - 1 - i] = w[i];
}
}
/// <summary>
/// Given alf, the parameter a of the Laguerre polynomials, this routine
/// returns arrays x[0..n - 1] and w[0..n - 1] containing the abscissas and weights
/// of the n-point Gauss-Laguerre quadrature formula.The smallest abscissa is
/// returned in x[0], the largest in x[n - 1].
/// </summary>
/// <param name="x"></param>
/// <param name="w"></param>
/// <param name="alf"></param>
/// <exception cref="Exception"></exception>
public static void gaulag(double[] x, double[] w, double alf)
{
const int MAXIT = 10;
const double EPS = 1.0e-14;
double p2 = 0.0;
double pp = 0.0;
double z = 0.0;
int n = x.Length;
for (int i = 0; i < n; i++)
{
if (i == 0)
{
z = (1.0 + alf) * (3.0 + 0.92 * alf) / (1.0 + 2.4 * n + 1.8 * alf);
}
else if (i == 1)
{
z += (15.0 + 6.25 * alf) / (1.0 + 0.9 * alf + 2.5 * n);
}
else
{
int ai = i - 1;
z += ((1.0 + 2.55 * ai) / (1.9 * ai) + 1.26 * ai * alf / (1.0 + 3.5 * ai)) * (z - x[i - 2]) / (1.0 + 0.3 * alf);
}
int its = 0;
for (; its < MAXIT; its++)
{
double p1 = 1.0;
p2 = 0.0;
for (int j = 0; j < n; j++)
{
double p3 = p2;
p2 = p1;
p1 = ((2 * j + 1 + alf - z) * p2 - (j + alf) * p3) / (j + 1);
}
pp = (n * p1 - (n + alf) * p2) / z;
double z1 = z;
z = z1 - p1 / pp;
if (Math.Abs(z - z1) <= EPS)
{
break;
}
}
if (its >= MAXIT)
{
throw new Exception("too many iterations in gaulag");
}
x[i] = z;
w[i] = -Math.Exp(Globals.gammln(alf + n) - Globals.gammln((double)n)) / (pp * n * p2);
}
}
/// <summary>
/// This routine returns arrays x[0..n - 1] and w[0..n - 1] containing the
/// abscissas and weights of the n-point Gauss-Hermite quadrature formula.
/// The largest abscissa is returned in x[0], the most negative in x[n - 1].
/// </summary>
/// <param name="x"></param>
/// <param name="w"></param>
/// <exception cref="Exception"></exception>
public static void gauher(double[] x, double[] w)
{
const double EPS = 1.0e-14;
const double PIM4 = 0.7511255444649425;
const int MAXIT = 10;
double p2 = 0.0;
double pp = 0.0;
double z = 0.0;
int n = x.Length;
int m = (n + 1) / 2;
for (int i = 0; i < m; i++)
{
if (i == 0)
{
z = Math.Sqrt((double)(2 * n + 1)) - 1.85575 * Math.Pow((double)(2 * n + 1), -0.16667);
}
else if (i == 1)
{
z -= 1.14 * Math.Pow((double)n, 0.426) / z;
}
else if (i == 2)
{
z = 1.86 * z - 0.86 * x[0];
}
else if (i == 3)
{
z = 1.91 * z - 0.91 * x[1];
}
else
{
z = 2.0 * z - x[i - 2];
}
int its = 0;
for (; its < MAXIT; its++)
{
double p1 = PIM4;
p2 = 0.0;
for (int j = 0; j < n; j++)
{
double p3 = p2;
p2 = p1;
p1 = z * Math.Sqrt(2.0 / (j + 1)) * p2 - Math.Sqrt((double)j / (j + 1)) * p3;
}
pp = Math.Sqrt((double)(2 * n)) * p2;
double z1 = z;
z = z1 - p1 / pp;
if (Math.Abs(z - z1) <= EPS)
{
break;
}
}
if (its >= MAXIT)
{
throw new Exception("too many iterations in gauher");
}
x[i] = z;
x[n - 1 - i] = -z;
w[i] = 2.0 / (pp * pp);
w[n - 1 - i] = w[i];
}
}
/// <summary>
/// Given alf and bet, the parameters a and b of the Jacobi polynomials, this
/// routine returns arrays x[0..n - 1] and w[0..n - 1] containing the abscissas and
/// weights of the n-point Gauss-Jacobi quadrature formula.The largest
/// abscissa is returned in x[0], the smallest in x[n - 1].
/// </summary>
/// <param name="x"></param>
/// <param name="w"></param>
/// <param name="alf"></param>
/// <param name="bet"></param>
/// <exception cref="Exception"></exception>
public static void gaujac(double[] x, double[] w, double alf, double bet)
{
const int MAXIT = 10;
const double EPS = 1.0e-14;
double p2 = 0.0;
double pp = 0.0;
double temp = 0.0;
double z = 0.0;
int n = x.Length;
for (int i = 0; i < n; i++)
{
if (i == 0)
{
double an = alf / n;
double bn = bet / n;
double r1 = (1.0 + alf) * (2.78 / (4.0 + n * n) + 0.768 * an / n);
double r2 = 1.0 + 1.48 * an + 0.96 * bn + 0.452 * an * an + 0.83 * an * bn;
z = 1.0 - r1 / r2;
}
else if (i == 1)
{
double r1 = (4.1 + alf) / ((1.0 + alf) * (1.0 + 0.156 * alf));
double r2 = 1.0 + 0.06 * (n - 8.0) * (1.0 + 0.12 * alf) / n;
double r3 = 1.0 + 0.012 * bet * (1.0 + 0.25 * Math.Abs(alf)) / n;
z -= (1.0 - z) * r1 * r2 * r3;
}
else if (i == 2)
{
double r1 = (1.67 + 0.28 * alf) / (1.0 + 0.37 * alf);
double r2 = 1.0 + 0.22 * (n - 8.0) / n;
double r3 = 1.0 + 8.0 * bet / ((6.28 + bet) * n * n);
z -= (x[0] - z) * r1 * r2 * r3;
}
else if (i == n - 2)
{
double r1 = (1.0 + 0.235 * bet) / (0.766 + 0.119 * bet);
double r2 = 1.0 / (1.0 + 0.639 * (n - 4.0) / (1.0 + 0.71 * (n - 4.0)));
double r3 = 1.0 / (1.0 + 20.0 * alf / ((7.5 + alf) * n * n));
z += (z - x[n - 4]) * r1 * r2 * r3;
}
else if (i == n - 1)
{
double r1 = (1.0 + 0.37 * bet) / (1.67 + 0.28 * bet);
double r2 = 1.0 / (1.0 + 0.22 * (n - 8.0) / n);
double r3 = 1.0 / (1.0 + 8.0 * alf / ((6.28 + alf) * n * n));
z += (z - x[n - 3]) * r1 * r2 * r3;
}
else
{
z = 3.0 * x[i - 1] - 3.0 * x[i - 2] + x[i - 3];
}
double alfbet = alf + bet;
int its = 1;
for (; its <= MAXIT; its++)
{
temp = 2.0 + alfbet;
double p1 = (alf - bet + temp * z) / 2.0;
p2 = 1.0;
for (int j = 2; j <= n; j++)
{
double p3 = p2;
p2 = p1;
temp = 2 * j + alfbet;
double a = 2 * j * (j + alfbet) * (temp - 2.0);
double b = (temp - 1.0) * (alf * alf - bet * bet + temp * (temp - 2.0) * z);
double c = 2.0 * (j - 1 + alf) * (j - 1 + bet) * temp;
p1 = (b * p2 - c * p3) / a;
}
pp = (n * (alf - bet - temp * z) * p1 + 2.0 * (n + alf) * (n + bet) * p2) / (temp * (1.0 - z * z));
double z1 = z;
z = z1 - p1 / pp;
if (Math.Abs(z - z1) <= EPS)
{
break;
}
}
if (its > MAXIT)
{
throw new Exception("too many iterations in gaujac");
}
x[i] = z;
w[i] = Math.Exp(Globals.gammln(alf + n) + Globals.gammln(bet + n) - Globals.gammln(n + 1.0) - Globals.gammln(n + alfbet + 1.0)) * temp * Math.Pow(2.0, alfbet) / (pp * p2);
}
}
/// <summary>
/// Computes the abscissas and weights for a Gaussian quadrature formula from
/// the Jacobi matrix.On input, a[0..n - 1] and b[0..n - 1] are the coefficients
/// of the recurrence relation for the set of monic orthogonal polynomials. The
/// quantity u0= S(a, b)W(x) dx is input as amu0.The abscissas x[0..n - 1] are
/// returned in descending order, with the corresponding weights in w[0..n - 1].
/// The arrays a and b are modified.Execution can be speeded up by modifying
/// tqli and eigsrt to compute only the zeroth component of each eigenvector.
/// </summary>
/// <param name="a"></param>
/// <param name="b"></param>
/// <param name="amu0"></param>
/// <param name="x"></param>
/// <param name="w"></param>
public static void gaucof(double[] a, double[] b, double amu0, double[] x, double[] w)
{
int n = a.Length;
for (int i = 0; i < n; i++)
{
if (i != 0)
{
b[i] = Math.Sqrt(b[i]);
}
}
Symmeig sym = new Symmeig(a, b);
for (int i = 0; i < n; i++)
{
x[i] = sym.d[i];
w[i] = amu0 * sym.z[0, i] * sym.z[0, i];
}
}
/// <summary>
/// Computes the abscissas and weights for a Gauss-Radau quadrature formula.On
/// input, a[0..n - 1] and b[0..n - 1] are the coefficients of the recurrence
/// relation for the set of monic orthogo- nal polynomials corresponding to the
/// weight function. (b[0] is not referenced.) The quantity u0=S(a, b)W(x)dx is
/// input as amu0.x1 is input as either endpoint of the interval.The
/// abscissas x[0..n - 1] are returned in descending order, with the
/// corresponding weights in w[0..n - 1]. The arrays a and b are modified.
/// </summary>
/// <param name="a"></param>
/// <param name="b"></param>
/// <param name="amu0"></param>
/// <param name="x1"></param>
/// <param name="x"></param>
/// <param name="w"></param>
public static void radau(double[] a, double[] b, double amu0, double x1, double[] x, double[] w)
{
int n = a.Length;
if (n == 1)
{
x[0] = x1;
w[0] = amu0;
}
else
{
double p = x1 - a[0];
double pm1 = 1.0;
double p1 = p;
for (int i = 1; i < n - 1; i++)
{
p = (x1 - a[i]) * p1 - b[i] * pm1;
pm1 = p1;
p1 = p;
}
a[n - 1] = x1 - b[n - 1] * pm1 / p;
gaucof(a, b, amu0, x, w);
}
}
/// <summary>
/// Computes the abscissas and weights for a Gauss-Lobatto quadrature formula.
/// On input, the vectors a[0..n - 1] and b[0..n - 1] are the coefficients of the
/// recurrence relation for the set of monic orthogonal polynomials
/// corresponding to the weight function. (b[0] is not referenced.) The
/// quantity u0=S(a, b)W(x)dx is input as amu0.x1 amd xn are input as the
/// endpoints of the interval.The abscissas x[0..n - 1] are returned in
/// descending order, with the corresponding weights in w[0..n - 1]. The arrays a
/// and b are modified.
/// </summary>
/// <param name="a"></param>
/// <param name="b"></param>
/// <param name="amu0"></param>
/// <param name="x1"></param>
/// <param name="xn"></param>
/// <param name="x"></param>
/// <param name="w"></param>
/// <exception cref="Exception"></exception>
public static void lobatto(double[] a, double[] b, double amu0, double x1, double xn, double[] x, double[] w)
{
int n = a.Length;
if (n <= 1)
{
throw new Exception("n must be bigger than 1 in lobatto");
}
double pl = x1 - a[0];
double pr = xn - a[0];
double pm1l = 1.0;
double pm1r = 1.0;
double p1l = pl;
double p1r = pr;
for (int i = 1; i < n - 1; i++)
{
pl = (x1 - a[i]) * p1l - b[i] * pm1l;
pr = (xn - a[i]) * p1r - b[i] * pm1r;
pm1l = p1l;
pm1r = p1r;
p1l = pl;
p1r = pr;
}
double det = pl * pm1r - pr * pm1l;
a[n - 1] = (x1 * pl * pm1r - xn * pr * pm1l) / det;
b[n - 1] = (xn - x1) * pl * pr / det;
gaucof(a, b, amu0, x, w);
}
}
}
