公式原理,请移步二阶回归曲线拟合
以下是我改编的4阶实现
从data.txt中,读入指定数量(x,y),进行拟合四阶曲线,打印输出各系数值。
改编成n阶,仅需要修改do{}while();循环即可。
#include <stdio.h>
#include<stdlib.h>
#include <math.h>
#define N 1e-13
#define T 4 //阶数
#define NUM 101 //数据数量
#define OUTPUT
int main() {
int i, j, k;
//初始化数组
long double x[NUM + 1] = { 0 };
long double y[NUM + 1] = { 0 };
//读入数据
FILE *fp;
fopen_s(&fp, "data.txt", "r");
int checkNum;
for (i = 0; i <= NUM; i++) {
if (!fscanf_s(fp, "%lf", &x[i]))break;
if (!fscanf_s(fp, "%lf", &y[i]))break;
}
fclose(fp);
#ifdef OUTPUT
for (i = 0; i < 10; i++) {
printf("%16.9lf ", x[i]);
printf("%16.9lf\n", y[i]);
}
printf("\n\n");
#endif
long double a[T + 1] = { 0 }; // 参数
long double m[T + 1] = { 0 };// 中间变量
long double z[T + 1] = { 0 };// 中间变量
long double sumXY[T * 2 + 1][2] = { 0 };//存储各幂值累加和,如sumxy[3][3] 存储 (x ^ 3 * y ^ 3)的累加和
for (j = 0; j <= T * 2 + 1; j++) {
for (k = 0; k <= 1; k++) {
for (i = 0; i < NUM; i++)
{
sumXY[j][k] += (pow(x[i], j) * pow(y[i], k));
}
}
}
#ifdef OUTPUT
for (i = 0; i < 5; i++) {
printf("%16.9lf ", sumXY[i][0]);
printf("%16.9lf\n", sumXY[i][1]);
}
printf("\n\n");
#endif
int countNum = 0;
do {
m[4] = a[4];
a[4] = (sumXY[4][1] - a[0] * sumXY[4][0] - a[1] * sumXY[5][0] - a[2] * sumXY[6][0] - a[3] * sumXY[7][0])/ sumXY[8][0];
z[4] = (a[4] - m[4])*(a[4] - m[4]);
m[3] = a[3];
a[3] = (sumXY[3][1] - a[0] * sumXY[3][0] - a[1] * sumXY[4][0] - a[2] * sumXY[5][0] - a[4] * sumXY[7][0]) / sumXY[6][0];
z[3] = (a[3] - m[3])*(a[3] - m[3]);
m[2] = a[2];
a[2] = (sumXY[2][1] - a[0] * sumXY[2][0] - a[1] * sumXY[3][0] - a[3] * sumXY[5][0] - a[4] * sumXY[6][0]) / sumXY[4][0];
z[2] = (a[2] - m[2])*(a[2] - m[2]);
m[1] = a[1];
a[1] = (sumXY[1][1] - a[0] * sumXY[1][0] - a[2] * sumXY[3][0] - a[3] * sumXY[4][0] - a[4] * sumXY[5][0]) / sumXY[2][0];
z[1] = (a[1] - m[1])*(a[1] - m[1]);
m[0] = a[0];
a[0] = (sumXY[0][1] - a[1] * sumXY[1][0] - a[2] * sumXY[2][0] - a[3] * sumXY[3][0] - a[4] * sumXY[4][0]) / NUM;
z[0] = (a[0] - m[0])*(a[0] - m[0]);
printf("%d\n", countNum++);//统计运算次数
} while ((z[0] > N || z[1] > N || z[2] > N || z[3] > N || z[4] > N ));
//打印输出
printf("a=%9.6f,\nb=%9.6f,\nc=%9.6f\nd=%9.6f\ne=%9.6f\n", a[4], a[3], a[2], a[1], a[0]);
printf("四阶拟合方程为 y=%9.6fx^4 + %9.6fx^3 + %9.6fx^2 + %9.6fx^1 + %9.6f ", a[4], a[3], a[2], a[1], a[0]);
getchar();
return 0;
}