最小二乘法的编程实现

python 最小二乘法求解超定

python 最小二乘法求解超定最小二乘法是一种优化技术，用于求解超定方程组，也就是方程的数量大于未知数的数量的方程组。

在Python中，我们可以使用NumPy库中的linalg.lstsq函数来实现最小二乘法。

首先，我们需要理解最小二乘法的基本原理。

最小二乘法的基本思想是通过最小化误差的平方和来找到最佳函数匹配。

在超定方程组的情况下，我们无法找到一个精确的解，因为方程的数量超过了未知数的数量。

但是，我们可以找到一个最佳近似解，这个解能使得所有方程的残差平方和最小。

在Python中使用最小二乘法求解超定方程组的基本步骤如下：导入NumPy库。

定义超定方程组的系数矩阵A和目标向量b。

使用numpy.linalg.lstsq(A, b)函数求解超定方程组。

以下是一个示例代码：pythonimport numpy as np# 定义超定方程组的系数矩阵A和目标向量bA = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9], [10, 11, 12]])b = np.array([1, 2, 3, 4])# 使用numpy.linalg.lstsq函数求解超定方程组x, residuals, rank, s = np.linalg.lstsq(A, b, rcond=None)print("解向量x：", x)print("残差：", residuals)print("矩阵A的秩：", rank)print("奇异值：", s)注意，numpy.linalg.lstsq函数返回四个值：解向量x，残差，矩阵A的秩，以及A的奇异值。

其中，解向量x就是我们要求的近似解。

以上就是Python中使用最小二乘法求解超定方程组的方法。

最小二乘法C语言的实现

实验三.最小二乘法C语言的实现1.实验目的：进一步熟悉曲线拟合的最小二乘法。

掌握编程语言字符处理程序的设计和调试技术。

2.实验要求：输入：已知点的数目以及各点坐标。

输出：根据最小二乘法原理以及各点坐标求出拟合曲线。

3.程序流程：（1）输入已知点的个数；（2）分别输入已知点的X坐标；（3）分别输入已知点的Y坐标；（4）通过调用函数，求出拟合曲线。

最小二乘法原理如下：根据一组给定的实验数据，求出自变量x与因变量y的函数关系，只要求在给定点上的误差的平方和最小.当时，即(4.4.1)这里是线性无关的函数族，假定在上给出一组数据，以及对应的一组权，这里为权系数，要求使最小，其中(4.4.2)(4.4.2)中实际上是关于的多元函数，求I的最小值就是求多元函数I的极值，由极值必要条件，可得(4.4.3)根据内积定义引入相应带权内积记号(4.4.4)则(4.4.3)可改写为这是关于参数的线性方程组，用矩阵表示为(4.4.5)(4.4.5)称为法方程.当线性无关，且在点集上至多只有n个不同零点，则称在X上满足Haar条件，此时(4.4.5)的解存在唯一。

记(4.4.5)的解为从而得到最小二乘拟合曲线(4.4.6) 可以证明对，有故(4.4.6)得到的即为所求的最小二乘解.它的平方误差为(4.4.7) 均方误差为在最小二乘逼近中，若取，则，表示为(4.4.8)此时关于系数的法方程(4.4.5)是病态方程，通常当n≥3时都不直接取作为基。

