Romberg算法的m程序

f为m文件定义的函数，从a到b积分，n为次数参数，求R(n,n)

最后结果以列表形式输出，同时将最后两行的结果存储在r中。

%Romberg Integration

function r = romberg(f,a,b,n)

h = b-a

r(1,1) = h*(feval(f,a)+feval(f,b))/2

for i =2:n

sum = 0

for k =1:2^(i-2)

sum = sum+feval(f,a+(k-0.5)*h)

end

r(2,1)=(r(1,1)+h*sum)/2

for j = 2:i

r(2,j)=r(2,j-1)+(r(2,j-1)-r(1,j-1))/(4^(j-1)-1)

end

r(2,:)

h = h/2

r(1,:)=r(2,:)

end

误差界eps%被积函数为f(x)=(x^3+sin(x))/x积分区间为[0.3,0.8]

#include<stdio.h>

#include<stdlib.h>

#include<math.h>

int main(void)

{

/* int i,j,n

float *d=(float *)malloc(sizeof(float)*N)

float *x=(float *)malloc(sizeof(float)*N)

float *y=(float *)malloc(sizeof(float)*N)

float *u=(float *)malloc(sizeof(float)*N)

free(a)free(b)free(c)free(d)free(x)free(y)free(l)

return 0

*/

double b=0.8

double a=0.3

double h=0.0//

double eps=1.0e-5//误差界eps

int kmax=20//最大递推次数

double T1=0.0,S1=0.0,C1=0.0,R1=0.0,T2=0.0,S2=0.0,C2=0.0,R2=0.0

double sum

double *x,*f(x)

int i

h=b-a

T1=h/2*((pow(a,3)+sin(a))/a+(pow(b,3)+sin(b))/b)

printf("T1:%13.12f\n",T1)

for(int k=0k<kmaxk++)

{

h=(b-a)/(pow(2,k+1))

x=(double *)malloc(sizeof(double)*int(pow(2,k)))

for(i=0i<pow(2,k)i++)

{

x[i]=a+(2*i+1)*h

}

fx=(double *)malloc(sizeof(double)*int(pow(2,k)))

sum=0.0

for(i=0i<pow(2,k)i++)

{

fx[i]=(pow(x[i],3)+sin(x[i]))/x[i]

sum+=fx[i]

}

T2=T1/2+sum*h

printf("T2:%13.12f\n",T2)

S2=T2+(T2-T1)/3

printf("S%d:%13.12f\n",int(pow(2,k)),S2)

if(k<2)

{

if(k==1)

{

C2=S2+(S2-S1)/15

printf("C1:%13.12f\n",C2)

}

}

else

{

C2=S2+(S2-S1)/15

printf("C%d:%13.12f\n",int(pow(2,k-1)),C2)

R2=C2+(C2-C1)/63

printf("R%d:%13.12f\n",int(pow(2,k-2)),R2)

if(fabs(R2-R1)<eps)

break

R1=R2

}

T1=T2S1=S2C1=C2

free(x)free(fx)

}

printf("所求积分I=%13.12f\n",R2)

return 0

}

很久之前写的了，改了一下，满足你的要求了。

这里a是积分的下限。b是积分上限。epsilon是一个精确度，如果越小的话，迭代次数越多，越精确。

这个程序会输出每次迭代的过程。

不要问我龙贝格的算法了。。我已经忘了，这个程序应该是我3年前写的。。

#include <iostream>

#include <cmath>

#include <iomanip>

using namespace std

double f(double x) //函数f(x)

{

    if(x == 0)

        return 1

   return sin(x) / x

}

int main()

{

    double a = 0, b = 1, epsilon = 0.00000001

    int m = 1, k = 1

    double h = (b - a) / 2.0

    double T0 = h * (f(a) + f(b)), T = 3

    double F = 0

    while(fabs(T - T0) >= 3 * epsilon)

    {

       if(m != 1)

            T0 = T

        F = 0

        k = pow(2., m - 1)

       for(int i = 1 i <= k  i++)

        {

            F += f(a + (2 * i - 1) * h)

        }

        T = T0 / 2.0 + h * F

        m += 1

        h /= 2.0

        cout << setprecision(16) << m << "次迭代" << " T = " << T << endl

    }

    return 0

}

---