【添加新的自适应程序】这个简单的辛普森算法对吗？

sslchi

写了一个简短的代码，用辛普森公式计算积分，一开始的时候觉得误差很大
后来发现是算法中的最后一步不对，文中有注释。当然这是最简单的一个，使用不方便。有没有更简洁、更容易调用的代码？
自适应格式需要用到递归(貌似递归很费时啊)。如果可以的话，怎么才使程序以一个函数为输入量？？
这样的话就可以很方便的调用了（类似于matlab里面的函数句柄）。
另外的一个问题是，我发现精度只能达到 e-8 左右，这又是什么原因呢？ 如何能够使精度达到最高呢(相同的节点个数)

[Fortran] 纯文本查看 复制代码

program main 
    !use mkl95_blas
    integer,parameter :: M=100
    real(kind=8)::N,f(M),d(M),res,s
    s=sin(0.0)!计算在a处的值
    N=M/1.0   !用于计算步长（步长之倒数）
    f= (/(i,i=1,M,1)/)/N !计算节点
    d=sin(f)             !节点处的函数值
    do counter=1,M,2     ! 计算奇数节点的累加值(0算第0个点)
        s=s+d(counter)*4.0/6.0
    end do
    do counter=2,M-2,2   ! 将计算偶数节点的值累加(不包括最后一个节点)
        s=s+d(counter)*2.0/6.0
    end do
    s=s+d(M)/6.0             !将最后一个节点加上去。。。就是这里忘了除以6
    s=2*s/N
    res=s-(cos(0.0)-cos(1.0)) !计算误差
    write(*,*) res
    read(*,*) 
end

今天又参照Moler(matlab作者)的matlab程序写了一个自适应格式的辛普森公式，但是还是无法避免把被积函数封装到module里面去了，
另外，由于目前对Fortran语法不是很熟悉，所以并没有设置报错和默认的误差上限(tol)，如果有人感兴趣，希望可以帮助改进，将被积函数
放到module外面，然后就可作为一个可以被其他程序随便调用的子程序存在(现在想要改变被积函数需要修改module)，欢迎将修改精简的代码
和我分享

[Fortran] 纯文本查看 复制代码

module Simpson
    !________________________________________________________
    ! Code by sslchi @ 2014
    ! 您可以随意使用、修改，同时欢迎将改进精简的代码和我分享
    !______________________________________________________
    contains
    
    
    recursive subroutine SimpRc(func,a,b,tol,fa,fc,fb,Q,fcount)
    !_________________________________________________________
    ! 递归子程序，用复化辛普森公式计算函数func在区间[a,b]上的
    ! 黎曼积分，但是该程序没有提供判断积分是否可积，仅能用于
    ! 计算常积分
    ! call SimRc(func,a,b,tol,fa,fc,fb,Q,fcount)
    ! func 是外部函数(用于计算被积函数在x出的值call func(f,x))，
    ! a,b分别代表区间[a,b]的左右端点
    ! tol 表示绝对误差的上限
    ! fa,fb,fc 表示func 在a,b,c(a/2+b/2)处的函数值
    ! Q表示积分值，fcount表示所用节点个数
    !_________________________________________________________
    implicit real*8(a-z)
    external func
    integer :: fcount,ka,kb
    h=b-a
    c=(b+a)/2d0
    call func(fd,(a+c)/2d0)
    call func(fe,(c+b)/2d0)
    Q1=h/6d0*(fa+4d0*fc+fb)
    Q2=h/12d0*(fa+4d0*fd+2d0*fc+4d0*fe+fb)
    del=dabs(Q1-Q2)
    if(del<tol) then 
        Q = Q2 + (Q2 - Q1)/15 !Richard 外推法
        fcount = 2;
    else
        call SimpRc(func, a, c, tol, fa, fd, fc, Qa,ka);
        call SimpRC(func, c, b, tol, fc, fe, fb, Qb,kb);
        Q = Qa + Qb;
        fcount = ka + kb + 2
    end if
    end subroutine SimpRc
    
    
    subroutine fun1(f,x)
    !___________________________________________
    !
    !fun1用于计算函数f在x处的值
    !call fun1(f,x)
    !____________________________________________
    implicit real*8(a-z)
    f=sin(x)
    end subroutine fun1
    
