标签：
现在发现每写一篇随笔，就要在前面添加些牢骚话，各位看客如果嫌烦，直接绕道吧。
近期重新拾起C语言，因为工作的需要。
图像这个行当，matlab可以作为测试，但是真正应用的话还得转成C，所以这就是我这段时间苦逼的开始。
因为需要用到多项式变换，其中的系数求解又牵涉到线性方程组和最小二乘法的求解，所以在此，单开小灶来讲解最小二乘法和列主元高斯消元法。
有关最小二乘法的详细介绍可以参考维基百科：
相信有点数学功底的人都能看懂，这里不加详解。
在此贴上C函数代码
1 int LeastSquare(double* a, double* d, int c, int n) 2 { 3 /* 4 最小二乘法求超线性方程组 5 A*C=B 6 a表示增广矩阵A，d表示输出结果C，c表示未知数C的个数，n表示线性方程组的个数 7 例： 8 线性方程组： 9 6x+7y+9z=19 10 4x+5y+7z=15 11 6x+2y+3z=11 12 7x+3y+2z=7 13 则参数说明为： 14 d=[x; 15 y; 16 z] 17 a=[6 7 9 11； 18 4 5 7 15； 19 6 2 3 11； 20 7 3 2 7] 21 */ 22 if (a == NULL || d == NULL) 23 return -1; 24 double* equ = (double*)malloc(sizeof(double)*c*(c + 1)); 25 for (int i = 0; i < c; i++){ 26 for (int j = 0; j < c + 1; j++){ 27 equ[j + i*(c + 1)] = 0.0; 28 } 29 } 30 31 double* col = (double*)malloc(sizeof(double)*n); 32 if (equ == NULL) 33 return -2; 34 for (int i = 0; i < c; i++){ 35 for (int k = 0; k < n; k++){ 36 col[k] = a[i + k*n]; 37 cout << col[k] << " "; 38 } 39 cout << endl; 40 for (int j = 0; j < c + 1; j++){ 41 for (int l = 0; l < n; l++){ 42 equ[j + i*(c + 1)] = equ[j + i*(c + 1)] + col[l] * a[j + l*(c + 1)]; 43 } 44 } 45 } 46 for (int i = 0; i < c; i++){ 47 for (int j = 0; j < c + 1; j++){ 48 cout << equ[j + i*(c + 1)] << " "; 49 } 50 cout << endl; 51 } 52 53 GaussianElimination(equ, d, c); 54 return 0; 55 }
有关列主元高斯消元法可以在任何一本数值分析的书中找到讲解，这里我参考维基百科进行编程。
代码如下
1 int GaussianElimination(double* a, double* b, int n) 2 { 3 /* 4 这是利用高斯消元法求解线性方程组的解 5 输入： 6 a：输入增广矩阵，维数n*(n+1) 7 n:未知数的个数 8 输出： 9 b:输出矩阵，维数n*1 10 */ 11 if (a == NULL || b == NULL) 12 return -1; 13 /*for (int i = 0; i < n; i++){ 14 for (int j = 0; j < n + 1; j++){ 15 cout << a[j + i*(n + 1)] << " "; 16 } 17 cout << endl; 18 } 19 cout << "-------------" << endl;*/ 20 // 转化三角矩阵 21 double* temp = (double*)malloc(sizeof(double)*(n + 1)); 22 double max = 0.0; //每一列的最大值 23 int maxj = 0; //每一列最大值所在行数 24 for (int i = 0; i < n; i++){ 25 max = a[i + i*(n + 1)]; 26 maxj = i; 27 for (int k = 0; k < n + 1; k++){ 28 temp[k] = a[k + i*(n + 1)]; 29 } 30 //寻找列最大值 31 for (int j = i + 1; j < n; j++){ 32 if (fabs(a[i + j*(n + 1)]) > max){ 33 max = a[i + j*(n + 1)]; 34 maxj = j; 35 } 36 } 37 //交换行 38 if (maxj != i){ 39 for (int k = 0; k < n + 1; k++){ 40 a[k + i*(n + 1)] = a[k + maxj*(n + 1)] / max; 41 a[k + maxj*(n + 1)] = temp[k]; 42 } 43 } 44 else{ 45 for (int k = 0; k < n + 1; k++){ 46 a[k + i*(n + 1)] = a[k + i*(n + 1)] / max; 47 } 48 } 49 50 //求三角 51 double div = 0.0; 52 for (int j = i + 1; j < n; j++){ 53 div = a[i + j*(n + 1)]; 54 for (int k = i; k < n + 1; k++){ 55 a[k + j*(n + 1)] = a[k + j*(n + 1)] - a[k + i*(n + 1)] * div; 56 } 57 } 58 /*for (int i = 0; i < n; i++){ 59 for (int j = 0; j < n + 1; j++){ 60 cout << a[j + i*(n + 1)] << " "; 61 } 62 cout << endl; 63 } 64 cout << "-------------" << endl;*/ 65 } 66 67 /*for (int i = 0; i < n; i++){ 68 for (int j = 0; j < n + 1; j++){ 69 cout << a[j + i*(n + 1)] << " "; 70 } 71 cout << endl; 72 } 73 cout << "-------------" << endl;*/ 74 //反迭代求解 75 for (int i = n - 1; i >-1; i--){ 76 b[i] = a[i*(n + 1) + n]; 77 for (int j = n - 1; j > i; j--){ 78 b[i] = b[i] - a[j + i*(n + 1)] * b[j]; 79 } 80 } 81 /*for (int i = 0; i < n; i++){ 82 cout << b[i] << endl; 83 }*/ 84 return 0; 85 }
这是我用来测试的整个代码，仅供参考：
1 #include <iostream> 2 #include <math.h> 3 using namespace std; 4 #define epsinon 0.001 5 6 int Jacobi(double* a,double* b,int n) 7 { 8 /* 9 这是雅可比迭代法求线性方程 10 a表示增广矩阵 11 b表示输出结果 12 n表示未知数的个数 13 */ 14 if (a == NULL || b==NULL){ 15 return -1; 16 } 17 for (int i = 0; i < n; i++){ 18 b[i] = 0.0; 19 cout << b[i] << " "; 20 } 21 cout << endl; 22 double tmp = b[n-1]; 23 double c = 0.0; 24 double d = 0.0; 25 do{ 26 tmp = b[n - 1]; 27 for (int i = 0; i < n; i++){ 28 c = 0.0; 29 for (int j = 0; j < n ; j++){ 30 if (j != i){ 31 c = c + a[j + i*(n + 1)] * b[j]; 32 } 33 } 34 b[i] = (-a[n + i*(n + 1)] - c) / a[i + i*(n + 1)]; 35 cout << b[i] << " "; 36 } 37 cout << endl; 38 d = fabs((tmp - b[n - 1])); 39 } while (d > epsinon); 40 41 for (int i = 0; i < n; i++){ 42 cout << b[i] << " "; 43 } 44 cout << endl; 45 return 0; 46 47 } 48 49 int GaussianElimination(double* a, double* b, int n) 50 { 51 /* 52 这是利用高斯消元法求解线性方程组的解 53 输入： 54 a：输入增广矩阵，维数n*(n+1) 55 n:未知数的个数 56 输出： 57 b:输出矩阵，维数n*1 58 */ 59 if (a == NULL || b == NULL) 60 return -1; 61 /*for (int i = 0; i < n; i++){ 62 for (int j = 0; j < n + 1; j++){ 63 cout << a[j + i*(n + 1)] << " "; 64 } 65 cout << endl; 66 } 67 cout << "-------------" << endl;*/ 68 // 转化三角矩阵 69 double* temp = (double*)malloc(sizeof(double)*(n + 1)); 70 double max = 0.0; //每一列的最大值 71 int maxj = 0; //每一列最大值所在行数 72 for (int i = 0; i < n; i++){ 73 max = a[i + i*(n + 1)]; 74 maxj = i; 75 for (int k = 0; k < n + 1; k++){ 76 temp[k] = a[k + i*(n + 1)]; 77 } 78 //寻找列最大值 79 for (int j = i + 1; j < n; j++){ 80 if (fabs(a[i + j*(n + 1)]) > max){ 81 max = a[i + j*(n + 1)]; 82 maxj = j; 83 } 84 } 85 //交换行 86 if (maxj != i){ 87 for (int k = 0; k < n + 1; k++){ 88 a[k + i*(n + 1)] = a[k + maxj*(n + 1)] / max; 89 a[k + maxj*(n + 1)] = temp[k]; 90 } 91 } 92 else{ 93 for (int k = 0; k < n + 1; k++){ 94 a[k + i*(n + 1)] = a[k + i*(n + 1)] / max; 95 } 96 } 97 98 //求三角 99 double div = 0.0; 100 for (int j = i + 1; j < n; j++){ 101 div = a[i + j*(n + 1)]; 102 for (int k = i; k < n + 1; k++){ 103 a[k + j*(n + 1)] = a[k + j*(n + 1)] - a[k + i*(n + 1)] * div; 104 } 105 } 106 /*for (int i = 0; i < n; i++){ 107 for (int j = 0; j < n + 1; j++){ 108 cout << a[j + i*(n + 1)] << " "; 109 } 110 cout << endl; 111 } 112 cout << "-------------" << endl;*/ 113 } 114 115 /*for (int i = 0; i < n; i++){ 116 for (int j = 0; j < n + 1; j++){ 117 cout << a[j + i*(n + 1)] << " "; 118 } 119 cout << endl; 120 } 121 cout << "-------------" << endl;*/ 122 //反迭代求解 123 for (int i = n - 1; i >-1; i--){ 124 b[i] = a[i*(n + 1) + n]; 125 for (int j = n - 1; j > i; j--){ 126 b[i] = b[i] - a[j + i*(n + 1)] * b[j]; 127 } 128 } 129 /*for (int i = 0; i < n; i++){ 130 cout << b[i] << endl; 131 }*/ 132 return 0; 133 } 134 135 int LeastSquare(double* a, double* d, int c, int n) 136 { 137 /* 138 最小二乘法求超线性方程组 139 A*C=B 140 a表示增广矩阵A，d表示输出结果C，c表示未知数C的个数，n表示线性方程组的个数 141 例： 142 线性方程组： 143 6x+7y+9z=19 144 4x+5y+7z=15 145 6x+2y+3z=11 146 7x+3y+2z=7 147 则参数说明为： 148 d=[x; 149 y; 150 z] 151 a=[6 7 9 11； 152 4 5 7 15； 153 6 2 3 11； 154 7 3 2 7] 155 */ 156 if (a == NULL || d == NULL) 157 return -1; 158 double* equ = (double*)malloc(sizeof(double)*c*(c + 1)); 159 for (int i = 0; i < c; i++){ 160 for (int j = 0; j < c + 1; j++){ 161 equ[j + i*(c + 1)] = 0.0; 162 } 163 } 164 165 double* col = (double*)malloc(sizeof(double)*n); 166 if (equ == NULL) 167 return -2; 168 for (int i = 0; i < c; i++){ 169 for (int k = 0; k < n; k++){ 170 col[k] = a[i + k*n]; 171 cout << col[k] << " "; 172 } 173 cout << endl; 174 for (int j = 0; j < c + 1; j++){ 175 for (int l = 0; l < n; l++){ 176 equ[j + i*(c + 1)] = equ[j + i*(c + 1)] + col[l] * a[j + l*(c + 1)]; 177 } 178 } 179 } 180 for (int i = 0; i < c; i++){ 181 for (int j = 0; j < c + 1; j++){ 182 cout << equ[j + i*(c + 1)] << " "; 183 } 184 cout << endl; 185 } 186 187 GaussianElimination(equ, d, c); 188 return 0; 189 } 190 191 192 193 int main() 194 { 195 double a[16] = {6,7,9,19,4,5,7,15,6,2,3,11,7,3,2,7}; 196 //double equ[12] = { 2,-4,6,3,4,-9,2,5,1,-1,3,4 }; 197 double* d = (double*)malloc(sizeof(double) * 3); 198 LeastSquare(a,d,3,4); 199 //GaussianElimination(equ, d, 3); 200 for (int i = 0; i < 3; i++){ 201 cout << d[i] << " "; 202 } 203 cout << endl; 204 return 0; 205 }
这几个代码以前也写了很多，只是长久以后，经不起岁月侵蚀，希望以后经常复习，以资鼓励。
标签：
原文地址：http://www.cnblogs.com/gucolin/p/4383863.html
