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

1、2021-9-261物理问题的求解，常常涉及到积分，因此本节介绍在物理问题的求解，常常涉及到积分，因此本节介绍在FortranFortran中进行中进行的方法。的方法。数值求该面积：数值求该面积：一般都是将积分区间一般都是将积分区间a,ba,b分成有限多个分成有限多个（每段宽（每段宽h=(b-a)/nh=(b-a)/n），然后），然后“”求出每小段的求出每小段的“面积面积”，再把它们加起来，以此作为积分的近似值。再把它们加起来，以此作为积分的近似值。badxxfI)(定积分：定积分：表示表示f(x)f(x)正下方与正下方与x x轴所围的轴所围的“面积面积”。ab)(xfyxO12ni) 1(1

2、hiaxi)(ihaxi2021-9-262这种方法把每小段用宽这种方法把每小段用宽h h、高、高f(a+ih)f(a+ih)的的“矩形矩形”来近似，即：来近似，即：)(ihafhsininiibaihafhsdxxfI11)()(据此，计算的函数子程据此，计算的函数子程序程序如：序程序如：！矩形积分子程序矩形积分子程序：function sjx(function sjx(f f,a,b,n),a,b,n)h=(b-a)/nh=(b-a)/nsjx=0sjx=0do i=1,ndo i=1,n sjx=sjx+h sjx=sjx+h* *f(a+if(a+i* *h)h)end doend d

3、oendend2021-9-263（2 2）梯形法）梯形法每小段用每小段用“梯形梯形”-上底上底f(xf(xi-1i-1) )、下底、下底f(xf(xi i) )、高、高h h，“近似近似”其面积，即：其面积，即：)()(21iiixfxfhs于是于是, ,得到计算的子程序：得到计算的子程序：! !梯形积分子程序梯形积分子程序：subroutine stx(subroutine stx(f f,a,b,n,s),a,b,n,s)h=(b-a)/nh=(b-a)/ns=0.5s=0.5* *h h(f(a)+f(b)(f(a)+f(b)do i=1,n-1do i=1,n-1 s=s+h s=s

4、+h* *f(a+if(a+i* *h)h)end doend doendend111211)()()(2)()()()()()(2)(niinbaxhfbfafhbfxfxfxfxfafhdxxfI注意注意：利用了：利用了x x0 0=a=a、x xn n=b=b，且且x xi i=a+ih=a+ih。2021-9-264（3 3）辛卜生（）辛卜生（SimpsonSimpson）法）法该四边形的面积为：该四边形的面积为：4)(2)()()(61)2131()()(22222222232222CxxBxxACBxAxCBxAxxxCxBxAxdxCBxAxdxxfsixiiiiiixixxxx

5、xxiiiiiii这种方法把这种方法把“每每2 2个小块个小块”用用“1 1个抛物线为顶的四边形个抛物线为顶的四边形”来来“近似近似”，如右图所示：，如右图所示：) 1( ii1ix2ixixxCBxAxy2设抛物线为ihaxhiaxhiaxiii,) 1(,)2(121222,2iiiiixxxhxx)()(4)(312iiiixfxfxfhs2021-9-265)()(2)()( 5 . 032 )(2)(4)()(3)()()( 2)()()( 4)()(3)()(4)()()(4)()()(4)(3)()(1212111212224212312122243221020niiininii

6、innnnnxxbaxfxfbfafhxfxfbfafhxfxfxfxfxfxfbfafhxfxfxfxfxfxfxfxfxfhdxxfdxxfIn因此，积分近似为：因此，积分近似为：积分区间要分成积分区间要分成偶数偶数(2n)(2n)小段，每段长度小段，每段长度h=(b-a)/(2n)h=(b-a)/(2n)，而，而且且f(xf(x0 0)=a)=a（即第一段的左边）、（即第一段的左边）、f(xf(x2n2n)=b)=b（最后一段（最后一段2n2n的右的右边），边），x xi i=a+ih=a+ih（即第（即第i i段的右边）。段的右边）。2021-9-266Subroutine simp(

7、a,b,f,n,s)Subroutine simp(a,b,f,n,s)h=(b-a)/(2h=(b-a)/(2* *n)n)s=0.5s=0.5* *(f(a)-f(b)(f(a)-f(b)do i=1,ndo i=1,n s1=f(a+(2 s1=f(a+(2* *i-1)i-1)* *h)h) s2=f(a+2 s2=f(a+2* *i i* *h)h) s=s+2 s=s+2* *s1+s2s1+s2end doend do s=(b-a) s=(b-a)* *s/(3s/(3* *n)n)endend根据上式，就可以编写相应的计算根据上式，就可以编写相应的计算FortranFortr

8、an子程序：子程序：例题例题：运用：运用3 3种基本数值积分方法计算种基本数值积分方法计算1021)1ln(dxxxS！矩形法：！矩形法：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=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)/

9、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.27218912021-9-267！辛卜生法：！辛卜生法：2n=1002n=100f(x)=log(1+x)/(1+xf(x)=log(1+x)/(1+x* * *2)2)write(write(* *, ,* *) input a,b,n=?) input a,b,n=?readread* *, a,b,n, a,b,nh=(b-a)/(2h=(b-a)/(2* *n)n)s=0.5s=0.5* *(f(a)-f(b)(f(a)

10、-f(b)do i=1,ndo i=1,n s1=f(a+(2 s1=f(a+(2* *i-1)i-1)* *h)h) s2=f(a+2 s2=f(a+2* *i i* *h)h) s=s+2 s=s+2* *s1+s2s1+s2end doend do s=(b-a) s=(b-a)* *s/(3s/(3* *n)n)write(write(* *, ,* *) 2) 2* *n,sn,sEndEnd运行结果：运行结果：100100，0.27219820.27219822021-9-2682021-9-269表控输出表控输出“Write(Write(* *, ,* *) S”) S”语句或语

11、句或“printprint* *,S”,S”语句，是把语句，是把计算结果计算结果“S”S”输出到显示器上。输出到显示器上。但是很多时候需要把计算数据保存到某个文件之中，这时需但是很多时候需要把计算数据保存到某个文件之中，这时需要利用要利用配套语句实现：配套语句实现：Open(Open(8 8，file=jfz.dat,status=unknown)file=jfz.dat,status=unknown)write(write(8 8, ,* *) x,s) x,s设备号设备号文件名文件名而且而且OpenOpen语句语句，。2021-9-2610 xyROdldq),(yxMr),(ba),(y

12、xV线电荷密度为线电荷密度为、半径为、半径为R R、圆心角为、圆心角为的均匀带电圆弧电势的的均匀带电圆弧电势的数值模拟，如图所示。要求：数值模拟，如图所示。要求：1 1、推出计算磁场、推出计算磁场 的公式，编写出的公式，编写出FortranFortran计算程序；计算程序；2 2、用、用OriginOrigin绘制电势的分布图；绘制电势的分布图；注意：注意：写上班级、学号、名字；写上班级、学号、名字；( (双面双面) )打印（控制在打印（控制在1 1张张A4A4纸）；纸）；由学习委员统一交；由学习委员统一交；可取：可取：),(yxV3,3,3, 1, 140yxR2021-9-26112222

13、0)sin()cos(4),(RyRxdRdVyxVM圆弧2200)sin()cos(44RyRxRdrdldVM22)()(byaxrRddlRbRa,sin,cos求求xyxy平面上任一点平面上任一点M(x,y)M(x,y)的电势的电势V(x,y)V(x,y)：计算计算V(x,y)V(x,y)的的FortranFortran程序：程序：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,3

14、00,1 Y=0.01*j Call Vsimp(t1,t2,x,y,Vm) Write(*,*) x,y,Vm Write(9,*) x,y,Vm End doEnd doEndFunction 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

15、=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 绘制绘制V(x,y)V(x,y)的分布：的分布：（1 1）画出）画出V(x,y)V(x,y)的分布：的分布：3 3维图；维图；-3x,y-3x,y3（2 2）画出）画出V(x,y)V(x,y)沿着沿着x x方向的分布：方向的分布：2 2维图；取维图；取y=0.4,0.6, y=0.4,0.6, 1.2,1.4;1.2,1.4;（3 3）画出）画出V(x,y)V(x,y)沿着沿着y y方向的分布：方向的分布：2 2维图；取维图；取x=-0.6,

16、0.6, x=-0.6,0.6, -0.8,0.8;-0.8,0.8;14xyROlId),(baMr),(yxB一段载流为一段载流为I I、半径为、半径为R R、圆心角为、圆心角为的圆弧导线磁场的数值模的圆弧导线磁场的数值模拟，如图所示。要求：拟，如图所示。要求：1 1、推出计算磁场、推出计算磁场 的公式，编写出的公式，编写出FortranFortran计算程序；计算程序；2 2、用、用OriginOrigin绘制磁场的分布图；绘制磁场的分布图；注意：注意：写上班级、学号、名字；写上班级、学号、名字；( (双面双面) )打印（控制在打印（控制在1 1张张A4A4纸）；纸）；由学习委员统一交；

17、由学习委员统一交；可取：可取：B3,3,3, 1, 140yxRI15222/3220)sin()cos()sincos(4),(dRbRabaRkIRbaBjdyidxl dybxarjybixar22)()()()(kdyxadxybybxaIybxadydxkjirIrl drIBd)()()()(40)()(0442/32203030dRdydRdxRyRxcos,sinsin,cos求求xyxy平面上任一点平面上任一点M(x,y)M(x,y)的磁场的磁场B(x,y)B(x,y)：Program mainDouble precision t1,t2,x,y,BmOpen(9,file=

18、 Bm.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 Bsimp(t1,t2,x,y,Bm) Write(*,*) x,y,Bm Write(9,*) x,y,Bm End doEnd doEnd16Function 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