C#，数值计算——用于积分函数与方法的Stiel类的计算方法与源程序
1 文本格式
using System;
namespace Legalsoft.Truffer
{
public class Stiel
{
public class pp : UniVarRealValueFun, RealValueFun
{
public Stiel st { get; set; } = null;
public pp()
{
}
public double funk(double[] x)
{
double pval = st.p(x[0]);
return pval * st.wt1(x[0], x[0]) * pval;
}
public double funk(double t)
{
double x = st.fx(t);
double pval = st.p(x);
return pval * st.wt2(x) * st.fdxdt(t) * pval;
}
}
public class ppx : UniVarRealValueFun, RealValueFun
{
public Stiel st { get; set; } = null;
public ppx()
{
}
public double funk(double[] x)
{
return st.ppfunc.funk(new double[] { x[0], x[1] }) * x[0];
}
public double funk(double t)
{
return st.ppfunc.funk(t) * st.fx(t);
}
}
public pp ppfunc { get; set; } = new pp();
public ppx ppxfunc { get; set; } = new ppx();
public int j { get; set; }
public int n { get; set; }
public double aa { get; set; }
public double bb { get; set; }
public double hmax { get; set; }
public double[] a { get; set; }
public double[] b;
public Quadrature s1 { get; set; }
public Quadrature s2 { get; set; }
public double p(double x)
{
double pval = 0.0;
double pj;
double pjm1;
if (j == 0)
{
return 1.0;
}
else
{
pjm1 = 0.0;
pj = 1.0;
for (int i = 0; i < j; i++)
{
pval = (x - a[i]) * pj - b[i] * pjm1;
pjm1 = pj;
pj = pval;
}
}
return pval;
}
public Stiel(int nn, double aaa, double bbb, double hmaxx)
{
this.n = nn;
this.aa = aaa;
this.bb = bbb;
this.hmax = hmaxx;
this.a = new double[nn];
this.b = new double[nn];
s1 = new DErule(ppfunc, aa, bb, hmax);
s2 = new DErule(ppxfunc, aa, bb, hmax);
}
public Stiel(int nn, double aaa, double bbb)
{
this.n = nn;
this.aa = aaa;
this.bb = bbb;
this.a = new double[nn];
this.b = new double[nn];
s1 = new Trapzd(ppfunc, aa, bb);
s2 = new Trapzd(ppxfunc, aa, bb);
}
public double quad(Quadrature s)
{
const double EPS = 3.0e-11;
double MACHEPS = float.Epsilon;
const int NMAX = 11;
double olds = 0.0;
double sum;
s.n = 0;
for (int i = 1; i <= NMAX; i++)
{
sum = s.next();
if (i > 3)
{
if (Math.Abs(sum - olds) <= EPS * Math.Abs(olds))
{
return sum;
}
}
if (i == NMAX)
{
if (Math.Abs(sum) <= MACHEPS && Math.Abs(olds) <= MACHEPS)
{
return 0.0;
}
}
olds = sum;
}
throw new Exception("no convergence in quad");
}
public void get_weights(double[] x, double[] w)
{
double amu0;
double c;
double oldc = 1.0;
if (n != x.Length)
{
throw new Exception("bad array size in Stiel");
}
for (int i = 0; i < n; i++)
{
j = i;
c = quad(s1);
b[i] = c / oldc;
a[i] = quad(s2) / c;
oldc = c;
}
amu0 = b[0];
GaussianWeights.gaucof(a, b, amu0, x, w);
}
public double wt1(double x, double del) { return -9999; }
public double wt2(double x) { return -9999; }
public double fx(double t) { return -9999; }
public double fdxdt(double t) { return -9999; }
}
}
2 代码格式
using System;namespace Legalsoft.Truffer
{public class Stiel{public class pp : UniVarRealValueFun, RealValueFun{public Stiel st { get; set; } = null;public pp(){}public double funk(double[] x){double pval = st.p(x[0]);return pval * st.wt1(x[0], x[0]) * pval;}public double funk(double t){double x = st.fx(t);double pval = st.p(x);return pval * st.wt2(x) * st.fdxdt(t) * pval;}}public class ppx : UniVarRealValueFun, RealValueFun{public Stiel st { get; set; } = null;public ppx(){}public double funk(double[] x){return st.ppfunc.funk(new double[] { x[0], x[1] }) * x[0];}public double funk(double t){return st.ppfunc.funk(t) * st.fx(t);}}public pp ppfunc { get; set; } = new pp();public ppx ppxfunc { get; set; } = new ppx();public int j { get; set; }public int n { get; set; }public double aa { get; set; }public double bb { get; set; }public double hmax { get; set; }public double[] a { get; set; }public double[] b;public Quadrature s1 { get; set; }public Quadrature s2 { get; set; }public double p(double x){double pval = 0.0;double pj;double pjm1;if (j == 0){return 1.0;}else{pjm1 = 0.0;pj = 1.0;for (int i = 0; i < j; i++){pval = (x - a[i]) * pj - b[i] * pjm1;pjm1 = pj;pj = pval;}}return pval;}public Stiel(int nn, double aaa, double bbb, double hmaxx){this.n = nn;this.aa = aaa;this.bb = bbb;this.hmax = hmaxx;this.a = new double[nn];this.b = new double[nn];s1 = new DErule(ppfunc, aa, bb, hmax);s2 = new DErule(ppxfunc, aa, bb, hmax);}public Stiel(int nn, double aaa, double bbb){this.n = nn;this.aa = aaa;this.bb = bbb;this.a = new double[nn];this.b = new double[nn];s1 = new Trapzd(ppfunc, aa, bb);s2 = new Trapzd(ppxfunc, aa, bb);}public double quad(Quadrature s){const double EPS = 3.0e-11;double MACHEPS = float.Epsilon;const int NMAX = 11;double olds = 0.0;double sum;s.n = 0;for (int i = 1; i <= NMAX; i++){sum = s.next();if (i > 3){if (Math.Abs(sum - olds) <= EPS * Math.Abs(olds)){return sum;}}if (i == NMAX){if (Math.Abs(sum) <= MACHEPS && Math.Abs(olds) <= MACHEPS){return 0.0;}}olds = sum;}throw new Exception("no convergence in quad");}public void get_weights(double[] x, double[] w){double amu0;double c;double oldc = 1.0;if (n != x.Length){throw new Exception("bad array size in Stiel");}for (int i = 0; i < n; i++){j = i;c = quad(s1);b[i] = c / oldc;a[i] = quad(s2) / c;oldc = c;}amu0 = b[0];GaussianWeights.gaucof(a, b, amu0, x, w);}public double wt1(double x, double del) { return -9999; }public double wt2(double x) { return -9999; }public double fx(double t) { return -9999; }public double fdxdt(double t) { return -9999; }}
}