最小二乘法用于多项式的拟合及C语言程序实现
改写自:
N元线性函数拟合的C++实现_c++连续分段线性函数拟合-CSDN博客
最小二乘法用于直线,多项式,圆,椭圆的拟合及程序实现_椭圆拟合系数的应用-CSDN博客
将C++代码改为C语言
注意要做防止除数为0的保护
1 #include <stdio.h> 2 #include "stdlib.h" 3 #include "math.h" 4 //#include <vector> 5 //using namespace std; 6 7 #define MAX_EXP 4 /* x最高次幂 */ 8 #define MATRIX_DIM (MAX_EXP + 1) /* 矩阵阶数 */ 9 #define SMPNUM 5 /* 采样个数 */ 10 11 struct point 12 { 13 double x; 14 double y; 15 }; 16 17 /* 采样点想, y */ 18 float sampleX[SMPNUM] = {0.0}; 19 float sampleY[SMPNUM] = {0.0}; 20 21 void getFile(char *File); //获取文件数据 22 void getCoeff(float sampleX[SMPNUM], float sampleY[SMPNUM], float coeff[MATRIX_DIM]); //获取系数 23 24 int main() 25 { 26 int i; 27 char *File = "XY.txt"; 28 //vector<point> sample; 29 //doubleVector coefficient; 30 //sample = getFile(File); 31 getFile(File); 32 printf("拟合多项式阶数n="); 33 //scanf("%d", &n); 34 // coefficient = getCoeff(sample, n); 35 //n = 3; 36 37 float coeff[MATRIX_DIM] = {0}; 38 getCoeff(sampleX, sampleY, coeff); 39 printf("\n拟合矩阵的系数为:\n"); 40 for (i = 0; i < MATRIX_DIM; i++) 41 { 42 printf("a%d = %lf\n", i, coeff[i]); 43 } 44 45 return 0; 46 } 47 //矩阵方程 48 void getCoeff(float sampleX[SMPNUM], float sampleY[SMPNUM], float coeff[MATRIX_DIM]) 49 { 50 double sum; 51 int i, j, k; 52 53 float matX[MATRIX_DIM][MATRIX_DIM]; //公式3左矩阵 54 float matY[MATRIX_DIM][MATRIX_DIM]; //公式3右矩阵 55 float temp2[MATRIX_DIM]; 56 //公式3左矩阵 57 for (i = 0; i < MATRIX_DIM; i++) 58 { 59 for (j = 0; j < MATRIX_DIM; j++) 60 { 61 sum = 0.0; 62 for (k = 0; k < SMPNUM; k++) 63 { 64 sum += pow(sampleX[k], j + i); 65 } 66 matX[i][j] = sum; 67 } 68 } 69 70 //打印matFunX矩阵 71 printf("matFunX矩阵:\n"); 72 for (i = 0; i < MATRIX_DIM; i++) 73 { 74 for (j = 0; j < MATRIX_DIM; j++) 75 { 76 printf("%f\t", matX[i][j]); 77 } 78 printf("\n"); 79 } 80 81 //公式3右矩阵 82 for (i = 0; i < MATRIX_DIM; i++) 83 { 84 //temp.clear(); 85 sum = 0; 86 for (k = 0; k < SMPNUM; k++) 87 { 88 sum += sampleY[k] * pow(sampleX[k], i); 89 } 90 matY[i][0] = sum; 91 } 92 //printf("matFunY.size=%d\n", matFunY.size()); 93 //打印matFunY的矩阵 94 printf("matFunY的矩阵:\n"); 95 for (i = 0; i < MATRIX_DIM; i++) 96 { 97 printf("%f\n", matY[i][0]); 98 } 99 //矩阵行列式变换,将matFunX矩阵变为下三角矩阵,将matFunY矩阵做怎样的变换呢? 100 //AX=Y在将X变换为下三角矩阵X'时,是左右两边同乘ratio 101 double num1, num2, ratio; 102 for (i = 0; i < MATRIX_DIM; i++) 103 { 104 num1 = matX[i][i]; 105 for (j = i + 1; j < MATRIX_DIM; j++) 106 { 107 num2 = matX[j][i]; 108 ratio = num2 / num1; 109 for (k = 0; k < MATRIX_DIM; k++) 110 { 111 matX[j][k] = matX[j][k] - matX[i][k] * ratio; 112 } 113 matY[j][0] = matY[j][0] - matY[i][0] * ratio; 114 } 115 } 116 //打印matFunX行列变化之后的矩阵 117 printf("matFunX行列变换之后的矩阵:\n"); 118 for (i = 0; i < MATRIX_DIM; i++) 119 { 120 for (j = 0; j < MATRIX_DIM; j++) 121 { 122 printf("%f\t", matX[i][j]); 123 } 124 printf("\n"); 125 } 126 //打印matFunY行列变换之后的矩阵 127 printf("matFunY行列变换之后的矩阵:\n"); 128 for (i = 0; i < MATRIX_DIM; i++) 129 { 130 printf("%f\n", matY[i][0]); 131 } 132 //计算拟合曲线的系数 133 //doubleVector coeff(n + 1, 0); 134 for (i = MATRIX_DIM - 1; i >= 0; i--) 135 { 136 if (i == MATRIX_DIM - 1) 137 { 138 coeff[i] = matY[i][0] / matX[i][i]; 139 } 140 else 141 { 142 for (j = i + 1; j < MATRIX_DIM; j++) 143 { 144 matY[i][0] = matY[i][0] - coeff[j] * matX[i][j]; 145 } 146 coeff[i] = matY[i][0] / matX[i][i]; 147 } 148 } 149 return; 150 } 151 152 //获取文件数据 153 void getFile(char *File) 154 { 155 int i = 0, j = 0, k = 0; 156 //vector<point> dst; 157 158 FILE *fp = fopen(File, "r"); 159 160 if (fp == NULL) 161 { 162 printf("Open file error!!!\n"); 163 exit(0); 164 } 165 166 point temp; 167 double num; 168 169 while (fscanf(fp, "%lf", &num) != EOF) 170 { 171 if (i % 2 == 0) 172 { 173 temp.x = num; 174 sampleX[j++] = num; 175 176 177 } 178 else 179 { 180 temp.y = num; 181 sampleY[k++] = num; 182 } 183 i++; 184 } 185 fclose(fp); 186 return; 187 //return dst; 188 }
【推荐】编程新体验,更懂你的AI,立即体验豆包MarsCode编程助手
【推荐】凌霞软件回馈社区,博客园 & 1Panel & Halo 联合会员上线
【推荐】抖音旗下AI助手豆包,你的智能百科全书,全免费不限次数
【推荐】博客园社区专享云产品让利特惠,阿里云新客6.5折上折
【推荐】轻量又高性能的 SSH 工具 IShell:AI 加持,快人一步
· 实操Deepseek接入个人知识库
· CSnakes vs Python.NET:高效嵌入与灵活互通的跨语言方案对比
· 【.NET】调用本地 Deepseek 模型
· Plotly.NET 一个为 .NET 打造的强大开源交互式图表库
· 上周热点回顾(2.17-2.23)