1、2018/4/3,1,一、基本的数值积分方法,数值求该面积：一般都是将积分区间a,b分成有限多个小段（每段宽h=(b-a)/n），然后“近似”求出每小段的“面积”，再把它们加起来，以此作为积分的近似值。,2018/4/3,2,主要有3种对小段面积的近似方法：,2018/4/3,3,注意：利用了x0=a、xn=b，且xi=a+ih。,2018/4/3,4,（3）辛卜生（Simpson）法,该四边形的面积为：,2018/4/3,5,因此，积分近似为：,注意：积分区间要分成偶数(2n)小段，每段长度h=(b-a)/(2n)，而且f(x0)=a（即第一段的左边）、f(x2n)=b（最后一段2n的右边）
2、，xi=a+ih（即第i段的右边）。,2018/4/3,6,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+s2end do s=(b-a)*s/(3*n)end,根据上式，就可以编写相应的计算Fortran子程序：,例题：运用3种基本数值积分方法计算,！矩形法：n=100f(x)=log(1+x)/(1+x*2)write(*,*) input a,b,n=?read*, a,b,nh=(b-a)/ns=0.0do i=1,n s1
3、=f(a+i*h) s=s+h*s1end dowrite(*,*) n,sEnd运行结果：100，0.2739220,！梯形法：n=100f(x)=log(1+x)/(1+x*2)write(*,*) input a,b,n=?read*, a,b,nh=(b-a)/ns=0.5*h*(f(a)+f(b)do i=1,n-1 s=s+h*f(a+i*h)end dowrite(*,*) n,sEnd运行结果：100，0.2721891,2018/4/3,7,！辛卜生法：2n=100f(x)=log(1+x)/(1+x*2)write(*,*) input a,b,n=?read*, a,b,
4、nh=(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+s2end do s=(b-a)*s/(3*n)write(*,*) 2*n,sEnd运行结果：100，0.2721982,2018/4/3,8,2018/4/3,9,表控输出“Write(*,*) S”语句或“print*,S”语句，是把计算结果“S”输出到显示器上。但是很多时候需要把计算数据保存到某个文件之中，这时需要利用Open语句+write语句配套语句实现：,二、计算数据保存到某个文件中的方法,而且Open语句放在类型说明之后
5、，数据文件默认放在相应的项目工作间的文件夹之中。,2018/4/3,10,计算物理课程考试题（1）：,2018/4/3,11,求xy平面上任一点M(x,y)的电势V(x,y)：,计算V(x,y)的Fortran程序：,Program mainDouble precision t1,t2,x,y,VmOpen(9,file=Vm.dat,status=unknown)T1=3.14/3T2=6.28/3Do i=-300,300,1X=0.01*i Do j=-300,300,1 Y=0.01*j Call Vsimp(t1,t2,x,y,Vm) Write(*,*) x,y,Vm Write(
6、9,*) x,y,Vm End doEnd doEnd,Function f(t,x,y) Double precision t,x,y,fF=1.0/dsqrt(x-dcos(t)*2+(y-dsin(t)*2)End,Subroutine Vsimp(a,b,x,y,s)Double precision a,b,x,y,s,f,h,s1,s2External fN=1000H=(b-a)/(2*n)S=0.5*(f(a,x,y)-f(b,x,y)Do i=1,n S1=f(a+(2*i-1)*h),x,y) S2=f(a+2*i*h),x,y) S=s+2*s1+s2End do S=(b
7、-a)*s/(3*n)End,绘制V(x,y)的分布：,（1）画出V(x,y)的分布：3维图；-3x,y3,（2）画出V(x,y)沿着x方向的分布：2维图；取y=0.4,0.6, 1.2,1.4;,（3）画出V(x,y)沿着y方向的分布：2维图；取x=-0.6,0.6, -0.8,0.8;,14,计算物理课程考试题（2）：,15,求xy平面上任一点M(x,y)的磁场B(x,y)：,Program mainDouble precision t1,t2,x,y,BmOpen(9,file=Bm.dat,status=unknown)T1=3.14/3T2=6.28/3Do i=-300,300,1
8、X=0.01*i Do j=-300,300,1 Y=0.01*j Call Bsimp(t1,t2,x,y,Bm) Write(*,*) x,y,Bm Write(9,*) x,y,Bm End doEnd doEnd,16,Function f(t,x,y) Double precision t,x,y,fF=(1.0-x*dcos(t)-y*dsin(t)/(dsqrt(x-dcos(t)*2+(y-dsin(t)*2)*3)End,Subroutine Bsimp(a,b,x,y,s)Double precision a,b,x,y,s,f,h,s1,s2External fN=100
9、0H=(b-a)/(2*n)S=0.5*(f(a,x,y)-f(b,x,y)Do i=1,n S1=f(a+(2*i-1)*h),x,y) S2=f(a+2*i*h),x,y) S=s+2*s1+s2End do S=(b-a)*s/(3*n)End,计算B(x,y)的Fortran程序：,绘制B(x,y)的分布：,（1）画出B(x,y)的分布：3维图；-3x,y3,（2）画出B(x,y)沿着x方向的分布：2维图；取y=0.4,0.6, 1.2,1.4;,（3）画出B(x,y)沿着y方向的分布：2维图；取x=-0.6,0.6, -0.8,0.8;,18,图2 磁场沿着水平方向的分布情况。,图3.竖直方向磁场的分布。,
