求积分的两点高斯公式

高斯求积公式是变步长数值积分的一种，基本形式是计算[-1,1]上的定积分。下面简单说明一下思想（仅仅是说明，而非证明）： 假设现在要求f(x)在[-1,1]上的积分值，只允许计算一次f(x)的值，你会怎么做呢？

显然我们会选取一点x0,计算出f(x0),然后用A=f(x0)2作为近似值。现在问题是怎样选取x0，使得结果尽可能精确呢？

直觉告诉我们选取区间中点最合适，这也就是所谓的中点公式，也就是1点高斯求积公式。

如果选取个点作为计算节点，同样可以按公式：A=k1f(x1)+k2f(x2)++knf(xn)来计算近似值，关键就是如何确定节点xi和系数ki（i=1,2,3,,n) 理论证明对于n个节点的上述求积公式，最高有2n-1次的代数精度，高斯公式就是使得上述公式具有2n-1次代数精度的积分公式。至于如何确定公式中的节点和系数，最常见的是利用勒让德多项式，具体的这里不方便说，你查查相关资料吧。

这里向你推荐一下克鲁特算法（其实就是对高斯列主元消元法进行优化，使之更适合于计算机编程），首先将矩阵A进行LU分解（将系数矩阵分解成一个上三角矩阵和一个下三角矩阵），分解的过程中用到了隐式的主元寻找法，同时利用克鲁特算法可以将两个nn矩阵压缩到一个nn矩阵中，大大节省了存储空间提高了计算速度。

方程可化为LUx=B，令Ux=y --->Ly=B

然后利用回代先求y，再利用y求x

因为该方法在求解过程中不涉及增广矩阵所以矩阵B几乎不参与什么运算，所以它的计算速度应该能够达到高斯列主元消元法的三倍，但原理与其基本一致。

而且我在程序中使用了动态数组方便你今后进行扩展。

以下程序按照《矩阵论第二版》和《C语言数值计算法方法大全》编写，LU分解部分程序主要参考了《C语言数值计算法方法大全》第二章的程序

如果你需要详细的理论讲解我可以将这两本书和源程序发给你，上面的论述相当详细足够你答辩用的了，我的邮箱hu_hu605@163com

计算结果：

A矩阵:

2 2 5

3 4 7

1 3 3

B矩阵:

5

6

5

解矩阵:

x 1=-7

x 2=0333333

x 3=366667

Press any key to continue

#include <cmath>

#include <iostream>

#include <iomanip>

#include <cstdlib>

#include <functional>

#include <vector>

#include <algorithm>

using namespace std;

#define TINY 10e-20 //A small number

#define N 3

void ludcmp(vector<vector<float> > &a, int n, vector<int> &indx, float &d);//对矩阵进行LU分解

void lubksb(vector<vector<float> > &a, int n, vector<int> &indx, vector<float> &b);//对矩阵进行前向和后向回代

void root(vector<vector<float> > &x,vector<float> &col);//解方程结果保存在y中

void iniv(vector<vector<float> > &x,vector<float> line,int n);//对二维动态数组进行初始化

void main()

{

int i,j,n=N;//输入矩阵的维数

float A[N][N]={{2,2,5},{3,4,7},{1,3,3}};//左边A矩阵

float B[N]={5,6,5};//右边B矩阵

vector<vector<float> > x;//建立动态二维数组存放A，保证你的程序进行扩展时只改A，B，N

vector<float> line;

vector<float> y(n);//建立动态数组存放B

iniv(x,line,n);

yclear();

for(i=0;i<n;i++)//将A赋给x，B赋给y

{

ypush_back(B[i]);

for(j=0;j<n;j++)

{

x[i]push_back(A[i][j]);

}

}

cout<<"A矩阵:"<<endl;

for(i=0;i<n;i++)

{

for(j=0;j<n;j++)

{

cout<<setw(2)<<setiosflags(ios::left)<<setw(2)<<x[i][j]<<" ";

}

cout<<endl;

}

cout<<"B矩阵:"<<endl;

for(i=0;i<n;i++)

{

cout<<setw(2)<<setiosflags(ios::left)<<setw(2)<<y[i]<<endl;

}

root(x,y);//求根

cout<<"解矩阵:"<<endl;

for(i=0;i<n;i++)

{

cout<<setw(2)<<setiosflags(ios::left)<<"x"<<i+1<<"="<<y[i]<<endl;

}

cout<<endl;

}

void root(vector<vector<float> > &x,vector<float> &col)

