用Romberg算法计算下列积分的C语言程序

误差界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

}

先用另外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

% introductory 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

---