程序流程图：↓↓程序：#include <math.h> #include <stdio.h> #include <stdlib.h> #include<malloc.h>float average(int n,float *x) {int i; float av; av=0;for(i=0;i<n;i++) av+=*(x+i);return(av);}//平方和float spfh(int n,float *x){int i;float a,b;a=0;for(i=0;i<n;i++)a+=(*(x+i))*(*(x+i));return(a);}//和平方float shpf(int n,float *x){int i;float a,b;a=0;for(i=0;i<n;i++)a=a+*(x+i);b=a*a/n;return(b);}//两数先相乘，再相加float dcj(int n,float *x,float *y) {int i;float a;a=0;for(i=0;i<n;i++)a+=(*(x+i))*(*(y+i));return(a);}//两数先相加，再相乘float djc(int n,float *x,float *y) {int i;float a=0,b=0;for(i=0;i<n;i++){a=a+*(x+i);b=b+*(y+i);}a=a*b/n;}//系数afloat xsa(int n,float *x,float *y){float a,b,c,d,e;a=spfh(n,x);b=shpf(n,x);c=dcj(n,x,y);d=djc(n,x,y);e=(c-d)/(a-b);//printf("%f %f %f %f",a,b,c,d);return(e);}float he(int n,float *y){int i;float a;a=0;for(i=0;i<n;i++)a=a+*(y+i);return(a);}float xsb(int n,float *x,float *y,float a){ float b,c,d;b=he(n,y);c=he(n,x);d=(b-a*c)/n;return(d);}void main(){ int n,i;float *x,*y,a,b;printf("请输入将要输入的有效数值组数n的值:"); scanf("%d",&n);x=(float*)calloc(n,sizeof(float));if(x==NULL){printf("内存分配失败");exit(1);}y=(float*)calloc(n,sizeof(float));if(y==NULL){printf("内存分配失败");exit(1);}printf("请输入x的值\n");for(i=0;i<n;i++)scanf("%f",x+i);printf("请输入y的值,请注意与x的值一一对应:\n"); for(i=0;i<n;i++)scanf("%f",y+i);for(i=0;i<n;i++){ printf("x[%d]=%3.2f ",i,*(x+i));printf("y[%d]=%3.2f\n",i,*(y+i));}a=xsa(n,x,y);b=xsb(n,x,y,a);printf("经最小二乘法拟合得到的一元线性方程为:\n"); printf("f(x)=%3.2fx+%3.2f\n",a,b);}。

最小二乘法C语言的实现

实验三.最小二乘法C语言的实现1.实验目的：进一步熟悉曲线拟合的最小二乘法。

掌握编程语言字符处理程序的设计和调试技术。

2.实验要求：输入：已知点的数目以及各点坐标。

输出：根据最小二乘法原理以及各点坐标求出拟合曲线。

3.程序流程：（1）输入已知点的个数；（2）分别输入已知点的X坐标；（3）分别输入已知点的Y坐标；（4）通过调用函数，求出拟合曲线。

最小二乘法原理如下：根据一组给定的实验数据，求出自变量x与因变量y的函数关系，只要求在给定点上的误差的平方和最小.当时，即(4.4.1)这里是线性无关的函数族，假定在上给出一组数据，以及对应的一组权，这里为权系数，要求使最小，其中(4.4.2)(4.4.2)中实际上是关于的多元函数，求I的最小值就是求多元函数I的极值，由极值必要条件，可得(4.4.3)根据内积定义引入相应带权内积记号(4.4.4)则(4.4.3)可改写为这是关于参数的线性方程组，用矩阵表示为(4.4.5) (4.4.5)称为法方程.当线性无关，且在点集上至多只有n个不同零点，则称在X上满足Haar条件，此时(4.4.5)的解存在唯一。

记(4.4.5)的解为从而得到最小二乘拟合曲线(4.4.6) 可以证明对，有故(4.4.6)得到的即为所求的最小二乘解.它的平方误差为(4.4.7) 均方误差为在最小二乘逼近中，若取，则，表示为(4.4.8)此时关于系数的法方程(4.4.5)是病态方程，通常当n≥3时都不直接取作为基。

