下载文档原格式
NET与Matlab结合——最⼩⼆乘法直线拟合（C#）NET与Matlab结合 —— 最⼩⼆乘法直线拟合（C#）⾸先是⼀个.m⽂件drawgraph.m，确保它能够在Matlab⾥运⾏。
我这⾥是最⼩⼆乘法直线拟合程序。
%最⼩⼆乘法直线拟合%Created by Safirst C. Ke 2007.8.29 Wed 14:51function drawgraph(coords)%传⼊的参数为两⾏向量，第⼀⾏为x坐标，第⼆⾏为坐标。
%axis ([0 100 0 100]);grid on;hold on;%显⽰欲拟合的点的位置plot(coords(1,:), coords(2,:), '*');%分解x，y坐标x = coords(1,:)y = coords(2,:)'b = size(coords);c = ones(1, b(2));MT = [c; x];M = MT';%f为直线函数，f = mx + b;f = inv(MT * M) * MT * y['y = ', num2str(f(2)), 'x + ', num2str(f(1))]%显⽰最终拟合的直线x = -max(x):max(x);y = f(1) + f(2) * x;plot(x, y);xlabel('X轴');ylabel('Y轴');title('最⼩⼆乘法直线拟合 by Safirst C. Ke');legend(['y = ', num2str(f(2)), 'x + ', num2str(f(1))]);然后将这个⽂件包含在.NET的类库⼯程中，并进⾏编译。
这⾥需要理解它的过程，毕竟.NET不能编译.m⽂件。
怎么做到的呢？通过设置这个⼯程的⽣成事件属性，添加为call PlotDemoBuild.bat然后在PlotDemoBuild.bat这个⽂件⾥⾯写好⽤Matlab编译器mcc编译的命令⾏，最重要的部分就是mcc -M -silentsetup -vg -B "dotnet:PlotDemoComp,Plotter,2.0,private" -d ../../src ../../drawgraph.m这样的话，点击⽣成，就会通过mcc产⽣dll，即我们需要的类库。
最小二乘法曲线拟合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语言实现最小二乘法，并且给出一个简单的示例程序。
首先，我们需要明确最小二乘法的原理。
其核心思想是寻找一个最佳的线性函数，使得该函数与给定的数据点之间的误差最小。
具体来说，我们可以定义误差函数为：E = Σ(yi - a - bx_i)^2其中，a和b是待求的系数，xi和yi代表第i个数据点的横纵坐标。
我们的目标是寻找a和b的取值，使得E最小。
为了求解这个问题，我们可以采用求导的方法来求得a和b的解析式。
具体来说，我们可以对E关于a和b分别求导，并令其为0，得到如下方程组：Σyi = na + bΣxiΣxiyi = aΣxi + bΣ(x_i)^2其中，n为数据点的数量。
通过解这个方程组，我们可以得到a 和b的解析式：a = (ΣyiΣ(x_i)^2 - ΣxiyiΣxi) / (nΣ(x_i)^2 - (Σxi)^2)b = (nΣxiyi - ΣxiΣyi) / (nΣ(x_i)^2 - (Σxi)^2)有了这个解析式，我们就可以使用C语言来实现最小二乘法。
下面是一个简单的示例程序：#include <stdio.h>int main(){int n = 5; // 数据点的数量double xi[] = {1, 2, 3, 4, 5}; // 横坐标double yi[] = {2.1, 3.9, 6.1, 8.2, 10.2}; // 纵坐标double sum_xi = 0, sum_yi = 0, sum_xi2 = 0, sum_xiyi = 0; for (int i = 0; i < n; i++) {sum_xi += xi[i];sum_yi += yi[i];sum_xi2 += xi[i] * xi[i];sum_xiyi += xi[i] * yi[i];}double a = (sum_yi * sum_xi2 - sum_xi * sum_xiyi) / (n * sum_xi2 - sum_xi * sum_xi);double b = (n * sum_xiyi - sum_xi * sum_yi) / (n * sum_xi2 - sum_xi * sum_xi);printf('a = %.2f, b = %.2f', a, b);return 0;}在这个程序中，我们假设有5个数据点，分别对应横坐标1到5，纵坐标2.1到10.2。
最小二乘法 c语言最小二乘法是一种常用的数学方法，用于通过已知数据点拟合出一条最佳拟合曲线。
在本文中，我们将讨论如何使用C语言实现最小二乘法。
我们需要明确最小二乘法的基本原理。
最小二乘法的目标是找到一条曲线，使得该曲线上的点到已知数据点的距离之和最小。
具体地，我们假设已知数据点的集合为{(x1, y1), (x2, y2), ..., (xn, yn)}，我们需要找到一条曲线y = f(x)，使得f(xi)与yi的差的平方和最小。
那么，如何在C语言中实现最小二乘法呢？首先，我们需要定义一个函数来计算拟合曲线上的点f(xi)。
在这个函数中，我们可以使用多项式的形式来表示拟合曲线。
例如，我们可以选择使用一次多项式y = ax + b来拟合数据。
然后，我们可以使用最小二乘法的公式来计算出最优的a和b的值。
接下来，我们需要编写一个函数来计算拟合曲线上每个点f(xi)与已知数据点yi的差的平方和。
通过遍历已知数据点的集合，并计算每个点的差的平方，最后将所有差的平方求和，即可得到拟合曲线的误差。
然后，我们可以使用梯度下降法来最小化误差函数。
梯度下降法是一种优化算法，通过不断迭代来找到误差函数的最小值。
在每次迭代中，我们通过计算误差函数对参数a和b的偏导数，来更新a和b的值。
通过多次迭代，最终可以找到最优的a和b的值，从而得到最佳拟合曲线。
我们可以编写一个主函数来调用以上的函数，并将最终的拟合曲线绘制出来。
在这个主函数中，我们可以读取已知数据点的集合，并调用最小二乘法函数来计算拟合曲线的参数。
然后，我们可以使用绘图库来绘制已知数据点和拟合曲线，并将结果输出到屏幕上。
通过以上的步骤，我们就可以使用C语言实现最小二乘法了。
当然，在实际应用中，我们可能会遇到更复杂的数据和更高阶的多项式拟合。
但是基本的原理和方法是相同的，只是需要做一些适当的调整。
总结一下，最小二乘法是一种常用的数学方法，用于通过已知数据点拟合出一条最佳拟合曲线。
c语言最小二乘法
C语言中的最小二乘法是一种常用的统计学工具，用于拟合一组数据的最佳直线。
它的原理是通过对这组数据进行线性回归，找到最能代表这组数据的一条直线，从而求出直线的斜率和截距，并计算出最小二乘法的误差平方和。
在C语言中，可以通过利用数组和循环语句来实现最小二乘法的计算，具有简单易学、方便快捷等优点。
一般来说，使用最小二乘法能够更准确地预测未知数据的趋势，并为数据的分析和应用提供依据。
- 1 -。
c语言最小二乘法拟合直线最小二乘法是一种拟合直线或曲线的方法，适用于给定一组实验数据的情况。
在C语言中，可以使用以下步骤实现最小二乘法拟合直线：1. 定义并初始化数据变量。
例如：double x[10] = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10}; double y[10] = {2.5,3.7,4.1, 4.8,5.5,6.3,7.2,8.4, 8.5,9.1};2. 求出数据变量的平均值。
例如：double sumx = 0, sumy = 0, avgx, avgy; for(int i=0; i<10; i++) { sumx += x[i]; sumy +=y[i]; } avgx = sumx / 10; avgy = sumy / 10;3. 求出数据变量的协方差和方差。
例如：double s_xx = 0,s_xy = 0, s_yy = 0; for(int i=0; i<10; i++) { s_xx += (x[i] - avgx) * (x[i] - avgx); s_xy += (x[i] - avgx) * (y[i] - avgy); s_yy += (y[i] - avgy) * (y[i] - avgy); }4. 求出拟合直线的斜率和截距。
例如：double slope = s_xy / s_xx; double intercept = avgy - slope * avgx;5. 输出结果。
例如：printf("拟合直线方程为 y = %.2fx+ %.2f", slope, intercept);通过以上步骤，就可以使用最小二乘法拟合给定数据的直线。
C++最小二乘法拟合直线有一堆的X, Y值需要拟合直线，如下X: x1, x2, ...... x nY: y1, y2, …… y n如何得到拟合成y = kx + b的等式呢？用最小二乘法可以实现。
经过计算的k和b的公式如下k=，b =-kC++ 代码实现如下/*最小二乘法y = kx + bX，Y数据分别存在数组中*/#include <iostream>using namespace std;typedef struct{double k; //拟合后的直线斜率double b; //拟合后的直线截距}Linekb;//传入X, Y数组，返回斜率k,截距bLinekb CalcLine(double srcX[], double srcY[], int nX){Linekb LKB;double sumX = 0, sumY = 0, s_xy = 0, s_xx = 0;for(int i = 0; i < nX; i++){sumX += srcX[i];sumY += srcY[i];s_xy += srcX[i] * srcY[i];s_xx += srcX[i] * srcX[i];}double _x = sumX / nX;double _y = sumY / nX;LKB.k = (s_xy - nX * _x * _y) / (s_xx - nX * _x * _x);LKB.b = _y - LKB.k * _x;return LKB;}int main(){//y = 3x + 5double nXCrd[10] = {1.0, 2.0, 3.3, 4.0, 5.4, 6.0, 7.0, 8.0, 9.0, 10.0};double nYValue[10] = {8.0, 11.0, 14.0, 17.3, 20.0, 23.7, 26.0, 29.0, 32.0, 35.0};Linekb lkb = CalcLine(nXCrd, nYValue, 10);cout << "该直线是：" << "y = " << lkb.k << "X + " << lkb.b << endl;cout << endl;system("pause");return 0;}运行结果如下与原直线y = 3x+5 相比误差很小。
最小二乘法直线拟合的C++编程//面向过程#include<iostream>#include<cmath>using namespace std;void linear_equation(double x[],double y[],int n,double&a,double&b) {int i;double sum_x=0, sum_y=0, sum_xy=0, sum_x2=0; for(i=0;i<n;i++)sum_x+=x[i];for(i=0;i<n;i++)sum_y+=y[i];for(i=0;i<n;i++)sum_xy+=x[i]*y[i];for(i=0;i<n;i++)sum_x2+=pow(x[i],2);a=(sum_xy*sum_x-sum_y*sum_x2)/(pow(sum_x,2)-n*sum_x2);b=(sum_x*sum_y-n*sum_xy)/(pow(sum_x,2)-n*sum_x2); }void main(){char independent_variable,dependent_variable;int i,n;double a,b;double x[100],y[100];cout<<"Please input dependent variable:\n"; cin>>dependent_variable;cout<<"Please input independent variable:\n";cin>>independent_variable;cout<<"Please input the number N:\n"; cin>>n;cout<<"Please input the number of independentvariable"<<"("<<independent_variable<<")"<<endl; for(i=0;i<n;i++) cin>>x[i];cout<<"Please input the number of dependentvariable"<<"("<<dependent_variable<<")"<<endl;for(i=0;i<n;i++)cin>>y[i];linear_equation(x,y,n,a,b);cout<<"a= "<<a<<endl;cout<<"b= "<<b<<endl;if(a==0){if(b==0)cout<<dependent_variable<<"=0"<<endl; else if(b==1)cout<<dependent_variable<<"="<<independent_variable<<endl;}else{if(b==1)cout<<dependent_variable<<"="<<a<<"+"<<independent_variable<<endl; elsecout<<dependent_variable<<"="<<a<<"+"<<b<<independent_variable<<endl ; }}//面向对象#include<iostream>#include<cmath>using namespace std;class linear_equation{public:void getnumber();void equation();void disp();private:char independent_variable,dependent_variable;int i,n;double a,b;double x[100],y[100];};void linear_equation::getnumber(){cout<<"Please input dependent variable:\n"; cin>>dependent_variable;cout<<"Please input independent variable:\n";cin>>independent_variable;cout<<"Please input the number N:\n"; cin>>n;cout<<"Please input the number of independentvariable"<<"("<<independent_variable<<")"<<endl; for(i=0;i<n;i++) cin>>x[i];cout<<"Please input the number of dependentvariable"<<"("<<dependent_variable<<")"<<endl;for(i=0;i<n;i++)cin>>y[i];}void linear_equation::equation(){double sum_x=0, sum_y=0, sum_xy=0, sum_x2=0; for(i=0;i<n;i++) sum_x+=x[i];for(i=0;i<n;i++)sum_y+=y[i];for(i=0;i<n;i++)sum_xy+=x[i]*y[i];for(i=0;i<n;i++)sum_x2+=pow(x[i],2);a=(sum_xy*sum_x-sum_y*sum_x2)/(pow(sum_x,2)-n*sum_x2);b=(sum_x*sum_y-n*sum_xy)/(pow(sum_x,2)-n*sum_x2); }void linear_equation::disp(){cout<<"a= "<<a<<endl;cout<<"b= "<<b<<endl;if(a==0){if(b==0)cout<<dependent_variable<<"=0"<<endl; else if(b==1)cout<<dependent_variable<<"="<<independent_variable<<endl;}else{if(b==1)cout<<dependent_variable<<"="<<a<<"+"<<independent_variable<<endl; elsecout<<dependent_variable<<"="<<a<<"+"<<b<<independent_variable<<endl; }}void main(){linear_equation s1;s1.getnumber();s1.equation();s1.disp();}。
c语言最小二乘法最小二乘法是一种经典的线性回归分析方法，广泛应用于数据拟合和趋势线的绘制。
在计算机科学领域中，使用C语言代码来实现最小二乘法，可以快速高效地处理大量数据。
本文将介绍在C语言中实现最小二乘法的步骤和基本原理。
1. 定义输入和输出参数在C语言中，我们需要先定义要输入和输出的参数。
输入参数包括X数组和Y数组，分别表示自变量和因变量的值。
为了便于计算，我们需要定义X的平方和、X乘Y的和、X的和和Y的和。
double X[100], Y[100];double SumX2, SumXY, SumX, SumY;int Count;2. 输入数据利用scanf函数输入数据，可以通过循环操作，一次性输入整个数组。
如果您觉得输入过程繁琐，可以使用文件读取函数读取数据。
printf("请输入X、Y、数据个数: ");scanf("%lf%lf%d", &X[i], &Y[i], &Count);3. 计算和最小二乘法的核心是计算和，利用输入的数据计算出SumX2、SumXY、SumX、SumY四个和数。
这里可以使用循环操作，遍历整个数组并计算值。
for (i = 0; i < Count; i++) {SumX += X[i];SumY += Y[i];SumX2 += (X[i] * X[i]);SumXY += (X[i] * Y[i]);}4. 计算斜率和截距利用计算出来的和，我们可以计算出最小二乘法的斜率和截距。
否则，在调用最小二乘法算法时将无法获得正确的结果。
double Slope = (Count * SumXY - SumX * SumY) / (Count * SumX2 - SumX * SumX);double Intercept = (SumY - Slope * SumX) / Count;5. 输出结果计算出来斜率和截距，就可以使用printf函数输出结果。
/*最小二乘法的曲线拟合*/#include<stdio.h>#include<math.h>#include<stdlib.h>#define max100void main(){int i,j,k,m,N,mi;float mx,temp;float X[max][max],Y[max],x[max],y[max],a[max];FILE*fp1;if((fp1=fopen("in1.txt","r"))==NULL)/*输入拟合曲线的次数m以及已知的数据组数N*/{printf("Can't open this file!\n");exit(0);}for(i=0;i<2;i++)fscanf(fp1,"%d%d",&m,&N);fclose(fp1);FILE*fp2;if((fp2=fopen("in2.txt","r"))==NULL)/*输入已知的N组数据坐标*/{printf("Can't open this file!\n");exit(0);}for(i=0;i<N;i++)fscanf(fp2,"%f%f",&x[i],&y[i]);fclose(fp2);for(i=0;i<=m;i++)/*由给定的点得系数矩阵*/for(j=i;j<=m;j++){temp=0;for(k=0;k<N;k++)temp=temp+pow(x[k],(i+j));X[i][j]=temp;X[j][i]=X[i][j];}for(i=0;i<=m;i++)/*由给定的点得右端矩阵列向量*/{temp=0;for(k=0;k<N;k++)temp=temp+y[k]*pow(x[k],i);Y[i]=temp;}for(j=0;j<m;j++)/*列主元素消去法解该线性方程组，得拟合曲线多项式的系数*/{for(i=j+1,mi=j,mx=fabs(X[j][j]);i<=m;i++)if(fabs(X[i][j])>mx){mi=i;mx=fabs(X[i][j]);}if(j<mi){temp=Y[j];Y[j]=Y[mi];Y[mi]=temp;for(k=j;k<=m;k++){temp=X[j][k];X[j][k]=X[mi][k];X[mi][k]=temp;}}for(i=j+1;i<=m;i++){temp=-X[i][j]/X[j][j];Y[i]+=Y[j]*temp;for(k=j;k<=m;k++)X[i][k]+=X[j][k]*temp;}}a[m]=Y[m]/X[m][m];for(i=m-1;i>=0;i--){a[i]=Y[i];for(j=i+1;j<=m;j++)a[i]-=X[i][j]*a[j];a[i]/=X[i][i];}FILE*fp3;/*输出拟合所得的曲线方程*/fp3=fopen("out.txt","w");fprintf(fp3,"P(x)=(%f)",a[0]); for(i=1;i<=m;i++)fprintf(fp3,"+(%f)*x^%d",a[i],i); fclose(fp3);}。
