最小二乘法拟合n阶曲线
#include <math.h>
#include <stdio.h>
#define N 4
void arryans(double a[][N], double b[N], double x[N], int n)
{
int i, j, k, mi;
double mx, tmp;
for (i = 0; i < n - 1; i++) {
for (j = i + 1, mi = i, mx = fabs(a[i][i]); j < n; j++)
if (fabs(a[j][i]) > mx) {
mi = j;
mx = fabs(a[j][i]);
}
if (i < mi) {
tmp = b[i];
b[i] = b[mi];
b[mi] = tmp;
for (j = i; j < n; j++) {
tmp = a[i][j];
a[i][j] = a[mi][j];
a[mi][j] = tmp;
}
}
for (j = i + 1; j < n; j++) {
tmp = -a[j][i] / a[i][i];
b[j] += b[i] * tmp;
for (k = i; k < n; k++)
a[j][k] += a[i][k] * tmp;
}
}
x[n - 1] = b[n - 1] / a[n - 1][n - 1];
for (i = n - 2; i >= 0; i--) {
x[i] = b[i];
for (j = i + 1; j < n; j++)
x[i] -= a[i][j] * x[j];
x[i] /= a[i][i];
}
}
void arrytrans(double a[][N], double b[][N], int m, int n)
{
int i, j;
for (i = 0; i < m; i++)
for (j = 0; j < n; j++)
b[j][i] = a[i][j];
}
void arrymul(double a[][N], double b[][N], double c[][N], int m, int l, int n)
{
int i, j, k;
for (i = 0; i < m; i++)
for (j = 0; j < n; j++) {
c[i][j] = 0;
for (k = 0; k < l; k++)
c[i][j] += a[i][k] * b[k][j];
}
}
void arrycon(double a[], double b[][N], int m, int n)
{
int i, j, k;
for (i = 0; i < m; i++)
for (j = 0; j < n; j++) {
b[i][j] = 1;
for (k = 0; k < j; k++)
b[i][j] *= a[i];
}
}
void arryzeng(double a[][N], double b[], double c[][N], int m, int n)
{
int i, j;
for (i = 0; i < m; i++)
for (j = 0; j <= n; j++) {
if (j == n)
c[i][j] = b[i];
else
c[i][j] = a[i][j];
}
}
void arrymuls(double a[][N], double b[], double c[], int m, int n) {
int i, j;
for (i = 0; i < m; i++) {
c[i] = 0;
for (j = 0; j < n; j++)
c[i] += a[i][j] * b[j];
}
}
double power(double a[], double x, int i, int j)
{
if (i >= j)
return a[i] * x;
else
return x * power(a, x, i + 1, j) + a[i];
}
void CalLS(double x[N], double y[N], int n, double k[N]) {
int m;
double a[N][N], b[N][N], c[N][N], d[N];
n++;
m = N;
arrycon(x, a, m, n);
arrytrans(a, b, m, n);
arrymul(b, a, c, n, m, n);
arrymuls(b, y, d, n, m);
arryans(c, d, k, n);
}