程序流程图：开始↓输入已知点个数n↓输入已知点的X坐标↓输入已知点的Y坐标输出结果程序：#include <math.h>#include <stdio.h>#include <stdlib.h>#include<malloc.h>float average(int n,float *x){int i;float av;av=0;for(i=0;i<n;i++)av+=*(x+i);return(av);}//平方和float spfh(int n,float *x){int i;float a,b;a=0;for(i=0;i<n;i++)a+=(*(x+i))*(*(x+i));return(a);}//和平方float shpf(int n,float *x){int i;float a,b;a=0;for(i=0;i<n;i++)a=a+*(x+i);b=a*a/n;return(b);}//两数先相乘，再相加float dcj(int n,float *x,float *y) {int i;float a;a=0;for(i=0;i<n;i++)a+=(*(x+i))*(*(y+i));return(a);}//两数先相加，再相乘float djc(int n,float *x,float *y) {int i;float a=0,b=0;for(i=0;i<n;i++){a=a+*(x+i);b=b+*(y+i);}a=a*b/n;}//系数afloat xsa(int n,float *x,float *y){float a,b,c,d,e;a=spfh(n,x);b=shpf(n,x);c=dcj(n,x,y);d=djc(n,x,y);e=(c-d)/(a-b);//printf("%f %f %f %f",a,b,c,d);return(e);}float he(int n,float *y){int i;float a;a=0;for(i=0;i<n;i++)a=a+*(y+i);return(a);}float xsb(int n,float *x,float *y,float a){ float b,c,d;b=he(n,y);c=he(n,x);d=(b-a*c)/n;return(d);}void main(){ int n,i;float *x,*y,a,b;printf("请输入将要输入的有效数值组数n的值:"); scanf("%d",&n);x=(float*)calloc(n,sizeof(float));if(x==NULL){printf("内存分配失败");exit(1);}y=(float*)calloc(n,sizeof(float));if(y==NULL){printf("内存分配失败");exit(1);}printf("请输入x的值\n");for(i=0;i<n;i++)scanf("%f",x+i);printf("请输入y的值,请注意与x的值一一对应:\n"); for(i=0;i<n;i++)scanf("%f",y+i);for(i=0;i<n;i++){ printf("x[%d]=%3.2f ",i,*(x+i));printf("y[%d]=%3.2f\n",i,*(y+i));}a=xsa(n,x,y);b=xsb(n,x,y,a);printf("经最小二乘法拟合得到的一元线性方程为:\n"); printf("f(x)=%3.2fx+%3.2f\n",a,b);}。

(完整word版)最小二乘法拟合圆公式推导及matlab实现

2009-01-17 |最小二乘法(least squares analysis) 是一种数学优化技术，它通过最小化误差的平方和找到一组数据的最佳函数匹配。

最小二乘法是用最简的方法求得一些绝对不可知的真值，而令误差平方之和为最小。

小二乘法通常用于曲线拟合(least squares fitti ng) 。

这里有拟合圆曲线的公式推导过程和vc实现。

最小二乘法拟會圆曲线；= (x- +R2 = +- 2By4-B2令a=-2J4b = -2Bc = J^ +矿-0可得圆曲线方程的另一个册式Ix2 -\-y3十切十u = 0只要求出参数就可以求得圆心半径的参教；d)样本集（禺＜并e （123…N）中点到圆心的距离为a：打=（禺・4）2+（E傢点（耳乙）到圆边嫌的距离的平方与和半径平方的差为：@=£2_衣=（圣.4）2+（込.8）2_氏2=血2+込2+込+&乙+卍令Q（a,b,c）为Q的平方和：Q（aM = Z^2=工【（*/ + §2 + 込+b 齐+C）］2求参数a f b,c使得Q（a,g的值最小值。

解・PTT •平方差Qgg大于0,因此函数存在大于或等于0的极小值，极大值为无穷大.F（a,M）对a,吐求偏导，令偏导等于0,得到极值点，比较所有极值点的函数值即可得到最小值.绘仏"疋）=工2窗 +里+込+埒+c）Xjda —=0 迤(a,bQ =匸2阳+貯+込+坷+训=0範仏上疋）=工2（禺2+乙2+込 +空+° = 0 d解这个方程组。

