数值方法
gauss列主元素消去法
在gauss消元过程中位于矩阵A(k)(k为1,2,3,4……n)的对角线(k,k)位置上的元素称为主元素,因为在计算解的分量x(k)的过程时,a(k,k)都作除数,在计算机中,为了减少误差,应避免用较小的数做除数。
具体过程如图一:
图一
利用该公式编程求解:
图二
具体代码如下:
#include<stdio.h>
#define N 4
int main(){
int gauss(int n, double *a, double b[], double x[]);//用指针a来引用二维数组,b数组用于存放常向量,x数组用于存放线性方程的解。
int i;
double x[N];
double a[N][N] ={ {1,2,1,-2} , {2,5,3,-2} , {-2,-2,3,5} , {1,3,2,5} };
double B[N] = {-1,3,15,9};
if( gauss(N, &a[0][0], B, x) != 0){//确定方程有解
for(i = 0; i<=N-1; i++){
printf("x[%d] = %8.2lf\n",i,x[i]);
}
}
else
printf("线性方程组无解");
return 0;
}
int gauss(int n, double *a, double b[], double x[]){//高斯列主元素消去法
int i,j,k;
double sum, temp;
for(i = 0; i<=n-1; i++){
for(j = 0; j<=n-1; j++){
printf("%8.2lf",a[i*n+j]);
}
printf("%8.2lf",b[i]);
printf("\n");
}//打印初始矩阵
printf("\n");
for(k = 0; k<=n-2; k++){
for(i = k+1; i<=n-1; i++){
if(a[k*n+k] == 0){//确定方程有解(待会下面有解释)
printf("calculater failed!");
return 0;
}
temp = -a[i*n+k] / a[k*n+k];
b[i] = b[i] + b[k] * temp;
for(j = k+1; j<=n-1; j++){//第一列可以不用消,对最后结果并不影响也可以写成 j = k,这种写法便于看消去过程;
a[i*n+j] = a[i*n+j] + a[k*n+j] * temp;
}
}//消元
for(i = 0; i<=n-1; i++){
for(j = 0; j<=n-1; j++){
printf("%8.2lf",a[i*n+j]);
}
printf("%8.2lf",b[i]);
printf("\n");
}//输出中间的状态
printf("\n");
}
for(k = n-1; k>=0; k--){
sum = 0;
for(j = k+1; j<=n-1; j++){
sum = sum+a[k*n+j]*x[j];
}
x[k] = (b[k] - sum) / a[k*n+k];
} //回代,求解
return 1;
}
确定线性方程有解解释:假设a(n,n) = 0,则方程中必有a(n,n)*x(n) =b(n),此时x(n)无解,导致后来的解向量也无解。
代码运行效果如下图所示:
j= k+1的运行结果如图三:
图三
j = k时的运行结果如图四所示(可以看消去过程):
图四