龍貝格演算法計算積分

Ⅰ 分別用牛頓-科特斯公式,復合積分公式,龍貝格演算法,高斯勒讓德公式計算下面的積分

摘要 計算定積分的方法很多，而高斯—勒讓德公式就是其中之一。 高斯積分法是精度最高的插值型數值積分，具有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++;