#用C语言求解N阶线性矩阵方程Axb的简单解法.docx

#用C语言求解N阶线性矩阵方程Axb的简单解法.docx

《#用C语言求解N阶线性矩阵方程Axb的简单解法.docx》由会员分享，可在线阅读，更多相关《#用C语言求解N阶线性矩阵方程Axb的简单解法.docx（10页珍藏版）》请在冰豆网上搜索。

#用C语言求解N阶线性矩阵方程Axb的简单解法

用C语言求解N阶线性矩阵方程Ax=b的简单解法

一、描述问题：

题目：

求解线性方程组Ax=b，写成函数。

其中，A为n×n的N阶矩阵，x为需要求解的n元未知数组成的未知矩阵，b为n个常数组成的常数矩阵。

即

运行程序时的具体实例为：

转化为矩阵形式（为检验程序的可靠性，特意选取初对角线元素为0的矩阵方程组）即为：

二、分析问题并找出解决问题的步骤：

由高等代数知识可知，解高阶线性方程组有逆矩阵求解法、增广矩阵求解法等，而在计算机C语言中，有高斯列主消元法、LU分解法、雅克比迭代法等解法。

为了和所学的高等代数知识相一致，选择使用“高斯简单迭代消元法”，和高等代数中的“增广矩阵求解法”相一致。

以下简述高斯消元法的原理：

算法基本原理：

首先，为了能够求解N阶线性方程组（N由用户输入），所以需要定义一个大于N维的数组a[dim+1][dim+1]（dim为设定的最大维数，防止计算量溢出），当用户输入的阶数N超过设定值时提示重启程序重新输入。

进而，要判断方程组是否有解，无解提示重启程序重新输入，有解的话要判断是有无数不定解还是只有唯一一组解，在计算中，只有当原方程组有且只有一组解时算法才有意义，而运用高等代数的知识，只有当系数矩阵对应的行列式|A|≠0时，原方程组才有唯一解，所以输入系数矩阵后要计算该系数矩阵的行列式|A|（定义了getresult（n）函数计算），当行列式|A|=0时同样应提示重启程序重新输入，|A|≠0时原方程组必然有且仅有唯一一组解。

判断出方程组有且仅有唯一一组解后，开始将系数矩阵和常数矩阵（合并即为增广矩阵）进行初等行变换（以 a11 为基元开始，将第j列上j行以下的所有元素化为0），使系数矩阵转化为上三角矩阵。

这里要考虑到一种特殊情况，即交换到第j-1列后，第j行第j列元素 ajj=0 ，那此时不能再以 ajj 为基元。

当变换到第j列时，从j行j列的元素 ajj 以下的各元素中选取第一个不为0的元素，通过第三类初等行变换即交换两行将其交换到 ajj 的位置上，然后再进行消元过程。

交换系数矩阵中的两行，相当于两个方程的位置交换了。

再由高斯消元法，将第j列元素除 ajj 外第j行以下的其他元素通过第二种初等行变换化为0，这样，就能使系数矩阵通过这样的行变换化为一个上三角矩阵，即

，

当系数矩阵A进行初等行变换时，常数矩阵也要进行对应的初等行变换，即此时 

那么有

接下来，进行“反代”，由 

可求出 

 ，再往上代入 

 即可求出 

 以此类推，即可从 xn推到 xn-1 ，再推到xn-2 直至 x1 。

至此，未知矩阵x的所有元素就全部求出，即求出了原方程组有且仅有的唯一一组解。

基本原理示意图：

三、编写程序

1.#include

2.#include

3.#include

4.#definedim10                    //定义最大的维数10，为防止计算值溢出

5.doublea[dim+1][dim+1],b[dim+1],x[dim+1];  //定义双精度数组

6.doubletemp;

7.doublegetarray（intn）;              //定义输入矩阵元素的函数

8.doubleshowarray（intn）;              //定义输出化简系数矩阵过程的函数

9.intn,i,j,k,p,q;

10.doublemain（）

