C语言实现龙贝格积分

《数值分析》实验习题2

用Romberg算法求下列积分，允许误差eps=0.00001

/*****************************************************************《数值分析》实验习题2机械工程现代制造技术教育部重点实验室 *****************************************************************/#include<stdio.h>#include<malloc.h>#include<math.h>#define PI 3.14159265typedef struct node  {int num;    double T;    double S;    double C;    double R;    struct node *next; }ParaList;double fun(double );double calT(ParaList* ,double ,double );double calS(ParaList*);double calC(ParaList*);double calR(ParaList*);int main(){double a = 0, b = 3, eps = 0.00001;//第二题参数 //double a = 1, b = 3, eps = 0.00001;//printf("请输入积分下限\n");//scanf("%lf", &a);//printf("请输入积分上限\n");//scanf("%lf", &b);//printf("请输入eps");//scanf("%lf", &eps);//算前四行的T，S，C的值并存入参数链表中 ParaList *phead, *p;int i;phead = (ParaList*)malloc(sizeof(ParaList));phead->num = 1;phead->T = (b - a) / 2 * (fun(a) + fun(b));phead->next = p = (ParaList*)malloc(sizeof(ParaList));p->num = 2;p->T = calT(phead, a, b);p->S = calS(phead);p->next = (ParaList*)malloc(sizeof(ParaList));p->next->num = 3;p->next->T = calT(p, a, b);p->next->S = calS(p);p->next->C = calC(p);p = p->next;p->next = (ParaList*)malloc(sizeof(ParaList));p->next->num = 4;p->next->T = calT(p, a, b);p->next->S = calS(p);p->next->C = calC(p);p->next->R = calR(p);p = p->next;for(i=5; ;i++){p->next = (ParaList*)malloc(sizeof(ParaList));p->next->num = i;p->next->T = calT(p, a, b);p->next->S = calS(p);p->next->C = calC(p);p->next->R = calR(p);if(abs(p->next->R - p->R) < eps){printf("计算结果为%.8lf，共进行了%d次迭代", p->next->R, p->next->num);break;}} }//double fun(double x){////第一题函数 //double fresult;//fresult = x * pow(1 + pow(x, 2), 0.5);//return fresult;//}double fun(double x){//第二题函数 double fresult;fresult = pow(3, x) * pow(x, 1.4) * (5* x + 7) * sin(pow(x, 2) * PI / 180);return fresult;}double calT(ParaList* s, double a, double b){//传入某一行的指针和积分上下限，用于计算下一行的T值 int i;double sum = 0, tresult;for(i=0; i<pow(2, s->num -1); i++)sum += fun(a + (i + 0.5) * (b - a) / pow(2, s->num -1));tresult = 0.5 * (s->T + (b - a) / pow(2, s->num -1) * sum);return tresult;}double calS(ParaList* t){//传入某一行的指针，用于计算下一行的S值 double sresult;sresult =(4 * t->next->T - t->T) / 3;return sresult; } double calC(ParaList* s){//传入某一行的指针，用于计算下一行的C值 double cresult;cresult = (16 * s->next->S - s->S) / 15;return cresult;}double calR(ParaList* c){//传入某一行的指针，用于计算下一行的R值 double rresult;rresult = (64 * c->next->C - c->C) / 63;return rresult;}