{

int n=xsize(),i=0,j=0;

vector<int> index(n);//用于记录寻找主元素过程中对矩阵的初等变换

indexclear();

float m=10;//记录变换方式，此程序中无用

ludcmp(x,n,index,m);//进行LU分解

lubksb(x,n,index,col);//根据分解结果进行回带

}

//以下程序按照《矩阵论第二版》和《C语言数值计算法方法大全》编写，LU分解部分程序主要参考了《C语言数值计算法方法大全》第二章的程序

//如果你需要详细的理论讲解我可以将这两本书和源程序发给你，我的邮箱hu_hu605@163com

void ludcmp(vector<vector<float> > &a, int n, vector<int> &indx, float &d)

{

int i,imax,j,k;

float big=0,dum=0,sum=0,temp=0;

vector<float> vv(n);

vvclear();

d=10;

for (i=0;i<n;i++)

{

big=00;

for (j=0;j<n;j++)

if ((temp=fabs(a[i][j])) > big)

big=temp;

vv[i]=10/big;

}

for (j=0;j<n;j++)

{

for (i=0;i<j;i++)

{

sum=a[i][j];

for (k=0;k<i;k++)

sum -= a[i][k]a[k][j];

a[i][j]=sum;

}

big=00;

for (i=j;i<n;i++)

{

sum=a[i][j];

for (k=0;k<j;k++)

sum -= a[i][k]a[k][j];

a[i][j]=sum;

if ( (dum=vv[i]fabs(sum)) >= big)

{

big=dum;

imax=i;

}

}

if (j != imax)

{

for (k=0;k<n;k++)

{

dum=a[imax][k];

a[imax][k]=a[j][k];

a[j][k]=dum;

}

d = -(d);

vv[imax]=vv[j];

}

indx[j]=imax;

if (a[j][j] == 00)

a[j][j]=TINY;

if (j != n)

{

dum=10/(a[j][j]);

for (i=j+1;i<n;i++)

a[i][j] = dum;

}

}

}

void lubksb(vector<vector<float> > &a, int n, vector<int> &indx, vector<float> &b)

{

int i,ii=0,ip,j;

float sum;

for(i=0;i<n;i++)//按LU分解时寻找主元所进行的初等变换进行反边变换。

{

ip=indx[i];

sum=b[ip];

b[ip]=b[i];

b[i]=sum;

}

sum=0;

for (i=1;i<n;i++)

{

sum=0;

for(j=0;j<i;j++)

{

sum+=a[i][j]b[j];

}

b[i]=b[i]-sum;

}

b[n-1]=b[n-1]/a[n-1][n-1];

for (i=n-2;i>=0;i--)

{

sum=0;

for(j=i+1;j<n;j++)

{

sum+=a[i][j]b[j];

}

b[i]=(b[i]-sum)/a[i][i];

}

}

void iniv(vector<vector<float> > &x,vector<float> line,int n)

{

int i,j;

for(i=0;i<n;i++)

{

xpush_back(line);

for(j=0;j<n;j++)

{

x[i]clear();

}

}

}

void gaussj(double a[], int n, double b[])

{

int i,j,k,l,ll,irow,icol;

double big,pivinv,dum;

int ipiv[50], indxr[50], indxc[50];

for (j=0;j<=n-1;j++)

{

ipiv[j]=0;

}

for (i=0;i<=n-1;i++)

{

big=00;

for (j=0;j<=n-1;j++)

{

if(ipiv[j]!=1)

{

for(k=0;k<=n-1;k++)

{

if(ipiv[k]==0)

{

if(fabs(a[jn+k])>=big)

{

big=fabs(a[jn+k]);

irow=j;

icol=k;

}

else if(ipiv[k]>1)

{

cout<<"singular matrix";

}

}

}

}

}

ipiv[icol]=ipiv[icol]+1;

if(irow!=icol)

{

for(l=0;l<=n-1;l++)

{

dum=(a[irown+l]);

a[irown+l]=a[icoln+l];

a[icoln+l]=dum;

}

dum=b[irow];

b[irow]=b[icol];

b[icol]=dum;

}

indxr[i]=irow;

indxc[i]=icol;

if(a[icoln+icol]==00)

{

cout<< "singular matrix";

}

pivinv=10/(a[icoln+icol]);

a[icoln+icol]=10;

for(l=0;l<=n-1;l++)

{

a[icoln+l]=a[icoln+l]pivinv;

}

b[icol]=b[icol]pivinv;

for(ll=0;ll<=n-1;ll++)

{

if(ll!=icol)

{

dum=a[lln+icol];

a[lln+icol]=00;

for(l=0;l<=n-1;l++)

{

a[lln+l]=a[lln+l]-a[icoln+l]dum;

}

b[ll]=b[ll]-b[icol]dum;

}

}

}

for(l=n-1;l<=0;l--)

{

if(indxr[l]!=indxc[l])

{

for(k=0;k<=n-1;k++)

{

dum=a[kn+indxr[l]];

a[kn+indxr[l]]=a[kn+indxc[l]];

a[kn+indxr[l]]=dum;

}

}

}

}