11.{

12.  

13.printf（"请输入系数矩阵的阶数n（n<10）:

"）;

14.scanf（"%d",&n）;

15.  /*判断矩阵阶数是否超过界定值*/

16.  if（n>dim）

17.  {

18.      printf（"错误：

元数超过初设定的值%d，请重启程序重新输入\n",dim）;

19.      exit（0）;

20.  }

21.  /*输入系数矩阵和常数矩阵（即增广矩阵）的元素*/

22.  getarray（n）;

23.   

24.  /*使对角线上的主元素不为0*/

25.  for（j=1;j<=n-1;j++）

26.  {

27.      if（a[j][j]==0）

28.        for（i=j+1;i<=n;i++）

29.        {

30.          if（a[i][j]!

=0）

31.          {

32.              /*交换增广矩阵的第i行和第j行的所有元素*/

33.              for（k=1;k<=n;k++）

34.              {

35.                a[i][k]+=a[j][k];

36.                a[j][k]=a[i][k]-a[j][k];

37.                a[i][k]-=a[j][k];

38.              }

39.              b[i]+=b[j];

40.              b[j]=b[i]-b[j];

41.              b[i]-=b[j];

42.          }

43.          continue;    //找到第j列第一个不为0的元素即跳回第一层循环

44.        }

45.  }

46.  /*开始用高斯简单迭代消元法进行求解计算*/

47.  for（j=1;j<=n-1;j++）

48.  {

49.      /*使系数矩阵转化为上三角矩阵，常数矩阵相应进行变换*/

50.      for（i=j+1;i<=n;i++）

51.      {

52.        temp=a[i][j]/a[j][j];

53.        b[i]=b[i]-temp*b[j];

54.        for（k=1;k<=n;k++）

55.          a[i][k]=a[i][k]-temp*a[j][k];

56.        printf（"\n通过初等行变换增广矩阵矩阵C化为：

\n"）;

57.        /*输出进行初等行变换的过程*/

58.        printf（"C="）;

59.        for（p=1;p<=n;p++）

60.        {

61.          for（q=1;q<=n;q++）

62.              printf（"\t%.3f",a[p][q]）;

63.          printf（"\t%.3f\n",b[p]）;

64.        }

65.        printf（"\n"）;

66.      }

67.  }

68.  /*输出最终的增广矩阵C*/

69.  showarray（n）;

70.  /*开始按顺序反代求解x[i]（i=n,n-1,n-2,…,2,1）*/

71.  x[n]=b[n]/a[n][n];

72.  for（j=n-1;j>=1;j--）

73.  {

74.      x[j]=b[j];

75.      for（k=n;k>=j+1;k--）

76.        x[j]=x[j]-x[k]*a[j][k];

77.      x[j]=x[j]/a[j][j];

78.  }

79.  printf（"\n原方程组的唯一一组实数解为：

\n"）;

80.  for（j=1;j<=n;j++）

81.      printf（"x[%d]=%.3f\n",j,x[j]）;

82.}

83./*定义矩阵输入函数getarray（n）并打印以作检查*/

84.doublegetarray（intn）

85.{

86.printf（"\n请输入该矩阵各行的实数（以空格隔开）\n"）;

87.for（i=1;i<=n;i++）

88.  {

89.printf（"\n第%d行:

\t",i）;

90.for（j=1;j<=n;j++）

91.      {

92.        scanf（"%lf",&a[i][j]）;

93.        printf（"a[%d][%d]=%.3f",i,j,a[i][j]）;

94.        printf（"\n"）;

95.      }

96.  }

97.  printf（"\nA="）;

98.for（i=1;i<=n;i++）

99.  {

100.for（j=1;j<=n;j++）

101.printf（"\t%.3f",a[i][j]）;

102.printf（"\n"）;

103.  }

104.  printf（"\n"）;

105.  /*输入常数矩阵的各个数*/

106.  for（i=1;i<=n;++i）

107.  {

108.      printf（"请输入常数b[%d]=",i）;

109.      scanf（"%lf",&b[i]）;

110.  }

111.}

112./*定义增广矩阵C输出函数showarray（n）*/

113.doubleshowarray（intn）

114.{

115.printf（"\n通过初等行变换最终增广矩阵矩阵C化为：

\n"）;

116.  printf（"C="）;

117.  for（i=1;i<=n;i++）

118.  {

119.      for（j=1;j<=n;j++）

120.        printf（"\t%.3f",a[i][j]）;

121.      printf（"\t%.3f",b[i]）;

122.      printf（"\n"）;

123.  }

124.  temp=1;

125.  for（i=1;i<=n;i++）

126.      temp*=a[i][i];

127.  printf（"\n矩阵的行列式|A|=%f\n",temp）;

128.  /*判断原线性方程组是否有唯一解*/

129.  if（temp==0）

130.  {

131.      printf（"\n该方程组无唯一解，请重新启动程序输入\n"）;

132.      exit（0）;

133.  }

134.}

复制代码

程序执行结果：

四、误差分析

由程序执行结果图可知，该C语言程序所求得的该N阶矩阵方程即N维线性方程组的解为

即

由于程序中所有变量除了增广矩阵的角标以外都定义为double型，而double型变量的精确度是16位，所以程序运行过程中变量的有效数字至多有15位，而为了程序执行时界面的清爽，将每个变量的有效数字只取了小数点后3位，就运行的具体程序来说，这小数点的后三维数字均为有效数字，所以本程序的误差至多为0.001即小数点后三位。

而在该具体的5维线性方程组中，用克拉默法则计算出系数行列式 |A|=665，其精确解为

所以各个解均在误差范围内。