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

1、用C语言求解N阶线性矩阵方程 Ax=b的简单解法、描述问题:题目求解线性方程组 Ax=b，写成函数。其中， A为nxn的N阶矩阵，x为需要求解的n 元未知数组成的未知矩阵，b为n个常数组成的常数矩阵。即1222如丁 IIan2* « - fl-nA J14运行程序时的具体实例为：仇2十5x3十=0+ 5 定4 + 6x5 = 04 丸4 + 5冷=15丸 1 + 6x2 = 1 、直1 + 5x2 + 6x3 = 00的矩阵方程组）即为:转化为矩阵形式（为检验程序的可靠性，特意选取初对角线元素为01560 0150 0 0 15 6 0 015 6 001尤irOi6%051010S

2、二、分析问题并找出解决问题的步骤：由高等代数知识可知，解高阶线性方程组有逆矩阵求解法、增广矩阵求解法等，而在计算机C语言中，有高斯列主消元法、LU分解法、雅克比迭代法等解法。为了与所学的高等代数知识相一致，选择使用“高斯简单迭代消元法”，与高等代数中的“增广矩阵求解法”相一致。以下简述高斯消元法的原理：算法基本原理：首先，为了能够求解 N阶线性方程组（N由用户输入），所以需要定义一个大于 N维的数组 adim+1dim+1 （ dim为设定的最大维数，防止计算量溢出），当用户输入的阶数N超过设定值时提示重启程序重新输入。进而，要判断方程组是否有解，无解提示重启程序重新输入，有解的话要判断是有无

3、数不定 解还是只有唯一一组解，在计算中，只有当原方程组有且只有一组解时算法才有意义，而运用高等代数的知识，只有当系数矩阵对应的行列式|A|工0时，原方程组才有唯一解，所以输入系数矩阵后要计算该系数矩阵的行列式|A| （定义了 getresult（n）函数计算），当行列式|A|=0 时同样应提示重启程序重新输入，|A|工0时原方程组必然有且仅有唯一一组解。判断出方程组有且仅有唯一一组解后，开始将系数矩阵和常数矩阵（合并即为增广矩阵）进行初等行变换（以 aii为基元开始，将第j列上j行以下的所有元素化为 0）,使系数矩阵转化为上三角矩阵。这里要考虑到一种特殊情况，即交换到第j-1列后，第j行第j列

4、元素ajj=0，那此时不能再以ajj为基元。当变换到第j列时，从j行j列的元素ajj以下的各元素中选取第一个不为o的元素，通过第三类初等行变换即交换两行将其交换到 的位置上，然后再进行消元过程。交换系数矩阵中的两行，相当于两个方程的位置交换了。再由高斯消元法，将第j列元素除ajj外第j行以下的其他元素通过第二种初等行变换化为o,这样，就能使系数矩阵通过这样的行变换化为一个上三角矩阵，即ana12Af =0a22切00当系数矩阵a进行初等行变换时，常数矩阵也要进行对应的初等行变换，即此那么有龙2X nJ接下来，进行“反代”，由略Z =兀1可求出”佩=耳十爲，再往上代入务-X + g-炕舅斗=由L

5、即可求出_i J以此类推，即可从 Xn推到Xn-1，再推到Xn-2直至X1。至此，未知矩阵 X 的所有元素就全部求出，即求出了原方程组有且仅有的唯一一组解。基本原理示意图：三、编写程序开始si按行»入系麴E阵A优简系数矩阵为上三角扌醉输出上三角師的讦列式的aiAi输入常数抵阵bF -1si用嵩斷消 元法求解1. #include<stdio.h>2. #include<stdlib.h>3. #include<math.h>j卄14-+矩阵元素超出设疽|正11方程无解或有础个認互换増广扌阱 中第1行和第J 行的对应元素是我到该列第f不九0的元素呵#

