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

sslchi · 发表于 2014-11-18 22:05:32

本帖最后由 sslchi 于 2014-11-19 22:24 编辑

写了一个简短的代码，用辛普森公式计算积分，一开始的时候觉得误差很大
后来发现是算法中的最后一步不对，文中有注释。当然这是最简单的一个，使用不方便。有没有更简洁、更容易调用的代码？
自适应格式需要用到递归(貌似递归很费时啊)。如果可以的话，怎么才使程序以一个函数为输入量？？
这样的话就可以很方便的调用了（类似于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

珊瑚虫 · 发表于 2014-11-18 22:15:36

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

module Simpson
!----------------------------------------module coment
!  Version     :  V1.0    
!  Coded by    :  syz 
!  Date        :  2010-5-29
!-----------------------------------------------------
!  Description :   复合Simpson法模块
!    
!-----------------------------------------------------
!  Contains    :
!      1.     方法函数
!      2.     测试函数
!-----------------------------------------------------
!  Post Script :
!      1.
!      2. 
!-----------------------------------------------------

contains 

subroutine solve(func,s,a,b,n)
!---------------------------------subroutine  comment
!  Version   :  V1.0    
!  Coded by  :  syz 
!  Date      :  2010-5-29
!-----------------------------------------------------
!  Purpose   :  复合Simpson公式计算数值积分
!    
!-----------------------------------------------------
!  Input  parameters  :
!       1.         func 为外部子程序
!       2.         
!       3.         a,b积分区间
!       4.         n 区间划分个数
!  Output parameters  :
!       1.         s 积分结果
!       2.
!  Common parameters  :
!
!----------------------------------------------------
!  Post Script :
!       1.
!       2.
!----------------------------------------------------
implicit real*8(a-z)
external func
integer::n,k

s=0d0
h=(b-a)/n/2d0

call func(f1,a)
call func(f2,b)

s=f1+f2

!k=0 情况
call func(f1,a+h)
s=s+4d0*f1

do k=1,n-1

t1=a+(2d0*k+1)*h

t2=a+2d0*k*h

call func(f3,t1)
call func(f4,t2)

s=s+f3*4d0+f4*2d0  

end do

s=s*h/3d0

end subroutine solve


subroutine fun1(f,x)
!---------------------------------subroutine  comment
!  Version   :  V1.0    
!  Coded by  :  syz 
!  Date      :  2010-5-29
!-----------------------------------------------------
!  Purpose   :  需要计算的函数
!    
!-----------------------------------------------------
!  Input  parameters  :
!       1.     x  自变量
!       2. 
!  Output parameters  :
!       1.     f  因变量
!       2.
!  Common parameters  :
!
!----------------------------------------------------
!  Post Script :
!       1.
!       2.
!----------------------------------------------------
implicit real*8(a-z)

f=x**2+dsin(x)
end subroutine fun1


end module Simpson



program main
!--------------------------------------program comment
!  Version   :  V1.0    
!  Coded by  :  syz 
!  Date      :  2010-5-29
!-----------------------------------------------------
!  Purpose   :  复合Simpson法计算数值积分主函数
!    
!-----------------------------------------------------
!  In put data  files :
!       1.
!       2.
!  Output data files  :
!       1.   result.txt计算结果
!       2.
!-----------------------------------------------------
!  Post Script :
!       1.  可以通过公式预先估计在一定的精度下需要划分
!           多少节点
!-----------------------------------------------------
use Simpson

implicit real*8(a-z)


open(unit=11,file='result.txt')

write(11,101)

call solve(fun1,s,-2d0,2d0,40)

write(11,102)s
call solve(fun1,s,-2d0,2d0,80)

write(11,103)s

call solve(fun1,s,-2d0,2d0,200)

write(11,104)s


101 format(/,T5,'复合Simpson法计算数值积分',/)
102 format(T5,'n=40',F15.8)
103 format(T5,'n=80',F15.8)
104 format(T5,'n=200',F14.8)

end program main

不知道这个对你有没有用，没看懂你要问什么

珊瑚虫 · 发表于 2014-11-18 22:17:00

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

module autoSimpson
!----------------------------------------module coment
!  Version     :  V1.0    
!  Coded by    :  syz 
!  Date        :  2010-5-29
!-----------------------------------------------------
!  Description :   自动变步长Simpson法模块
!    
!-----------------------------------------------------
!  Contains    :
!      1.     方法函数
!      2.     测试函数
!-----------------------------------------------------
!  Post Script :
!      1.
!      2. 
!-----------------------------------------------------