高斯求积公式是变步长数值积分的一种，基本形式是计算［-1，1］上的定积分。

假设现在要求f（x）在［-1，1］上的积分值，只允许计算一次f（x）的值，会选取一点x0，计算出f（x0），用A=f（x0）2作为近似值。现在问题是怎样选取x0，使得结果尽可能精确，直觉告诉我们选取区间中点最合适，这也就是所谓的中点公式，也就是1点高斯求积公式。

含义

由于磁力线总是闭合曲线，因此任何一条进入一个闭合曲面的磁力线必定会从曲面内部出来，否则这条磁力线就不会闭合起来了。如果对于一个闭合曲面，定义向外为正法线的指向，则进入曲面的磁通量为负，出来的磁通量为正，那么就可以得到通过一个闭合曲面的总磁通量为0。这个规律类似于电场中的高斯定理，因此也称为高斯定理。

百度百科-高斯公式

% 高斯-勒让德数值积分

fun=@(x)sin(x)x^3; % fun：积分表达式

% a,b：积分上下限

a = 0;

b = pi;

tol=1e-8; % tol：积分精度，默认1e-6

% 计算求积节点

syms x

p=sym2poly(diff((x^2-1)^(n+1),n+1))/(2^nfactorial(n));

tk=roots(p); % 求积节点

% 计算求积系数

Ak=zeros(n+1,1);

for i=1:n+1

xkt=tk;

xkt(i)=[];

pn=poly(xkt);

fp=@(x)polyval(pn,x)/polyval(pn,tk(i));

Ak(i)=quadl(fp,-1,1,tol); % 求积系数

end

% 积分变量代换，将[a,b]变换到[-1,1]

xk=(b-a)/2tk+(b+a)/2;

% 检验积分函数fun有效性

fun=fcnchk(fun,'vectorize');

% 计算变量代换之后积分函数的值

fx=fun(xk)(b-a)/2;

% 计算积分值

ql=sum(Akfx)

quadl(fun,0,pi) % 调用MATLAB内部积分函数检验

//TurboC 20太落后了，建议使用VC++60。

#include"stdioh"

#include"mathh"

//最大49阶

#define N 50

void Gauss(float U[N][N],int n);

void main()

{

int n,i,j;

float U[N][N];

printf("------------特殊说明---------------\n");

printf("当输出的数据含有<-1#IND>时,表示在计算过程中数据已经出现溢出!\n");

printf("-----------------------------------\n");

printf("输入对应方程的阶数:");

scanf("%d",&n);

for(i=0;i<N;i++)

for(j=0;j<N;j++)

U[i][j]=0;

printf("输入方程组的增广矩阵:\n");

for(i=0;i<n;i++)

for(j=0;j<=n;j++)

scanf("%f",&U[i][j]);

Gauss(U,n);

}

//高斯选列主元消去法

void Gauss(float U[N][N],int n)

{

int i,j,m,row;

float max,t,sum;

float result[50];

for(m=0;m<n-1;m++)

{

//选取主元

max=U[m][m];

for(i=m;i<n;i++)

{

if(fabs(max)<fabs(U[i][m]))

{

max=U[i][m];

row=i;

}

}

if(fabs(max)<001)

{

printf("主元接近于零,方法失效!\n");

return;

}

else

{

if(max!=U[m][m])

{

for(j=m;j<=n;j++)

{

t=U[m][j];

U[m][j]=U[row][j];

U[row][j]=t;

}

}

}

//消元

for(i=m+1;i<n;i++)

{

float t1,t2;

t1=U[i][m];

t2=U[m][m];

U[i][m]=0;

for(j=m+1;j<=n;j++)

U[i][j]=U[i][j]t2-U[m][j]t1;

}

}

//回代求解

for(i=n-1;i>=0;i--)

{

if(i==n-1) result[i]=U[i][i+1]/U[i][i];

else

{

sum=0;

for(j=i+1;j<n;j++)

sum=U[i][j]result[j]+sum;

result[i]=(U[i][n]-sum)/U[i][i];

}

}

//输出根

printf("高斯选列主元消去法求得的解为:\n");

for(i=0;i<n;i++)

printf("%33f ",result[i]);

printf("\n");

}