(2)(3)(4)di（诵先消去c(2) W ⑷*工扎得：Ng 代'+Y-+aX +bY + c)X -工莎‘ +严 +aX +bY+c)x^X = 0 N^(X 2 +Y : +bY)X -^(X : +Y : +aX +bY)x^X =0("工禺2_工兀工兀)a + (“Y*占一工禺工齐仏(*+ + M 工*必2 -工牡丁 +去2)工禺=0(3) *N_⑷*工£得：N 工(X’ + y' + oZ +bY+c)Y-^(X 2 +Y- +aX +bY + c)x^Y =Q 吧(/+护 +aX +bY)Y +Y : +aX +dK)xVy =o （N'X 必一工禺工齐归+ （“丫呼一工§工齐）3 +“Y+N 工厅一 g af +严）三齐=o C =〔NgQ -gX 二X)D = (N 工尤F -工龙三卩)E-N^X 、+N^XY -工疔+丫‘)工XG = (NM 旷-三丫工丫)H =NW X'Y 七NT H -工 2’ +K-)YK可解得：|G? + Db + 5 = 0Da+Gb + H = 0HD-EG a = r CG-D 、v HC- ED o =D' _GC 工（疔+齐2）+幺工兀+c ―― ---------------------------------------------- N得A 、B 、R 的估计拟合值：R= - Ja‘ +2?' -牡 2(6)matlab 实现：function [R,A,B]=circ(x,y,N)x1 = 0;x2 = 0;x3 = 0;y1 = 0;y2 = 0;y3 = 0;x1y1 = 0;x1y2 = 0;x2y1 = 0;for i = 1 : Nx1 = x1 + x(i);x2 = x2 + x(i)*x(i);x3 = x3 + x(i)*x(i)*x(i);y1 = y1 + y(i);y2 = y2 + y(i)*y(i);y3 = y3 + y(i)*y(i)*y(i); x1y1 = x1y1 + x(i)*y(i); x1y2 = x1y2 +x(i)*y(i)*y(i); x2y1 = x2y1 + x(i)*x(i)*y(i); endC = N * x2 - x1 * x1;D = N * x1y1 - x1 * y1;E = N * x3 + N * x1y2 - (x2 + y2) * x1;G = N * y2 - y1 * y1;H = N * x2y1 + N * y3 - (x2 + y2) * y1;a = (H * D - E * G)/(C * G - D * D);b = (H * C - E * D)/(D * D - G * C);c = -(a * x1 + b * y1 + x2 + y2)/N;A = a/(-2); %x 坐标B = b/(-2); %y 坐标R = sqrt(a * a + b * b - 4 * c)/2;void CViewActionImageTool::LeastSquaresFitting(){if (m_nNum<3){ return; } int i=0;double X1=0;double Y1=0;double X2=0;double Y2=0;double X3=0;double Y3=0;double X1Y1=0;double X1Y2=0;double X2Y1=0;for (i=0;i<m_nNum;i++){X1 = X1 + m_points[i].x;Y1 = Y1 + m_points[i].y;X2 = X2 + m_points[i].x*m_points[i].x;Y2 = Y2 + m_points[i].y*m_points[i].y;X3 = X3 + m_points[i].x*m_points[i].x*m_points[i].x;Y3 = Y3 + m_points[i].y*m_points[i].y*m_points[i].y;X1Y1 = X1Y1 + m_points[i].x*m_points[i].y;X1Y2 = X1Y2 + m_points[i].x*m_points[i].y*m_points[i].y;X2Y1 = X2Y1 + m_points[i].x*m_points[i].x*m_points[i].y; } double C,D,E,G ,H,N;double a,b,c;N = m_nNum;C = N*X2 - X1*X1;D = N*X1Y1 - X1*Y1;E = N*X3 + N*X1Y2 - (X2+Y2)*X1;G = N*Y2 - Y1*Y1;H = N*X2Y1 + N*Y3 - (X2+Y2)*Y1;a = (H*D-E*G)/(C*G-D*D);b = (H*C-E*D)/(D*D-G*C);c = -(a*X1 + b*Y1 + X2 + Y2)/N;double A,B,R;A = a/(-2);B = b/(-2);R = sqrt(a*a+b*b-4*c)/2; m_fCenterX = A; m_fCenterY = B;m_fRadius = R; return;}。

matlab加权最小二乘法拟合编程

