第七讲物理学中定积分是数值计算方法课题.ppt

* * 第七讲 Fortran中的数值积分 物理问题的求解，常常涉及到积分，因此本节介绍在Fortran中进行数值积分的方法。 一、基本的数值积分方法 数值求该面积：一般都是将积分区间[a,b]分成有限多个小段（每段宽h=(b-a)/n），然后“近似”求出每小段的“面积”，再把它们加起来，以此作为积分的近似值。 定积分： 表示f(x)正下方与x轴所围的“面积”。 * * 主要有3种对小段面积的近似方法： （1）矩形法 这种方法把每小段用宽h、高f(a+ih)的“矩形”来近似，即： 据此，计算的函数子程序程序如： ！矩形积分子程序： 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 end * * （2）梯形法 每小段用“梯形”---上底f(xi-1)、下底f(xi)、高h，“近似”其面积，即： 于是,得到计算的子程序： !梯形积分子程序： 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 注意：利用了x0=a、xn=b，且xi=a+ih。 * * （3）辛卜生（Simpson）法 该四边形的面积为： 这种方法把“每2个小块”用“1个抛物线为顶的四边形”来“近似”，如右图所示： * * 因此，积分近似为： 注意： 积分区间要分成偶数(2n)小段，每段长度h=(b-a)/(2n)，而且f(x0)=a（即第一段的左边）、f(x2n)=b（最后一段2n的右边），xi=a+ih（即第i段的右边）。 * * Subroutine simp(a,b,f,n,s) h=(b-a)/(2*n) s=0.5*(f(a)-f(b)) do i=1,n 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 根据上式，就可以编写相应的计算Fortran子程序： 例题：运用3种基本数值积分方法计算 ！矩形法：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 s1=f(a+i*h) s=s+h*s1 end do write(*,*) n,s End 运行结果：100，0.2739220 ！梯形法：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 s=s+h*f(a+i*h) end do write(*,*) n,s End 运行结果：100，0.2721891 * * ！辛卜生法：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 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 * * * * 表控输出“Write(*,*) S”语句或“print*,S”语句，是把计算结果“S”输出到显示器上。 但是很多时候需要把计算数据保存到某个文件之中，这时需要利用Open语句+write语句配套语句实现： 二、计算数据保存到某个文件中的方法 Open(8，file=‘jfz.dat’,status=‘unknown’) write(8,*) x,s ┆ 设备号 文件名 而且Open语句放在类型说明之后，数据文件默认放在相应的项目工作间的文件夹之中。 * * 《计算物理》课程考试题（1）： 线电荷密度为λ、半径为R、圆心角为θ的均匀带电圆弧电势的数值模拟，如图所示。要求： 1、推出计算磁场 的公式，编写出Fortran计算程序； 2、用Origin绘制电势的分布图； 注意： 1、写上班级、学号、名字；(双面)打印（控制在1张A4纸）；由学习委员统一交； 2、可取： * * 求xy平面上任一点M(x,y)的电势V(x,y)： 计算V(x,y)的Fortran程序： Program main Double precision t1,t2,x,y,