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

C#，数值计算——函数计算，Ratfn的计算方法与源程序

news 2025/7/27 21:42:19

1 文本格式

using System;

namespace Legalsoft.Truffer
{
public class Ratfn
{
private double[] cofs { get; set; }
private int nn { get; set; }
private int dd { get; set; }

public Ratfn(double[] num, double[] den)
{
this.cofs = new double[num.Length + den.Length - 1];
this.nn = num.Length;
this.dd = den.Length;

for (int j = 0; j < nn; j++)
{
cofs[j] = num[j] / den[0];
}
for (int j = 1; j < dd; j++)
{
cofs[j + nn - 1] = den[j] / den[0];
}
}

public Ratfn(double[] coffs, int n, int d)
{
this.cofs = Globals.CopyFrom(coffs);
this.nn = n;
this.dd = d;
}

public double get(double x)
{
double sumn = 0.0;
double sumd = 0.0;
for (int j = nn - 1; j >= 0; j--)
{
sumn = sumn * x + cofs[j];
}
for (int j = nn + dd - 2; j >= nn; j--)
{
sumd = sumd * x + cofs[j];
}
return sumn / (1.0 + x * sumd);
}

public static Ratfn pade(double[] cof)
{
int n = (cof.Length - 1) / 2;
double[,] q = new double[n, n];
double[,] qlu = new double[n, n];
double[] x = new double[n];
double[] y = new double[n];
for (int j = 0; j < n; j++)
{
y[j] = cof[n + j + 1];
for (int k = 0; k < n; k++)
{
q[j, k] = cof[j - k + n];
}
}

LUdcmp lu = new LUdcmp(q);
lu.solve( y, x);

for (int j = 0; j < 4; j++)
{
lu.mprove(y, x);
}
for (int k = 0; k < n; k++)
{
double sum = cof[k + 1];
for (int j = 0; j <= k; j++)
{
sum -= x[j] * cof[k - j];
}
y[k] = sum;
}

double[] num = new double[n + 1];
double[] denom = new double[n + 1];
num[0] = cof[0];
denom[0] = 1.0;
for (int j = 0; j < n; j++)
{
num[j + 1] = y[j];
denom[j + 1] = -x[j];
}
return new Ratfn(num, denom);
}

public static Ratfn ratlsq(UniVarRealValueFun fn, double a, double b, int mm, int kk, ref double dev)
{
const int NPFAC = 8;
const int MAXIT = 5;
const double BIG = 1.0e99;
const double PIO2 = 1.570796326794896619;

int ncof = mm + kk + 1;
int npt = NPFAC * ncof;
double[] bb = new double[npt];
double[] coff = new double[ncof];
double[] ee = new double[npt];
double[] fs = new double[npt];
double[] wt = new double[npt];
double[] xs = new double[npt];
double[,] u = new double[npt, ncof];
Ratfn ratbest = new Ratfn(coff, mm + 1, kk + 1);

dev = BIG;
for (int i = 0; i < npt; i++)
{
if (i < (npt / 2) - 1)
{
double hth = PIO2 * i / (npt - 1.0);
xs[i] = a + (b - a) * Globals.SQR(Math.Sin(hth));
}
else
{
double hth = PIO2 * (npt - i) / (npt - 1.0);
xs[i] = b - (b - a) * Globals.SQR(Math.Sin(hth));
}
fs[i] = fn.funk(xs[i]);
wt[i] = 1.0;
ee[i] = 1.0;
}
double e = 0.0;
for (int it = 0; it < MAXIT; it++)
{
for (int i = 0; i < npt; i++)
{
double power = wt[i];
bb[i] = power * (fs[i] + Globals.SIGN(e, ee[i]));
for (int j = 0; j < mm + 1; j++)
{
u[i, j] = power;
power *= xs[i];
}
power = -bb[i];
for (int j = mm + 1; j < ncof; j++)
{
power *= xs[i];
u[i, j] = power;
}
}
SVD svd = new SVD(u);
svd.solve(bb, coff);
double devmax = 0.0;
double sum = 0.0;
Ratfn rat = new Ratfn(coff, mm + 1, kk + 1);
for (int j = 0; j < npt; j++)
{
ee[j] = rat.get(xs[j]) - fs[j];
wt[j] = Math.Abs(ee[j]);
sum += wt[j];
if (wt[j] > devmax)
{
devmax = wt[j];
}
}
e = sum / npt;
if (devmax <= dev)
{
ratbest = rat;
//ratbest.CopyFrom(rat);
dev = devmax;
}
//Console.Write(" ratlsq iteration= ");
//Console.Write(it);
//Console.Write(" max error= ");
//Console.Write("{0,10}", devmax);
//Console.Write("{0}", "\n");
}
return ratbest;
}

}
}

