龙贝格算法计算积分

Ⅰ 分别用牛顿-科特斯公式,复合积分公式,龙贝格算法,高斯勒让德公式计算下面的积分

摘要 计算定积分的方法很多，而高斯—勒让德公式就是其中之一。 高斯积分法是精度最高的插值型数值积分，具有2n+1阶精度，并且高斯积分总是稳定。而高斯求积系数，可以由Lagrange多项式插值系数进行积分得到。 高斯—勒让德求积公式是构造高精度差值积分的最好方法之一。他是通过让节点和积分系数待定让函数f(x)以此取i=0,1,2....n次多项式使其尽可能多的能够精确成立来求出积分节点和积分系数。高斯积分的代数精度是2n-1，而且是最高的。通常运用的是（-1，1）的积分节点和积分系数，其他积分域是通过变换x=(b-a)t/2 +(a+b)/2 变换到-1到1之间积分。

Ⅱ 龙贝格积分公式的Fortran程序代码

C++的行么？
Romberg.h：文件

// Romberg.h: interface for the Romberg class.
//
//////////////////////////////////////////////////////////////////////

#if !defined(AFX_ROMBERG_H__3DC3390F_E15E_4BB7_98E5_64C7538BF4DD__INCLUDED_)
#define AFX_ROMBERG_H__3DC3390F_E15E_4BB7_98E5_64C7538BF4DD__INCLUDED_

#if _MSC_VER > 1000
#pragma once
#endif // _MSC_VER > 1000

#include "afxtempl.h" //包含数组CArray模板
#include "math.h"

typedef double (*Fun)(double);

//Romberg算法
//方法取自《计算方法引论》(第二版)徐萃薇、孙绳武着 高等教育出版社 第119页
//但有改动，只使用了两个一维数组
//fun--被积函数的地址，[a，b]--积分区间，error--误差；
//count--计算次数，达到次数即使不满足精度也返回积分值，默认count=0--以误差为准

class CRomberg
{
public:
static double Romberg(Fun fun,double a,double b,double error,int count=0);
CRomberg();
virtual ~CRomberg();

};

#endif // !defined(AFX_ROMBERG_H__3DC3390F_E15E_4BB7_98E5_64C7538BF4DD__INCLUDED_)

Romberg.cpp：文件

// Romberg.cpp: implementation of the Romberg class.
//
//////////////////////////////////////////////////////////////////////

#include "stdafx.h"
//#include "NumericalMethods.h"
#include "Romberg.h"

#ifdef _DEBUG
#undef THIS_FILE
static char THIS_FILE[]=__FILE__;
#define new DEBUG_NEW
#endif

typedef CArray<double,double> CDoubleArray;

//在C++中可以直接使用fun(a)来表示C中的(*fun)(a)
//而不必定义下面的这个宏，这里只是明确一下:)
#define fun(a) (*fun)(a)

#define array1(a) (*array1)[a]
#define array0(a) (*array0)[a]

//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////

CRomberg::CRomberg()
{

}

CRomberg::~CRomberg()
{

}

//方法取自《计算方法引论》(第二版)徐萃薇、孙绳武着 高等教育出版社 第119页
//但有改动，只使用了两个一维数组
//fun--被积函数的地址，[a，b]--积分区间，error--误差；
//count--计算次数，达到次数即使不满足精度也返回积分值，默认count<=0--以误差为准
double CRomberg::Romberg(Fun fun, double a, double b, double error,int count)
{
double h=(b-a)/2,F=0.0f;
int i,w,m=0,n=1;
CDoubleArray * array0=new CDoubleArray();
CDoubleArray * array1=new CDoubleArray();
CDoubleArray * arraytemp;

array0->Add(h*(fun(a)+fun(b)));
do
{
F=0.0f,m++;//注意m初值为0
for(i=1;i<=n;i++)
F+=fun(a+(2*i-1)*h);
array1->SetAtGrow(0,h*F+array0(0)/2);
w=4;
for(i=0;i<m;i++)
{
array1->SetAtGrow(i+1,(w*array1(i)-array0(i))/(w-1));
w<<=1;
}
h/=2,n<<=1;
arraytemp=array0,array0=array1,array1=arraytemp;
F=fabs(array0(m)-array1(m-1));

if(F<error) break;
if (count>0)
if(m>count) break;

}while(true);

F=array0(m);
delete array0;
delete array1;

return F;
}

Ⅲ 求龙贝格积分公式fortran程序，不要其他的。

龙贝格积分公式Sn=(4T2n-Tn)/3,
Cn=(4^2S2n-Sn)/(4^2-1),
Rn=(4^3C2n-Cn)/(4^3-1),
其中,2n,n都是下标.
同理,依次类推.
这是在变步长求积过程中的三个加速公式,将粗糙的积分近似值迅速加工成精度较高的积分近似值的求积方法为龙为个求积算法.
这个地方输入法问题,麻烦你的眼睛辨认了！

