多项式拟合

《数值计算》实验报告学院：软件学院专业：软件工程班级：12级1班实验名称多项式拟合姓名杜倩学号1402120110 成绩实验报告内容要求：实验三：编写多项式拟合程序。

并用该程序解决下列问题：假定某天的气温变化记录如下表，试用最小二乘方法找出这一天的气温变化规律。

ht/ 1 2 3 4 5 6 7 8 9 10 11 12 13 CT︒/14 14 14 14 15 16 18 20 22 23 25 28 31 ht/14 15 16 17 18 19 20 21 22 23 24 CT︒/32 31 29 27 25 24 22 20 18 17 16考虑下列类型函数，计算误差平方和，并作图比较效果。

1.二次函数2.三次函数3.四次函数4.函数))((2ctbaeC--=（提高：非线性拟合问题）一.实验目的通过比较不同次数的多项式拟合效果，了解多项式拟合原理。

二.实验原理最小二乘法三.实验环境 vc6.0四.实验过程（编写的程序）#include "stdafx.h"#include "stdlib.h"#include "string.h"#include "math.h"static double * qr_fraction(double *a, int m, int n, double **q);//QR分解static double * up_tria_inv_n_ord(double *t,int n);//n阶上三角矩阵求逆static double * array_mut(double *A, double *x, int m, int n);//m * n型,矩阵乘以向量static double * array_trans_mut(double *A, double *x, int m, int n);//m * n型矩阵转置乘以向量static double * polyfit(double *x, double *y, int len,int order);//根据待拟合向量,给出拟合的order阶多项式系数static double polyval(double *coef, int order, double x);//多项式带入数值进行插值运算static double *matrix_trans_mul(double *A, double *B,int order);static double *inv(double *A, int order);// order阶矩阵求逆int main(){double x[10] = {1,2,3,4,5,6,7,8,9,10};double y[10] = {1,3,7,8,5,3.5,1,-2.5,-8,-16};double *coef = NULL;int i = 0;/* 拟合多项式f(t) = coef[0]*t^3 + coef[1]*t^2 + coef[2]*t^1 + coef[3] */coef = polyfit(x, y, sizeof(x)/sizeof(x[0]), 5);//5阶多项式拟合//以下结果与matlab拟合的结果一致printf("coef[0-5]=%lf%lf%lf%lf%lf%lf\n",coef[0],coef[1],coef[2],coef[3],coef[4],coef[5]); double yy = polyval(coef, 5, 11);free(coef);//polyfit 返回的是堆内存double A[4][4] = {{1,2,3,4},{2,3,4,7},{1,1,0,11},{-1,7,0,11}};double *INV_A = inv(&A[0][0], 4);for(i = 0; i < 4; i++){printf("%lf %lf %lf %lf\n",INV_A[4*i+0],INV_A[4*i+1],INV_A[4*i+2],INV_A[4*i+3]); }return 0;}static double * qr_fraction(double *a, int m, int n, double **q)//m >= n,列满秩{int i = 0;int j = 0;int k = 0;//调用malloc()必须自己free()，不然多次调用就会内存泄露的double *Q = (double *)malloc(sizeof(double)*m*n); //Qm*ndouble *R = (double *)malloc(sizeof(double)*n*n); //Rn*nmemset(R,0,sizeof(double)*n*n);for(i = 0; i < n; i++) //A矩阵共n列{double tmp = 0;for(j = 0; j < m; j++)//求A矩阵各列的模,m个元素平方和，再开方{tmp += a[n*j + i]*a[n*j + i]; // i = 0时，a[0] a[n] a[2*n] ... a[(m-1)*n]}tmp = sqrt(tmp);//得到矩阵列的模R[i*n + i] = tmp;//R[i][i] 即R矩阵的对角元///////////////////////////////第一列的模得到后，就可以得到归一化的Q1(Q的第一列)for(j = 0; j < m; j++)//求A矩阵各列的模,m个元素平方和，再开方{Q[n*j + i] = a[n*j + i] / tmp; // i = 0时，a[0] a[n] a[2*n] ... a[(m-1)*n]}for(j = i + 1; j < n; j++)//{tmp = 0;for(k = 0; k < m; k++) //R[i][j] = <qi, aj> qi与aj内积Q的j列，A的k列{tmp += Q[n*k + i] * a[n*k + j];}R[n*i + j] = tmp; //得到R[i][j]for(k = 0; k < m; k++){ a[n*k + j] = a[n*k + j] - Q[n*k + i] * R[n*i + j];//a[k][j] = a[k][j] - Q[k][i] * R[i][j] }}}*q = Q;return R;}static double * up_tria_inv_n_ord(double *t,int n){int i = 0;int j = 0;int k = 0;int m = 0;double tmp = 0;//调用malloc()必须自己free()，不然多次调用就会内存泄露的double *a = (double *)malloc(sizeof(double)*n*n); //an*ndouble *inv = (double *)malloc(sizeof(double)*n*n); //Rn*nmemset(a,0,sizeof(double)*n*n);memset(inv,0,sizeof(double)*n*n);for(i = 0; i < n; i++){inv[i*n + i] = 1/t[i*n + i];//即inv[i][i] = 1/t[i][i];}for(i = 0; i < n - 1; i++)//n - 1项,ai(i+1) (i = 1，n - 1){a[i*n + i + 1] = -inv[(i+1)*n + i + 1] * t[i*n + i + 1];//即inv[i][i] = 1/a[i][i]; }m = 2;while(m < n)//这部分需要认真考虑，先算什么，再算什么{for(i = 0,j = i + m; i < n - 2; i++)//n - 2项{tmp = 0;for(k = i + 1; k < j; k++){tmp += a[i*n + k] * t[k*n + j];}a[i*n + j] = -inv[j*n + j] * (t[i*n + j] + tmp);//即inv[i][i] = 1/a[i][i];}m++;}for(i = 0; i < n - 1; i++){for(j = i + 1; j < n; j++){inv[i*n + j] = inv[i*n + i] * a[i*n + j];}}free(a);return inv;}/*********************************************************************** ******double * array_mut(double *A, double *x, int len)参数：A指向m*n矩阵，x指向向量，len为向量长度返回：指针，指向向量y = Ax************************************************************************ ****///static double * array_mut(double *A, double *x, int m, int n)//m * n型//{// static double y[POLYFIT_NUM] = {0}; //不会超过1000个点拟合一个点的// for(int i = 0; i < m; i++)// {// y[i] = 0; //这条不能少！多次调用就有问题，static的后遗症// for(int j = 0; j < n; j++)// {// y[i] += A[n*i + j] * x[j];// }// }//// return y;//}static double * array_mut(double *A, double *x, int m, int n)//m * n型{double *y = (double *)malloc(sizeof(double)*m); //不会超过1000个点拟合一个点的for(int i = 0; i < m; i++){y[i] = 0;for(int j = 0; j < n; j++){y[i] += A[n*i + j] * x[j];}}return y;}/*********************************************************************** ******double * array_trans_mut(double *A, double *x, int len)参数：A指向m*n矩阵，x指向向量，len为向量长度返回：指针，指向向量y = A^T * x 即A的转置乘以x向量************************************************************************ ****///static double * array_trans_mut(double *A, double *x, int m, int n)//m * n型//{// static double y[POLYFIT_NUM] = {0}; //// for(int i = 0; i < n; i++) //n*m * n// {// y[i] = 0;// for(int j = 0; j < m; j++)// {// y[i] += A[n*j + i] * x[j]; //多次调用就有问题，static的后遗症// }// }//// return y;//}static double * array_trans_mut(double *A, double *x, int m, int n)//m * n型{double *y = (double *)malloc(sizeof(double)*n);for(int i = 0; i < n; i++) //n*m * n{y[i] = 0;for(int j = 0; j < m; j++){y[i] += A[n*j + i] * x[j];}}return y;}/*********************************************************************** ******void polyfit(double *x, double *y, int len, int order)参数：x,y为待拟合的向量，len为向量长度,order为多项式阶（次）数返回：coef为order阶多项式系数************************************************************************ ****/static double * polyfit(double *x, double *y, int len, int order){int i = 0;int j = 0;int colum = order + 1;double *coef = NULL;double *A = (double *)malloc(sizeof(double)*len*colum); //A len 行order + 1列，polyfit 里的V矩阵double *Q = NULL;double *R = NULL;double *inv_R = NULL;//以下代码生成A矩阵,A其实为范德蒙行列式，最后一列为1for(i = 0; i < len; i++) //i对应行{A[i*colum + colum - 1] = 1.0;}for(j = colum - 2 ; j > -1; j--) //j对应列0 -- 2 列{for(i = 0; i < len; i++) //i对应行{A[i*colum + j] = x[i] * A[i*colum + j + 1];}}//至此,A矩阵生成完毕R = qr_fraction(A, len, colum, &Q);inv_R = up_tria_inv_n_ord(R,colum); //返回的是堆内存//free(R)就出错啦！//参数p计算公式为：p = R^(-1) * Q^T * y//故需先计算Q^T * ydouble *yy = array_trans_mut(Q, y, len, colum);coef = array_mut(inv_R, yy, colum, colum);//coef = array_mut(inv_R, y, colum, colum);free(yy);free(inv_R);free(Q);free(R);free(A);return coef;}static double polyval(double *coef, int order, double x){int i = 0;double sum = 0;for(i = 0; i < order + 1; i++){sum += coef[i] * pow(x, order - i);}return sum;}/*****************************************************************************A*B的转置************************************************************************ ****/static double *matrix_trans_mul(double *A, double *B,int order){double *result = (double *)malloc(sizeof(double) * order * order);int i = 0;int j = 0;int k = 0;double tmp = 0;for(i = 0; i < order; i++){for(j = 0; j < order; j++){tmp = 0;for(k = 0; k < order; k++){//*(result + i*order + j) = (*(A + i*order + k)) * (*(B + j*order + k));tmp += (*(A + i*order + k)) * (*(B + j*order + k));}*(result + i*order + j) = tmp;}}return result;}/*********************************************************************** ******n阶矩阵求逆************************************************************************ ****/static double *inv(double *A, int order){//double *INV = (double *)malloc(sizeof(double) * order * order);double *Q = NULL;double *R = NULL;double *INV_R = NULL;double *INV_A = NULL;R = qr_fraction(A, order, order, &Q);INV_R = up_tria_inv_n_ord(R, order);INV_A = matrix_trans_mul(INV_R, Q, order);free(Q);free(R);free(INV_R);return INV_A;}五、实验结果及分析六.实验反思对用c语言编程完成该题，在编写过程查找相关资料，更深入了解多项拟合，相比之前学的方法稍微简单一些，变成却复杂。