用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

}

---