6、define dim 10/定义最大的维数10，为防止计算4.5.6.7.8.9.10.11.12.13.14.15.16.17.18.19.20.21.22.23.24.25.26.27.28.29.30.31.32.33.值溢岀定义双精度数组定义输入矩阵元素的函数定义输岀化简系数矩阵过程的函double adim+1dim+1,bdim+1,xdim+1; /double temp;double getarray(int n);/double showarray(int n);/数int n,i,j,k,p,q;double main()printf(" 请输入系数矩阵的阶数 n

7、(n<10):"); scanf("%d",&n);/* 判断矩阵阶数是否超过界定值 */if(n>dim)printf("错误：元数超过初设定的值d,请重启程序重新输入n",dim);exit(0);/*输入系数矩阵和常数矩阵(即增广矩阵)的元素 */getarray(n);/*使对角线上的主元素不为0*/for(j=1;jv=n_1;j+)if(ajj=0)for(i=j+1;i<=n;i+)if(aij!=O)#define dim 10/定义最大的维数10，为防止计算#define dim 10/定义最大的维数

8、10，为防止计算/*交换增广矩阵的第i行与第j行的所有元素*/34.35.36.37.38.39.40.41.42.43.44.45.46.47.48.49.50.51.52.53.54.55.56.57.58.59.60.61.62.63.64.for(k=1;k<=n;k+) a【i【k+=aj【k;aj【k=a【i【k-aj【k;aik-=ajk;b【i+=bjbj=b【i-bj;bi-=bj;continue; /找到第j列第一个不为0的元素即跳回第一层循环/*开始用高斯简单迭代消元法进行求解计算*/for(j=1;jv=n_1;j+)/*使系数矩阵转化为上三角矩阵，常数矩阵相应

9、进行变换*/for(i=j+1;i<=n;i+)temp=a【ij/ajj;bi=bi-temp*bj;for(k=1;k<=n;k+)aik=aik-temp*ajk;printf("n通过初等行变换增广矩阵矩阵C化为：n");/*输岀进行初等行变换的过程*/printf("C=");for(p=1;p<=n;p+)for(q=1;qv=n;q+)printf("t%.3f",apq);printf("t%.3fn",bp);65.66.67.68.69.70.71.72.73.74.75.76

10、.77.78.79.80.81.82.83.84.85.86.87.88.89.90.91.92.93.94.95.96.printf("n");/*输岀最终的增广矩阵C*/showarray(n);/*开始按顺序反代求解xi(i=n,n-1,n-2,，2,1)*/xn=bn/ann;for(j=n-1;j>=1;j-)xj=bj;for(k=n;k>=j+1;k-)xj=xj-xk*ajk;xj=xj/ajj;printf("n原方程组的唯一一组实数解为：n");for(j=1;j<=n;j+)printf("x%d= %.

11、3fn",j,xj);/*定义矩阵输入函数getarray(n)并打印以作检查*/double getarray(int n)printf("n请输入该矩阵各行的实数(以空格隔开) n");for(i=1;i<=n;i+)printf("n 第 %d 行:t",i);for(j=1;j<=n;j+)scanf("%lf", &aij);97. printf("a%d%d= %.3f",i,j,aij);98. prin tf("n");99. 100. 101. p

12、rintf("nA=");102. for(i=1;i<=n;i+)103. 104. for(j=1;j<=n;j+)105. printf("t%.3f",aij);106. printf("n");107. 108. printf("n");109. /*输入常数矩阵的各个数*/110. for(i=1;i<=n;+i)111. 112. printf("请输入常数 b%d = ",i);113. scanf("%lf",&bi);114. 1

13、15. 116.116. /*定义增广矩阵 C输出函数showarray(n)*/117. double showarray(int n)118. 119. printf("n通过初等行变换最终增广矩阵矩阵C化为：n");120. printf("C=");121. for(i=1;i<=n;i+)122. 123. for(j=1;jv=n;j+)124. printf("t%.3f",aij);125. printf("t%.3f",bi);126. printf("n");129.1

14、30. temp=1;131. for(i=1;i<=n;i+)132. temp*=aii;133. printf("n矩阵的行列式 |A|=%fn",temp);134. /*判断原线性方程组是否有唯一解*/135. if(temp=0)136. 137. printf("n该方程组无唯一解，请重新启动程序输入n");138. exit(0);139. 140. 复制代码程序执行结果:S3 C: Nindov.ssystem32cmd.exe请输入系数矩阵的阶数介5<10：5请输入该矩阵各行的实数（以空格隔开）/str y第仆atlHll

15、 aCiC2a(ll3 alH41 aClJCSl0 15 6=0.000=1.000=5.000 6.000=0.000第2行三a2Hll a2K21 a2C3J a(2K4 a2E5J0 0 15 =0.000 =0.000=1.000-5.000=6.000第3仃：a3ClJaC3H2 a3H3J a3JC4J aE3)C50 0 0 1 =0.000 =0.000 0.000 1.000» 5.000第4仃： a4JLl a【4J【4JC3J a(4J41 a4JC55 6 0 0=5.000=6.000=0.000 0.000=0.000第5行：a5Cl a5JC2 a5K

16、3J a5J41 a5JC5J15 6 0=1.000=5.000=6.000-0.000=0.0000.0001.0005.0006.0000.0000.0001.0005.0000.0000.0000.0001.0005.0006.0000.0000.0001.0005.0006.0000.0000 0006.0005.0000.0000.00012 3 4 5 入入入入入 QSrQ&r J 厶 主冃主冃主冃主冃主冃 Ik TP 1 Tk Tk TP通过初等行变换増广矩阵矩阵C化为:1.0005.0006.0000.0000.000-19.000-30.0000.0000.0001

17、.0005.0006.0000.0000.0000.0001.0000.0000.0001.0005 0000.000 0 000 0.0005.0006.0000.0001.0000.0001.0000 000通过初等行变换増广矩阵矩阵C化为X1 0005.0加6 000P R00-19.000-30.0000.0000.0003.4210.0000.0000.0000.0000.000) 0000 6 15加加加UM0 00 0 00.0BD 000 0.0531.0000.000通过初等行变换最终增广矩阵矩阵C化为:四、误差分析由程序执行结果图可知，该C语言程序所求得的该 N阶矩阵方程即N维线性方程组的解为r 2.2661 -1.722X =1.057 -0. 594 L 0.319Xj = L. Z66 = -1.722 寓3 1 r 057jr4 = -0.594 xs = 0.319由于程序中所有变量除了增广矩阵的角标以外都定义为double型，而double型变量的精确度是16位，所以程序运行过程中变量的有效数字至多有15位，而为了程序执行时界面的清爽，将每个变量的有效数字只取了小数点后3位，就运行的具体程