线性方程组求解的几种常用方法-c++代码实现
来源:互联网 发布:手机淘宝怎么秒杀东西 编辑:程序博客网 时间:2024/06/05 20:49
写在博客之前
- 代码可能比较杂乱无章,初学者,大神勿喷。
- 代码后会附上算法伪码,可以对照看看。
Doolittle直接分解算法
算法伪码:
A[n][n]—系数矩阵、b[n] —常数项矩阵
L[n][n]、U[n][n]、y[n]、x[n]—L、U矩阵与解向量
Step1 L赋值:对角线元素赋值1,上三角元素赋值0;
Step2 U赋值:下三角元素赋值为0;
Step3 For k=1 To n Do Step4、Step5 /开始分解 /
Step4
For j=k To n Do /求 U中第 k 行元素/
u[k][j] = a[k][j]
For r = 1 To k-1 Do
u[k][j] = u[k][j] - l[k][r]*u[r][j]
EndFor r
EndFor j
Step5
For i=k+1 To n Do /求 L 中第 k 列元素/
l[i][k] = a[i][k];
For r = 1 To k-1 Do
l[i][k] = l[i][k] - l[i][r]*u[r][k];
EndFor r
l[i][k] = l[i][k] / u[k][k]
EndFor i
EndFor k /结束分解 /
Step6
For i = 1 To n Do /由 Ly=b 求出y[i] /
y[i] = b[i]
For j = 1 To i-1 Do
y[i] = y[i] - l[i][j]*y[j];
EndFor j
EndFor i
Step7
For i = n Downto 1 Do /由 Ux=y 求出 x[i] /
x[i] = y[i];
For j = i+1 To n Do
y[i] = y[i] - u[i][j]*x[j]
EndFor j
x[i] = y[i]/u[i][i]
EndFor i
void Doolittle(int n,double *A,double *b)//n为阶数 A为系数矩阵 b为常数矩阵{ double *L = new double[n*n];//开辟L矩阵空间 double *U = new double[n*n];//开辟U矩阵空间 double *y = new double[n];//开辟y矩阵空间 for (int i = 0; i < n; i++) { for (int j = 0; j < n; j++) { *(U + i*n + j) = 0;//暂时全部赋值为0 if (i==j) { *(L + i*n + j) = 1;//对角线赋值为1 } else { *(L + i*n + j) = 0;//其他暂时赋值为0 } } } for (int k = 0; k < n; k++)//计算u和l矩阵的值 { for (int j = k; j < n; j++) { *(U + k*n + j) = *(A + k*n + j);//第一行 for (int r = 0; r < k; r++)//接下来由L的前一列算u的下一行 { *(U + k*n + j) = *(U + k*n + j) - (*(L + k*n + r)*(*(U + r*n + j))); } } for (int i = k+1; i < n; i++)//计算L的列 { *(L + i*n + k) = *(A + i*n + k); for (int r = 0; r < k; r++) { *(L + i*n + k) = *(L + i*n + k) - (*(L + i*n + r)*(*(U + r*n + k))); } *(L + i*n + k) = *(L + i*n + k) / (*(U + k*n + k)); } } for (int i = 0; i < n; i++)//由Ly=b算y { *(y + i) = *(b + i); for (int j = 0; j < i; j++) { *(y + i) = *(y + i) - *(L + i*n + j)*(*(y + j)); } } for (int i = n-1; i >= 0; i--)//由Ux=y算x { *(x + i) = *(y + i); for (int j = i+1; j < n; j++) { *(y + i) = *(y + i) - *(U + i*n + j)*(*(x + j)); } *(x + i) = *(y + i) / (*(U + i*n + i)); } cout << "解:\n";//得出解 for (int i = 0; i < n; i++) { cout <<"x"<<i+1<<":"<< *(x + i) << endl; } delete[]L;//释放空间 delete[]U; delete[]y;}
追赶法——三对角线方程组求解
注:三对角线方程组形如
算法伪码:
数据说明:
定义一维数组 a[n],b[n],c[n],l[n],u[n],d[n],x[n],y[n]
Step1
输入数组a、b、c、d各元素
Step2
a[1]=0 ; l[1]=b[1];
u[1]=c[1]/l[1] ; y[1]=d[1]/l[1];
Step3
For i = 1 To n Do /* 追的过程,求 yi */
l[i] = b[i]-a[i]*u[i-1]
u[i] = c[i]/l[i]
y[i] = (d[i]-a[i]*y[i-1])/l[i]
EndFor i
Step4
x[n] = y[n] /* 赶的过程,求 xi */
For i = n-1 Downto 1 Do
x[i] = y[i]-u[i]*x[i+1]
EndFor i
void Catch_Solution(int n,double *a,double *b,double *c,double *d)//参数依次为:阶数 系数a 系数b 系数c 常数列d{ double *l = new double[n]; double *u = new double[n]; double *x = new double[n]; double *y = new double[n];//开空间存储数据 a[0] = 0; l[0] = b[0]; u[0] = c[0] / l[0]; y[0] = d[0] / l[0]; for (int i = 1; i < n; i++) { l[i] = b[i] - a[i] * u[i - 1]; u[i] = c[i] / l[i]; y[i] = (d[i] - a[i] * y[i - 1]) / l[i]; } x[n - 1] = y[n - 1]; for (int i = n-1; i >= 0; i--)//求x { x[i] = y[i] - u[i] * x[i + 1]; } for (int i = 0; i < n; i++) { cout << x[i] << endl; } delete[]x; delete[]y; delete[]l; delete[]u;}
LDLT分解法
注:要求系数矩阵为对称正定矩阵
算法伪码:
数据说明
a[n][n] ——存放系数矩阵A的系数;
b[n] ——存放方程组右端常数项;
x[n],y[n],z[n] ——即LTx=y,Dy=z,Lz=b;
L[n][n] ——存放L矩阵中的元素;
D[n] ——存放D矩阵中的元素;
Step1
输入系数矩阵、常数项矩阵元素a[i][j],b[i]
Step2
按下面的公式计算(k=1,…,n):
/* L对角线元素赋值为1*/
Step2-1
For i = 1 To n Do
l[i][i]=1
EndFor i
Step2-2
For k = 1 To n Do
d[k] = a[k][k]
For j=1 To k-1 Do
d[k] = d[k]-l[k][j]*l[k][j]*d[j]
EndFor j
For i = k+1 To n Do
l[i][k]=a[i][k]
For j = 1 To k-1 Do
l[i][k] = l[i][k]-l[i][j]*l[k][j]*d[j]
EndFor j
l[i][k] = l[i][k]/d[k]
EndFor i
EndFor k
Step3
For i = 1 To n Do /* 求zi */
z[i] = b[i]
For j = 1 To i-1 Do
z[i] = z[i]-l[i][j]*z[j]
EndFor j
EndFor i
Step4
For i = 1 To n Do /* 求 yi */
y[i]=z[i]/d[i]
EndFor i
Step5
For i = n Downto 1 Do /* 求 xi */
x[i] = y[i];
For j = i+1 To n Do
x[i] = x[i]-l[j][i]*x[j]
EndFor j
EndFor i
void LDL(int n, double *a, double *b)//LDL法,参数依次:阶数 系数矩阵a 常数矩阵b { double *U = new double[n*n]; double *y = new double[n]; double *z = new double[n]; double *L = new double[n*n]; double *D = new double[n];for (int i = 0; i < n; i++)//用LU先算出L U { for (int j = 0; j < n; j++) { *(U + i*n + j) = 0;//暂时全部赋值为0 if (i == j) { *(L + i*n + j) = 1;//对角线赋值为1 } else { *(L + i*n + j) = 0;//其他暂时赋值为0 } } } for (int k = 0; k < n; k++)//计算u和l矩阵的值 { for (int j = k; j < n; j++) { *(U + k*n + j) = *(a + k*n + j);//第一行 for (int r = 0; r < k; r++)//接下来由L的前一列算u的下一行 { *(U + k*n + j) = *(U + k*n + j) - (*(L + k*n + r)*(*(U + r*n + j))); } } for (int i = k + 1; i < n; i++)//计算L的列 { *(L + i*n + k) = *(a + i*n + k); for (int r = 0; r < k; r++) { *(L + i*n + k) = *(L + i*n + k) - (*(L + i*n + r)*(*(U + r*n + k))); } *(L + i*n + k) = *(L + i*n + k) / (*(U + k*n + k)); } } for (int i = 0; i < n; i++)//把D赋值 { *(D + i) = *(U + i*n + i); } for (int i = 0; i < n; i++)//由Lz=b算z { *(z + i) = *(b + i); for (int j = 0; j < i; j++) { *(z + i) = *(z + i) - *(L + i*n + j)*(*(z + j)); } } for (int i = 0; i < n; i++)//算y { *(y + i) = *(z + i) / (*(D + i)); } double *temp = new double[n*n]; for (int i = 0; i < n; i++)//这里实现对L的转置 { for (int j = 0; j < n; j++) { *(temp + i*n + j) = *(L + j*n + i); } } for (int i = 0; i < n; i++) { for (int j = 0; j < n; j++) { *(L + i*n + j) = *(temp + i*n + j); } } delete[]temp;//释放 for (int i = n-1; i >= 0; i--)//最后算x { *(x + i) = *(y + i); for (int j = i+1; j <n; j++) { *(x + i) = *(x + i) - *(L + i*n + j)*(*(x + j)); } } for (int i = 0; i < n; i++) { cout << "解为:\n"; cout << *(x + i) << endl; } delete[]U; delete[]y; delete[]z; delete[]L; delete[]D;}
高斯列主元法
算法伪码
数据说明:
N—所能求解方程组的最大阶数;
A[N][N]—表示系数矩阵;b[N]—表示常数项矩阵;
x[N]—表示解向量;n—方程组实际阶数;
Aug[N][N+1]—表示由A和b构成的增广矩阵;
Mtemp[N+1]—用于交换行;r—主元所在行;Pe—主元
Step1
由A[n][n]和b[n]生成增广矩阵Aug[n][n+1]
Step2 For k=1 To n-1 Do
/* 开始消元 */
Step3
r=k ; Pe=fabs(Aug[k][k]) /* 寻找主元 */
Step4
For i=k+1 To n Do
If Pe<
fabs(Aug[i][k]) Then
Pe=fabs(Aug[i][k]) ; r=i
EndIf
EndFor
Step5
If r< >k Then
Swap(r , k) /交换r,k两行 /
EndIf
Step6
If Aug[k][k]=0 Then /主元为0,失败/
Output(“Method Failed!”) ; Stop
EndIf
/* 将第 k 列中从 Aug[k+1][k] 开始的元素消为 0 */
Step7
For i=k+1 To n Do
m=Aug[i][k]/Aug[k][k] /* 取倍数m */
For j=k To n+1 Do
Aug[i][j]=Aug[i][j]-m*Aug[k][j]
EndFor j
EndFor i
EndFor k /* 消元完成*/
Step8
从 x[n] 开始逐步回代,求出解向量 x
x[n]=Aug[n][n+1]/Aug[n][n]
For i=n-1 Downto 1 Do
sum=0.0
For j=i+1 To n Do
sum=sum+Aug[i][j]*x[j]
EndFor j
x[i]=(Aug[i][n]-sum)/Aug[i][i]
EndFor i
void Swap(int a,int b,double *Aug,int n)//交换矩阵行 参数依次:a,b行(交换行) Aug交换的矩阵 n阶数{ double *temp = new double[n + 1]; for (int i = 0; i <n+1; i++) { *(temp + i) = *(Aug + a*(n + 1) + i); } for (int i = 0; i <n + 1; i++) { *(Aug + a*(n + 1) + i) = *(Aug + b*(n + 1) + i); } for (int i = 0; i <n + 1; i++) { *(Aug + b*(n + 1) + i) = *(temp + i); } delete[]temp;}void Guss(int n,double *Aug,double *x)//高斯列主元法,n为阶数,Aug为增广矩阵(由系数矩阵和常数矩阵合并),x为解向量{ for (int k = 0; k < n-1; k++) { double pe; int r; r = k; pe = fabs(*(Aug + k*(n + 1) + k));//假定k行k列为最大主元 for (int i = k+1; i < n; i++)//寻找最大主元 { if (pe<fabs(*(Aug + i*(n + 1) + k))) { pe = abs(*(Aug + i*(n + 1) + k)); r = i; } } if (r!=k)//后面的行为最大主元时,交换,使主元始终在前 { Swap(r, k,Aug,n); } if (pe==0) { cout << "方法失败!"; break; } for (int i = k+1; i < n; i++)//开始化简矩阵 { double m; m = *(Aug + i*(n + 1) + k) / pe; for (int j = k; j < n+1; j++) { *(Aug + i*(n + 1) + j) = *(Aug + i*(n + 1) + j) - m*(*(Aug + k*(n + 1) + j)); } } } *(x + n) = *(Aug + (n-1)*(n + 1) + n + 1) / (*(Aug + (n-1)*(n + 1) + n));//计算出最后一个未知数的值,再往回带 for (int i = n-1; i > 0; i--)//回带公式的实现 { double sum = 0; for (int j = i+1; j <= n; j++) { double m = (*(x + 1)); sum +=( *(Aug + (i - 1)*(n + 1) + j-1))*(*(x + j - 1)); } *(x + i-1) = (*(Aug + (i-1)*(n + 1) + n) - sum) / (*(Aug + (i-1)*(n + 1) + (i-1))); } cout << "增广矩阵阵为:\n"; for (int i = 0; i < n; i++) { for (int j = 0; j < n+1; j++) { cout << *(Aug + i*(n + 1) + j) << '\t'; } cout << endl; } cout << "解向量为:\n"; for (int i = 0; i < n; i++) { cout << "x"<<i+1<<":"<<*(x + i) << endl; }}
后记
第一次写博客,大神轻喷。这个算法的实现留给有用的人,可以去优化吧。
- 线性方程组求解的几种常用方法-c++代码实现
- Matlab求解线性方程组Ax=b的几种常见方法Matlab求解线性方程组Ax=b的几种常见方法
- 线性方程组的求解(C++)
- 各种线性方程组求解算法的C++实现
- 龙贝格法求解线性方程组(c++)
- C语言求解线性方程组
- 数值分析课实验-求解线性方程组的直接法C代码
- 【线性代数】线性方程组的求解
- lapack 线性方程组求解 C接口
- 行主元素消去法求解线性方程组的Pascal代码
- 数值计算线性方程组求解实现
- 矩阵的LU分解求解线性方程组(C++实现)
- c语言实现,用最小二乘法求解方程数多于未知变量的线性方程组的最适解(即矛盾方程组)
- c语言实现,用最小二乘法求解方程数多于未知变量的线性方程组的最适解(即矛盾方程组)
- 求解逆序对数的几种方法
- 有关线性方程组求解的理解
- 一种求解线性方程组的技巧
- 优化C代码常用的几招
- XAMPP 使用教程
- HTML基础知识
- 结构体如何分配内存
- 简单实现静态/动态顺序表
- VIM中的保存和退出
- 线性方程组求解的几种常用方法-c++代码实现
- C++类型推断
- 【微服务】——理论篇
- 高精度小数
- 使用git reset回退git add操作
- Action获取表单提交数据+向jsp传递数据【重要】
- 网络流24题 05圆桌聚餐
- 以zybpub新建csdn帐号,命名来源于我注册的域名zyb.pub。
- Docker实践2