一、概述最小二乘法（Least Squares Method）是一种常用的数学优化方法，通过最小化残差的平方和来拟合实际数据与理论模型之间的关系。

在实际应用中，我们常常需要对数据进行加权处理，以提高拟合效果和准确度。

而Matlab作为一种强大的数学建模和仿真软件，提供了丰富的函数和工具来实现加权最小二乘法的拟合编程。

二、加权最小二乘法原理1. 最小二乘法原理最小二乘法是一种常用的拟合方法，通过最小化实际观测值和理论值之间的误差来寻找最佳拟合曲线或曲面。

其数学表达为：minimize ||Ax - b||^2其中A为设计矩阵，x为拟合参数，b为观测值向量。

最小二乘法可以看作是一种优化问题，通过求解参数x的最优值来实现最佳拟合。

2. 加权最小二乘法原理在实际情况下，我们往往会遇到观测值有不同的权重或方差的情况，此时可以使用加权最小二乘法来提高拟合效果。

加权最小二乘法的数学表达为：minimize ||W^(1/2)(Ax - b)||^2其中W为权重矩阵，将不同观测值的权重考虑在内，通过加权的方式来优化拟合效果。

三、Matlab实现加权最小二乘法1. 数据准备在进行加权最小二乘法的拟合编程前，首先需要准备实际观测数据和设计矩阵A。

还需要考虑观测值的权重矩阵W，根据实际情况来确定不同观测值的权重。

2. 加权最小二乘法函数Matlab提供了丰富的函数和工具来实现加权最小二乘法的拟合。

其中，可以使用lsqcurvefit或lsqnonlin等函数来进行加权最小二乘法的拟合计算。

通过传入设计矩阵A、观测值向量b和权重矩阵W，以及拟合参数的初始值，来实现加权最小二乘法的拟合计算。

3. 拟合结果评估完成加权最小二乘法的拟合计算后，我们需要对拟合结果进行评估。

主要包括残差分析、拟合效果的可视化等方面。

通过分析残差的分布和拟合曲线与实际观测值的符合程度，来评估拟合效果的优劣。

四、实例分析1. 示例一：线性模型拟合假设我们有一组线性关系的实际观测数据，且各观测值具有不同的权重。

最小二乘法C语言编程

