课题八曲线拟合的最小二乘法实验目标：在某冶炼过程中，通过实验检测得到含碳量与时间关系的数据如下，试求含碳量y与时间t#include<iostream>#include<>using namespace std;int Array(double ***Arr, int n){double **p;int i;p=(double **)malloc(n*sizeof(double *));if(!p)return 0;for(i=0;i<n;i++){p[i]=(double *)malloc(n*sizeof(double));if(!p[i])return 0;}*Arr=p;return 1;}void main(){int n,i,j,k;double **A,*B,*S;S=new double[3];B=new double[3];cout<<"共2有®DN个节¨²点Ì：êo";cin>>n;cout<<"请输º入¨节¨²点Ì值¦Ì（ꡧXi）ê：êo"<<endl;double *X;X = new double[n];for(i=0;i<n;i++){cin>>X[i];}cout<<"请输º入¨节¨²点Ì函¡¥数ºy值¦Ì（ꡧYi）ê：êo"<<endl;double *Y;Y = new double[n];for(i=0;i<n;i++){cin>>Y[i];}if(!Array(&A,3))cout<<"内¨²存ä分¤配失º¡ì败㨹！ê"; else{for(i=0;i<3;i++){for(j=0;j<3;j++){A[i][j]=0;}}for(i=0;i<n;i++){A[0][0]+=1;// cout<<A[0][0];A[0][1]+=X[i];A[0][2]+=X[i]*X[i];A[1][0]=A[0][1];A[1][1]=A[0][2];A[1][2]+=X[i]*X[i]*X[i];A[2][0]=A[0][2];A[2][1]=A[1][2];A[2][2]+=X[i]*X[i]*X[i]*X[i];}}for(i=0;i<3;i++){B[i]=0;}for(i=0;i<n;i++){B[0]+=Y[i];B[1]+=X[i]*Y[i];B[2]+=X[i]*X[i]*Y[i];}for( k = 0; k < 2 ;k++){for(i = k+1; i<3; i++){double T = A[i][k]/A[k][k];B[i] = B[i] - T * B[k];for ( j = k+1 ; j < 3 ; j++ ){A[i][j] = A[i][j] - T * A[k][j];}}}S[2] = B[2]/A[2][2];for (i = 1; i >=0 ; i--){double Temp = 0;for (int j = i+1; j<2 ;j++)Temp = Temp + A[i][j] * S[j];S[i] = (B[i] - Temp) /A[i][i];}cout<<"拟a合曲¨²线为a：êo"<<endl;cout<<"y="<<S[0]<<" + "<<S[1]<<" x + "<<S[2]<<"x^2 "<<endl<<endl<<"误¨®差为a：êo"<<endl;for(i=0;i<n;i++){cout<<"Y"<<i<<"-Y(t"<<i<<")="<<S[0]+S[i]*X[i]+S[2]*X[i ]*X[i]-Y[i]<<endl;}}。
最小二乘法曲线拟合C语言实现简单思路如下：1，采用目标函数对多项式系数求偏导，得到最优值条件，组成一个方程组；2，方程组的解法采用行列式变换（两次变换：普通行列式——三角行列式——对角行列式——求解），行列式的求解算法上优化过一次了，目前还没有更好的思路再优化运算方法，限幅和精度准备再修改修改目前存在的问题：1，代码还是太粗糙2，数学原理可行，但是计算机运算有幅度溢出和精度问题，这方面欠考虑，导致高阶大数据可能拟合失败下面附简单注释版代码：#include "stdafx.h"#include "stdlib.h"#include "math.h"//把二维的数组与一维数组的转换，也可以直接用二维数组，只是我的习惯是不用二维数组#define ParaBuffer(Buffer,Row,Col) (*(Buffer + (Row) * (SizeSrc + 1) + (Col)))///************************************************************** *********************从txt文件里读取double型的X，Y数据txt文件里的存储格式为X1 Y1X2 Y2X3 Y3X4 Y4X5 Y5X6 Y6X7 Y7X8 Y8函数返回X，Y，以及数据的数目（以组为单位）*************************************************************** ********************/static int GetXY(const char* FileName, double* X, double* Y, int* Amount){FILE* File = fopen(FileName, "r");if (!File)return -1;for (*Amount = 0; !feof(File); X++, Y++, (*Amount)++)if (2 != fscanf(File, (const char*)"%lf %lf", X, Y))break;fclose(File);return 0;}/************************************************************** *********************打印系数矩阵，只用于调试，不具备运算功能对于一个N阶拟合，它的系数矩阵大小是（N + 1）行（N + 2）列double* Para：系数矩阵存储地址int SizeSrc：系数矩阵大小（SizeSrc）行（SizeSrc + 1）列***********************************************************************************/static int PrintPara(double* Para, int SizeSrc){int i, j;for (i = 0; i < SizeSrc; i++){for (j = 0; j <= SizeSrc; j++)printf("%10.6lf ", ParaBuffer(Para, i, j));printf("\r\n");}printf("\r\n");return 0;}/************************************************************** *********************系数矩阵的限幅处理，防止它溢出，目前这个函数很不完善，并不能很好地解决这个问题原理：矩阵解行列式，同一行乘以一个系数，行列式的解不变当然，相对溢出问题，还有一个精度问题，也是同样的思路，现在对于这两块的处理很不完善，有待优化以行为单位处理*************************************************************** ********************/static int ParalimitRow(double* Para, int SizeSrc, int Row){int i;double Max, Min, Temp;for (Max = abs(ParaBuffer(Para, Row, 0)), Min = Max, i = SizeSrc; i; i--){Temp = abs(ParaBuffer(Para, Row, i));if (Max < Temp)Max = Temp;if (Min > Temp)Min = Temp;}Max = (Max + Min) * 0.000005;for (i = SizeSrc; i >= 0; i--)ParaBuffer(Para, Row, i) /= Max;return 0;}/************************************************************** *********************同上，以矩阵为单位处理*************************************************************** ********************/static int Paralimit(double* Para, int SizeSrc){int i;for (i = 0; i < SizeSrc; i++)if (ParalimitRow(Para, SizeSrc, i))return -1;return 0;}/************************************************************** *********************系数矩阵行列式变换*************************************************************** ********************/static int ParaPreDealA(double* Para, int SizeSrc, int Size){int i, j;for (Size -= 1, i = 0; i < Size; i++){for (j = 0; j < Size; j++)ParaBuffer(Para, i, j) = ParaBuffer(Para, i, j) * ParaBuffer(Para, Size, Size) - ParaBuffer(Para, Size, j) * ParaBuffer(Para, i, Size);ParaBuffer(Para, i, SizeSrc) = ParaBuffer(Para, i, SizeSrc) * ParaBuffer(Para, Size, Size) - ParaBuffer(Para, Size, SizeSrc) * ParaBuffer(Para, i, Size);ParaBuffer(Para, i, Size) = 0;ParalimitRow(Para, SizeSrc, i);}return 0;}/************************************************************** *********************系数矩阵行列式变换，与ParaPreDealA配合完成第一次变换，变换成三角矩阵*************************************************************** ********************/static int ParaDealA(double* Para, int SizeSrc){int i;for (i = SizeSrc; i; i--)if (ParaPreDealA(Para, SizeSrc, i))return -1;return 0;}/************************************************************** *********************系数矩阵行列式变换*************************************************************** ********************/static int ParaPreDealB(double* Para, int SizeSrc, int OffSet) {int i, j;for (i = OffSet + 1; i < SizeSrc; i++){for (j = OffSet + 1; j <= i; j++)ParaBuffer(Para, i, j) *= ParaBuffer(Para, OffSet, OffSet);ParaBuffer(Para, i, SizeSrc) = ParaBuffer(Para, i, SizeSrc) * ParaBuffer(Para, OffSet, OffSet) - ParaBuffer(Para, i, OffSet) * ParaBuffer(Para, OffSet, SizeSrc);ParaBuffer(Para, i, OffSet) = 0;ParalimitRow(Para, SizeSrc, i);}return 0;}/************************************************************** *********************系数矩阵行列式变换，与ParaPreDealB配合完成第一次变换，变换成对角矩阵，变换完毕***********************************************************************************/static int ParaDealB(double* Para, int SizeSrc){int i;for (i = 0; i < SizeSrc; i++)if (ParaPreDealB(Para, SizeSrc, i))return -1;for (i = 0; i < SizeSrc; i++){if (ParaBuffer(Para, i, i)){ParaBuffer(Para, i, SizeSrc) /= ParaBuffer(Para, i, i);ParaBuffer(Para, i, i) = 1.0;}}return 0;}/************************************************************** *********************系数矩阵变换*************************************************************** ********************/static int ParaDeal(double* Para, int SizeSrc){PrintPara(Para, SizeSrc);Paralimit(Para, SizeSrc);PrintPara(Para, SizeSrc);if (ParaDealA(Para, SizeSrc))return -1;PrintPara(Para, SizeSrc);if (ParaDealB(Para, SizeSrc))return -1;return 0;}/************************************************************** *********************最小二乘法的第一步就是从XY数据里面获取系数矩阵double* Para：系数矩阵地址const double* X：X数据地址const double* Y：Y数据地址int Amount：XY数据组数int SizeSrc：系数矩阵大小（SizeSrc）行（SizeSrc + 1）列*************************************************************** ********************/static int GetParaBuffer(double* Para, const double* X, const double* Y, int Amount, int SizeSrc){int i, j;for (i = 0; i < SizeSrc; i++)for (ParaBuffer(Para, 0, i) = 0, j = 0; j < Amount; j++)ParaBuffer(Para, 0, i) += pow(*(X + j), 2 * (SizeSrc - 1) - i);for (i = 1; i < SizeSrc; i++)for (ParaBuffer(Para, i, SizeSrc - 1) = 0, j = 0; j < Amount; j++) ParaBuffer(Para, i, SizeSrc - 1) += pow(*(X + j), SizeSrc - 1 - i);for (i = 0; i < SizeSrc; i++)for (ParaBuffer(Para, i, SizeSrc) = 0, j = 0; j < Amount; j++)ParaBuffer(Para, i, SizeSrc) += (*(Y + j)) * pow(*(X + j), SizeSrc - 1 - i);for (i = 1; i < SizeSrc; i++)for (j = 0; j < SizeSrc - 1; j++)ParaBuffer(Para, i, j) = ParaBuffer(Para, i - 1, j + 1);return 0;}/************************************************************** *********************整个计算过程*************************************************************** ********************/int Cal(const double* BufferX, const double* BufferY, int Amount, int SizeSrc, double* ParaResK){double* ParaK = (double*)malloc(SizeSrc * (SizeSrc + 1) * sizeof(double));GetParaBuffer(ParaK, BufferX, BufferY, Amount, SizeSrc);ParaDeal(ParaK, SizeSrc);for (Amount = 0; Amount < SizeSrc; Amount++, ParaResK++) *ParaResK = ParaBuffer(ParaK, Amount, SizeSrc);free(ParaK);return 0;}/************************************************************** *********************测试main函数数据组数：20阶数：5***********************************************************************************/int main(int argc, char* argv[]){//数据组数int Amount;//XY缓存，系数矩阵缓存double BufferX[1024], BufferY[1024], ParaK[6];if (GetXY((const char*)"test.txt", (double*)BufferX, (double*)BufferY, &Amount))return 0;//运算Cal((const double*)BufferX, (const double*)BufferY, Amount, sizeof(ParaK) / sizeof(double), (double*)ParaK);//输出系数for (Amount = 0; Amount < sizeof(ParaK) / sizeof(double); Amount++)printf("ParaK[%d] = %lf!\r\n", Amount, ParaK[Amount]);//屏幕暂留system("pause");return 0;}。
c语言最小二乘法拟合曲线C语言中，可以使用最小二乘法来拟合曲线。
最小二乘法是一种常用的数学优化方法，用于找到一条曲线，使得曲线和实际数据之间的误差最小。
下面是一个简单的示例代码，使用最小二乘法来拟合一条直线的曲线。
c#include <stdio.h>// 最小二乘法拟合直线void leastSquareFit(int n, double x[], double y[], double* slope, double* intercept) {// 计算 x 和 y 的平均值double sumX = 0, sumY = 0;for (int i = 0; i < n; i++) {sumX += x[i];sumY += y[i];}double meanX = sumX / n;double meanY = sumY / n;// 计算直线的斜率double numerator = 0, denominator = 0;for (int i = 0; i < n; i++) {numerator += (x[i] - meanX) * (y[i] - meanY);denominator += (x[i] - meanX) * (x[i] - meanX);}*slope = numerator / denominator;// 计算直线的截距*intercept = meanY - (*slope) * meanX;}int main() {// 原始数据double x[5] = {1, 2, 3, 4, 5};double y[5] = {2, 4, 6, 8, 10};// 拟合结果double slope, intercept;leastSquareFit(5, x, y, &slope, &intercept);printf("拟合直线方程：y = %.2fx + %.2f\n", slope, intercept);return 0;}运行以上代码，将得到拟合直线方程为：`y = 2.00x + 0.00`。
实验名称：曲线拟合的最小二乘法 实验目的了解曲线拟合的最小二乘法实验类型设计型实验环境Windows XP TC实验内容相关知识:已知C ［a,b ］中函数f （x ）的一组实验数据（x i ，y i ）（i=0，1，…，m),其中y i =f （x i ）。
设);,,1,0)((m n n j x j <= ϕ是C ［a ，b]上线性无关函数族。
在)}(,),(),({10x x x span n ϕϕϕφ =中找函数f(x) 曲线拟合的最小二乘解∑==nj j j x a x S 0*)()(ϕ，其法方程(组）为：),,1,0(),(0n k d a nj k j j k==∑=ϕϕ其中,∑==mi i k i j i k j x x x 0)()()(),(ϕϕωϕϕ∑=≡=mi k i k i i k d x x f x f 0)()()(),(ϕωϕ k=0，1，…，n特别是，求函数f(x ） 曲线拟合的线性最小二乘解b ax x S +=)(*的计算公式为:∑∑∑∑∑∑======-+-=mi i m i i mi i i mi i mi i mi i x x m y x x y x b 02202)()1())(())((∑∑∑∑∑=====-+-+=mi i m i i mi i m i i m i i i x x m y x y x m a 0220)()1())(()1(数据结构：两个一维数组或一个二维数组 算法设计:（略）实验用例： 已知函数y=f （x)的一张表：处的近似值，并画出实验数据和直线。
编写代码:＃include〈stdio.h>＃include<stdlib.h>＃include<graphics.h〉double qiuhe1(double a［10]［2]，int p){int i；double y;y=0;for（i=0;i〈10;i++)y=y+a[i］[p］;return y;｝double qiuhe2(double a［10］[2]，int p){int i；double y=0;for(i=0；i<10；i++）y=y+a［i][0］＊a[i]［p］；return y;｝double nihe（double a［10][2］，double x）{double a1，b，y；a1=（10*qiuhe2(a，1）-qiuhe1(a,0)*qiuhe1(a，1)）/（10*qiuhe2（a,0)-qiuhe1(a，0）*qiuhe1（a,0))；b=（qiuhe2(a,0）＊qiuhe1(a,1)—qiuhe1(a,0）*qiuhe2(a，1）)/（10*qiuhe2(a,0)—qiuhe1(a，0)*qiuhe1(a,0)）;y=a1＊x+b;return y；｝int main（）{double a[10］[2］={0,68，10，67.1,20,66.4,30，65.6,40,64.6，50，61.8， 60，61.0，70,60.8,80，60.4,90，60｝；double x,x1,q=1;c har c［12];i nt i;long n;int arw[6］=｛515，235，520，240,515，245｝；int arw1［6］=｛315,45,320，40，325,45｝;int gdriver=IBM8514;int gmode=IBM8514HI；initgraph（&gdriver， &gmode, ”c：\\TC20\\BGI");cleardevice(）；printf(”input x:\n"）;scanf("%lf"，＆x）;printf（"%f\n",nihe(a,x））；n=nihe（a，x）＊1000000+1;c[0］='y’;c[1]=’=';c［4］='。
