【添加新的自适应程序】这个简单的辛普森算法对吗？
本帖最后由 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
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
不知道这个对你有没有用，没看懂你要问什么 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
自适应的 上述代码摘自《fortran95-2003科学计算与工程光盘各章代码》 本帖最后由 sslchi 于 2014-11-18 22:37 编辑
珊瑚虫 发表于 2014-11-18 22:17
上述代码摘自《fortran95-2003科学计算与工程光盘各章代码》
这个在 调用的时候是不是要重新写 fun1这个函数啊？我的意思是有木有类似于matlab里面函数句柄的东西来作为程序里面的输入参数。这个是将被积分函数封装在module里面了，用起来不方便啊。 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 嗯，这个程序的目的是计算积分，但是如果要是求积只作为一个程序的子程序，这么做是否有点不合适呢？总不能每次调用还要去把这个module找到把？ sslchi 发表于 2014-11-18 22:34
这个在 调用的时候是不是要重新写 fun1这个函数啊？我的意思是有木有类似于matlab里面函数句柄的东西来作 ...
fun1里面写你打算积分的函数，你说的把matlab混合编程的问题，没搞过，等待高手解答 fun1 也可以作为外部函数，不一定非在 module 里。
你说的每次，是指每次使用，还是每次编译？？换一个函数的时候，你是不想重新编译？ 楼主说的应该是做成.exe文件还能实现改变函数的形式
页:
[1]
2