|
说明： 代码较为简洁没有过多的说明，如有不明白之处可查阅相关最小二乘法计算步骤资料和求解线性方程组的资料。另外该方法只能实现二元N次拟合，多元方程不适用。 以下是最小二乘法类的实现： public classMatrixEquation { private double[,] gaussMatrix; private int coe; public MatrixEquation() { } public MatrixEquation(double[] arrX,double[] arrY, int n) { coe = n; gaussMatrix=GetGauss(GetXPowSum(arrX, n), GetXPowYSum(arrX, arrY, n), n); } public double[,] GetGaussMatrix() { return gaussMatrix; } public double[] GetResult() { return ComputeGauss(gaussMatrix,coe); } /// <summary> /// 计算获取x散点的幂次和数组 /// </summary> /// <paramname="arrX">x散点序列</param> /// <param name="n">函数拟合次数</param> /// <returns></returns> protected double[] GetXPowSum(double[]arrX, int n) { int m = arrX.Length;//X散点的个数 double[] xPow = new double[2 * n +1]; //存储X散点的幂次值 for (int i = 0; i < xPow.Length;i++) { if (i == 0) { xPow = m; } else { //计算x的i次方和 double max = 0; for (int j = 0; j < m;j++) { if (arrX[j] == 0) max = max + 1; else max = max +Math.Pow(arrX[j], i); } xPow =Math.Round(max,4); } } return xPow; } /// <summary> /// 计算获取xy的幂次和序列 /// </summary> /// <paramname="arrX">x散点序列</param> /// <paramname="arrY">y散点序列</param> /// <param name="n">拟合曲线次数</param> /// <returns></returns> protected double[] GetXPowYSum(double[]arrX, double[] arrY, int n) { int m = arrX.Length;//X散点的个数 double[] xyPow = new double[n + 1];//仓储X散点的幂次值 for (int i = 0; i <xyPow.Length; i++) { //计算xy的i次方和 double max = 0; for (int j = 0; j < m; j++) { if (arrX[j] == 0) max = max + 1; else max = max +Math.Pow(arrX[j], i) * arrY[j]; } xyPow =Math.Round( max,4); } return xyPow; } /// <summary> /// 获取高斯矩阵(增广矩阵) /// </summary> /// <paramname="arrX">X的幂次和</param> /// <paramname="arrXY">XY的幂次和</param> /// <param name="n">拟合曲线次数</param> /// <returns></returns> protected double[,] GetGauss(double[]arrX, double[] arrXY, int n) { double[,] gauss = new double[n+1, n+ 2]; for (int i = 0; i < n + 1; i++) { int j; int m = i; for (j = 0; j < n + 1; j++) { gauss[i, j] = arrX[m]; m++; } gauss[i,j] = arrXY; } return gauss; } /// <summary> /// 求解拟合曲线的系数 /// </summary> /// <paramname="gauss">线性方程的增广矩阵</param> /// <param name="n">方程次数</param> /// <returns></returns> protected double[]ComputeGauss(double[,] gauss, int n) { double[] a = new double[n + 1]; double s; int matrixLine = n + 1; for (int i = 0; i < n + 1; i++) a = 0; //循环每列 for (int j = 0; j < matrixLine;j++) { //每列J行以后的绝对值最大值 double max = 0; int k = j; for (int i = j; i <matrixLine; i++) { if (Math.Abs(gauss[i, j])> max) { max = gauss[i, j]; k = i; } } //判断j行否为最大值行 若不是将j行调换为最大值行 if (k != j) { double temp; for (int m = j; m <matrixLine+1; m++) { temp = gauss[j, m]; gauss[j, m] = gauss[k,m]; gauss[k, m] = temp; } } if (max == 0) { //奇异矩阵无解 return a; } //进行初等行变换得到上三角矩阵 for (int i = j + 1; i <matrixLine; i++) { s = gauss[i, j]; for (int m = j; m <matrixLine + 1; m++) { gauss[i, m]=Math.Round( gauss[i, m] - gauss[j, m] * s / gauss[j, j],6); } } } //根据倒推方式一次计算现行方程的解 for (int i = matrixLine - 1; i>= 0; i--) { s = 0; for (int j = i + 1; j <matrixLine; j++) { s += gauss[i, j] * a[j]; } a =Math.Round( (gauss[i,matrixLine] - s) / gauss[i, i],6); } //返回方程的解 即拟合曲线的系数 return a; } }
| 
IP卡
狗仔卡
发表于 2015-12-3 21:32:00
提升卡
置顶卡
喧嚣卡
变色卡
千斤顶
显身卡