經典問題用高斯約當演算法求解線性方程組。這裡要求對任意形式的線性方程組都能夠妥善處理,不能只適用於方程個數和未知量數目相等的特殊情形。 先用迴圈結構將增廣矩陣轉換為階梯形矩陣,迴圈結束時得到階梯型矩陣非零行行數,同時得到一個鏈表其中存放有各非零行主元的列標,列標在鏈表中按從左到右的順序依次遞減。然後根 ...
經典問題用高斯約當演算法求解線性方程組。這裡要求對任意形式的線性方程組都能夠妥善處理,不能只適用於方程個數和未知量數目相等的特殊情形。
先用迴圈結構將增廣矩陣轉換為階梯形矩陣,迴圈結束時得到階梯型矩陣非零行行數,同時得到一個鏈表其中存放有各非零行主元的列標,列標在鏈表中按從左到右的順序依次遞減。然後根據線性代數中線性方程組的解的情況及判別準則判斷方程是否有解,有多少個解。當線性方程組有解時,需要用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 }