C语言求解线性方程组
程序员文章站
2022-03-20 21:53:42
经典问题用高斯约当算法求解线性方程组。这里要求对任意形式的线性方程组都能够妥善处理,不能只适用于方程个数和未知量数目相等的特殊情形。 先用循环结构将增广矩阵转换为阶梯形矩阵,循环结束时得到阶梯型矩阵非零行行数,同时得到一个链表其中存放有各非零行主元的列标,列标在链表中按从左到右的顺序依次递减。然后根 ......
经典问题用高斯约当算法求解线性方程组。这里要求对任意形式的线性方程组都能够妥善处理,不能只适用于方程个数和未知量数目相等的特殊情形。
先用循环结构将增广矩阵转换为阶梯形矩阵,循环结束时得到阶梯型矩阵非零行行数,同时得到一个链表其中存放有各非零行主元的列标,列标在链表中按从左到右的顺序依次递减。然后根据线性代数中线性方程组的解的情况及判别准则判断方程是否有解,有多少个解。当线性方程组有解时,需要用convert函数将其转换为简化行阶梯型矩阵,然后输出唯一解或一般解
C语言代码如下:
1 #include <stdio.h> 2 #include <malloc.h> 3 #include <math.h> 4 #define N 5 //增广矩阵列数 5 #define M 3 //增广矩阵行数 6 struct maincol 7 { 8 int col; //存放各主元下标的结构体类型 9 struct maincol *next; 10 }; 11 typedef struct maincol mc1; 12 13 int test(int s, int t, float a[][N]); //判断增广矩阵的s行至M行与t列至N列相交形成的子矩阵是否为零矩阵,若是返回0,若不是返回第一个不为零的列的列标 14 void add(mc1 *head, int col, mc1 **tail); //函数,用于新建一节点,其中保存有主元col列标,然后按递减顺序将其插入maincol类型的链表 15 void convert(float a[][N], int row, mc1 *head); //函数,用于将阶梯型矩阵转化为简化行阶梯形矩阵 16 17 void main() 18 { 19 float a[M][N]; //增广矩阵 20 char str[N+1]; 21 int i, j; 22 int s, t; //子矩阵左上角元素行列下标 23 int row, col; //row用于存放阶梯形矩阵非零行行数 24 float swap; 25 mc1 *head, *tail, *psnew; 26 27 for (i=0; i<M; i++) //输入并初始化增广矩阵 28 { 29 printf("请输入增广矩阵第%d行\n", i+1); 30 scanf("%s", str); 31 for (j=0; j<N; j++) 32 a[i][j]=str[j]-48; 33 } 34 35 head=(mc1 *) malloc(sizeof(mc1)); 36 head->next=NULL; 37 tail=head; 38 s=t=1; //子矩阵即为增广矩阵本身,用增广矩阵左上角元素行列标初始化s,t 39 40 while((col=test(s, t, a))!=0) //子矩阵不为零矩阵 41 { 42 if (s==M) //增广矩阵已化为阶梯形矩阵 43 { 44 row=s; //记录非零行行数 45 add(head, col, &tail); //最后一个非零行主元列标放入maincol类型链表 46 break; //结束循环 47 } 48 else 49 { 50 j=s-1; 51 for (i=s; i<M; i++) 52 { 53 if (fabs(a[j][col-1])<fabs(a[i][col-1])) //列选主元 54 j=i; 55 } 56 57 if (s-1!=j) 58 { 59 for (i=col-1; i<N; i++) 60 { 61 swap=a[j][i]; 62 a[j][i]=a[s-1][i]; //列选主元 63 a[s-1][i]=swap; 64 } 65 } 66 67 if (col==N) //增广矩阵已经化为阶梯形矩阵 68 { 69 row=s; //记录非零行行数 70 add(head, col, &tail); //最后一个非零行主元列标放入maincol类型链表 71 break; //结束循环 72 } 73 74 for (i=s; i<M; i++) 75 a[i][col-1]=-(a[i][col-1]/a[s-1][col-1]); 76 77 for (i=col; i<N; i++) //消元 78 { 79 for (j=s; j<M; j++) 80 a[j][i]=a[j][col-1]*a[s-1][i]+a[j][i]; 81 } 82 83 add(head, col, &tail); //将消元后得到的主元列标col放入maincol类型链表 84 s++; 85 t=col+1; //更新s,t,使s,t成为消元后得到的新的子矩阵的左上角元素行列标,为test函数操作新的子矩阵作准备 86 continue; //开始新一轮循环 87 } 88 } 89 if (col==0) //从循环控制条件退出循环 90 row=s-1; //此时增广矩阵已成为阶梯形矩阵,非零行函数就是s-1 91 92 if (row==0) //利用线性方程组解的判别准则判断是否有解,有多少个解 93 { 94 printf("线性方程组有无穷多组解\n"); //增广矩阵为零矩阵,无穷多组解 95 printf("一般解:\n"); 96 for (i=1; i<N; i++) 97 printf("x%d=t%d\n", i, i); //输出解 98 } 99 else 100 { 101 psnew=head->next; 102 if (psnew->col==N) //阶梯形矩阵最后一主元在最后一列,无解 103 printf("线性方程组无解\n"); 104 else 105 { 106 convert(a, row, head); //用convert函数将阶梯形矩阵进一步化为简化行阶梯形矩阵 107 if (row==N-1) //非零行行数等于未知量个数,有唯一解 108 { 109 printf("线性方程组有唯一解:\n"); 110 for (i=1; i<=row; i++) //输出唯一解 111 printf("x%d=%f\n", i, a[i-1][N-1]); 112 } 113 else //非零行行数小于未知量个数,有无穷多组解 114 { 115 int *m; 116 m=(int *) malloc((N-1-row)*sizeof(int)); 117 118 i=N-1-row; 119 for (j=N-1; j>=1; j--) 120 { 121 if (j!=psnew->col) 122 { 123 m[--i]=j; //从所有未知量标号中筛选出*未知量标号 124 if (i==0) 125 break; 126 } 127 else 128 { 129 if (psnew->next!=NULL) 130 psnew=psnew->next; 131 } 132 } 133 printf("线性方程组有无穷多组解\n"); 134 printf("一般解:\n"); 135 i=row; 136 for (psnew=head->next; psnew!=NULL; psnew=psnew->next) 137 { 138 printf("x%d=%f", psnew->col, a[i-1][N-1]); //输出一般解 139 for (j=0; j<N-1-row; j++) 140 { 141 if(m[j]<psnew->col) 142 { 143 printf("-%dx%d", 0, m[j]); 144 } 145 else 146 { 147 printf("-%fx%d", a[i-1][m[j]-1], m[j]); 148 } 149 } 150 printf("\n"); 151 i--; 152 } 153 } 154 } 155 } 156 } 157 158 int test(int s, int t, float a[][N]) //判断增广矩阵的s行至M行与t列至N列相交形成的子矩阵是否为零矩阵,若是返回0,若不是返回第一个不为零的列的列标 159 { 160 int i, j; 161 162 for (j=t-1; j<N; j++) 163 { 164 for (i=s-1; i<M; i++) 165 { 166 if (a[i][j]!=0) 167 return (j+1); 168 } 169 } 170 return (0); 171 } 172 173 void add(mc1 *head, int col, mc1 **tail) //函数,用于新建一节点,其中保存有主元col列标,然后按递减顺序将其插入maincol类型的链表 174 { 175 mc1* psnew; 176 177 psnew=(mc1 *) malloc(sizeof(mc1)); 178 psnew->col=col; 179 180 if(head->next==NULL) 181 { 182 psnew->next=NULL; 183 head->next=psnew; 184 *tail=psnew; 185 } 186 else 187 { 188 psnew->next=head->next; 189 head->next=psnew; 190 } 191 } 192 193 void convert(float a[][N], int row, mc1 *head) //函数,用于将阶梯型矩阵转化为简化行阶梯形矩阵 194 { 195 mc1 *psnew, *pq; 196 int i, j, k, m; 197 198 psnew=head->next; 199 for (i=row-1; i>=0; i--) 200 { 201 if (a[i][psnew->col-1]!=1) //各非零行主元化为1 202 { 203 for (j=psnew->col; j<N; j++) 204 a[i][j]=a[i][j]/a[i][psnew->col-1]; 205 } 206 207 psnew=psnew->next; 208 } 209 210 psnew=head->next; //向上消元把除第一个主元外其余主元所在列中在主元上方的部分变为零 211 for (i=row-1; i>=1; i--) 212 { 213 m=N-psnew->col-(row-i); //获取未知量标号1,2,--,N-1中位于i+1非零行主元列标号右边的*未知量标号个数 214 for (j=i-1; j>=0; j--) 215 { 216 pq=head->next; //pq指向存放最后一个非零行主元列标号的节点 217 for (k=N; k>psnew->col; k--) 218 { 219 if (k!=pq->col) 220 { 221 a[j][k-1]=-(a[i][k-1]*a[j][psnew->col-1])+a[j][k-1]; //从右向左作初等行变换直至i+1行主元所在列右边的列位置,期间跳过i+2----row行主元所在的列 222 m--; 223 if (m==0) 224 break; 225 } 226 else 227 { 228 if (pq->next!=psnew) 229 pq=pq->next; 230 } 231 } 232 } 233 psnew=psnew->next; //递进至上一行主元,为新一轮向上消元作准备 234 } 235 }
下一篇: C 中数字数据类型在不同机器上所占字节数