Ⅳ 龙贝格求积，c语言代码

#include<iostream.h>
#include<math.h>
# define Precision 0.000001//积分精度要求
# define e 2.71828183
#define MAXRepeat 10 //最大允许重复

double function(double x)//被积函数
{
double s;
s=pow(e,x)*cos(x);
return s;
}
double Romberg(double a,double b,double f(double x))
{
int m,n,k;
double y[MAXRepeat],h,ep,p,xk,s,q;
h=b-a;
y[0]=h*(f(a)+f(b))/2.0;//计算T`1`(h)=1/2(b-a)(f(a)+f(b));
m=1;
n=1;
ep=Precision+1;
while((ep>=Precision)&&(m<MAXRepeat))
{
p=0.0;
for(k=0;k<n;k++)
{
xk=a+(k+0.5)*h; // n-1
p=p+f(xk); //计算∑f(xk+h/2),T
} // k=0
p=(y[0]+h*p)/2.0; //T`m`(h/2),变步长梯形求积公式
s=1.0;
for(k=1;k<=m;k++)
{
s=4.0*s;// pow(4,m)
q=(s*p-y[k-1])/(s-1.0);//[pow(4,m)T`m`(h/2)-T`m`(h)]/[pow(4,m)-1],2m阶牛顿柯斯特公式，即龙贝格公式
y[k-1]=p;
p=q;
}
ep=fabs(q-y[m-1]);//前后两步计算结果比较求精度
m=m+1;
y[m-1]=q;
n=n+n; // 2 4 8 16
h=h/2.0;//二倍分割区间
}
return q;
}

main()
{
double a,b,Result;
cout<<"请输入积分下限:"<<endl;
cin>>a;
cout<<"请输入积分上限:"<<endl;
cin>>b;
Result=Romberg( a, b, function);
cout<<"龙贝格积分结果："<<Result<<endl;
return 0;
}

本文来自CSDN博客，转载请标明出处：http://blog.csdn.net/liuhongyi666/archive/2008/12/23/3589002.aspx

Ⅳ 龙贝格求积公式的简介

龙贝格求积公式也称为逐次分半加速法。它是在梯形公式、辛普森公式和柯特斯公式之间的关系的基础上，构造出一种加速计算积分的方法。 作为一种外推算法, 它在不增加计算量的前提下提高了误差的精度.
在等距基点的情况下，用计算机计算积分值通常都采用把区间逐次分半的方法进行。这样，前一次分割得到的函数值在分半以后仍可被利用,且易于编程。

Ⅵ 关于matlab，龙贝格法求积分，求高手看一下错在哪里

%正解地
function y=Romberg(a,b,e)
n=10;%初始量，大致设一下
h=b-a;
t=zeros(n+1,n+1);
t(1,1)=h/4*(f(a)+2*f((a+b)/2)+f(b));
for i=1:n
t(1,i+1)=1/2*t(1,i);
for x=a+h/(2^i):h/2^(i-1):b-h/(2^i)
t(1,i+1)=t(1,i+1)+h/(2^i)*f(x);
end
end

for i=1:2
t(2,i)=4/3*t(1,i+1)-1/3*t(1,i);
end
i=3;
while(abs(t(2,i-1)-t(2,i-2))>e)
t(2,i)=4/3*t(1,i+1)-1/3*t(1,i);
y=t(2,i);
i=i+1;
if i>n break;
end
end

for i=1:2
t(3,i)=16/15*t(2,i+1)-1/15*t(2,i);
end
i=3;
while(abs(t(3,i-1)-t(3,i-2))>e)
t(3,i)=16/15*t(2,i+1)-1/15*t(2,i);
y=t(3,i);
i=i+1;
if i>n-1 break;
end
end

Ⅶ 用龙贝格法求积分

先用另外2种方法。

format long
%【1】精确值。符号积分
it=int('(2/sqrt(pi))*exp(-x)',0,1)
Accurate=eval(it)
y=inline('(2/sqrt(pi))*exp(-x)')
%【2】Simpson方法
Simpson=quad(y,0,1)
delta=Simpson-Accurate

结果：
Accurate = 0.713271669674918

y = Inline function:
y(x) = (2/sqrt(pi))*exp(-x)

Simpson = 0.713271671228492

delta = 1.553574158208448e-009

