龙贝格算法就积分的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 = 0k <nk++)

{

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 = 1k <= mk++)

{

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 = 0i <mi++)

{

int j

if (!(i % j))

{

cout<<T[i]<<endl

}

else

{

cout<<T[i]<<" "

}

j++

下面是使用龙贝格算法求积分的matlab程序代码

clear

clc

format

long

f='4/(1+x^2)'

%这是被积函数

x='x'

%这是弯消拆被积自变量

a=0

%这是积分下限

b=1

%这是积分上限

e=1e-5

%这是积分误差埋枣限制

%以下是龙贝格积分算法，是目前最为成熟的积分算法，具有收敛速度快，精度可以自定义的优点

%

I为积分的估计值

%

n为迭代次数，2^(n-1)是桥隐等分区间的份数

T(1,1)=(b-a)/2*(subs(f,x,a)+subs(f,x,b))

T=double(T)

n=2

h=b-a

T(2,1)=T(1,1)/2+h/2*double(subs(f,x,a+h/2))

T(2,2)=4/3*T(2,1)-1/3*T(1,1)

d=T(2,2)-T(1,1)

while

d>e

n=n+1

h=h/2

T(n,1)=T(n-1,1)/2

for

i=1:2^(n-2)

T(n,1)=T(n,1)+h/2*double(subs(f,x,a+(i-1/2)*h))

end

for

i=2:n

k=4^(i-1)

T(n,i)=k/(k-1)*T(n,i-1)-1/(k-1)*T(n-1,i-1)

end

d=abs(T(n,n)-T(n-1,n-1))

end

I=T(n,n)

%输出计算值

望采纳！谢谢！

// Romberg.cpp : Defines the entry point for the console application.

//

#include "stdafx.h"

const double e = 1.0E-8

const double end = 1.0E-6

//积分上下限

const double a = 0.0

const double b = 1.0

//被积函数

double function(long double x)

{

if(abs(x)<e) return 1.0

else return sin(x)/x

}

int min(int x, int y)

{

return x>y?y:x

}

int _tmain(int argc, _TCHAR* argv[])

{

double * temp

double * t = new double[1]

int m = 1

int j = 0

double h

t[0] = (b-a)/2*(function(a)+function(b))//t_0[0]

while(true)

{

temp = new double[m]//创建数组temp[m]

int k = sizeof(double)*min(4,m)

::memcpy(temp,t,k)//将t数组的数拷到temp数组

delete []t

t = new double[m+1]//当前行比上一行多一个数

h = (b-a)/pow(2.0,m)//步长

t[0] = temp[0]/2

for(int i = 1i<=pow(2,m-1)i++)

{

t[0] += h*function(a+(2*i-1)*h)//相当于t_k[0]

}

for(int i = 1i<=min(3,m)i++)

{

t[i] = (pow(4.0,i)*t[i-1]-temp[i-1])/(pow(4.0,i)-1)//t_m[i]

}

//产生4行数后判断精度是否满足要求

if(m>3)

{

if(abs(t[3]-temp[3])/t[3]<end)//判断第4列（八阶收敛）

{

//满足精度要求，输出结果，跳出循环

printf("%d\n",m)//循环次数

printf("%2.15f\n",t[3])//积分结果

//std::cout<<t[3]<<std::endl//这是C++的输出函数。

break

}

}

delete []temp

m++

}

delete []t

return 0

}

//由于第k+1行的数可以由第k行的数求出，与1，2,...,k-1行没有直接关系，可以用两个数组桐卖大存储计算结果，前面算过的数不符合要求就

//直接删掉了。这样可以节约内存。这两个数组的长度是变化的，随着数表配掘行局竖数逐渐加长。

---