第三章物理学中定积分的数值计算方法.docx
- 文档编号:7066621
- 上传时间:2023-01-16
- 格式:DOCX
- 页数:15
- 大小:137.64KB
第三章物理学中定积分的数值计算方法.docx
《第三章物理学中定积分的数值计算方法.docx》由会员分享,可在线阅读,更多相关《第三章物理学中定积分的数值计算方法.docx(15页珍藏版)》请在冰豆网上搜索。
第三章物理学中定积分的数值计算方法
第三章物理学中定积分的数值计算方法
一、填空题
1、库仑常数k等于9×109mV/C,真空中的介电常数ε0等于8.85×10-12F/m。
2、对于电量为Q的点电荷,在距离r处产生的电场强度为
。
3、已知定积分
,被积分函数为
,积分区间为
。
将该区间N等分,步长
,用曲线下的虚矩形面积和近似替代积分值,该方法称为矩形法。
积分近似计算公式为
。
4、毕奥—萨伐尔定律所描述的公式为
。
5、玻尔兹曼常数是k=1.38×1023J/K。
6、麦克斯韦速率分布律公式
。
7、在计算物理中求解定积分的方法有辛普森法、龙贝格法、高斯求积法等。
二、简答
1、写出库仑常数、真空中的介电常数和玻尔兹曼常数的值。
答:
库仑常数k=9×109mV/C,真空中的介电常数ε0=8.85×10-12F/m,玻尔兹曼常数是k=1.38×1023J/K。
2、什么是矩形法?
答:
已知定积分
,被积分函数为
,积分区间为
。
将该区间N等分,步长
,用曲线下的虚矩形面积和近似替代积分值,该方法称为矩形法。
积分近似计算公式为
。
3、毕奥—萨伐尔定律和麦克斯韦速率分布律公式。
答:
毕奥—萨伐尔定律所描述的公式为
。
麦克斯韦速率分布律公式
。
三、计算、编程
1、将区间(0,π/2)二等分,分别利用矩形法、梯形法、抛物线法计算积分
。
解:
将区间(0,π/2)二等分。
矩形法计算结果:
梯形法计算结果:
抛物线法计算结果:
若将区间10等分,有G1=1.0764,G2=0.9974,G3=1.000003。
计算程序:
open(1,file=’int.dat’)
write(*,*)’inputa,b,N=?
’
read(*,*)a,b,N
!
method1:
矩形法
y1=0.0
do10j=0,N-1
x1=a
x1=x1+float(j)*(b-a)/float(N)
10y1=y1+f(x1)*(b-a)/float(N)
write(1,*)N,y1
write(*,*)N,y1
!
method2:
梯形法
y2=0.0
do20j=0,N
x2=a
x2=x2+float(j)*(b-a)/float(N)
if(j.eq.0.or.j.eq.N)then
y2=y2+0.5*f(x2)*(b-a)/float(N)
else
y2=y2+f(x2)*(b-a)/float(N)
endif
20continue
write(1,*)N,y2
write(*,*)N,y2
!
method3:
抛物线法(N:
偶数)
y3=0.0
do30j=0,N
x3=a
x3=x3+float(j)*(b-a)/float(N)
if(j.eq.0.or.j.eq.N)then
y3=y3+1./3.*f(x3)*(b-a)/float(N)
else
k=j-2*int(j/2)
if(k.eq.0)then
y3=y3+2./3.*f(x3)*(b-a)/float(N)
else
y3=y3+4./3.*f(x3)*(b-a)/float(N)
endif
endif
30continue
write(1,*)N,y3
write(*,*)N,y3
end
functionf(x)
f=cos(x)
end
2、设有一长直导线均匀带电,线电荷密度为λ,长度为2l。
求空间任意一点P的电势,试利用变步长辛普森求积分法求解。
(取,a=l,b=-l,x0=0.125,y0=0,λ=5×10-9,l=0.05,ε0=8.85×10-12F/m)解:
在导线上取一小段dx,视为点电荷,其电量为λdx,它在P点产生的电势为:
因此有:
利用解析法可得以上定积分结果为:
采用数值法进行计算,则有
计算程序:
subroutinesimp(a,b,ep,s2,f)
h=b-a
t1=h/2.*(f(a)+f(b))
n=1
5s=0.
dok=0,n-1
s=s+f(a+(k+0.5)*h)
enddo
t2=t1/2.+h/2.*s
s2=t2+(t2-t1)/3.
if(n/=1)goto20
15n=n+n
h=h/2.
t1=t2
s1=s2
goto5
20if(abs(s2)<=1.)d=abs(s2-s1)
if(abs(s2)>1.)d=abs((s2-s1)/s2)
if(d>=ep)goto15
return
end
programmain
externalf
callsimp(0.075,-0.075,1e-5,s2,f)
write(*,"(20x,'s=',f10.4)")s2
end
functionf(x)
f=5*10**-9/4*3.14*8.85*10**-12*1./(0.125-x)
end
3、求带电量为Q的均匀带电圆环(半径为r)在轴线上x点的电势(线电荷密度为λ),试利用变步长辛普森求积分法求解(取λ=1,r=1,ε=1×10-5,ε0=8.85×10-12F/m)。
解:
圆环上任一点电荷元λds在轴线上x点的电势大小为
其中
为圆环上点电荷元到x的距离,则
计算程序:
subroutinesimp(a,b,ep,s2,f)
h=b-a
t1=h/2.*(f(a)+f(b))
n=1
5s=0.
dok=0,n-1
s=s+f(a+(k+0.5)*h)
enddo
t2=t1/2.+h/2.*s
s2=t2+(t2-t1)/3.
if(n/=1)goto20
15n=n+n
h=h/2.
t1=t2
s1=s2
goto5
20if(abs(s2)<=1.)d=abs(s2-s1)
if(abs(s2)>1.)d=abs((s2-s1)/s2)
if(d>=ep)goto15
return
end
programmain
externalf
callsimp(0.,6.282,1e-5,s2,f)
write(*,"(20x,'s=',f10.4)")s2
end
functionf(x)
f=1./4*3.14*8.85*10**-12*(x**+1)**1/2
end
4、求带电量为Q的均匀带电圆环(半径为r)在轴线上x点的电场强度(线电荷密度λ),试利用变步长辛普森求积分法求解(取λ=1,r=1,ε=1×10-5,l=3,ε0=8.85×10-12F/m)。
解:
如图所示,圆环上任一点电荷λds在轴线上x点的电场强度dE大小为
根据电荷分布对称性可知,场强沿着与x轴垂直的方向为0,场强仅沿着x方向有值,而dE沿着x方向的投影为
其中
。
对上式做积分可得解析解
其中
。
计算程序:
subroutinesimp(a,b,ep,s2,f)
h=b-a
t1=h/2.*(f(a)+f(b))
n=1
5s=0.
dok=0,n-1
s=s+f(a+(k+0.5)*h)
enddo
t2=t1/2.+h/2.*s
s2=t2+(t2-t1)/3.
if(n/=1)goto20
15n=n+n
h=h/2.
t1=t2
s1=s2
goto5
20if(abs(s2)<=1.)d=abs(s2-s1)
if(abs(s2)>1.)d=abs((s2-s1)/s2)
if(d>=ep)goto15
return
end
programmain
externalf
callsimp(0.,6.282,1e-5,s2,f)
write(*,"(20x,'s=',f10.4)")s2
end
functionf(x)
f=1./4*3.14*8.85*10**-12*x/3**3
end
5、我国第一颗人造地球卫星轨道是一个椭圆,椭圆周长的计算公式为:
其中a是椭圆的半长轴,c是地球的中心与轨道中心(椭圆中心)的距离,地球半径R=6371km,近地点距离h=439km,远地点距离H=2384km,a=7782.5km,c=972.5km。
计算程序:
subroutinesimp(a,b,ep,s2,f)
h=b-a
t1=h/2.*(f(a)+f(b))
n=1
5s=0.
dok=0,n-1
s=s+f(a+(k+0.5)*h)
enddo
t2=t1/2.+h/2.*s
s2=t2+(t2-t1)/3.
if(n/=1)goto20
15n=n+n
h=h/2.
t1=t2
s1=s2
goto5
20if(abs(s2)<=1.)d=abs(s2-s1)
if(abs(s2)>1.)d=abs((s2-s1)/s2)
if(d>=ep)goto15
return
end
programmain
externalf
callsimp(0.,1.5708,1e-5,s2,f)
write(*,"(20x,'s=',f10.4)")s2
end
functionf(x)
f=4*7782.5*sqrt(1.-(972.5/7782.5)**2*sin(x)**2)
end
6、利用两点高斯公式求积分
的近似值。
解:
由于a=0,b=1,因此作变换
计算程序:
externalf
real*8a,b,G
a=0.
b=1.
callgauss(a,b,G,f)
write(*,*)G
end
subroutineGauss(a,b,G,f)
real*8t(3),W(3),a,b,G,x
datat/0.,0.774597,-0.774597/
dataW/0.888889,0.555556,0.555556/
G=0.0
do10i=1,3
x=0.5*((b+a)+(b-a)*t(i))
10G=G+w(i)*f(x)
G=0.5*(b-a)*G
return
end
functionf(x)
real*8f,x
f=sqrt(1.+x**2)
end
7、用4点高斯求积公式编程计算积分
计算程序:
externalf
real*8a1,b1,a2,b2,G
a1=1.
b1=1.5
a2=1.4
b2=2.0
callgauss(a1,b1,a2,b2,G,f)
write(*,*)G
end
subroutineGauss(a1,b1,a2,b2,G,f)
real*8t(4),W(4),a1,b1,a2,b2,G,x,y
datat/0.3399810,-0.3399810,0.8611363,-0.8611363/
dataW/0.6521452,0.6521452,0.3478548,0.348548/
G=0.0
doi=1,3
x=0.5*((b2+a2)+(b2-a2)*t(i))
doj=1,3
y=0.5*((b1+a1)+(b1-a1)*t(j))
G=G+w(i)*w(j)*f(x,y)
enddo
enddo
G=0.25*(b2-a2)*(b1-a1)*G
return
end
functionf(x,y)
real*8f,x,y
f=dlog(x+2.*y)
end
8、利用龙贝格法计算单摆的振荡周期
其中,
计算程序:
subroutineromb(a,b,ep,r2,f)
h=b-a
t1=h/2.*(f(a)+f(b))
n=1
5s=0.0
dok=0,n-1
s=s+f(a+(k+0.5)*h)
enddo
t2=t1/2.+h/2.*s
s2=t2+(t2-t1)/3.
if(n/=1)goto20
15n=n+n
h=h/2
t1=t2
s1=s2
goto5
20c2=s2+(s2-s1)/15.
if(n/=2)goto40
30c1=c2
goto15
40r2=c2+(c2-c1)/63.
if(n/=4)goto60
50r1=r2
goto30
60if(abs(r2)-1.)70,70,80
70if(abs(r2-r1)>=ep)goto50
return
80if(abs((r2-r1)/r2)>=ep)goto50
return
end
programmain
externalf
callromb(0.,1.5708,1e-5,r2,f)
write(*,"(3x,'B=',e12.6)")r2
end
functionf(x)
f=4*sqrt(5/9.8)/(1-sin(10/3.14159)**2*sin(x)**2)**(1/2)
end
9、应用龙贝格法计算积分
要求误差不超过10-5。
计算程序:
subroutineromb(a,b,ep,r2,f)
h=b-a
t1=h/2.*(f(a)+f(b))
n=1
5s=0.0
dok=0,n-1
s=s+f(a+(k+0.5)*h)
enddo
t2=t1/2.+h/2.*s
s2=t2+(t2-t1)/3.
if(n/=1)goto20
15n=n+n
h=h/2
t1=t2
s1=s2
goto5
20c2=s2+(s2-s1)/15.
if(n/=2)goto40
30c1=c2
goto15
40r2=c2+(c2-c1)/63.
if(n/=4)goto60
50r1=r2
goto30
60if(abs(r2)-1.)70,70,80
70if(abs(r2-r1)>=ep)goto50
return
80if(abs((r2-r1)/r2)>=ep)goto50
return
end
programmain
externalf
callromb(0.1,0.6,1e-5,r2,f)
write(*,"(3x,'B=',e20.6)")r2
end
functionf(x)
f=0.02792*(2-x)/(1.449*x+1)**0.8*(1-x)**1.2*x
end
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 第三 物理学 积分 数值 计算方法