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

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

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

写了一个简短的代码,用辛普森公式计算积分,一开始的时候觉得误差很大
后来发现是算法中的最后一步不对,文中有注释。当然这是最简单的一个,使用不方便。有没有更简洁、更容易调用的代码?
自适应格式需要用到递归(貌似递归很费时啊)。如果可以的话,怎么才使程序以一个函数为输入量??
这样的话就可以很方便的调用了(类似于matlab里面的函数句柄)。
另外的一个问题是,我发现精度只能达到 e-8 左右,这又是什么原因呢? 如何能够使精度达到最高呢(相同的节点个数)
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),欢迎将修改精简的代码
和我分享module Simpson
    !________________________________________________________
    ! Code by sslchi @ 2014
    ! 您可以随意使用、修改,同时欢迎将改进精简的代码和我分享
    !______________________________________________________
    contains
   
   
    recursive subroutine SimpRc(func,a,b,tol,fa,fc,fb,Q,fcount)
    !_________________________________________________________
    ! 递归子程序,用复化辛普森公式计算函数func在区间上的
    ! 黎曼积分,但是该程序没有提供判断积分是否可积,仅能用于
    ! 计算常积分
    ! call SimRc(func,a,b,tol,fa,fc,fb,Q,fcount)
    ! func 是外部函数(用于计算被积函数在x出的值call func(f,x)),
    ! 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

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)
!---------------------------------subroutinecomment
!Version   :V1.0   
!Coded by:syz
!Date      :2010-5-29
!-----------------------------------------------------
!Purpose   :复合Simpson公式计算数值积分
!   
!-----------------------------------------------------
!Inputparameters:
!       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)
!---------------------------------subroutinecomment
!Version   :V1.0   
!Coded by:syz
!Date      :2010-5-29
!-----------------------------------------------------
!Purpose   :需要计算的函数
!   
!-----------------------------------------------------
!Inputparameters:
!       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 datafiles :
!       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

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)
!---------------------------------subroutinecomment
!Version   :V1.0   
!Coded by:syz
!Date      :2010-5-29
!-----------------------------------------------------
!Purpose   :自动变步长Simpson积分方法函数
!   
!-----------------------------------------------------
!Inputparameters:
!       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)
!---------------------------------subroutinecomment
!Version   :V1.0   
!Coded by:syz
!Date      :2010-5-29
!-----------------------------------------------------
!Purpose   :复合Simpson公式计算数值积分
!   
!-----------------------------------------------------
!Inputparameters:
!       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)
!---------------------------------subroutinecomment
!Version   :V1.0   
!Coded by:syz
!Date      :2010-5-29
!-----------------------------------------------------
!Purpose   :需要计算的函数
!   
!-----------------------------------------------------
!Inputparameters:
!       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 datafiles :
!       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 级

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文件还能实现改变函数的形式
页: [1] 2
查看完整版本: 【添加新的自适应程序】这个简单的辛普森算法对吗?