第七讲 物理学中定积分是数值计算方法资料

且f(x0)=a（即第一段的左边）、f(x2n)=b（最后一段2n的右 边），xi=a+ih（即第i段的右边）。
2020/9/15
5
根据上式，就可以编写相应的计算Fortran子程序：
Subroutine simp(a,b,f,n,s) h=(b-a)/(2*n) s=0.5*(f(a)-f(b)) do i=1,n
s=s+h*f(a+i*h) end do write(*,*) n,s End 运行结果：100，0.2721891
7
！辛卜生法：2n=100 f(x)=log(1+x)/(1+x**2) write(*,*) 'input a,b,n=?' read*, a,b,n h=(b-a)/(2*n) s=0.5*(f(a)-f(b)) do i=1,n
其面积，即：
si
h[ f 2
(xi1)
f
(xi )]
I
b
a
f
(x)dx
h {[ 2
f
(a)
f
(x1)] [
f
( x1 )
f
(x2 )] [
f
( xn1 )
f
(b)]}
h
n1
[ 2
f
(a)
f
(b)]
i 1
hf
(xi )
注意：利用了x0=a、xn=b， !梯形积分子程序：
且xi=a+ih。
I
b
f (x)dx
x2n f (x)dx
a
x0
h {[ f 3
(x0 ) 4 f
(x1)
f
(x2 )] [ f
(x2 ) 4 f
(x3 )
f
(x4 )]
[ f (x2n2 ) 4 f (x2n1) f (x2n )]}
h[ 3
f
(a)
f
(b) 4[
f
(x1)
f
(x3)
s1=f(a+(2*i-1)*h) s2=f(a+2*i*h) s=s+2*s1+s2 end do s=(b-a)*s/(3*n) write(*,*) 2*n,s End 运行结果：100，0.2721982
2020/9/15
8
二、计算数据保存到某个文件中的方法
表控输出“Write(*,*) S”语句或“print*,S”语句，是把计 算结果“S”输出到显示器上。 但是很多时候需要把计算数据保存到某个文件之中，这时需 要利用Open语句+write语句配套语句实现：
s1=f(a+i*h) s=s+h*s1 end do write(*,*) n,s End 运行结果：100，0.2739220
2020/9/15
！梯形法：n=100 f(x)=log(1+x)/(1+x**2) write(*,*) 'input a,b,n=?' read*, a,b,n h=(b-a)/n s=0.5*h*(f(a)+f(b)) do i=1,n-1
“近似”，如右图所示： 该四边形的面积为：
设抛物线为y Ax2 Bx C
si
xi f (x)dx
xi2
xi ( Ax2 Bx C)dx
xi2
(1 Ax3 1 Bx2 Cx) xi
(i 1)
i
3
2
xi2
xi2
xi1
பைடு நூலகம்
1 6
( xi
xx2 )[(
Axi22
Bxi2
C)
( Axi2
Bxi
s1=f(a+(2*i-1)*h) s2=f(a+2*i*h) s=s+2*s1+s2 end do s=(b-a)*s/(3*n) end
2020/9/15
6
例题：运用3种基本数值积分方法计算
S
1
0
ln(1 x) 1 x2
dx
！矩形法：n=100 f(x)=log(1+x)/(1+x**2) write(*,*) 'input a,b,n=?' read*, a,b,n h=(b-a)/n s=0.0 do i=1,n
设备号
文件名
Open(8，file=‘jfz.dat’,status=‘unknown’)
┆ write(8,*) x,s
而且Open语句放在类型说明之后，数据文件默认放在相应的 项目工作间的文件夹之中。
f
( x2 n1 )]
2[ f (x2 ) f (x4 ) f (x2n2 )]
h {[ f (a) 3
n
f (b)] 4
i 1
n1
f (x2i1) 2
i 1
f (x2i )}
2h {0.5[ f 3
(a)
f
(b)]
n
[2 f
i 1
(x2i1)
f
(x2i )]}
注意：
积分区间要分成偶数(2n)小段，每段长度h=(b-a)/(2n)，而
C)
A(xi2 xi )2 2B(xx2 xi ) 4C]
x xi
xi2 a (i 2)h, xi1 a (i 1)h, xi a ih
xi xi2 2h, xi xi2 2xi1
h
2020/9/15
si
[ 3
f
(xi2 ) 4 f
(xi1)
f
(xi )]
4
因此，积分近似为：
第七讲 Fortran中的数值积分
物理问题的求解，常常涉及到积分，因此本节介绍在Fortran 中进行数值积分的方法。
一、基本的数值积分方法
y
定积分：
f (x)
b
I a f (x)dx
12 i n
表示f(x)正下方与x轴所围的“面积”。O a
bx
xi1 [a (i 1)h] xi (a ih)
subroutine stx(f,a,b,n,s)
h=(b-a)/n
于是,得到计算的子程序：
s=0.5*h(f(a)+f(b)) do i=1,n-1
s=s+h*f(a+i*h)
end do
end
2020/9/15
3
（3）辛卜生（Simpson）法
这种方法把“每2个小块”用“1个抛物线为顶的四边形”来
数值求该面积：一般都是将积分区间[a,b]分成有限多个小段 （每段宽h=(b-a)/n），然后“近似”求出每小段的“面积”， 再把它们加起来，以此作为积分的近似值。
2020/9/15
1
主要有3种对小段面积的近似方法：
（1）矩形法 这种方法把每小段用宽h、高f(a+ih)的“矩形”来近似，即：
si h f (a ih)
I
b
a
f
(x)dx
n
si
n
[h
f
(a
ih)]
i 1
i 1
据此，计算的函数子程 序程序如：
！矩形积分子程序： function sjx(f,a,b,n) h=(b-a)/n
sjx=0
do i=1,n
sjx=sjx+h*f(a+i*h)
end do
2020/9/15
end
2
（2）梯形法
每小段用“梯形”---上底f(xi-1)、下底f(xi)、高h，“近似”