n阶多项式拟合与n阶矩阵求逆的C语言实现

最新推荐文章于 2024-05-30 09:42:44 发布
最新推荐文章于 2024-05-30 09:42:44 发布
阅读量3.1k 收藏 9
点赞数 1
分类专栏： C/C++ 文章标签： c++

/********************************************************************************* 
该程序实现对n维输入x，n维输出y，给出n阶多项式拟合系数，功能和matlab中polyfit(x,y,n)类似，   
涉及知识有QR分解、AX=Y的最小二乘问题、上三角矩阵求逆问题。
其中QR分解 采用施密特正交法， 
AX=Y的最小二乘问题以下类似这样一个问题：3维空间中一个向量，和x、y平面中的哪个向量的
欧式距离最小？这个问题就很形象了，答案也很简单，即该向量在  x、y平面的投影向量。所以AX
 =Y的最小二乘解与此类似，在此不多说了。
 上三角矩阵求逆这个是有公式的，下面提到的一篇论文里给出了，不过我推导了一下，发现几处打
印错误。
*********************************************************************************/

#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;
}

/**************************************************************************/
/*****************************************************************************
double * qr_fraction(double *a, int m, int n, double **q)
参数： a为待分解m*n矩阵，q用来回传Q阵
返回：R阵   满足a = QR
说明：算法采用施密特正交法，参考《高等数值算法与应用（三）》
****************************************************************************/
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矩阵的对角元
/
//第一列的模得到后，就可以得到归一化的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)
参数： t指向n * n上三角阵
返回： 指向逆矩阵的指针
说明：上三角阵求逆具体算法见论文《三角矩阵求逆的一种方法》,不过论文里面有一个地方印刷错误，具体地方是：
      a[i][j] = -1/t[j][j] * (t[j][j] + sum() ,  应该为a[i][j] = -1/t[j][j] * (t[i][j] + sum() ),
****************************************************************************/
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*n
    double *inv = (double *)malloc(sizeof(double)*n*n); //Rn*n

memset(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其实为范德蒙行列式，最后一列为 1
for(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 * y
double *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;
}
Qrui
关注
  • 1
    点赞
  • 踩
  • 9
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
专栏目录
c++ 三次多项式拟合_多项式回归（Polynomial Curve Fitting）
weixin_39966465的博客
11-22 1739
PRML这序章写的太好了，用代码实现一下切身体会一下知识点.假设有N个观察样例x，我们的目标函数是 (2 )加上一个高斯分布的噪声.这与我们实际工作中遇到的数据集其实也很相似：数据有一定的规律性，但每个数据点也都存在一定的噪声。噪声的产生可能是由于数据的随机性造成的，但更多的情况下是因为有我们未观察到的规律。我们的目标是利用多项式去探索隐含在 (2 )下的方程式。更蛋疼的是我们观察到的...
n次多项式拟合
热门推荐
jayloncheng的博客
05-13 3万+
n次多项式拟合一、用最小二乘法进行n次拟合原理给定一组数据作n次拟合曲线，转化为均方误差的极小问题。得到法方程：，其中解法方程得到a0，a1,...,an.注意：当拟合多项式阶数n&gt;5时，法方程的系数矩阵时病态的。二、算例数据集：x = -3,-2,-1,0,1,2,3y = 4,2,3,0,-1,-2,-5程序实现：#n次拟合曲线：y=a0+a1x+a2x^2+...+anx^n impo...
纯C++做多项式拟合
最新发布
qq_37033647的博客
05-30 326
纯C++的最小二乘法拟合多项式曲线
C++——多项式拟合
GeekZW的博客
07-20 1万+
                                           C++——多项式拟合 目标：利用C++对txt或者xml中的数据，进行高阶或低阶多项式拟合         为方便以后查找，代码以及详细资料已打包，并上传至云盘（链接：https://pan.baidu.com/s/1bvUBIoxv7Avxeq_Cz6xOZQ 密码：u9qe） 打包的内容如下：  ...
c++多项式拟合
06-12
#ifndef FUNCTION_H_ #define FUNCTION_H_ #include #include #include "polyfit.h" #include using namespace std; dxs::dxs() { ifstream fin("多项式拟合.txt"); fin>>n; x=new float[n]; y=new float[n]; for(int i=0;i>x[i]; } for(i=0;i>y[i]; } cout<>nn; m=nn+1; u=new float*[m]; for(i=0;i<m;i++) { u[i]=new float[m+1]; }//创建m行，m+1列数组 } void dxs::dfine() { for(int i=0;i<m;i++) { for(int j=0;j<m+1;j++) { u[i][j]=0; } } for(i=0;i<m;i++) { for(int j=0;j<m;j++) { for(int k=0;k<n;k++) { u[i][j]=u[i][j]+pow(x[k],j+i); } } } for(i=0;i<m;i++) { for(int k=0;k<n;k++) { u[i][m]=u[i][m]+pow(x[k],i)*y[k]; } } } void dxs::show() { for(int i=0;i<m;i++) { for(int j=0;j<m+1;j++) { cout<<u[i][j]<<" ";//<<endl; } cout<<endl; } ////显示具有m行m+1列u数组的各元素值 } void dxs::select_main(int k,float **p,int m) { double d; d=*(*(p+k)+k); //cout<<d; int l=k; int i=k+1; for(;i fabs(d)) { d=*(*(p+i)+k); l=i; } else continue; } if(d==0) cout<<"错误"; else { if(k!=l) { for(int j=k;j<m+1;j++) { double t; t=*(*(p+l)+j); *(*(p+l)+j)=*(*(p+k)+j); *(*(p+k)+j)=t; } } } } void dxs::gaosi() { for(int k=0;k<m;k++) { select_main(k,u,m);//调用列主元函数 for(int i=1+k;i<m;i++) { // *(*(p+i)+k)=(float) *(*(p+i)+k) / *(*(p+k)+k); u[i][k]=(float) u[i][k] / u[k][k]; } for(i=k+1;i<m;i++) { for(int j=k+1;j=0;i--) { float a=0; for(int j=i+1;j<m;j++) { //a=a + (*(*(p+i)+j) * *(*(p+j)+m)); a=a+u[i][j] * u[j][m]; } //*(*(p+i)+n-1)= (*(*(p+i)+n-1) - a) / *(*(p+i)+i); u[i][m]= (u[i][m] -a) / u[i][i]; } cout<<"方程组的解为："<<endl; for(i=0;i<m;i++) { cout<<"a"<<i+1<<"="; cout<<u[i][m]<<endl; // l[i]=*(*(p+i)+n-1); } cout<<"y="<<u[0][m]; for(i=1;i<m;i++) { cout<<showpos<<u[i][m]<<"x"; if(i!=1)cout<<"^"<<noshowpos<<i; } cout<<endl; } dxs::~dxs() { delete[]x,y; delete []*u; } #endif
多项式拟合c++程序
01-18
对数据进行多项式拟合，可得拟合多项式的系数，法方程的系数矩阵，平方误差。
c++ 多项式拟合算法
Kyle
01-18 3504
#ifndef CZY_MATH_FIT #define CZY_MATH_FIT #include <vector> /* 多项式拟合 */ namespace czy{ /// /// \brief 曲线拟合类 /// class Fit{ std::vector<double> factor; ///<拟合后的方程系数 double ssr; ...
N*N阶方阵求逆矩阵C代码
07-28
N*N阶方阵求逆矩阵C代码
vc++实现的多项式拟合程序
09-05
vc++实现的多项式拟合程序，可绘出拟合后的曲线
多项式拟合的C实现
07-27
多项式拟合的C实现,Linux下使用Gcc或者g++编译
C++实现最小二乘法一元回归和多项式拟合
03-23
此类的具体使用方法见：http://blog.csdn.net/czyt1988/article/details/21743595 在进行曲线拟合时用的最多的是最小二乘法，其中以一元函数（线性）和多元函数（多项式）居多，下面这个类专门用于进行多项式拟合，可以根据用户输入的阶次进行多项式拟合，算法来自于网上，和GSL的拟合算法对比过，没有问题。此类在拟合完后还能计算拟合之后的误差：SSE(剩余平方和)，SSR(回归平方和)，RMSE(均方根误差)，R-square(确定系数)。
多项式拟合
04-28
多项式拟合,源代码，完整工程文件，VS2010
离散点的C++多项式曲线拟合函数
01-24
其中包含部分原始数据，直接运行即可。
C程序——三次样条插值算法、最小二乘法多项式拟合等
10-13
三次样条插值（自然边界条件）算法、最小二乘法曲线拟合、多项式相关系数的计算等C语言源程序代码
计算机数值方法之最小二乘法拟合多项式C语言
sinat_41892454的博客
03-17 3113
给定数据点（xi ,yi），用最小二乘法拟合数据的多项式，并求平方误差。 有了之前多个程序的磨练，这次程序非常简单，用公式算出矩阵中的待求数值，再用高斯消元法求出a0和a1，写出拟合方程，再带入x,y值求出误差平方和即可。 代码 #include&lt;stdio.h&gt; #include&lt;math.h&gt; int i,j,k; float x[10000]; float...
用最小二乘法进行任意阶数的曲线拟合——C语言实现
weixin_41559101的博客
02-19 3899
/********************************************* 利用公式 ax=y–A’Ax=A’y 拟合曲线 * 用最小二乘法进行曲线拟合 * *********************************************/ #include &amp;amp;amp;amp;amp;amp;amp;amp;amp;amp;lt;stdio.h&amp;amp;amp;amp;amp;amp;amp;amp;am
C++实现多项式拟合求系数
YYSTINTERNET的博客
02-20 790
【代码】C++实现多项式拟合求系数。
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值