多项式拟合

《数值计算》实验报告

学院：软件学院专业：软件工程班级：12级1班

实验名称

多项式拟合

姓名杜倩学号1402120110 成绩

实验报告内容要求：

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

h

t/ 1 2 3 4 5 6 7 8 9 10 11 12 13 C

T︒/14 14 14 14 15 16 18 20 22 23 25 28 31 h

t/14 15 16 17 18 19 20 21 22 23 24 C

T︒/32 31 29 27 25 24 22 20 18 17 16

考虑下列类型函数，计算误差平方和，并作图比较效果。

1.二次函数

2.三次函数

3.四次函数

4.函数

)

)

(

(2

c

t

b

ae

C-

-

=（提高：非线性拟合问题）

一.实验目的通过比较不同次数的多项式拟合效果，了解多项式拟合原理。

二.实验原理最小二乘法

三.实验环境 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*n

double *R = (double *)malloc(sizeof(double)*n*n); //Rn*n

memset(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矩阵的对角元