我有VB的，自己很多年前写的，一直用，但是正算－>反算－>正算后，Y坐标与原来的差了05－07mm，不知道怎么回事，这两年工作忙也没有时间再深究，但是这样的计算精度做控制足够了，如果楼主或是者是哪位同仁见此贴能顺便把这个问题解决了，咱们就一起进步了！代码如下：

'高斯坐标正算

Private Sub DadiZs()

Dim t As Double, Itp As Double, X0 As Double, N As Double, L0 As Double

Dim V As Double, ll As Double, W As Double, M As Double

Lat = Radian(Lat)

Lon = Radian(Lon)

L0 = Radian(Lo)

If Tq = 0 Then

a = 6378245 '54椭球参数

b = 635686301877305

ep = 0006693421622966

ep1 = 0006738525414683

f = (a - b) / a

c = a ^ 2 / b

d = b ^ 2 / a

X0 = 1111348611 (Lat 180# / Pi) - (320057799 Sin(Lat) + 1339238 (Sin(Lat)) ^ 3 + 06973 (Sin(Lat)) ^ 5 + 00039 (Sin(Lat)) ^ 7) Cos(Lat)

'X0 = 1111348611 (Lat 180# / Pi) - (320057798 Sin(Lat) + 1339238 (Sin(Lat)) ^ 3 + 06972 (Sin(Lat)) ^ 5 + 00039 (Sin(Lat)) ^ 7) Cos(Lat)

Else

a = 6378140 '75椭球参数

b = 635675528815753

ep = 0006694384999588

ep1 = 0006739501819473

f = (a - b) / a

c = a ^ 2 / b

d = b ^ 2 / a

X0 = 1111330047 (Lat 180 / Pi) - (320098575 Sin(Lat) + 1339602 (Sin(Lat)) ^ 3 + 06976 (Sin(Lat)) ^ 5 + 00039 (Sin(Lat)) ^ 7) Cos(Lat)

End If

ll = Lon - L0

t = Tan(Lat)

Itp = ep1 Cos(Lat) ^ 2

W = Sqr(1 - ep Sin(Lat) ^ 2)

V = Sqr(1 + ep1 Cos(Lat) ^ 2)

M = c / V ^ 3

N = a / W

'x = X0 + N t (Cos(Lat)) ^ 2 ll ^ 2 / 2 + N t (5 - t t + 9 Itp + 4 Itp Itp) (Cos(Lat)) ^ 4 ll ^ 4 / 24 + N t (61 - 58 t ^ 2 + t ^ 4 + 270 Itp - 330 t ^ 2 Itp) (Cos(Lat)) ^ 6 ll ^ 6 / 720 + N t (1385 - 3111 t ^ 2 + 543 t ^ 4 - t ^ 6) Cos(Lat) ^ 8 ll ^ 8 / 40320

x = X0 + N t (Cos(Lat)) ^ 2 ll ^ 2 / 2 + N t (5 - t t + 9 Itp ^ 2 + 4 Itp ^ 4) (Cos(Lat)) ^ 4 ll ^ 4 / 24 + N t (61 - 58 t ^ 2 + t ^ 4 + 270 Itp ^ 2 - 330 t ^ 2 Itp ^ 2) (Cos(Lat)) ^ 6 ll ^ 6 / 720 + N t (1385 - 3111 t ^ 2 + 543 t ^ 4 - t ^ 6) Cos(Lat) ^ 8 ll ^ 8 / 40320

y = N Cos(Lat) ll + N (1 - t t + Itp) (Cos(Lat)) ^ 3 ll ^ 3 / 6 + N (5 - 18 t t + t ^ 4 + 14 Itp - 58 Itp t t) (Cos(Lat)) ^ 5 ll ^ 5 / 120 + N (61 - 479 t ^ 2 + 179 t ^ 4 - t ^ 6) Cos(Lat) ^ 7 ll ^ 7 / 5040

r = Sin(Lat) ll + Sin(Lat) (Cos(Lat)) ^ 2 ll ^ 3 (1 + 3 Itp + 2 Itp ^ 2) / 3 + Sin(Lat) (Cos(Lat)) ^ 4 ll ^ 5 (2 - t t) / 15

r = Degree(r)

y = y + 500000#

End Sub

'高斯反算

Private Sub DadiFs()

Dim t As Double, Itp As Double, X0 As Double, Bf As Double, N As Double

Dim v As Double, ll As Double, W As Double, M As Double, L0 As Double

L0 = Radian(Lo)

X0 = x 0000001

y = y - 500000#

If Tq = 0 Then

a = 6378245 '54椭球参数

b = 635686301877305

ep = 0006693421622966

ep1 = 0006738525414683

f = (a - b) / a

c = a ^ 2 / b

d = b ^ 2 / a

If X0 < 3 Then

Bf = 904353301294 X0 - 000000049604 X0 ^ 2 - 000075310733 X0 ^ 3 - 000000084307 X0 ^ 4 - 000000426055 X0 ^ 5 - 000000010148 X0 ^ 6

ElseIf X0 < 6 Then

Bf = 2711115372595 + 902468257083 (X0 - 3) - 000579740442 (X0 - 3) ^ 2 - 000043532572 (X0 - 3) ^ 3 + 000004857285 (X0 - 3) ^ 4 + 000000215727 (X0 - 3) ^ 5 - 000000019399 (X0 - 3) ^ 6

End If

Else

a = 6378140 '75椭球参数

b = 635675528815753

ep = 0006694384999588

ep1 = 0006739501819473

f = (a - b) / a

c = a ^ 2 / b

d = b ^ 2 / a

If X0 < 3 Then

Bf = 904369066313 X0 - 000000049618 X0 ^ 2 - 000075325505 X0 ^ 3 - 00000008433 X0 ^ 4 - 000000426157 X0 ^ 5 - 00000001015 X0 ^ 6

ElseIf X0 < 6 Then

Bf = 2711162289465 + 902483657729 (X0 - 3) - 000579850656 (X0 - 3) ^ 2 - 000043540029 (X0 - 3) ^ 3 + 000004858357 (X0 - 3) ^ 4 + 000000215769 (X0 - 3) ^ 5 - 000000019404 (X0 - 3) ^ 6

End If

End If

Bf = Bf Pi / 180#

t = Tan(Bf)

Itp = ep1 Cos(Bf) ^ 2

W = Sqr(1 - ep Sin(Bf) ^ 2)

v = Sqr(1 + ep1 Cos(Bf) ^ 2)

M = c / v ^ 3

N = a / W

Lat = Bf - 05 v ^ 2 t ((y / N) ^ 2 - (5 + 3 t t + Itp - 9 Itp t t) (y / N) ^ 4 / 12 + (61 + 90 t t + 45 t ^ 4) (y / N) ^ 6 / 360)

ll = ((y / N) - (1 + 2 t t + Itp) (y / N) ^ 3 / 6 + (5 + 28 t t + 24 t ^ 4 + 6 Itp + 8 Itp t t) (y / N) ^ 5 / 120) / Cos(Bf)

r = y t / N - y ^ 3 t (1 + t t - Itp) / (3 N ^ 3) + y ^ 5 t (2 + 5 t t + 3 t ^ 4) / (15 N ^ 5)

Lat = Degree(Lat)

Lon = Degree(L0 + ll)

r = Degree(r)

End Sub

有了正反算，换带也就完成了！

用到的子程序：

Public Const Pi = 314159265358979, p = 206264806

Public Cktq As String

'角度化弧度

Public Function Radian(a As Double) As Double

Dim Ro As Double

Dim c As Double

Dim Fs As Double

Dim Ib As Integer

Dim Ic As Integer

If a < 0 Then a = -a: t = 1

Ro = Pi / 180#

Ib = Int(a)

c = (a - Ib) 100#

Ic = Int(c + 0000000000001)

Fs = (c - Ic) 100#

If t = 1 Then Radian = -(Ib + Ic / 60# + Fs / 3600#) Ro Else Radian = (Ib + Ic / 60# + Fs / 3600#) Ro

End Function

'弧度化角度

Public Function Degree(a As Double) As Double

Dim Bo As Double

Dim Fs As Double

Dim Im As Integer

Dim Id As Integer

If a < 0 Then a = -a: t = 1

Bo = a

Call DMS(Bo, Id, Im, Fs)

If t = 1 Then Degree = -(Id + Im / 100# + Fs / 10000#) Else Degree = Id + Im / 100# + Fs / 10000#

End Function

Public Sub DMS(a As Double, Id As Integer, Im As Integer, Fs As Double)

Dim Bo As Double

Dim c As Double

c = a

c = 180# / Pi c

Id = Int(c)

Bo = (c - Id) 60

Im = Int(Bo)

Fs = (Bo - Im) 60

End Sub

'取位计算

Public Function Qw(a As Double, Ws As Integer) As Double

Qw = Int(a 10 ^ Ws + 05) / 10 ^ Ws

End Function

以上就是关于求积分的两点高斯公式全部的内容，包括:求积分的两点高斯公式、求C语言课程设计：用高斯列主元消元法解线性方程组、【编程求助】用c语言或者c++编程，实现用高斯消元法求解线性方程组Ax=b。等相关内容解答，如果想了解更多相关内容，可以关注我们，你们的支持是我们更新的动力！