【3】从网上找到一个，存为romberg.m
%=================
function R = romberg(f, a, b, n)
format long
% ROMBERG -- Compute Romberg table integral approximation.
%
% SYNOPSIS:
% R = romberg(f, a, b, n)
%
% DESCRIPTION:
% Computes the complete Romberg table approximation to the integral
%
% / b
% I = | f(x) dx
% / a
%
% PARAMETERS:
% f - The integrand. Assumed to be a function callable as
% y = f(x)
% with `x' in [a, b].
% a - Left integration interval endpoint.
% b - Right integration interval endpoint.
% n - Maximum level in Romberg table.
%
% RETURNS:
% R - Romberg table. Represented as an (n+1)-by-(n+1) lower
% triangular matrix of integral approximations.
%
% SEE ALSO:
% TRAPZ, QUAD, QUADL.

% NOTE: all indices adjusted for MATLAB's one-based indexing scheme.

% Pre-allocate the Romberg table. Avoids subsequent re-allocation which
% is often very costly.
R = zeros([n + 1, n + 1]);

% Initial approximation. Single interval trapezoidal rule.
R(0+1, 0+1) = (b - a) / 2 * (feval(f, a) + feval(f, b));

% First column of Romberg table. Increasingly accurate trapezoidal
% approximations.
for i = 1 : n,
h = (b - a) / 2^i;

s = 0;
for k = 1 : 2^(i-1),
s = s + feval(f, a + (2*k - 1)*h);
end

R(i+1, 0+1) = R(i-1+1, 0+1)/2 + h*s;
end

% Richardson extrapolation gives remainder of Romberg table.
%
% Note: The table is computed by columns rather than the more
% traditional row version. The reason is that this prevents frequent
% (and needless) re-computation of the `fac' quantity.
%
% Moreover, MATLAB matrices internally use ``column major'' ordering so
% this version is less harmful to computer memory cache systems. This
% reason is an implementational detail, though, and less important in
% introctory courses such as MA2501.
for j = 1 : n,
fac = 1 / (4^j - 1);
for m = j : n,
R(m+1, j+1) = R(m+1, j-1+1) + fac*(R(m+1, j-1+1) - R(m-1+1, j-1+1));
end
end

function ff=f(x)
ff=2/sqrt(pi)*exp(-x);
%=================

运行：
>> R=romberg('f', 0, 1, 5)

R =

0.771743332258054 0 0 0 0 0
0.728069946441243 0.713512151168973 0 0 0 0
0.716982762290904 0.713287034240791 0.713272026445579 0 0 0
0.714200167058928 0.713272635314936 0.713271675386546 0.713271669814180 0 0
0.713503839348432 0.713271730111600 0.713271669764711 0.713271669675476 0.713271669674932 0
0.713329714927254 0.713271673453528 0.713271669676323 0.713271669674920 0.713271669674918 0.713271669674918

Ⅷ MATLAB用龙贝格积分算法计算f=x^2在0到1上的积分

f = @(x) x.*x % f是被积函数
a = 0;b = 1; % a,b分别是积分的上下限
n = 5; % n+1是T数表的列数
delta = 10^(-8); % 允许误差

M=1;
h=b-a;
err=1;
J=0;
R=zeros(4,4);
R(1,1)=h*(feval(f,a)+feval(f,b))/2;
while ((err>delta)&&(J<n))||(J<4)
J=J+1;
h=h/2;
s=0;
for p=1:M
x=a+h*(2*p-1);
s=s+feval(f,x);
end
R(J+1,1)=R(J,1)/2+h*s;
M=2*M;
for K=1:J
R(J+1,K+1)=R(J+1,K)+(R(J+1,K)-R(J,K))/(4^K-1);
end
err=abs(R(J,J)-R(J+1,K+1));
end

quad11=R(J+1,J+1);

quad11 % 所求积分值
R % R是T数表

Ⅸ 龙贝格算法就积分的c语言程序

C++实现如下：
#include<iostream>
#include<cmath>

using namespace std;

const int MAXRepeat = 100; //最大允许重复

double function(double x)//被积函数,根据自己的需要手工输入
{
double s;
s = 1.0 / (1 + x);
return s;
}
void Romberg(double a, double b, double epsion, double f(double x))
{
int m = 1;
int n = 1;
int k;
double h;
double ep;
double p;
double xk;
double s;
double q;
double T[MAXRepeat];

h = b - a;
T[0] = 0.5 * h * (f(a) + f(b));
ep = epsion + 1.0;
while ((ep >= epsion) && (m < MAXRepeat))
{
p = 0.0;
for (k = 0; k < n; k++)
{
xk = a + (k + 0.5) * h; // n-1
p = p + f(xk); //计算∑f(xk+h/2),T
} // k=0
p = (T[0] + h * p) / 2.0; //T`m`(h/2),变步长梯形求积公式
s = 1.0;
for (k = 1; k <= m; k++)
{
s = 4.0 * s; //[pwww.hbbz08.com ow(4,m)T`m`(h/2)-T`m`(h)]/[pow(4,m)-1],2m阶牛顿柯斯特公式，即龙贝格公式
q = (s * p - T[k - 1]) / (s - 1.0);
T[k-1] = p;
p = q;
}
ep = fabs(q - T[m - 1]);
m++;
T[m - 1] = q;
n++; // 2 4 8 16
h /= 2.0;
}

for (int i = 0; i < m; i++)
{
int j;
if (!(i % j))
{
cout<<T[i]<<endl;
}
else
{
cout<<T[i]<<" ";
}
j++;