用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

}

%龙贝格求积算法function I=romberg(a,b)h=b-aT(1)=h/返埋2*(fun(a)+fun(b))m=1while 1h=h/2 S(1)=1/2*T(1)+h*sumf(2^(m-1),a,h) for j=1:mS(j+1)=S(j)+(S(j)-T(j))/(4^j-1) endif abs(S(m+1)-T(m))<1e-6break endT=Sm=m+1endI=S(m+1)end

function f=sumf(m,a,h)for j=1:my(j)=fun(a+(2*j-1)*h)endf=sum(y)end

function f=fun(x)f=x/态祥(4+x^2)end结果：

>>漏闭蚂 I=romberg(0,1)

I =

0.111571775646293

>>

首先解决怎么算，计算机肯定不会积分，所以我开始想用sinx的泰勒展开式，然后选3-4次作为近似，然后积分。听你说梯形法，是数值计算态悄的内容，刚好这学期在学，就把我调试的程序发一个给你吧这是romberg算法,把a 换为0，b换为桐兄pi就好了吧。附上书上的代码。

#include<stdio.h>

#include<math.h>

#define f(x) (sin(x))

#define N_H  20

#define MAXREPT 10

#define a  1.0

#define b  2.0

#define epsilon 0.00001

double computeT(double aa,double bb,long int n)

{

 int idouble sum=0.0double h=(bb-aa)/n

 for(i=1i<ni++)

 {

    sum+=f(aa+i*h)

 }

 sum+=(f(aa)+f(bb))/2

 return (h*sum)

}

void main()

{

 int i

 long int n=N_H,m=0

 double T[MAXREPT+1][2]

 T[0][1]=computeT(a,b,n)

 n*=2

 for (m=1m<MAXREPTm++)

 {

    for (i=0i<mi++)

   局闭袭 {

T[i][0]=T[i][1]

    }

    T[0][1]=computeT(a,b,n)

    n*=2

    for (i=1i<=mi++)

    {

 T[i][1]=T[i-1][1]+(T[i-1][1]-T[i-1][0])/(pow(2,2*m)-1)

 if((T[m-1][1]<T[m][1]+epsilon)&&(T[m-1][1]>T[m][1]-epsilon))

 {

    printf("the integrate is %lf\n",T[m][1])

    return

 }

    }

 }

 printf("return no solved...\n")

}

---