第1个回答 2008-04-25
解方程组有好几种方法:
一:高斯消去法,
程序段:
#define N 3
#include <iostream.h>
#include <math.h>
void main()
{
double a[N][N+1] = {0.101,2.304,3.555,1.183,
-1.347,3.712,4.623,2.137,
-2.835,1.072,5.643,3.035}; //增广矩阵
double x[N]; //存储解向量
double m; // 中间变量
int i, j, k;
//消元过程
for(k=0;k<N-1;k++){ // k<N-1 ?
if(a[k][k]==0) {
cout << "求解失败!";
break;
}
for(i=k+1;i<N;i++){
m = a[i][k] / a[k][k];
a[i][k] = 0;
for(j=k+1;j<=N;j++) // j<=N?
a[i][j] = a[i][j] - a[k][j] * m;
}
}
// 回代过程
x[N-1] = a[N-1][N] / a[N-1][N-1];
for(k=N-2;k>=0;k--){
for(j=k+1;j<N;j++){
a[k][N] = a[k][N] - a[k][j] * x[j];
}
x[k] = a[k][N] / a[k][k];
}
// 输出结果
for(k=0;k<N;k++)
cout << x[k]<<endl;
}
二:列主元高斯消去算法,
程序段:
#define N 3
#include <iostream.h>
#include <math.h>
void main()
{
double a[N][N+1] = {0.101,2.304,3.555,1.183,
-1.347,3.712,4.623,2.137,
-2.835,1.072,5.643,3.035}; //增广矩阵
double x[N]; //存储解向量
double m; // 中间变量
int i, j, k,s;
//消元过程
for(k=0;k<N-1;k++)
{ // k<N-1 ?
s=k;
for(i=k+1;i<N;i++)
{
if(fabs(a[s][k])<fabs(a[i][k]))
s=i;
}
if(s!=k)
{
for(j=0;j<=N;j++)
{
m=a[k][j];
a[k][j]=a[s][j];
a[s][j]=m;
}
}
if(a[k][k]==0)
{
cout << "求解失败!";
break;
}
for(i=k+1;i<N;i++)
{
m = a[i][k] / a[k][k];
a[i][k] = 0;
for(j=k+1;j<=N;j++) // j<=N?
a[i][j] = a[i][j] - a[k][j] * m;
}
}
// 回代过程
x[N-1] = a[N-1][N] / a[N-1][N-1];
for(k=N-2;k>=0;k--)
{
for(j=k+1;j<N;j++)
{
a[k][N] = a[k][N] - a[k][j] * x[j];
}
x[k] = a[k][N] / a[k][k];
}
// 输出结果
for(k=0;k<N;k++)
cout << x[k]<<endl;
}
三:顺序高斯消去算法相对应的LU分解方法:
程序段:
#include <iostream.h>
#include <math.h>
#define N 4
void main()
{
double a[N][N]={0.3e-15,59.14,3,1,
5.291,-6.13,-1,2,
11.2,9,5,2,
1,2,1,1 }; //存储矩阵A
double x[N],b[N]={59.17,46.78,1,2};
int i,j,k;
//****** 开始LU分解 *************
for(k=0;k<N-1;k++)
{
if(a[k][k]==0)
{
cout << "主元素为零!" << endl;
break;
}
for(i=k+1;i<N;i++)
{
a[i][k] = a[i][k] / a[k][k]; // 矩阵A的严格下三角部分存储L矩阵的严格下三角部分
for(j=k+1;j<N;j++)
{
a[i][j] = a[i][j] - a[i][k] * a[k][j];
} //计算U矩阵
}
}
//********* LU 分解完成 **************
for(i=0;i<N;i++)
{
for(j=0;j<N;j++)
{
cout << a[i][j] << " ";
}
cout << endl;
} // 输出 LU 分解
x[0]=b[0];
for(k=1;k<N;k++)
{
for(j=0;j<k;j++)
{
b[k]=b[k]-a[k][j]*x[j];
}
x[k]=b[k];
}
x[N-1]=x[N-1]/a[N-1][N-1];
for(k=N-2;k>=0;k--)
{
for(j=k+1;j<N;j++)
{
x[k] = x[k] - a[k][j] * x[j];
}
x[k] = x[k] / a[k][k];
}
for(k=0;k<N;k++)
cout << x[k]<<endl;
}
四:列主元高斯消去算法所对应的LU分解方法:
程序段:
#include <iostream.h>
#include <math.h>
#define N 4
void main()
{
double a[N][N]={0.3e-15,59.14,3,1,
5.291,-6.13,-1,2,
11.2,9,5,2,
1,2,1,1 }; //存储矩阵A
double x[N],b[N]={59.17,46.78,1,2};
int p[N]; // 记录分解过程中的行交换!
int i,j,k;
for(k=0;k<N;k++)
p[k] = k+1; // 初始化,相当于单位阵
//****** 开始LU分解 *************
for(k=0;k<N-1;k++)
{
int s,t; //s存储当前列绝对值最大的元素所在的行号,t为中间变量
s=k;
for(i=k+1;i<N;i++)
{
if(fabs(a[i][k])>fabs(a[s][k]))
s=i;
} // 选主元素
if(s!=k)
{
double m;
for(j=0;j<N;j++)// 注意这里要交换k、s两行所有的元素
{
m = a[k][j];
a[k][j] = a[s][j];
a[s][j] = m;
} // 交换行
t= p[s];
p[s] = p[k];
p[k] = t; //记录当前行交换
}
if(a[k][k]==0)
{
cout << "存在为零的主元素!" << endl;
continue;
}// 判定主元素为零,矩阵奇异,方程组解不唯一
for(i=k+1;i<N;i++)
{
a[i][k] = a[i][k] / a[k][k]; // 矩阵A的严格下三角部分存储L矩阵的严格下三角部分
for(j=k+1;j<N;j++)
{
a[i][j] = a[i][j] - a[i][k] * a[k][j];
} //计算U矩阵
}
}
//********* LU 分解完成 **************
for(i=0;i<N;i++)
{
for(j=0;j<N;j++)
{
cout << a[i][j] << " ";
}
cout << endl;
} // 输出 LU 分解
x[0]=b[0];
for(k=1;k<N;k++)
{
for(j=0;j<k;j++)
{
b[k]=b[k]-a[k][j]*x[j];
}
x[k]=b[k];
}
x[N-1]=x[N-1]/a[N-1][N-1];
for(k=N-2;k>=0;k--)
{
for(j=k+1;j<N;j++)
{
x[k] = x[k] - a[k][j] * x[j];
}
x[k] = x[k] / a[k][k];
}
for(k=0;k<N;k++)
cout << x[k]<<endl;
}
还有一些迭代法;
一:雅可比迭代
#define N 4 // 线性方程组的阶数
#include <iostream.h>
#include <math.h>
void main()
{
double a[N][N]={-4,1,1,1,1,-4,1,1,1,1,-4,1,1,1,1,-4}, //系数矩阵
b[N]={1,1,1,1}; //右端常数向量
double x0[N]={1,1,1,1},x[N]; // 迭代初始向量和迭代向量
double e=1e-5; // 精度要求
int M=5000; //最大迭代次数
int i,j,c_M=0;
double sum,current_e;
do{
current_e=0;
for(i=0;i<N;i++)
{
sum = 0;
for(j=0;j<N;j++)
{
if(j!=i)
{
sum = sum + a[i][j] * x0[j];
}
}
x[i] = (b[i] - sum) / a[i][i];
}// 更新迭代向量
c_M++; //迭代次数加1
for(i=0;i<N;i++)
{
if(fabs(x[i]-x0[i])>current_e)
current_e = fabs(x[i]-x0[i]);
} //计算当前误差
for(i=0;i<N;i++)
x0[i] = x[i]; //更新初始向量
}while(current_e>e&&c_M<M); //判断是否仍未达到精度要求且未达到最大迭代次数
for(i=0;i<N;i++)
cout << x[i] << endl; //输出结果
cout << c_M << endl; //输出迭代次数
}
二:高斯-塞德尔(Gauss-Seidel)迭代算法
程序
#define N 3 // 线性方程组的阶数
#include <iostream.h>
#include <math.h>
void main()
{
double a[N][N]={5,2,1,
2,8,-3,
1,-3,-6}, //系数矩阵
b[N]={8,21,1}; //右端常数向量
double x0[N]={1,1,1},x[N]; // 迭代初始向量和迭代向量
double e=1e-5; // 精度要求
int M=5000; //最大迭代次数
int i,j,c_M=0;
double sum,current_e;
do{
current_e=0;
for(i=0;i<N;i++)
{
sum = 0;
for(j=0;j<N;j++)
{
if(j<i)
{
sum = sum + a[i][j] * x[j];
}
if(j>i)
{
sum = sum + a[i][j] * x0[j];
}
}
x[i] = (b[i] - sum) / a[i][i];
}// 更新迭代向量
c_M++; //迭代次数加1
for(i=0;i<N;i++)
{
if(fabs(x[i]-x0[i])>current_e)
current_e = fabs(x[i]-x0[i]);
} //计算当前误差
for(i=0;i<N;i++)
x0[i] = x[i]; //更新初始向量
}while(current_e>e&&c_M<M); //判断是否仍未达到精度要求且未达到最大迭代次数
for(i=0;i<N;i++)
cout << x[i] << endl; //输出结果
cout << c_M << endl; //输出迭代次数
}
三:SOR迭代算法
程序
#define N 4 // 线性方程组的阶数
#include <iostream.h>
#include <math.h>
void main()
{
double a[N][N]={-4,1,1,1,
1,-4,1,1,
1,1,-4,1,
1,1,1,-4}, //系数矩阵
b[N]={1,1,1,1}; //右端常数向量
double x0[N]={1,1,1,1},x[N]; // 迭代初始向量和迭代向量
double e=1e-5,w=0.5; // 精度要求
int M=5000; //最大迭代次数
int i,j,c_M=0;
double sum,current_e;
do{
current_e=0;
for(i=0;i<N;i++)
{
sum = 0;
for(j=0;j<N;j++)
{
if(j<i)
{
sum = sum + a[i][j] * x[j];
}
if(j>i)
{
sum = sum + a[i][j] * x0[j];
}
}
x[i] = (1-w)*x0[i] + w*(b[i] - sum) / a[i][i];
}// 更新迭代向量
c_M++; //迭代次数加1
for(i=0;i<N;i++)
{
if(fabs(x[i]-x0[i])>current_e)
current_e = fabs(x[i]-x0[i]);
} //计算当前误差
for(i=0;i<N;i++)
x0[i] = x[i]; //更新初始向量
}while(current_e>e&&c_M<M); //判断是否仍未达到精度要求且未达到最大迭代次数
for(i=0;i<N;i++)
cout << x[i] << endl; //输出结果
cout << c_M << endl; //输出迭代次数
}