C#数值计算算法编程线性方程组类LEquations的设计.doc

C#数值计算算法编程线性方程组类LEquations的设计namespace CSharpAlgorithm.Algorithm/* * 求解线性方程组的类LEquations * author 周长发 * version 1.0 */public class LEquations private MatrixmtxLECoef;/ 系数矩阵private Matrix mtxLEConst;/ 常数矩阵/* * 基本构造函数 */public LEquations()/* * 指定系数和常数构造函数 * * param mtxCoef - 指定的系数矩阵 * param mtxConst - 指定的常数矩阵 */public LEquations(Matrix mtxCoef, Matrix mtxConst)Init(mtxCoef, mtxConst);/* * 初始化函数 * * param mtxCoef - 指定的系数矩阵 * param mtxConst - 指定的常数矩阵 * return bool 型，初始化是否成功 */public bool Init(Matrix mtxCoef, Matrix mtxConst)if (mtxCoef.GetNumRows() != mtxConst.GetNumRows()return false;mtxLECoef = new Matrix(mtxCoef);mtxLEConst = new Matrix(mtxConst);return true;/* * 获取系数矩阵 * * return Matrix 型，返回系数矩阵 */public Matrix GetCoefMatrix()return mtxLECoef;/* * 获取常数矩阵 * * return Matrix 型，返回系数矩阵 */public Matrix GetConstMatrix()return mtxLEConst;/* * 获取方程个数 * * return int 型，返回方程组方程的个数 */public int GetNumEquations()return GetCoefMatrix().GetNumRows();/* * 获取未知数个数 * * return int 型，返回方程组未知数的个数 */public int GetNumUnknowns()return GetCoefMatrix().GetNumColumns();/* * 全选主元高斯消去法 * * param mtxResult - Matrix对象，返回方程组的解 * return bool 型，方程组求解是否成功 */public bool GetRootsetGauss(Matrix mtxResult) int l,k,i,j,nIs=0,p,q;double d,t;/ 方程组的属性，将常数矩阵赋给解矩阵mtxResult.SetValue(mtxLEConst);double pDataCoef = mtxLECoef.GetData();double pDataConst = mtxResult.GetData();int n = GetNumUnknowns();/ 临时缓冲区，存放列数int pnJs = new intn;/ 消元l=1;for (k=0;k=n-2;k+) d=0.0;for (i=k;i=n-1;i+)for (j=k;jd) d=t; pnJsk=j; nIs=i;if (d = 0.0) l=0;else if (pnJsk!=k)for (i=0;i=n-1;i+) p=i*n+k; q=i*n+pnJsk;t=pDataCoefp; pDataCoefp=pDataCoefq; pDataCoefq=t;if (nIs!=k) for (j=k;j=n-1;j+) p=k*n+j; q=nIs*n+j;t=pDataCoefp; pDataCoefp=pDataCoefq; pDataCoefq=t; t=pDataConstk; pDataConstk=pDataConstnIs; pDataConstnIs=t; / 求解失败if (l=0) return false; d=pDataCoefk*n+k;for (j=k+1;j=n-1;j+) p=k*n+j; pDataCoefp=pDataCoefp/d; pDataConstk=pDataConstk/d;for (i=k+1;i=n-1;i+) for (j=k+1;j=0;i-) t=0.0;for (j=i+1;j=0;k-)if (pnJsk!=k) t=pDataConstk; pDataConstk=pDataConstpnJsk; pDataConstpnJsk=t;return true;/* * 全选主元高斯约当消去法 * * param mtxResult - Matrix对象，返回方程组的解 * return bool 型，方程组求解是否成功 */public bool GetRootsetGaussJordan(Matrix mtxResult) int l,k,i,j,nIs=0,p,q;double d,t;/ 方程组的属性，将常数矩阵赋给解矩阵mtxResult.SetValue(mtxLEConst);double pDataCoef = mtxLECoef.GetData();double pDataConst = mtxResult.GetData();int n = GetNumUnknowns();int m = mtxLEConst.GetNumColumns();/ 临时缓冲区，存放变换的列数int pnJs = new intn;/ 消元l=1;for (k=0;k=n-1;k+) d=0.0;for (i=k;i=n-1;i+)for (j=k;jd) d=t; pnJsk=j; nIs=i;if (d+1.0=1.0) l=0;else if (pnJsk!=k)for (i=0;i=n-1;i+) p=i*n+k; q=i*n+pnJsk;t=pDataCoefp; pDataCoefp=pDataCoefq; pDataCoefq=t;if (nIs!=k) for (j=k;j=n-1;j+) p=k*n+j; q=nIs*n+j;t=pDataCoefp; pDataCoefp=pDataCoefq; pDataCoefq=t; for (j=0;j=m-1;j+) p=k*m+j; q=nIs*m+j;t=pDataConstp; pDataConstp=pDataConstq; pDataConstq=t; / 求解失败if (l=0) return false; d=pDataCoefk*n+k;for (j=k+1;j=n-1;j+) p=k*n+j; pDataCoefp=pDataCoefp/d; for (j=0;j=m-1;j+) p=k*m+j; pDataConstp=pDataConstp/d; for (j=k+1;j=n-1;j+)for (i=0;i=n-1;i+) p=i*n+j;if (i!=k)pDataCoefp=pDataCoefp-pDataCoefi*n+k*pDataCoefk*n+j;for (j=0;j=m-1;j+)for (i=0;i=0;k-)if (pnJsk!=k)for (j=0;j=m-1;j+) p=k*m+j; q=pnJsk*m+j;t=pDataConstp; pDataConstp=pDataConstq; pDataConstq=t;return true;/* * 复系数方程组的全选主元高斯消去法 * * param mtxCoefImag - 系数矩阵的虚部矩阵 * param mtxConstImag - 常数矩阵的虚部矩阵 * param mtxResult - Matrix对象，返回方程组解矩阵的实部矩阵 * param mtxResultImag - Matrix对象，返回方程组解矩阵的虚部矩阵 * return bool 型，方程组求解是否成功 */public bool GetRootsetGauss(Matrix mtxCoefImag, Matrix mtxConstImag, Matrix mtxResult, Matrix mtxResultImag) int l,k,i,j,nIs=0,u,v;double p,q,s,d;/ 方程组的属性，将常数矩阵赋给解矩阵mtxResult.SetValue(mtxLEConst);mtxResultImag.SetValue(mtxConstImag);double pDataCoef = mtxLECoef.GetData();double pDataConst = mtxResult.GetData();double pDataCoefImag = mtxCoefImag.GetData();double pDataConstImag = mtxResultImag.GetData();int n = GetNumUnknowns();int m = mtxLEConst.GetNumColumns();/ 临时缓冲区，存放变换的列数int pnJs = new intn; / 消元for (k=0;k=n-2;k+) d=0.0;for (i=k;i=n-1;i+)for (j=k;jd) d=p;pnJsk=j;nIs=i; / 求解失败if (d = 0.0)return false; if (nIs!=k) for (j=k;j=n-1;j+) u=k*n+j; v=nIs*n+j;p=pDataCoefu; pDataCoefu=pDataCoefv; pDataCoefv=p;p=pDataCoefImagu; pDataCoefImagu=pDataCoefImagv; pDataCoefImagv=p; p=pDataConstk; pDataConstk=pDataConstnIs; pDataConstnIs=p;p=pDataConstImagk; pDataConstImagk=pDataConstImagnIs; pDataConstImagnIs=p; if (pnJsk!=k)for (i=0;i=n-1;i+) u=i*n+k; v=i*n+pnJsk;p=pDataCoefu; pDataCoefu=pDataCoefv; pDataCoefv=p;p=pDataCoefImagu; pDataCoefImagu=pDataCoefImagv; pDataCoefImagv=p;v=k*n+k;for (j=k+1;j=n-1;j+) u=k*n+j;p=pDataCoefu*pDataCoefv; q=-pDataCoefImagu*pDataCoefImagv;s=(pDataCoefv-pDataCoefImagv)*(pDataCoefu+pDataCoefImagu);pDataCoefu=(p-q)/d; pDataCoefImagu=(s-p-q)/d; p=pDataConstk*pDataCoefv; q=-pDataConstImagk*pDataCoefImagv;s=(pDataCoefv-pDataCoefImagv)*(pDataConstk+pDataConstImagk);pDataConstk=(p-q)/d; pDataConstImagk=(s-p-q)/d;for (i=k+1;i=n-1;i+) u=i*n+k;for (j=k+1;j=0;i-)for (j=i+1;j=0;k-)if (pnJsk!=k) p=pDataConstk; pDataConstk=pDataConstpnJsk; pDataConstpnJsk=p;p=pDataConstImagk; pDataConstImagk=pDataConstImagpnJsk; pDataConstImagpnJsk=p;return true;/* * 复系数方程组的全选主元高斯约当消去法 * * param mtxCoefImag - 系数矩阵的虚部矩阵 * param mtxConstImag - 常数矩阵的虚部矩阵 * param mtxResult - Matrix对象，返回方程组解矩阵的实部矩阵 * param mtxResultImag - Matrix对象，返回方程组解矩阵的虚部矩阵 * return bool 型，方程组求解是否成功 */public bool GetRootsetGaussJordan(Matrix mtxCoefImag, Matrix mtxConstImag, Matrix mtxResult, Matrix mtxResultImag)int l,k,i,j,nIs=0,u,v;double p,q,s,d;/ 方程组的属性，将常数矩阵赋给解矩阵mtxResult.SetValue(mtxLEConst);mtxResultImag.SetValue(mtxConstImag);double pDataCoef = mtxLECoef.GetData();double pDataConst = mtxResult.GetData();double pDataCoefImag = mtxCoefImag.GetData();double pDataConstImag = mtxResultImag.GetData();int n = GetNumUnknowns();int m = mtxLEConst.GetNumColumns();/ 临时缓冲区，存放变换的列数int pnJs = new intn;/ 消元for (k=0;k=n-1;k+) d=0.0;for (i=k;i=n-1;i+)for (j=k;jd) d=p;pnJsk=j;nIs=i; / 求解失败if (d = 0.0)return false; if (nIs!=k) for (j=k;j=n-1;j+) u=k*n+j; v=nIs*n+j;p=pDataCoefu; pDataCoefu=pDataCoefv; pDataCoefv=p;p=pDataCoefImagu; pDataCoefImagu=pDataCoefImagv; pDataCoefImagv=p; for (j=0;j=m-1;j+) u=k*m+j; v=nIs*m+j;p=pDataConstu; pDataConstu=pDataConstv; pDataConstv=p;p=pDataConstImagu; pDataConstImagu=pDataConstImagv; pDataConstImagv=p; if (pnJsk!=k)for (i=0;i=n-1;i+) u=i*n+k; v=i*n+pnJsk;p=pDataCoefu; pDataCoefu=pDataCoefv; pDataCoefv=p;p=pDataCoefImagu; pDataCoefImagu=pDataCoefImagv; pDataCoefImagv=p;v=k*n+k;for (j=k+1;j=n-1;j+) u=k*n+j;p=pDataCoefu*pDataCoefv; q=-pDataCoefImagu*pDataCoefImagv;s=(pDataCoefv-pDataCoefImagv)*(pDataCoefu+pDataCoefImagu);pDataCoefu=(p-q)/d; pDataCoefImagu=(s-p-q)/d; for (j=0;j=m-1;j+) u=k*m+j;p=pDataConstu*pDataCoefv; q=-pDataConstImagu*pDataCoefImagv;s=(pDataCoefv-pDataCoefImagv)*(pDataConstu+pDataConstImagu);pDataConstu=(p-q)/d; pDataConstImagu=(s-p-q)/d; for (i=0;i=n-1;i+)if (i!=k) u=i*n+k;for (j=k+1;j=n-1;j+) v=k*n+j; l=i*n+j;p=pDataCoefu*pDataCoefv; q=pDataCoefImagu*pDataCoefImagv;s=(pDataCoefu+pDataCoefImagu)*(pDataCoefv+pDataCoefImagv);pDataCoefl=pDataCoefl-p+q;pDataCoefImagl=pDataCoefImagl-s+p+q; for (j=0;j=0;k-)if (pnJsk!=k)for (j=0;j=m-1;j+) u=k*m+j;v=pnJsk*m+j;p=pDataConstu; pDataConstu=pDataConstv; pDataConstv=p;p=pDataConstImagu; pDataConstImagu=pDataConstImagv; pDataConstImagv=p;return true;/* * 求解三对角线方程组的追赶法 * * param mtxResult - Matrix对象，返回方程组解矩阵 * return bool 型，方程组求解是否成功 */public bool GetRootsetTriDiagonal(Matrix mtxResult) int k,j;double s; / 将常数矩阵赋给解矩阵mtxResult.SetValue(mtxLEConst);double pDataConst = mtxResult.GetData();int n = GetNumUnknowns();if (mtxLECoef.GetNumRows() != n)return false;/ 为系数矩阵对角线数组分配内存double pDiagData = new double3*n-2;/ 构造系数矩阵对角线元素数组k = j = 0;if (k = 0)pDiagDataj+ = mtxLECoef.GetElement(k, k);pDiagDataj+ = mtxLECoef.GetElement(k, k+1);for (k=1; kn-1; +k)pDiagDataj+ = mtxLECoef.GetElement(k, k-1);pDiagDataj+ = mtxLECoef.GetElement(k, k);pDiagDataj+ = mtxLECoef.GetElement(k, k+1);if (k = n-1)pDiagDataj+ = mtxLECoef.GetElement(k, k-1);pDiagDataj+ = mtxLECoef.GetElement(k, k);/ 求解for (k=0;k=0;k-)pDataConstk=pDataConstk-pDiagData3*k+1*pDataConstk+1; return true;/* * 一般带型方程组的求解 * * param nBandWidth - 带宽 * param mtxResult - Matrix对象，返回方程组解矩阵 * return bool 型，方程组求解是否成功 */public bool GetRootsetBand(int nBandWidth, Matrix mtxResult) int ls,k,i,j,nis=0,u,v;double p,t; / 带宽必须为奇数if (nBandWidth-1)%2 != 0)return false;/ 将常数矩阵赋给解矩阵mtxResult.SetValue(mtxLEConst);double pDataConst = mtxResult.GetData();/ 方程组属性int m = mtxLEConst.GetNumColumns();int n = GetNumUnknowns();if (mtxLECoef.GetNumRows() != n)return false;/ 带宽数组：带型矩阵的有效数据double pBandData = new doublenBandWidth*n;/ 半带宽ls = (nBandWidth-1)/2;/ 构造带宽数组for (i=0; in; +i)j = 0;for (k=Math.Max(0, i-ls); kMath.Max(0, i-ls)+nBandWidth; +k)if (k n)pBandDatai*nBandWidth+j+ = mtxLECoef.GetElement(i, k);elsepBandDatai*nBandWidth+j+ = 0;/ 求解for (k=0;k=n-2;k+) p=0.0;for (i=k;ip) p=t; nis=i; if (p = 0.0)return false;for (j=0;j=m-1;j+) u=k*m+j; v=nis*m+j;t=pDataConstu; pDataConstu=pDataConstv; pDataConstv=t; for (j=0;j=nBandWidth-1;j+) u=k*nBandWidth+j; v=nis*nBandWidth+j;t=pBandDatau; pBandDatau=pBandDatav; pBandDatav=t; for (j=0;j=m-1;j+) u=k*m+j; pDataConstu=pDataConstu/pBandDatak*nBandWidth; for (j=1;j=nBandWidth-1;j+) u=k*nBandWidth+j; pBandDatau=pBandDatau/pBandDatak*nBandWidth; for (i=k+1;i=ls;i+) t=pBandDatai*nBandWidth;for (j=0;j=m-1;j+) u=i*m+j; v=k*m+j;pDataConstu=pDataConstu-t*pDataConstv; for (j=1;j=nBandWidth-1;j+) u=i*nBandWidth+j; v=k*nBandWidth+j;pBandDatau-1=pBandDatau-t*pBandDatav; u=i*nBandWidth+nBandWidth-1; pBandDatau=0.0; if (ls!=(n-1) ls=