最小二乘法c语言编程最小二乘法c语言程序c语言最小二乘法最小二乘法c语言代码c语言编程九九乘法表最小二乘法最小二乘法公式最小二乘法曲线拟合最小二乘法原理偏最小二乘法
#include <hidef.h> /* common defines and macros */
#include "derivative.h" /* derivative-specific definitions */
x2=x*x*x*x;
x3+=x2;
}
*(pp++)=x3;
}
}
pp=*(b+1);
p=a9;
ppp=a8;
x3=0;
for(i=0;i<10;i++) {
x=*(p++);
x2=*(ppp++);
x3+=x*x2;
}
*pp=x3;
pp=*(b+2);
p=a9;
ppp=a8;
x3=0;
for(i=0;i<10;i++) {
*(p++)=(x1*x5-x2*x4);
p=*(b1+0);
pp=*(b4+0);
for (i=0;i<3;i++) {
for(j=0;j<3;j++){
x1=*(p++);
x=x1*c;
*(pp++)=x;
}
}
//////////////////////////求两矩阵的商//////////////////////////////////////

c++最小二乘法二次多项式拟合

c++最小二乘法二次多项式拟合最小二乘法是一种常用的曲线拟合方法，可以用于拟合二次多项式。

以下是用C++实现最小二乘法二次多项式拟合的示例代码：#include <iostream>#include <vector>#include <cmath>using namespace std;// 定义二次多项式拟合的函数void quadraticLeastSquaresFit(const vector<double>& x, const vector<double>& y, double& a, double& b, double& c){int n = x.size();// 定义矩阵A和向量Bdouble A[3][3] = {0};double B[3] = {0};// 构造矩阵A和向量Bfor(int i=0; i<n; i++){double xi = x[i];double yi = y[i];A[0][0] += xi * xi;A[0][1] += xi;A[0][2] += 1.0;A[1][1] += xi;A[1][2] += 1.0;A[2][2] += 1.0;B[0] += xi * yi;B[1] += yi;B[2] += 1.0;}// 求解线性方程组Ax = Bdouble detA = A[0][0] * (A[1][1] * A[2][2] - A[2][1] * A[1][2])- A[0][1] * (A[1][0] * A[2][2] - A[2][0] * A[1][2])+ A[0][2] * (A[1][0] * A[2][1] - A[2][0] * A[1][1]);double invA[3][3] = {0};// 计算矩阵A的逆矩阵invA[0][0] = (A[1][1] * A[2][2] - A[1][2] * A[2][1]) / detA;invA[0][1] = -(A[0][1] * A[2][2] - A[0][2] * A[2][1]) / detA;invA[0][2] = (A[0][1] * A[1][2] - A[0][2] * A[1][1]) / detA;invA[1][0] = -(A[1][0] * A[2][2] - A[1][2] * A[2][0]) / detA;invA[1][1] = (A[0][0] * A[2][2] - A[0][2] * A[2][0]) / detA;invA[1][2] = -(A[0][0] * A[1][2] - A[0][2] * A[1][0]) / detA;invA[2][0] = (A[1][0] * A[2][1] - A[1][1] * A[2][0]) / detA;invA[2][1] = -(A[0][0] * A[2][1] - A[0][1] * A[2][0]) /detA;invA[2][2] = (A[0][0] * A[1][1] - A[0][1] * A[1][0]) / detA;// 计算拟合的系数a、b和ca = invA[0][0] * B[0] + invA[0][1] * B[1] + invA[0][2] * B[2];b = invA[1][0] * B[0] + invA[1][1] * B[1] + invA[1][2] * B[2];c = invA[2][0] * B[0] + invA[2][1] * B[1] + invA[2][2] * B[2];}int main(){// 定义数据点的x和y值vector<double> x = {1.0, 2.0, 3.0, 4.0, 5.0};vector<double> y = {1.8, 2.9, 3.9, 5.1, 6.2};// 定义拟合的系数a、b和cdouble a, b, c;// 进行二次多项式拟合quadraticLeastSquaresFit(x, y, a, b, c);// 输出拟合结果cout << "拟合的二次多项式为：y = " << a << "x^2 + " << b << "x + " << c << endl;return 0;}此代码将给定的数据点进行二次多项式拟合，通过最小二乘法求解拟合的系数a、b和c，并输出拟合结果。

最小二乘法的编程实现

1、最小二乘法：1）（用1TAA方法计算逆矩阵）#include <stdio.h>#include <math.h>#include <stdlib.h>#include<string.h>#include <iostream.h>#define N 200#define n 9void Getdata(double sun[N])//从txt文档中读取数据（小数）{char data;char sunpot[10]={0000000000};//为防止结果出现‘烫’字int i=0,j=0;double d;FILE *fp=fopen("新建文本文档.txt","r");if(!fp){printf("can't open file\n");}while(!feof(fp)){data=fgetc(fp);if(data!='\n'){sunpot[i]=data;i++;}else if(data=='\n'){sunpot[i]='\0';//给定结束符d=atof(sunpot);//将字符串转换成浮点数sun[j]=d;j++;i=0;//将i复位}}}void Normal(double sun[N],double sun1[N])//将数据进行标准化{double mean,temp=0,variance=0;int i;for(i=0;i<N;i++)temp=temp+sun[i];mean=temp/N;for(i=0;i<N;i++){sun1[i]=sun[i]-mean;}for(i=0;i<N;i++){variance=variance+sun1[i]*sun1[i];}variance=variance/N;variance=sqrt(variance);for(i=0;i<N;i++){sun1[i]=sun[i]/variance;//减去均值除以均方差进行归一化}}void Matrix(double sun[N],double matrixl[n][n],double matrixr[n])//组建方程的左右两个矩阵{double matrix1[N-n][n];//方程左边系数矩阵double matrix_trans[n][N-n];double matrix2[N-n];//方程右边系数矩阵int i,j,k;double temp;for(i=0;i<N-n;i++)//组建N-nxn维方程矩阵{for(j=0;j<n;j++){matrix1[i][j]=sun[n+i-j-1];//此处与matlab不同，matlab是从1开始的}}for(i=0;i<N-n;i++)matrix2[i]=sun[i+n];for(i=0;i<n;i++)//将N-nxn维矩阵转置{for(j=0;j<N-n;j++){matrix_trans[i][j]=matrix1[j][i];}}for(i=0;i<n;i++)//组建nxn维方阵(A'*A){for(k=0;k<n;k++){temp=0;for(j=0;j<N-n;j++){temp+=matrix_trans[i][j]*matrix1[j][k];}matrixl[i][k]=temp;}}for(i=0;i<n;i++){temp=0;for(j=0;j<N-n;j++){temp+=matrix_trans[i][j]*matrix2[j];}matrixr[i]=temp;}}void ExchangeLine(double a[n][n], int i, int j, int p)//交换矩阵第i行和第j行{int k;double temp;for (k=0;k<=p-1;k++){temp=a[j][k];a[j][k]=a[i][k];a[i][k]=temp;}}void Gauss(double a[n][n],int k,int p)//高斯消元,用第k行消去其他行第一个非零元素{int i,j;double m;for(i=k+1;i<p;i++){if(a[i][k]==0)i++;m = a[i][k] / a[k][k];for (j = k; j < p; j ++){a[i][j] = a[i][j] - m * a[k][j];}}}double Determinant(double a[n][n],int p)//求a的行列式{int k,i,j;double result=1.0;for(k=0;k<p-1;k++)//逐一取对角元素以下矩阵进行消元，n-1是因为最后一个对角元素无需消元{i=k;j=k;while(a[i][j]==0)i++;if((i<p)&&i>k){ExchangeLine(a,k,i,p);result=(-1)*result;//换行之后行列式变号一次}else if(i>=p)//即找不到可以交换的行{result=0;break;}Gauss(a,k,p);//对该列进行消元}for(i=0;i<p;i++){result=result*a[i][i];}return result;}void Adjoint_matrix(double a[n][n],double b[n][n],int p,int i,int j)//求a[i][j]元素的伴随矩阵{double temp[n][n]={{0,0,0},{0,0,0},{0,0,0}};//实际上只要n-1xn-1维，但为了调用方便，此处写为nxn 维int k,m,l,r;b[j][i]=1;for (k = 0, l = 0; k <p-1; k++, l ++)//将原矩阵去除i行j列后复制到temp中{if (l == i) l ++;for (m = 0, r = 0; m <p-1; m ++, r ++){if (r == j) r ++;temp[k][m] = a[l][r];}}b[j][i]=b[j][i]*Determinant(temp,n-1);if ((i + j) % 2 == 1)//(-1)^(i+j)b[j][i]=-b[j][i];}void Inverse(double a[n][n],double b,double c[n][n])//求逆最后一步运算{int i,j;for(i=0;i<n;i++){for(j=0;j<n;j++){c[i][j]=a[i][j]/b;}}}void main(){int p=n,i,j;double sunpot[N],sunpot_normal[N];double fi[n];double matrix[n][n];//要换成返回的方阵double matrixr[n];double adjoint[n][n],inverse[n][n];double determinant,temp;Getdata(sunpot);//从txt文档中读取每年的太阳黑子数据Normal(sunpot,sunpot_normal);//数据归一化Matrix(sunpot_normal,matrix,matrixr);//组建方程for(i=0;i<p;i++){for(j=0;j<p;j++){Adjoint_matrix(matrix,adjoint,p,i,j);}}determinant=Determinant(matrix,p);if (determinant== 0)printf("该矩阵不存在逆矩阵！\n");printf("行列式：\n");printf("%10lf\n",determinant);Inverse(adjoint,determinant,inverse);printf("逆矩阵：\n");for(i=0;i<p;i++){for(j=0;j<p;j++){printf("%12lf",inverse[i][j]);}printf("\n");}printf("fi：\n");for(i=0;i<n;i++){temp=0;for(j=0;j<n;j++){temp+=inverse[i][j]*matrixr[j];}fi[i]=temp;printf("%10lf",fi[i]);}printf("\n");}计算结果：fi：1.327258 -0.622683 0.018571 0.129074 -0.141278 0.103665 -0.046903 0.060403 0.154991下面用matlab进行拟合度计算：clearM=load('太阳黑子数.txt');sunpot=(M(:,2));N=200;n=9;[Y,mu,sigma]=zscore(sunpot);B(1:N-n)=Y(n+1:N);B=B';fi=[1.327258 -0.622683 0.018571 0.129074 -0.141278 0.103665 -0.046903 0.060403 0.154991]; %由C语言计算出的fi值fi=[1;-fi'];ts2=idpoly(fi',[]);B=Y(1:N-n);plot(B);compare(B, ts2, 1);2）用构造矩阵1A I I A-→的方法求逆：（列主元高斯消去法）[;][;]部分主要程序：void BuildMatrix(double a[n][M],double b[n][n]){int i,j;for (i=0;i<n;i++){for (j=0;j<n;j++){a[i][j]=b[i][j];}a[i][i+n]=1.0;}}void ExchangeLine(double a[n][M],int i,int j){int k;double temp;for (k=0;k<M;k++){temp=a[i][k];a[i][k]=a[j][k];a[j][k]=temp;}}void Judge(double a[n][M],int i){int j,m=0;//m=0来表示是否需要换行double temp=a[i][i];for(j=i+1;j<n;j++){if(fabs(a[j][i])>fabs(temp)){temp=a[j][i];m=j;}}if(m)ExchangeLine(a,i,m);}void Gauss(double a[n][M],int k) //高斯消元法,用第k行消去其他行{ { int i,j;double m;for(i=k+1;i<n;i++){if(a[i][k]==0)i++;m = a[i][k] / a[k][k];for (j = k; j < M; j ++){a[i][j] = a[i][j] - m * a[k][j];}}}void Normal(double a[n][M]){int i,j;double temp;for(i=0;i<n;i++){temp=a[i][i];if(temp!=1.0){for(j=i;j<M;j++){a[i][j]=a[i][j]/temp;}}}}void Inverse(double a[n][M]){int i,j,k;double temp;for(i=n-1;i>=0;i--){for(j=i-1;j>=0;j--){temp=a[j][i];for(k=i;k<M;k++){a[j][k]=a[j][k]-temp*a[i][k];}}}}void GetInverse(double a[n][M],double b[n][n]) {int i,j;for(i=0;i<n;i++){for(j=0;j<n;j++){b[i][j]=a[i][j+n];}}}void main(){int i,j;double temp;double sunpot[N];double sunpot_normal[N];double matrixr[n];double matrix1[n][n]={0};//方阵double matrix[n][M]={0};//构造矩阵double inverse[n][n]={0};double fi[n];Getdata(sunpot);//从txt文档中读取每年的太阳黑子数据Normal(sunpot,sunpot_normal);//数据归一化Matrix(sunpot_normal,matrix1,matrixr);//组建方程BuildMatrix(matrix,matrix1);for(i=0;i<n-1;i++){Judge(matrix,i);Gauss(matrix,i);}for(i=0;i<n;i++)//判断是否可逆{if(matrix[i][i]==0){printf("矩阵不可逆！\n");exit(-1);}}Normal(matrix);Inverse(matrix);GetInverse(matrix,inverse);printf("逆矩阵：\n");for(i=0;i<n;i++){for(j=0;j<n;j++){printf("%10lf",inverse[i][j]);}printf("\n");}printf("fi：\n");for(i=0;i<n;i++){temp=0;for(j=0;j<n;j++){temp+=inverse[i][j]*matrixr[j];}fi[i]=temp;printf("%10lf",fi[i]);}}。