加载中…
个人资料
  • 博客等级:
  • 博客积分:
  • 博客访问:
  • 关注人气:
  • 获赠金笔:0支
  • 赠出金笔:0支
  • 荣誉徽章:
正文 字体大小:

编程实现数值积分的几种方法 c语言

(2010-04-08 13:31:01)
标签:

杂谈

分类: 考试系统

复化梯形公式,复化抛物线公式和Romberg求积法的算法程序:

以下程序均定义误差限为1*10^-5;

1)复化梯形公式:

#include <stdio.h>
#include <math.h>
#define e 1e-5
#define a 0                          //积分下限a
#define b 1                          //积分上限b
#define f(x) (4/(1+(x*x)))           //被积函数f(x)
int main()
{
    int i,n;
    double h,t0,t,g;
    n=1;                             //赋初值
    h=(double)(b-a)/2;
    t=h*(f(a)+f(b));
    do                             
    {
       t0=t;           
       g=0;
       for (i=1;i<=n;i++)
           g+=f((a+(2*i-1)*h));
       t=(t0/2)+(h*g);                //复化梯形公式
       n*=2;
       h/=2;
    }
    while (fabs(t-t0)>e);             //自定义误差限e
    printf("%.8lf",t);                //输出积分的近似值
  
return 0;
}

2)复化抛物线公式:

#include <stdio.h>
#include <math.h>
#define e 1e-5
#define a 0                          //积分下限a
#define b 1                          //积分上限b
#define f(x) (4/(1+(x*x)))           //被积函数f(x)
int main()
{
    int i,n;
    double f1,f2,f3,h,s0,s;
    f1=f(a)+f(b);                     //赋初值
    f2=f(((double)(b+a)/2));
    f3=0;
    s=((double)(b-a)/6)*(f1+4*f2);
    n=2;
    h=(double)(b-a)/4;
    do                                //复化抛物线算法
    {
                      f2+=f3;
                      s0=s;
                      f3=0;
                      for (i=1;i<=n;i++)
                          f3+=f((a+(2*i-1)*h));
                      s=(h/3)*(f1+2*f2+4*f3);
                      n*=2;
                      h/=2;
  
    }
    while (fabs(s-s0)>e);               //自定义误差限
    printf("%.8lf",s);
  
    return 0;
}

3)Romberg求积法:

#include <stdio.h>
#include <math.h>
#define e 1e-5
#define a 0                          //积分下限a
#define b 1                          //积分上限b
#define f(x) (4/(1+(x*x)))           //被积函数f(x)
double t[100][100];
int main()
{
    int n,k,i,m;
    double h,g,p;
    h=(double)(b-a)/2;
    t[0][0]=h*(f(a)+f(b));
    k=1;
    n=1;
    do                                //Romberg算法
    {
        g=0;
        for (i=1;i<=n;i++)
            g+=f((a+((2*i-1)*h)));
        t[k][0]=(t[k-1][0]/2)+(h*g);
        for (m=1;m<=k;m++)
        {
            p=pow(4,(double)(m));
            t[k-m][m]=(p*t[k-m+1][m-1]-t[k-m][m-1])/(p-1);
        }
        m-=1;
        h/=2;
        n*=2;
        k+=1;
    }
    while (fabs(t[0][m]-t[0][m-1])>e);      //自定义误差限e
    printf("%.8lf",t[0][m]);
  
    return 0;
}

给定精度,定义误差限为1*10^-5,分别求出步长的先验估计值:用复化梯形公式计算,要求h<0. 007746。用复化抛物线公式计算,要求h<0.131607。而实际上,在运用后验估计的程序中,以相同的精度,得到的复化梯形步长为h= 0.003906,得到的复化抛物线步长为h=0.0625,它们大致上分别为先验估计值的一半,符合要求。

在实践中比较上体三种算法的计算量,当取相同精度1*10^-5时,复化梯形调用函数f257次,孵化抛物线公式调用函数f17次,Romberg算法调用函数 f17次,从计算量上看,后两者较小。

0

阅读 收藏 喜欢 打印举报/Report
  

新浪BLOG意见反馈留言板 欢迎批评指正

新浪简介 | About Sina | 广告服务 | 联系我们 | 招聘信息 | 网站律师 | SINA English | 产品答疑

新浪公司 版权所有