contains 

subroutine solve(func,s,a,b,tol,n)
!---------------------------------subroutine  comment
!  Version   :  V1.0    
!  Coded by  :  syz 
!  Date      :  2010-5-29
!-----------------------------------------------------
!  Purpose   :  自动变步长Simpson积分方法函数
!    
!-----------------------------------------------------
!  Input  parameters  :
!       1.  func 外部函数
!       2. 
!       3.  a,b积分区间
!       4.  tol  积分误差容限
!
!  Output parameters  :
!       1.   s  积分结果
!       2.   n  实际区间划分个数
!  Common parameters  :
!
!----------------------------------------------------
!  Post Script :
!       1.  需要调用复合Simpson公式
!       2.
!----------------------------------------------------
implicit real*8(a-z)
external func
integer::n,i

!初始划分40个子区间
n=40

!最大允许划分20次
do i=1,20
call simp(func,s,a,b,n)

n=n*2
call simp(func,s1,a,b,n)

del=dabs(s-s1)

!满足精度后就停止循环
if(del<tol)  exit
end do

s=s1

end subroutine solve

subroutine simp(func,s,a,b,n)
!---------------------------------subroutine  comment
!  Version   :  V1.0    
!  Coded by  :  syz 
!  Date      :  2010-5-29
!-----------------------------------------------------
!  Purpose   :  复合Simpson公式计算数值积分
!    
!-----------------------------------------------------
!  Input  parameters  :
!       1.         func 为外部子程序
!       2.         
!       3.         a,b积分区间
!       4.         n 区间划分个数
!  Output parameters  :
!       1.         s 积分结果
!       2.
!  Common parameters  :
!
!----------------------------------------------------
!  Post Script :
!       1.
!       2.
!----------------------------------------------------
implicit real*8(a-z)
external func
integer::n,k

s=0d0
h=(b-a)/n/2d0

call func(f1,a)
call func(f2,b)

s=f1+f2

!k=0 情况
call func(f1,a+h)
s=s+4d0*f1

do k=1,n-1

t1=a+(2d0*k+1)*h

t2=a+2d0*k*h

call func(f3,t1)
call func(f4,t2)

s=s+f3*4d0+f4*2d0  

end do

s=s*h/3d0

end subroutine simp


subroutine fun1(f,x)
!---------------------------------subroutine  comment
!  Version   :  V1.0    
!  Coded by  :  syz 
!  Date      :  2010-5-29
!-----------------------------------------------------
!  Purpose   :  需要计算的函数
!    
!-----------------------------------------------------
!  Input  parameters  :
!       1.     x  自变量
!       2. 
!  Output parameters  :
!       1.     f  因变量
!       2.
!  Common parameters  :
!
!----------------------------------------------------
!  Post Script :
!       1.
!       2.
!----------------------------------------------------
implicit real*8(a-z)

f=1d0/(x**3-2*x-5)
end subroutine fun1


end module autoSimpson



program main
!--------------------------------------program comment
!  Version   :  V1.0    
!  Coded by  :  syz 
!  Date      :  2010-5-29
!-----------------------------------------------------
!  Purpose   :  自动变步长Simpson法计算数值积分主函数
!    
!-----------------------------------------------------
!  In put data  files :
!       1.
!       2.
!  Output data files  :
!       1.   result.txt计算结果
!       2.
!-----------------------------------------------------
!  Post Script :
!       1.  
!-----------------------------------------------------
use autoSimpson

implicit real*8(a-z)
integer::n

open(unit=11,file='result.txt')

write(11,101)

call solve(fun1,s,0d0,2d0,1d-7,n)

write(11,102)n,s




101 format(/,T5,'自动变步长Simpson法计算数值积分',/)
102 format(T5,'区间划分等份为：',I5,/,T5,'积分结果：',F15.8)


end program main

自适应的

珊瑚虫 · 发表于 2014-11-18 22:17:45

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

sslchi · 发表于 2014-11-18 22:34:38

本帖最后由 sslchi 于 2014-11-18 22:37 编辑

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

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

楚香饭 · 发表于 2014-11-18 22:47:47

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 23:27:22

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

珊瑚虫 · 发表于 2014-11-18 23:27:27

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

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

楚香饭 · 发表于 2014-11-19 08:32:11

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

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

爷且无辜 · 发表于 2014-11-19 16:44:53

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

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

爱心勋章

管理勋章

新人勋章

水王勋章

元老勋章

热心勋章