雅克比迭代法c语言
A. 雅可比迭代法的工作原理
雅克比迭代法就是众多迭代法中比较早且较简单的一种,其命名也是为纪念普鲁士着名数学家雅可比。雅克比迭代法的计算公式简单,每迭代一次只需计算一次矩阵和向量的乘法,且计算过程中原始矩阵A始终不变,比较容易并行计算。
概念
考虑线性方程组Ax =b时,一般当A为低阶稠密矩阵时,用主元消去法解此方程组是有效方法。但是,对于由工程技术中产生的大型稀疏矩阵方程组(A的阶数很高,但零元素较多,例如求某些偏微分方程数值解所产生的线性方程组),利用迭代法求解此方程组就是合适的,在计算机内存和运算两方面,迭代法通常都可利用A中有大量零元素的特点。雅克比迭代法就是众多迭代法中比较早且较简单的一种,其命名也是为纪念普鲁士着名数学家雅可比。
折叠编辑本段迭代过程
迭代过程
首先将方程组中的系数矩阵A分解成三部分,即:A = L+D+U,如图1所示,其中D为对角阵,L为下三角矩阵,U为上三角矩阵。
之后确定迭代格式,X^(k+1) =B*X^(k) +f,(这里^表示的是上标,括号内数字即迭代次数),如图2所示,其中B称为迭代矩阵,雅克比迭代法中一般记为J。(k = 0,1,......)
再选取初始迭代向量X^(0),开始逐次迭代。
折叠编辑本段收敛性
设Ax= b,其中A=D+L+U为非奇异矩阵,且对角阵D也非奇异,则当迭代矩阵J的谱半径ρ(J)<1时,雅克比迭代法收敛。
折叠编辑本段优缺点
雅克比迭代法的优点明显,计算公式简单,每迭代一次只需计算一次矩阵和向量的乘法,且计算过程中原始矩阵A始终不变,比较容易并行计算。然而这种迭代方式收敛速度较慢,而且占据的存储空间较大,所以工程中一般不直接用雅克比迭代法,而用其改进方法。
折叠编辑本段程序实现示例
#include<stdio.h>
#include<math.h>
#include <stdlib.h>
main(){
float e=0.001,z,m,a[3][3]={5,2,1,-1,4,2,2,-3,10},b[3]={-12,20,3},x[3]={0,0,0},y[3];
int n=3,j,i,k=1;
while(1) {
for(i=0;i<3;i++) {
for(j=0;j<3;j++)
m=m+a[i][j]*x[j];
m=m-x[i]*a[i][i];
y[i]=(b[i]-m)/a[i][i];
m=0;
}
i=0;
while(i<3) {
z=fabs(x[i]-y[i]);
if(z>e)
break;
i++;
}
if(i!=3) {
for(i=0;i<3;i++)
x[i]=y[i];
k++;
}
else if(i==3)
break;
}
printf("%f %f %f ",y[0],y[1],y[2]);
}
B. 雅克比迭代法求解线性方程组的C语言程序
void Solve ( double dCoef [] , double dY [] , unsigned int iOrder , double dErr)
{//用Jacobi迭代法解方程组, dCoef[]系数阵, Y[]向量, iOrder给出方程阶数, dErr给出精度
double res [Max]; //方程解
double res2[Max]; //保存上一阶方程解
if ( Max < iOrder )
{
printf ("最多支持%d阶方程组.", Max);
return;
}
for ( unsigned int i = 0 ; i < iOrder ; res2 [i++] = 0.0 ); //初始解向量 (0,0...)
while ( true )
{
bool bStopIterative = true;
for (unsigned int i = 0 ; i < iOrder ; ++i)
{
double dSum2 = 0;
for (unsigned int j = 0 ; j < iOrder ; j++)
{//求第二项
if ( j == i ) continue;
dSum2 += dCoef [i * iOrder + j] * res2 [j];
}
res[i] = 1/dCoef[i * iOrder + i] * ( dY[i] - dSum2 );
if ( abs ( res2[i] - res [i] ) > dErr )
bStopIterative = false;
}
if ( bStopIterative )
break;
for (unsigned int i = 0 ; i < iOrder ; i++ )
res2[ i ] = res[ i ];
}
//输出结果
for (unsigned int i = 0 ; i < iOrder ; i++)
printf ("x%d = %lf\n", i+1 , res[i]);
}
int main(int argc, char* argv[])
{
double a[] =
{
8,-3,2,
4,11,-1,
2,1,4
};
double b[3] =
{
20,33,12
};
Solve ( a , b , 15 , 1e-10);
getchar();
return 0;
}