2 代码格式

using System;namespace Legalsoft.Truffer
{public class Ratfn{private double[] cofs { get; set; }private int nn { get; set; }private int dd { get; set; }public Ratfn(double[] num, double[] den){this.cofs = new double[num.Length + den.Length - 1];this.nn = num.Length;this.dd = den.Length;for (int j = 0; j < nn; j++){cofs[j] = num[j] / den[0];}for (int j = 1; j < dd; j++){cofs[j + nn - 1] = den[j] / den[0];}}public Ratfn(double[] coffs, int n, int d){this.cofs = Globals.CopyFrom(coffs);this.nn = n;this.dd = d;}public double get(double x){double sumn = 0.0;double sumd = 0.0;for (int j = nn - 1; j >= 0; j--){sumn = sumn * x + cofs[j];}for (int j = nn + dd - 2; j >= nn; j--){sumd = sumd * x + cofs[j];}return sumn / (1.0 + x * sumd);}public static Ratfn pade(double[] cof){int n = (cof.Length - 1) / 2;double[,] q = new double[n, n];double[,] qlu = new double[n, n];double[] x = new double[n];double[] y = new double[n];for (int j = 0; j < n; j++){y[j] = cof[n + j + 1];for (int k = 0; k < n; k++){q[j, k] = cof[j - k + n];}}LUdcmp lu = new LUdcmp(q);lu.solve( y,  x);for (int j = 0; j < 4; j++){lu.mprove(y,  x);}for (int k = 0; k < n; k++){double sum = cof[k + 1];for (int j = 0; j <= k; j++){sum -= x[j] * cof[k - j];}y[k] = sum;}double[] num = new double[n + 1];double[] denom = new double[n + 1];num[0] = cof[0];denom[0] = 1.0;for (int j = 0; j < n; j++){num[j + 1] = y[j];denom[j + 1] = -x[j];}return new Ratfn(num, denom);}public static Ratfn ratlsq(UniVarRealValueFun fn, double a, double b, int mm, int kk, ref double dev){const int NPFAC = 8;const int MAXIT = 5;const double BIG = 1.0e99;const double PIO2 = 1.570796326794896619;int ncof = mm + kk + 1;int npt = NPFAC * ncof;double[] bb = new double[npt];double[] coff = new double[ncof];double[] ee = new double[npt];double[] fs = new double[npt];double[] wt = new double[npt];double[] xs = new double[npt];double[,] u = new double[npt, ncof];Ratfn ratbest = new Ratfn(coff, mm + 1, kk + 1);dev = BIG;for (int i = 0; i < npt; i++){if (i < (npt / 2) - 1){double hth = PIO2 * i / (npt - 1.0);xs[i] = a + (b - a) * Globals.SQR(Math.Sin(hth));}else{double hth = PIO2 * (npt - i) / (npt - 1.0);xs[i] = b - (b - a) * Globals.SQR(Math.Sin(hth));}fs[i] = fn.funk(xs[i]);wt[i] = 1.0;ee[i] = 1.0;}double e = 0.0;for (int it = 0; it < MAXIT; it++){for (int i = 0; i < npt; i++){double power = wt[i];bb[i] = power * (fs[i] + Globals.SIGN(e, ee[i]));for (int j = 0; j < mm + 1; j++){u[i, j] = power;power *= xs[i];}power = -bb[i];for (int j = mm + 1; j < ncof; j++){power *= xs[i];u[i, j] = power;}}SVD svd = new SVD(u);svd.solve(bb,  coff);double devmax = 0.0;double sum = 0.0;Ratfn rat = new Ratfn(coff, mm + 1, kk + 1);for (int j = 0; j < npt; j++){ee[j] = rat.get(xs[j]) - fs[j];wt[j] = Math.Abs(ee[j]);sum += wt[j];if (wt[j] > devmax){devmax = wt[j];}}e = sum / npt;if (devmax <= dev){ratbest = rat;//ratbest.CopyFrom(rat);dev = devmax;}//Console.Write(" ratlsq iteration= ");//Console.Write(it);//Console.Write("  max error= ");//Console.Write("{0,10}", devmax);//Console.Write("{0}", "\n");}return ratbest;}}
}

查看全文

http://www.lryc.cn/news/230321.html