#include "stdafx.h"
#include <stdio.h>
#include <conio.h>
#include <stdlib.h>
#include <cmath>
#include <iostream>
using namespace std
void polyfit(int n, double x[], double y[], int poly_n, double a[])
{
int i, j
double *tempx, *tempy, *sumxx, *sumxy, *ata
void gauss_solve(int n, double A[], double x[], double b[])
tempx = new double[n]
sumxx = new double[poly_n * 2 + 1]
tempy = new double[n]
sumxy = new double[poly_n + 1]
ata = new double[(poly_n + 1)*(poly_n + 1)]
for (i = 0i<ni++)
{
tempx[i] = 1
tempy[i] = y[i]
}
for (i = 0i<2 * poly_n + 1i++)
for (sumxx[i] = 0, j = 0j<nj++)
{
sumxx[i] += tempx[j]
tempx[j] *= x[j]
}
for (i = 0i<poly_n + 1i++)
for (sumxy[i] = 0, j = 0j<nj++)
{
sumxy[i] += tempy[j]
tempy[j] *= x[j]
}
for (i = 0i<poly_n + 1i++)
for (j = 0j<poly_n + 1j++)
ata[i*(poly_n + 1) + j] = sumxx[i + j]
gauss_solve(poly_n + 1, ata, a, sumxy)
delete [] tempx
tempx = NULL
delete [] sumxx
sumxx = NULL
delete [] tempy
tempy = NULL
delete [] sumxy
sumxy = NULL
delete [] ata
ata = NULL
}
void gauss_solve(int n, double A[], double x[], double b[])
{
int i, j, k, r
double max
for (k = 0k<n - 1k++)
{
max = fabs(A[k*n + k])/*find maxmum*/
r = k
for (i = k + 1i<n - 1i++)
if (max<fabs(A[i*n + i]))
{
max = fabs(A[i*n + i])
r = i
}
if (r != k)
for (i = 0i<ni++) /*change array:A[k]&A[r] */
{
max = A[k*n + i]
A[k*n + i] = A[r*n + i]
A[r*n + i] = max
}
max = b[k] /*change array:b[k]&b[r] */
b[k] = b[r]
b[r] = max
for (i = k + 1i<ni++)
{
for (j = k + 1j<nj++)
A[i*n + j] -= A[i*n + k] * A[k*n + j] / A[k*n + k]
b[i] -= A[i*n + k] * b[k] / A[k*n + k]
}
}
for (i = n - 1i >= 0x[i] /= A[i*n + i], i--)
for (j = i + 1, x[i] = b[i]j<nj++)
x[i] -= A[i*n + j] * x[j]
}
/*==================polyfit(n,x,y,poly_n,a)===================*/
/*=======拟合y=a0+a1*x+a2*x^2+……+apoly_n*x^poly_n========*/
/*=====n是数据个数 x y是数据值 poly_n是多项式的项数======*/
/*===返回a0,a1,a2,……a[poly_n]，系数比项数多一（常数项）=====*/
void main()
{
int n = 9, poly_n = 2
//double x[20] = { 1 ,2 , 3 , 4 , 7 ,8 ,9,11,12, 13,15,15,17 ,17, 19 ,18 ,20,19,20, 20 },
//y[20] = { 1,3,6 ,10 ,17,25,34,45,57,70,85,100,117,134,153,171,191,210 , 230 , 250 }
double x[9]={1,3,4,5,6,7,8,9,10},
y[9]={10,5,4,2,1,1,2,3,4}
double a[50]
polyfit(n, x, y, poly_n, a)
for (int i = 0i <poly_n + 1i++)/*这里是升序排列，Matlab是降序排列*/
cout <<"a[" <<i <<"]=" <<a[i] <<endl
}
运行结果，我是拟合的2次的，你可以拟合多次。
方程式：
0.267571*x*x-3.60531*x+13.4597=0
这个2次多项式的最低点还用我说吗，在那个区间上，自己代入；
#include <stdio.h>#include <math.h>
#define epsilon 1e-6
void nihe1(int n,int m,float sum_x,float sum_y,float sum_xy,float x2)
void nihe2(int n,int m,float sum_x,float sum_y,float sum_xy,float x2,float x2y,float x3,float x4)
int main(){
float x[100]={0.0}
float y[100]={0.0}
int n,i,flag=1
float sum_y=0.0,sum_x=0,x2=0,sum_xy=0.0,x3=0,x4=0,x2y=0.0
printf("请你输入需要测试的数据(先输入x[]，后输入y[])的个数:")
scanf("%d",&n)
for(i = 0i <ni++){
scanf("%f",&x[i])}
for(i = 0i <ni++){
scanf("%f",&y[i])}
for(i = 0i <ni++){
sum_x += x[i]
sum_y += y[i]
sum_xy += x[i]*y[i]
x2 += x[i]*x[i]
x2y += x[i]*x[i]*y[i]
x3 += x[i]*x[i]*x[i]
x4 += x[i]*x[i]*x[i]*x[i]}
printf("---------------请你输入的要拟合的函数------------------\n")
printf(" 1、拟合一次函数\n")
printf(" 2、拟合二次函数\n")
scanf("%d",&flag)
switch(flag){
case 1:
nihe1(n,flag+1,sum_x,sum_y,sum_xy,x2)break
case 2:
nihe2(n,flag+1,sum_x,sum_y,sum_xy,x2,x2y,x3,x4)break
default:
printf("ERROR\n")}
return 0}
void nihe1(int n,int m,float sum_x,float sum_y,float sum_xy,float x2){
int i,k,j
float t,s=0
float a[2][3] = {{(float)n,sum_x,sum_y},{sum_x,x2,sum_xy}}
n=m
//if(m == 3)
// a[3][4] = {{n,sum_x,sum_y},{sum_x,x2,x3,sum_xy},{x2,x3,x4,x2y}}
for(k=0k<n-1k++) {
for(i=k+1i<ni++)
if( abs((int)a[i][k]) >abs((int)a[k][k]) )
for(j=kj<n+1j++) {
t=a[k][j]
a[k][j]=a[i][j]
a[i][j]=t}
if( abs((int)a[k][k]) <epsilon) {
printf("\nError,主元消去法 cann't be durable,break at %d!\n",k+1)
return}
for(i=k+1i<ni++){
a[i][k]=a[i][k] / a[k][k]
for(j=k+1j<n+1j++)
a[i][j]=a[i][j]-a[i][k] * a[k][j]}}
a[n-1][n]=a[n-1][n] / a[n-1][n-1]
for(k=n-2k>=0k--) {
s=0
for(j=k+1j<nj++)
s+=a[k][j]*a[j][n]
a[k][n]=( a[k][n]-s ) / a[k][k]}
printf("\n*****The Result*****\n")
for(i=0i<ni++)
printf(" x[%d]=%.4f\n",i+1,a[i][n])
printf("函数为：p(x) = %.4f + (%.4f)*x\n",a[0][n],a[1][n])
getchar()}
void nihe2(int n,int m,float sum_x,float sum_y,float sum_xy,float x2,float x2y,float x3,float x4){
int i,k,j
float t,s=0
float a[3][4]=
{{(float)n,sum_x,x2,sum_y},{sum_x,x2,x3,sum_xy},{x2,x3,x4,x2y}}
n=m
for(k=0k<n-1k++) {
for(i=k+1i<ni++)
if( abs((int)a[i][k]) >abs((int)a[k][k]) )
for(j=kj<n+1j++) {
t=a[k][j]
a[k][j]=a[i][j]
a[i][j]=t}
if( abs((int)a[k][k]) <epsilon) {
printf("\nError,主元消去法 cann't be durable,break at %d!\n",k+1)
return}
for(i=k+1i<ni++){
a[i][k]=a[i][k] / a[k][k]
for(j=k+1j<n+1j++)
a[i][j]=a[i][j]-a[i][k] * a[k][j]} }
a[n-1][n]=a[n-1][n] / a[n-1][n-1]
for(k=n-2k>=0k--) {
s=0
for(j=k+1j<nj++)
s+=a[k][j]*a[j][n]
a[k][n]=( a[k][n]-s ) / a[k][k]}
printf("\n*****The Result*****\n")
for(i=0i<ni++)
printf(" x[%d]=%.4f\n",i+1,a[i][n])
printf("函数为：p(x) = %.4f + (%.4f)*x + (%.4f)*x*x\n",a[0][n],a[1][n],a[2][n])
getchar()}