end module Simpson
    
    
program main
    use Simpson
    integer :: fcount
    real(kind=8) a,b,c,fa,fb,fc,tol,Q
    a=0d0
    b=1d0
    tol=1d-9
    c = (a + b)/2d0
    call fun1(fa,a)
    call fun1(fc,c)
    call fun1(fb,b)
    call SimpRc(fun1, a, b, tol, fa, fc, fb,Q,k)
    fcount = k + 3
    res=Q-(cos(0d0)-cos(1d0))
    write(*,*)res,fcount
end program main

sslchi

楚香饭 发表于 2014-11-18 22:47
fun1 是需要重写的，fortran 是编译型语言，所有代码都需要编译。而 matlab 是解释型的。

封装在 module  ...

原来Fortran里面常数的双精度也有特别的表示

sslchi

珊瑚虫 发表于 2014-11-18 22:17
[mw_shl_code=fortran,true]module autoSimpson
!----------------------------------------module coment
...

这个不是严格的自适应，只是自动变步长，自适应是自动在某些区间上加密网格，所有的网格都加密显然是没有效率的。

爷且无辜

楼主说的应该是做成.exe文件还能实现改变函数的形式

楚香饭

fun1 也可以作为外部函数，不一定非在 module 里。

你说的每次，是指每次使用，还是每次编译？？换一个函数的时候，你是不想重新编译？

珊瑚虫

sslchi 发表于 2014-11-18 22:34
这个在 调用的时候是不是要重新写 fun1这个函数啊？我的意思是有木有类似于matlab里面函数句柄的东西来作 ...

fun1里面写你打算积分的函数，你说的把matlab混合编程的问题，没搞过，等待高手解答

sslchi

嗯，这个程序的目的是计算积分，但是如果要是求积只作为一个程序的子程序，这么做是否有点不合适呢？总不能每次调用还要去把这个module找到把？

楚香饭

fun1 是需要重写的，fortran 是编译型语言，所有代码都需要编译。而 matlab 是解释型的。

封装在 module 里应该更容易使用才对。

http://fcode.cn/code_gen-25-1.html 这个模块或许是你想要的。

PS：你提到的精度问题，应该是受单精度影响。把常数都改成双精度会提供到 E-11 级

[Fortran] 纯文本查看 复制代码

program main 
    !use mkl95_blas
    integer,parameter :: M=100
    real(kind=8)::N,f(M),d(M),res,s
    s=sin(0.0D0)!计算在a处的值
    N=M/1.0D0   !用于计算步长（步长之倒数）
    f= 1.0D0*(/(i,i=1,M,1)/)/N !计算节点
    d=sin(f)              !节点处的函数值
    do counter=1,M,2     ! 计算奇数节点的累加值(0算第0个点)
        s=s+d(counter)*4.0D0/6.0D0
    end do
    do counter=2,M-2,2   ! 将计算偶数节点的值累加(不包括最后一个节点)
        s=s+d(counter)*2.0D0/6.0D0
    end do
    s=s+d(M)/6.0D0             !将最后一个节点加上去。。。就是这里忘了除以6
    s=2.D0*s/N
    res=s-(cos(0.0D0)-cos(1.0D0)) !计算误差
    write(*,*) res
    read(*,*) 
end

sslchi

珊瑚虫 发表于 2014-11-18 22:17
上述代码摘自《fortran95-2003科学计算与工程光盘各章代码》

这个在 调用的时候是不是要重新写 fun1这个函数啊？我的意思是有木有类似于matlab里面函数句柄的东西来作为程序里面的输入参数。这个是将被积分函数封装在module里面了，用起来不方便啊。

珊瑚虫

上述代码摘自《fortran95-2003科学计算与工程光盘各章代码》